[QE-users] Negative phonon frequency at qpoints away from Gamma

Jie Peng pengjiesjtu at gmail.com
Fri Apr 22 22:39:57 CEST 2022


Dear all:

I know it is a question that has been asked so many times and I went over lots of them in the PW forum archive. But I still could not fix the problem.

I am calculating phonon dispersion of alpha-Uranium. First, the unit cell is relaxed to equilibrium followed by a SCF calculation. The input for the SCF calculation is copied below:
&control
    calculation='scf',
    restart_mode='from_scratch',
    pseudo_dir = '/home/peng276/pseudopotential',
    prefix='Uranium',
    outdir='.',
    forc_conv_thr =1E-5
    etot_conv_thr = 1E-6
    tstress=.true.
    tprnfor=.true.
    nstep=100
    !verbosity="high"
 /
 &system
    ibrav = 0,
    celldm(1) = 5.4051
    nat = 2,
    ntyp = 1,
    ecutwfc = 100,
    occupations='smearing',
    degauss=0.005
 /
 &electrons
    conv_thr = 1.0e-8
    electron_maxstep = 100,
    !mixing_beta=0.3
    diagonalization='cg'
 /
 &ions
    ion_dynamics = 'bfgs'
 /
 &cell
    cell_dynamics='bfgs',
    press=0.0,
/
ATOMIC_SPECIES
 U 238.029 U.pbe-spfn-rrkjus_psl.1.0.0.UPF

CELL_PARAMETERS (alat)
   0.490195165  -1.021315393   0.000000000
   0.490195165   1.021315393   0.000000000
   0.000000000   0.000000000   1.711961721

ATOMIC_POSITIONS (crystal)
U             0.9017856414        0.0982143586        0.2500000000
U             0.0982143586        0.9017856414        0.7500000000

K_POINTS {automatic}
  12 12 12 0 0 0

Then I start a phonon calculation on a qpoint mesh 4*4*4. The input ph.in is copied below:
phonons of alpha Uranium
&inputph
prefix = 'Uranium',
amass(1) = 238.029,
outdir = './',
fildvscf='Udv',
fildyn = 'Uranium.dyn',
ldisp=.true.
nq1=4,
nq2=4,
nq3=4,
/

The obtained dispersion is plotted in the figure below, in which I see negative frequency of the first acoustic phonon mode near X. Then I tried to improve computational accuracy by tuning the calculation parameters used in both Pw.x and Ph.x, particularly looking at the phonon frequencies of one of the negative frequency phonons at (0.5,0,0). An example output is copied below:
    Diagonalizing the dynamical matrix

     q = (    0.500000000   0.000000000   0.000000000 ) 

 **************************************************************************
     freq (    1) =      -0.939207 [THz] =     -31.328581 [cm-1]
     freq (    2) =       1.029847 [THz] =      34.352008 [cm-1]
     freq (    3) =       2.289485 [THz] =      76.368998 [cm-1]
     freq (    4) =       2.490969 [THz] =      83.089795 [cm-1]
     freq (    5) =       2.862652 [THz] =      95.487797 [cm-1]
     freq (    6) =       3.345818 [THz] =     111.604479 [cm-1]
 **************************************************************************.
The parameters I tried to tune include:
In SCF calculations:
1) varying Ecutwfc from 80 Ry to 140 Ry
2) varying kpoint grid from 8*8*8 to 16*16*16.
3) varying conv_thr from 1e-6 to 1e-10
4) varying degauss from 0.005 to 0.02 (since alpha-U is a metal, the gauss smearing might affect the energy and therefore force constant calculations)
In Ph.x calculations:
1) varying tr2_ph from 1e-12 to 1e-14
2) varying alpha_mix(1) from 0.3 to 0.7 (In principle this should not matter too much since it is the proportion of old charge density used in updating to the new charge density, no matter what value is used at the physical quantities should be converged in the end)

However, in all these calculations with different combinations of the parameter sets, I am consistently getting negative frequencies at (0.5,0,0). The possible reasons:
1) Unstable structure. It is unlikely this is the cause, the frequency at Gamma even before I impose the acoustic summation rule is very close to 0. Also the forces obtained from SCF calculation on the equilibrium structure (see below) is quite small.
     Forces acting on atoms (cartesian axes, Ry/au):

     atom    1 type  1   force =     0.00000000    0.00001190    0.00000000
     atom    2 type  1   force =    -0.00000000   -0.00001190    0.00000000

     Total force =     0.000017     Total SCF correction =     0.000003
     SCF correction compared to forces is large: reduce conv_thr to get better values


     Computing stress (Cartesian axis) and pressure

          total   stress  (Ry/bohr**3)                   (kbar)     P=        0.05
   0.00000066   0.00000000   0.00000000            0.10        0.00        0.00
   0.00000000  -0.00000037   0.00000000            0.00       -0.05        0.00
   0.00000000   0.00000000   0.00000077            0.00        0.00        0.11

2) Incommensurate qpoint grid with kpoint grid. I have seen some discussions on the PW forum that in metallic systems, the incommensurate qpoint and kpoint grid might lead to issues. However, I have tried using 12*12*12 and 16*16*16 kpoint grid in SCF calculations to do the phonon calculation on a 4*4*4 grid, but still get negative frequency at (0.5,0,0). Moreover, in the case of single qpoint phonon calculation, does it matter what kpoint mesh grid is used in the SCF calculation?

3) Unstable structure of the supercell displaced according to the phonon pattern. I think that something like soft modes could cause this negative frequency. However, the published dispersions of alpha-U all show similar shape to mine, with a sudden drop in the acoustic branch near X along the Gamma-X (see attached figure from a literature), but no negative frequency as in my case
The figure is taken from: Manley, M. E., et al. "Phonon dispersion in uranium measured using inelastic x-ray scattering." Physical Review B 67.5 (2003): 052302.

4) Spin-orbit coupling. In 5f electron systems, spin-orbital coupling might be important to include in calculations. However, there is still debate on whether it is necessary to include SOC in DFT calculations on alpha-U or not. Moreover, I have calculated dispersion of alpha-U in VASP without SOC and agrees with other published results well. So I would not go too deep into this, plus SOC calculations are expensive.


With all the attempts and thoughts I have made, I doo not know how to resolve this issue. I strongly believe it is a numerical issue and currently trying another pseudopotential (Ulatrsolft PBE instead of the PAW PBE used in all above calculations). I would appreciate if anyone could give me suggestions or advice.

Thank you very much!

Best
Jie



————————————————————————————————————————————————
Jie Peng, Postdoctoral research fellow
School of Materials Engineering, Purdue University, West Lafayette, Indiana
Email: pengjiesjtu at gmail.com <mailto:pengjiesjtu at gmail.com>
Cell: (240)-495-9445

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20220422/6128ce1e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: qpoint444.tif
Type: image/tiff
Size: 64459 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20220422/6128ce1e/attachment.tif>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screen Shot 2022-04-22 at 15.23.27.png
Type: image/png
Size: 154332 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20220422/6128ce1e/attachment.png>


More information about the users mailing list