ABINIT Discussion Forums

Dear Pascal and prof. Gonze,

I have tried your suggestions (changed tolvrs in DSs 4-17 and use anaddb input from tests/v2/Input/t22.ddb.in.gz), but with no success. The 2nd derivations do not change significantly (ie. the phonon frequencies) with increased tolvrs.

I am dried out of ideas how to reproduce the results from the PRL article. If needed, I can send (for any case) again (in the attachment) the inputs I used (and xmgrace plot). The psudos are the same and I used even better ecut and k-point mesh. I don't receive any warnings in the output file during these calcs, except

Total localisation tensor (bohr^2) in cartesian coordinates
WARNING : still subject to testing - especially symmetries.
direction matrix element
alpha beta real part imaginary part
1 1 1.0624706248 0.0000000000
1 2 0.5312353124 0.0000000000
1 3 0.5312353124 0.0000000000
2 1 0.5312353124 0.0000000000
2 2 1.0624706248 0.0000000000
2 3 0.5312353124 0.0000000000
3 1 0.5312353124 0.0000000000
3 2 0.5312353124 0.0000000000
3 3 1.0624706248 0.0000000000

WARNING : Localization tensor in reciprocal space incomplete,
transformation to cartesian coordinates may be wrong.

in DS 2.


I am very grateful for any response.

Thank you!

Yours

Igor Lukacevic

As posted to the mailing lists:

1) try running the rf abinit calculations with rfasr=1
(http://www.abinit.org/documentation/hel ... html#rfasr)
For some extreme pseudopotentials (with large breaking of the asr)
this is crucial: correct asr both in abinit _and_ anaddb.

2) are you really using the _same_ pseudopotential as in the PRL?
Otherwise, try each of your pseudos for ASR breaking.

3) Your phonon BS indeed looks a quite bad. Are you sure that the
structure is correct, and that your reciprocal space lines are aligned
correctly (ie using the primitive unit cell and standard
orientations)? Of course this will not solve your imaginary
frequencies.

4) Look through the _abinit_ outputs for the phonon frequencies: which
of them are negative? Send in your abinit output file.

g'luck

Matthieu

Dear prof. Verstraete,

first I would like to thank You for replying.

mverstra wrote:As posted to the mailing lists:

1) try running the rf abinit calculations with rfasr=1
(http://www.abinit.org/documentation/hel ... html#rfasr)
For some extreme pseudopotentials (with large breaking of the asr)
this is crucial: correct asr both in abinit _and_ anaddb.

I have submitted a rf phonon job with rfasr 1. I'll adapt the anaddb option also.

mverstra wrote:As posted to the mailing lists:

2) are you really using the _same_ pseudopotential as in the PRL?
Otherwise, try each of your pseudos for ASR breaking.

I am 99.99% sure. As prof. Gonze said "If I remember correctly, the one that had been used in the paper you mention
is the one in the table http://www.abinit.org/downloads/psp-lin ... nks/lda_ex". I checked these and they are the same as the ones I use (they're both LDA, extended norm-cons., they have the same number of valence electrons: Zr -> 12, O -> 6). The only thing that comes in mind to me here is that the internal parameters used to create these psps may be different (PRL is from 1998), in the case that it has significance. But how could I check that? I only know the ones I use

40zr.971106.mod

Zr, exnc+self, no core, rcs=1.75, rcp=1.55, rcd=1.7, ecut 25/34 (Ar+3d10)+4s2 4p
40.00000 12.00000 971106 zatom,zion,pspdat
4 3 2 2 2001 0 pspcod,pspxc,lmax,lloc,mmax,r2well
0 0 0 2 1.75484409801703478 l,e99.0,e99.9,nproj,rcpsp
.000 .000 .000 .000 rms,ekb1,ekb2,epsatm
1 0 0 2 1.54863439661046565 l,e99.0,e99.9,nproj,rcpsp
.000 .000 .000 .000 rms,ekb1,ekb2,epsatm
2 0 0 0 1.69654886360351576 l,e99.0,e99.9,nproj,rcpsp
.000 .000 .000 .000 rms,ekb1,ekb2,epsatm
.000 .000 .000 rchrg,fchrg,qchrg

mverstra wrote:As posted to the mailing lists:

3) Your phonon BS indeed looks a quite bad. Are you sure that the
structure is correct, and that your reciprocal space lines are aligned
correctly (ie using the primitive unit cell and standard
orientations)? Of course this will not solve your imaginary
frequencies.

