[gpaw-users] Inaccuracies in phonon dispersion calculation

[Holden Parks]

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?

-Holden

On Fri, Nov 3, 2017 at 2:37 AM, Jens Jørgen Mortensen <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
>> 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/20171106/b2c63c26/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: asap_si_dispersion.png
Type: image/png
Size: 57451 bytes
Desc: not available
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20171106/b2c63c26/attachment-0001.png>