[gpaw-users] Inaccuracies in phonon dispersion calculation

Holden Parks hparks at andrew.cmu.edu
Mon Nov 13 16:09:34 CET 2017


Hi Jens and Kristen,

Thank you both for your suggestions. Kristen, though I believe my structure
does converge to ~1meV, I was not specifically setting this as a
convergence criterion but will do so from now on. I also agree that using a
cutoff for setting up the dynamical matrix is correct - I had just ignored
it recently because I believe my problems in getting good dispersions were
caused by a more serious problem in the original force calculations.

Jens, would you please send me your script and/or dispersion curve? I would
like to run this calculation as well. Thank you!

-Holden

On Fri, Nov 10, 2017 at 2:48 AM, Kristen Kaasbjerg <kkaa at nanotech.dtu.dk>
wrote:

> 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 <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
> <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/20171113/80482522/attachment.html>


More information about the gpaw-users mailing list