I use the xreds from the Parlinskis' PRL (which the above mentioned PRL addresses to)

#Definition of the unit cell
acell 3*9.46
angdeg 90 90 90
brvltt -1
spgroup 225
#Definition of the atom types
ntypat 2
znucl 40 8 8
#Definition of the atoms
natom 3
typat 1 2 2
xred 0.00 0.00 0.00
0.25 0.25 0.25
0.75 0.75 0.75
#Gives the number of band, explicitely (do not take the default)
nband 14
#Definition of the k-point grid
nshiftk 4
shiftk 0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
ngkpt 3*6

diemac 4

The GS calcs gave me the lattice constant within 2% of the experimental one. Thus, I presumed that the structure is ok.

I'm not sure that I understand what You mean by: "reciprocal space lines are aligned correctly (ie using the primitive unit cell and standard orientations)".

mverstra wrote:As posted to the mailing lists:

4) Look through the _abinit_ outputs for the phonon frequencies: which
of them are negative? Send in your abinit output file.

g'luck

Matthieu

Unfortunately, I cannot attach the ouput, nor any other file. I always get a message: "The extension is not allowed."

The negative ones are at
0.16667 0.16667 0.00000
0.33333 0.33333 0.00000
0.50000 0.50000 0.00000
0.50000 0.33333 0.16667
-0.33333 0.50000 0.16667

Here's the list of q-points used
qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt4 1.66666667E-01 0.00000000E+00 0.00000000E+00
qpt5 3.33333333E-01 0.00000000E+00 0.00000000E+00
qpt6 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt7 1.66666667E-01 1.66666667E-01 0.00000000E+00
qpt8 3.33333333E-01 1.66666667E-01 0.00000000E+00
qpt9 5.00000000E-01 1.66666667E-01 0.00000000E+00
qpt10 -3.33333333E-01 1.66666667E-01 0.00000000E+00
qpt11 -1.66666667E-01 1.66666667E-01 0.00000000E+00
qpt12 3.33333333E-01 3.33333333E-01 0.00000000E+00
qpt13 5.00000000E-01 3.33333333E-01 0.00000000E+00
qpt14 -3.33333333E-01 3.33333333E-01 0.00000000E+00
qpt15 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt16 5.00000000E-01 3.33333333E-01 1.66666667E-01
qpt17 -3.33333333E-01 3.33333333E-01 1.66666667E-01
qpt18 -3.33333333E-01 5.00000000E-01 1.66666667E-01


I hope this will help You to help me further.

Cheers!

Igor Lukacevic

I checked the output file more and searched for the differences in variational and non-variational 2DEtotal. I found this for 2 q points and Gamma:

UNSTABLE q(1/8,1/8,0)
=====================

Thirteen components of 2nd-order total energy (hartree) are
1,2,3: 0th-order hamiltonian combined with 1st-order wavefunctions
kin0= 1.65000102E+03 eigvalue= 2.00776545E+02 local= -9.29626310E+02
4,5,6: 1st-order hamiltonian combined with 1st and 0th-order wfs
loc psp = -2.86633025E+03 Hartree= 6.39822440E+02 xc= -1.29343160E+02
note that "loc psp" includes a xc core correction that could be resolved
7,8,9: eventually, occupation + non-local contributions
edocc= 1.31009188E+02 enl0= 2.50930477E+02 enl1= -7.60810169E+02
1-9 gives the relaxation energy (to be shifted if some occ is /=2.0)
erelax= -1.81357021E+03
10,11,12 Non-relaxation contributions : frozen-wavefunctions and Ewald
fr.local= 1.12997278E+03 fr.nonlo= 3.83220247E+02 Ewald= 3.11475726E+02
13,14 Frozen wf xc core corrections (1) and (2)
frxc 1 = 0.00000000E+00 frxc 2 = 0.00000000E+00
Resulting in :
2DEtotal= 0.1109854717E+02 Ha. Also 2DEtotal= 0.302006827206E+03 eV
(2DErelax= -1.8135702103E+03 Ha. 2DEnonrelax= 1.8246687575E+03 Ha)
( non-var. 2DEtotal : 1.1098549462E+01 Ha)

--------
2DEtotal - non-var. 2DEtotal = -2.292e-06



