[gpaw-users] Inaccuracies in phonon dispersion calculation

[Jens Jørgen Mortensen]

On 11/07/2017 03:20 AM, Holden Parks wrote:
> Hi Jens,
>
> Thank you for the reply. I've tried running the plane-wave 
> calculations with your suggested cutoff but haven't seen an 
> improvement. I've confirmed that inaccuracies in the forces are the 
> cause of the problem in the dispersions, as I've run the same 
> calculation with the Brennan potential (part of the ASAP calculator) 
> and the problematic features aren't present (I've attached it here). 
> Do you have any thoughts on how I should move forward?

I got a nice looking phonon dispersion with a 7x7x7 cell and these 
parameters:

calc = GPAW(symmetry={'point_group': False},
             mode='lcao',
             basis='szp(dzp)',
             convergence={'density': 1e-7},
             xc='PBE',
             occupations=FermiDirac(0.01))

the important thing being a more tight convergence (I think).

Jens Jørgen

>
> -Holden
>
> On Fri, Nov 3, 2017 at 2:37 AM, Jens Jørgen Mortensen <jjmo at dtu.dk 
> <mailto:jjmo at dtu.dk>> wrote:
>
>     On 10/31/2017 11:19 PM, Holden Parks via gpaw-users wrote:
>
>         Dear GPAW users,
>
>         I am having a problem calculating the phonon dispersion for
>         bulk diamond silicon. This calculation is performed in ase,
>         and I at first suspected that something was wrong with the
>         implementation. However, the example with aluminum using EMT
>         on the ase manual works perfectly, so I believe my problem is
>         that I am incorrectly calculating the force constants with my
>         calculator.
>
>         I've been calculating phonon dispersion curves for silicon
>         with a variety of supercell and k-point sizes, grid spacing
>         values, and displacement (delta) values. I've also tried a
>         using GGA xc functionals rather than LDA, and turned
>         point-group symmetry off. However, there are nonzero (and
>         imaginary frequencies) at the gamma point, and the acoustic
>         dispersion curves don't degenerate where they're expected
>         (between the gamma and W points and between the L and gamma
>         points). I've attached an example dispersion plot and the code
>         I used to generate it.
>
>         Any suggestions at this point would be extremely helpful.
>         Thank you!
>
>
>     Try running your calculation in plane-wave mode.  Instead of
>     h=0.16, use something like:
>
>         calc = GPAW(..., mode=PW(300), ...)
>
>     and remember to do "from gpaw import PW".  This should give you
>     better forces.
>
>     Jens Jørgen
>
>         -Holden
>
>         ---------------------
>
>         from ase.build import bulk
>         from ase.dft.kpoints import ibz_points, bandpath
>         from ase.phonons import Phonons
>         from gpaw import GPAW, FermiDirac
>         import numpy as np
>
>         # Setup crystal and EMT calculator
>         atoms = bulk('Si', 'diamond', a=5.477)
>         calc = GPAW(kpts=(2, 2, 2),
>                         symmetry={'point_group': False},
>                         h=0.16,
>                         xc='PBE',
>                         occupations=FermiDirac(0.01))
>
>         ph = Phonons(atoms, calc, supercell=(4, 4, 4), delta=0.05)
>         ph.run()
>
>         # Read forces and assemble the dynamical matrix
>         ph.read(acoustic=True)
>
>         High-symmetry points in the Brillouin zone
>         G = [0, 0, 0]
>         X = [0.5, 0, 0.5]
>         W = [0.5,0.25,0.75]
>         L = [0.5, 0.5, 0.5]
>
>         point_names = ['$\Gamma$', 'X', 'W', 'L', '$\Gamma$']
>         path = [G, X, W, L, G]
>
>         # Band structure in meV
>         freq_conv = 2.41798902*(10**14)
>         path_kc, q, Q = bandpath(path, atoms.cell, 100)
>         omega_kn = freq_conv*ph.band_structure(path_kc)/(10**12)
>
>
>         import matplotlib.pyplot as plt
>         plt.figure(1, (8, 6))
>         plt.axes([.1, .07, .67, .85])
>
>         for n in range(len(omega_kn[0])):
>             omega_n = omega_kn[:, n]
>             plt.plot(q, omega_n, 'k-', lw=2)
>
>         plt.xticks(Q, point_names, fontsize=18)
>         plt.yticks(fontsize=18)
>         plt.xlim(q[0], q[-1])
>         plt.ylabel("Frequency ($\mathrm{THz}$)", fontsize=22)
>         plt.grid('on')
>
>         plt.show()
>
>
>         _______________________________________________
>         gpaw-users mailing list
>         gpaw-users at listserv.fysik.dtu.dk
>         <mailto:gpaw-users at listserv.fysik.dtu.dk>
>         https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users
>         <https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users>
>
>
>