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

[mkondrin at hppi.troitsk.ru]

Dear Dr. Jie Peng,

Recently I have encountered the similar problem but for pure hypothetical 
structure which I also do not know how to resolve.

Still just as a suggestions I believe that in you case your option 3) is 
the most probable. Pay attention that the experimentally observed structure 
of alpha uranium was investigated at room temperature but your calculation 
was performed at zero temperature. Is not it possible that some kind of 
displacive phase transition take place in between? So, alpha-uranium is 
just a high symmetry disordered phase (aristotype) which at low temperature 
would transform into lower symmetry more energetically favourable phase.

This hypothesis can be checked if you apply cryatl distortions 
corresponding to phonone X-mode to positions of atoms in the unit cell (it 
will result in its doubling) and try to relax it. Would not it lead to 
lowering of energy per atom in comparison to initial (smaller) structure?

Sincerely yours,
M. V. Kondrin

On Apr 23 2022, Jie Peng wrote:

>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
>
>