[gpaw-users] Inaccuracies in phonon dispersion calculation
Holden Parks
hparks at andrew.cmu.edu
Tue Nov 14 22:46:39 CET 2017
Dear Jens and Kristen,
Thank you! I am running this calculation as we speak. I should also note
that I tried re-plotting results from old calculations, and it seems as
though the biggest culprit was the force constant cutoff. Setting
'cutoff=5' in the Phonons.read(), as you did in your script, does the trick
wonderfully even for smaller systems (I just checked for 4x4x4 unit cells).
Thank you both so much for your help!
-Holden
On Tue, Nov 14, 2017 at 4:54 AM, Jens Jørgen Mortensen <jjmo at dtu.dk> wrote:
> On 11/13/2017 04:09 PM, Holden Parks wrote:
>
>> 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!
>>
>
> Attached.
>
> Jens Jørgen
>
>
>> -Holden
>>
>> On Fri, Nov 10, 2017 at 2:48 AM, Kristen Kaasbjerg <kkaa at nanotech.dtu.dk
>> <mailto: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
>> <mailto: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>
>> > <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>
>> > <mailto: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>
>> >
>> <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/20171114/706a04bc/attachment.html>
More information about the gpaw-users
mailing list