[gpaw-users] Inaccuracies in phonon dispersion calculation

[Jens Jørgen Mortensen]

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