STABLE q(1/2,0,0)
=================

Thirteen components of 2nd-order total energy (hartree) are
1,2,3: 0th-order hamiltonian combined with 1st-order wavefunctions
kin0= 3.02571359E+03 eigvalue= -1.19647927E+02 local= -1.67271751E+03
4,5,6: 1st-order hamiltonian combined with 1st and 0th-order wfs
loc psp = -2.48405700E+03 Hartree= 4.43026844E+02 xc= -1.04235126E+02
note that "loc psp" includes a xc core correction that could be resolved
7,8,9: eventually, occupation + non-local contributions
edocc= 8.62257748E+00 enl0= 5.72072210E+02 enl1= -1.82161232E+03
1-9 gives the relaxation energy (to be shifted if some occ is /=2.0)
erelax= -2.15283466E+03
10,11,12 Non-relaxation contributions : frozen-wavefunctions and Ewald
fr.local= 1.13170071E+03 fr.nonlo= 9.27281465E+02 Ewald= 1.00008560E+02
13,14 Frozen wf xc core corrections (1) and (2)
frxc 1 = 0.00000000E+00 frxc 2 = 0.00000000E+00
Resulting in :
2DEtotal= 0.6156069747E+01 Ha. Also 2DEtotal= 0.167515176944E+03 eV
(2DErelax= -2.1528346602E+03 Ha. 2DEnonrelax= 2.1589907299E+03 Ha)
( non-var. 2DEtotal : 6.1560705977E+00 Ha)

--------
2DEtotal - non-var. 2DEtotal = -8.507e-07



q=Gamma
=========

Seven components of 2nd-order total energy (hartree) are
1,2,3: 0th-order hamiltonian combined with 1st-order wavefunctions
kin0= 3.12538521E+02 eigvalue= -4.32992225E+01 local= -2.36545710E+02
4,5,6: 1st-order hamiltonian combined with 1st and 0th-order wfs
dotwf= -2.13296829E+02 Hartree= 7.96991653E+00 xc= -4.52776738E+00
7,8,9: eventually, occupation + non-local contributions
edocc= 4.37750396E+01 enl0= 2.67376382E+01 enl1= 0.00000000E+00
1-9 gives the relaxation energy (to be shifted if some occ is /=2.0)
erelax= -1.06648414E+02
No Ewald or frozen-wf contrib.: the relaxation energy is the total one
2DEtotal= -0.1066484144E+03 Ha. Also 2DEtotal= -0.290205094264E+04 eV
( non-var. 2DEtotal : -1.0664841454E+02 Ha)

--------
2DEtotal - non-var. 2DEtotal = 1.4e-07


