[gpaw-users] Inaccuracies in phonon dispersion calculation

[Kristen Kaasbjerg]

Good point Jens, it is indeed crucial to crank up the convergence criteria in phonon calculations. As the force is the derivative of the total energy, I usually converge that one to ~< 1 meV. In this respect, you should keep in mind that the criteria you specify for the total energy in GPAW is per electron, i.e. you have to convert to the total energy of all your supercell electrons.


Let me also emphasize that it is important to apply the cutoff to the force constant in real space once setting up the dynamical matrix. If you don't do that, the interatomic forces will have longer range in some directions compared to others, which is of course undesirable.


Best,
Kristen




________________________________
From: Jens Jørgen Mortensen
Sent: 10 November 2017 07:18
To: Holden Parks
Cc: gpaw-users at listserv.fysik.dtu.dk; Kristen Kaasbjerg
Subject: Re: [gpaw-users] Inaccuracies in phonon dispersion calculation

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>
>
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20171110/984ea2db/attachment-0001.html>