Are these differences too large (I don't have enough experience so that I can compare with other cases), so that I should increase the convergence parameters?

I also notice that for the unstable q the eigenvalue of the 0th-order hamiltonian combined with 1st-order wavefunctions is positive
2.00776545E+02
while for the stable q it is negative
-1.19647927E+02
I'm not sure if this has significance, i.e. induces the imaginary frequencies?
On the other hand, the frequencies of the acoustic modes at Gamma are very close to 0 (about -0.9 cm^-1), which should show that the convergence parameters are ok.


I appreciate any remarks and help.

Yours

Igor Lukacevic

We have changed the permissions of the Forum so that you should be able to attach ABINIT output and input files in your posts. If you still have problem, let me know the extension of the file you tried to attach.

OK. Thanks.

I'm attaching now the rf phonon, mrgddb, anaddb inputs.

Cheers!

Igor Lukacevic

Dear prof. Verstraete,

mverstra wrote:As posted to the mailing lists:

1) try running the rf abinit calculations with rfasr=1
(http://www.abinit.org/documentation/hel ... html#rfasr)
For some extreme pseudopotentials (with large breaking of the asr)
this is crucial: correct asr both in abinit _and_ anaddb.

g'luck

Matthieu

I 've tried the rfasr =1 option in the rf phonon calc. and set asr = 2 (and for any case asr = 1 also). Unfortunately, the results are the same. The acoustic modes are wildly imaginary still. Degeneracy of modes is also bad throughout.

I am without clues on how to continue.

Yours,

Igor Lukacevic

Dear colleagues,

I've restarted the calculations on ZrO2 from scratch. Now I optimized the lattice first with ionmov/optcell (previously I used Etot vs acell method).
First I again calculated the phonons at Gamma:

2.465697E-01 2.707347E-01 3.813680E-01 2.730116E+02 2.730117E+02 5.925013E+02 5.925013E+02 5.925013E+02 6.726585E+02

while the above mentioned PRL has: 0 0 0 267 267 583 583 583 667. The difference is in the ecut value (my is 100 and their 30), otherwise I would say a good agreement (pseudos and k-point mesh are the same).
Z* and diel. constant are also in a good agreement:
My Their
Z* 5.74 -2.86
e 5.74
Their
Z* 5.75 -2.86
e 5.75
For variational and non-variational 2DEtotal I get
2DEtotal= -0.1057232230E+03 Ha
non-var. 2DEtotal : -1.0572322287E+02 Ha
I suppose there are approximately enough the same.


Then I tried the phonons at X point (q 0.0 0.5 0.5). Here the results differ heavily:
My
-3.672780E+02 -2.828415E+02 1.439990E+02 3.115149E+02 3.291936E+02 3.455205E+02 5.859111E+02 6.967360E+02 7.379747E+02

Their (from the lowest to the highest phonon freq. and some are degenerate - I'm not 100% sure which)
-193 136 323 365 567 700

Here for variational and non-variational 2DEtotal I get
2DEtotal= 0.6067278158E+01 Ha
non-var. 2DEtotal : 6.0672827780E+00 Ha
but I cannot say if these are close enough. The cell is relaxed with the same parameters as in the rf calculation.

How come that this happens? I appreciate any comments.
One more question: when abinit stacks the phonon frequencies, does it always stack them from the lowest to the highest, or it first gives 3 acoustic ones and then optic ones? I.e. if there is an unstable optic mode (it has imaginary freq.), will abinit put its freq. in front of acoustic ones (because it starts with a minus sign) or it will put it between the optical modes and leave the acoustic ones in front?

Yours

Igor Lukacevic


Here's my q=X input file (increasing both tolvrs' to 10^(-18) changes nothing)

# Crystalline ZrO2 - fluorite : computation of the response to atomic displacements, at q=X

ndtset 2

#Ground state calculation
kptopt1 1 # Automatic generation of k points, taking
# into account the symmetry
tolvrs1 1.0d-15 # SCF stopping criterion
iscf1 5 # Self-consistent calculation, using algorithm 5

#Response Function calculation : electric field perturbation and phonons
rfphon2 1 # Activate the calculation of the atomic dispacement perturbations
rfatpol2 1 3 # All the atoms will be displaced
rfdir2 1 1 1 # All directions are selected. However, symmetries will be used to decrease
# the number of perturbations, so only the x electric field is needed
# (and this explains why in the second dataset, rfdir was set to 1 0 0).

nqpt2 1
qpt2 0.0 0.5 0.5 # This is a calculation at X point

getwfk2 1 # Uses as input wfs the output wfs of the dataset 1

kptopt2 3 # Automatic generation of k points,
# no use of symmetries to decrease
# the size of the k point set.
tolvrs2 1.0d-8
iscf2 5 # Self-consistent calculation, using algorithm 5


#######################################################################
#Common input variables

#Definition of the unit cell
acell 3*9.47118842767
rprim 0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

#Definition of the atom types
ntypat 2
znucl 40 8 8

#Definition of the atoms
natom 3
typat 1 2 2
xred 0.00 0.00 0.00
0.25 0.25 0.25
0.75 0.75 0.75

#Gives the number of band, explicitely (do not take the default)
nband 14

#Exchange-correlation functional
ixc 3

#Definition of the planewave basis set
ecut 100
ecutsm 0.5

#Definition of the k-point grid
ngkpt 3*4
nshiftk 4
shiftk 0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5

#Definition of the SCF procedure
nstep 100
diemac 4.0

I happened to calculate the phonon dispersion of cubic ZrO2 before and the results are consistent with earlier literatures. The attached are the input files for both the geometrical optimization and phonon dispersion calculation.

(Hope this late reply could be useful for the future searcher.)

Sincerely,
Guangfu Luo

P.S. I employed Abinit-6.0.3 in this calculation.

Dear Guangfu,

just how much before was that? Actually, with the help from prof. Gonze, we figured that there was a bug in the abinit version I used, which did not exist until some version (can't really remember which one). But its fixed now, and everything works well.

Thanks anyway for the files. They might sometimes prove useful.

Igor L.