PAW+phonons

[ksmirnov]

As for your results, something is troubling me. You said that you analyzed the frequencies with anaddb, but anaddb is supposed to impose the acoustic sum rule, therefore you should have the first three frequencies set to zero. Plus, what you are quoting in your posts look like the abinit results printed in the output file, are you sure you're looking at the output of anaddb?
Boris

It troubles me too. I'm perfectly sure that I cite the anaddb results. The input file for the code is

# Phonons at Gamma
#**********************
rfmeth 1
enunit 0
eivec 1
asr 1 # Acoustic Sum Rule (ASR) for interatomic force constants
chneut 2 # Acoustic Sum Rule (ASR) for effective charges

# Wavevector list
**************************
nph1l 1 # number of wavevectors in the phonon list
qph1l 0.0 0.0 0.0 1.0

I used the SAME file for calculations with NCPP and everything went as it should.

Konstantin.

[Boris]

That's strange.

anaddb should give you a file with the extension ".freq". Do you have this file? If you do, could you please provide what is inside?

I had already observed that the phonon frequencies can be wrong depending on the cutoff energy, that's why I asked you to change your ecut. Now I just found out that it's actually related to the value of boxcut(ratio) printed in the output file. It seems that if boxcut(ratio)<2, the second derivative of the total energy is wrong. This in turn yields wrong frequencies.

Could you please grep "boxcut(ratio)" out of your output file (from abinit) and check if it's greater than 2?

Thanks

[ksmirnov]

anaddb should give you a file with the extension ".freq". Do you have this file? If you do, could you please provide what is inside?

anaddb produces a file with the extension PHFRQ. The file is attached (I've changed the extension to get the file attached).

Could you please grep "boxcut(ratio)" out of your output file (from abinit) and check if it's greater than 2?

The quote is for the output file of the RF calculation:

$ grep 'boxcut(ratio)' phonons.test3.out
ecut(hartree)= 25.000 => boxcut(ratio)= 2.05924
ecut(hartree)= 40.000 => boxcut(ratio)= 2.00783
ecut(hartree)= 25.000 => boxcut(ratio)= 2.05924
ecut(hartree)= 40.000 => boxcut(ratio)= 2.00783

Obviously, the ratio > 2.

[Boris]

Hi Konstantin

Maybe your issue has something to do with the fact that the acoustic sum rule is still not satisfied, even with anaddb. I have run a series of calculations on different materials and had the same issue, i.e. when the first three acoustic phonon frequencies are exactly zero, I obtain a good phonon spectrum, but when they're not, the frequencies are wrong.

I don't know why anaddb is not imposing the acoustic sum rule. I don't think it's related to anaddb because when I set rfasr=1 in abinit, the acoustic sum rule is not satisfied either.

I'm going to investigate this and I'll let you know.

Also, always check that your boxcut(ratio) is greater than 2. If it's not, change your cutoff until it is. This has been reported and should be corrected in a future release.

Boris

[ksmirnov]

Dear Boris,

Thank you for the hints. The weird thing is that the ASR is obviously satisfied with NCPPs, while it does not with PAW. It seems to be strange because anaddb operates with a precomputed Hessian. I'll appreciate any feedback concerning the issue. Have a good holidays.
Best regards,
Konstantin.

[Boris]

Yes it is strange

What's puzzling me even more is that I calculated phonon frequencies at Gamma in different materials and the acoustic sum rule is satisfied when there is only one atom per unit cell. When I build a supercell with two atoms per unit cell, the asr is not satisfied anymore and the frequencies go nuts.

I first thought that the asr was not satisfied when natom>1, but now I just tried in UO2, with 3 atoms per unit cell, and the asr is satisfied with reasonable phonon frequencies at Gamma (compared to small displacements).

And as you said it seems to occur only with PAW.

Could you please post your input file for the NCPP calculation?

Thanks

[ksmirnov]

The NCPP input file for the phonons at Gamma attached; PSPs were taken from FHI database. At the end of the file you find extracts from abinit and anaddb outputs with the values of the phonon frequencies. As you see the ASR is not satisfied in the abinit output, while when anaddb has done its job the three low-frequencies are almost zero (although I used a rather soft convergence criterion for the geometry tolmxf=2.0d-5 in that calculation).

I also noticed that the strange behaviour of PAW depends on the number of atoms in unit cell. Thus, everything went well for Li2O, three large negative frequencies were obtained for a-quartz, four/five for a-cristobalite, and a lot for a-V2O5. Again, no problem was encountered for either system in NCPP calculations.

Best regards,
Konstantin.

[ksmirnov]

Hello everybody.

The thread was inactive for quite a long time. Any movements?

Very recently, I've tested the DFPT phonon calculations for alpha-quartz at the LDA/PAW level with v7.8.2 and have obtained the same weird result: large negative frequencies. Obviously, the issue, which is the root of the thread, was
not fixed (the release announcement of the lastest 7.10.2 version does not mention any changes in this direction).

On the other hand, LDA/PAW DFPT calculations of the electronic dielectric tensor seem to work. At least, the results for alpha-quartz are as follows:

LDA/NCPP
2.492037 0.000000 -0.000000
0.000000 2.492037 0.000000
-0.000000 -0.000000 2.522437

LDA/PAW
2.524434 0.000001 -0.000000
0.000001 2.524436 -0.000000
0.000000 -0.000000 2.557826

Hope this information will help tracing the source of the problem.

Best regards,

Konstantin.

[raul_l]

Same thing here. I tried this type of calculation with CdWO4 (12 atoms per unit cell) and got

Code: Select all

 Phonon frequencies in cm-1    :
- -3.771587E+02 -3.144719E+02 -2.516969E+02 -2.178755E+02 -9.971816E+01
- -9.688309E+01 -9.625479E+01  4.964200E+01  7.695275E+01  1.006255E+02
-  1.311926E+02  1.770384E+02  2.100342E+02  2.420567E+02  2.577156E+02
-  2.687643E+02  2.992373E+02  3.331531E+02  3.486120E+02  3.609707E+02
-  3.812165E+02  4.065609E+02  4.213171E+02  4.506210E+02  4.722604E+02
-  4.774440E+02  5.419136E+02  5.523574E+02  6.136546E+02  6.715841E+02
-  7.080649E+02  7.121983E+02  7.772886E+02  8.311704E+02  9.049231E+02
-  9.649267E+02

for qpt 0.5 0.5 0.0 (similar with other qpoints). I'm using v7.10.2.

[sheng]

Just to report that 7.10.4 has the same problem. Anaddb only manages to bring one acoustic mode to zero whereas the other two modes have frequencies of about negative one hundred. There is no such error for ncpp.

[Boris]

Hi everybody

We've been working a lot on DFPT+PAW and we have corrected a number of things in the code, in particular issues related to symmetries. The 7.10.4 version should work well on many systems, so for those of you who have not tried it yet, I suggest you do so and let us know about the results.

Some additional information: looking at the abinit frequencies, it is normal that the acoustic sum rule is more satisfied in NCPP than in PAW. This is because of the PAW formalism that calculates some integral quantities on the fine PAW grid. However, after going through anaddb, the acoustic sum rule should be perfectly satisfied. If it's not, it means that the calculation ran into an issue. In a "good" calculation, the first three abinit frequencies (not the anaddb frequencies) should be close to zero, but not exactly. From my experience, a "good" calculation has the first three frequencies approximately equal to 0.3-0.7 THz. It can be positive or negative. If you have larger values, then something probably went wrong.

In addition, we have noticed that in supercells with 50 or more atoms, the acoustic sum rule is not perfectly satisfied even after going through anaddb. It is currently under investigation.

Next production version of abinit will be further improved related to the DFPT+PAW part. Stay tuned!

[roginovicci]

Dear Dr. Boris Dorado,
Thanks a lot for working on DFPT+PAW problem! I may assume the problem partly (or fully) resolved. Talking about subject structure using LDA PAW datasets taken from JTH project and using acoustic sum rule the negative frequencies vanishes in 7.10.5 version. But the value of acoustic frequencies in ABINIT output still rather high. Its of order 40-20cm-1. Actually one is -40 while other two approximately 20. The stress tensor seems to be homogeneous:

Code: Select all

 Cartesian components of stress tensor (hartree/bohr^3)
  sigma(1 1)=  1.06276443E-06  sigma(3 2)=  0.00000000E+00
  sigma(2 2)=  1.06276443E-06  sigma(3 1)=  0.00000000E+00
  sigma(3 3)=  9.90957041E-07  sigma(2 1)=  0.00000000E+00

The plain wave basis is definitely full:

Code: Select all

ecut 40
pawecutdg 80

But I should say the system is extremely compact, as the closes distance between Si and O atoms is of order 3Bohr while the sum of paw radius is bigger:

Code: Select all

 chkpawovlp : WARNING -
  PAW SPHERES ARE OVERLAPPING !
   Distance between atoms   1 and   8 is  :   3.03025
   PAW radius of the sphere around atom   1 is:   1.90945
   PAW radius of the sphere around atom   8 is:   1.41465
   This leads to a (voluminal) overlap ratio of  1.80 %

Reducing paw radius in PAW datasets often leads to increasing of the plane-wave basis size, and I not sure it worth to do.

So, my question is it normal to have -40cm-1 of acoustic mode frequency? Or maybe PAW dataset should be enhanced some way or whatever?
The other question is when DFPT calculations with GGA functionals is going to be implemented witin PAW approach?

Bien sincèrement, Eugene.

[Boris]

Hello Eugene,

40 cm^-1 makes approximately 1.2 THz. Usually my frequencies at Gamma in the abinit output are around 0.5 or 0.6 so it's the same order of magnitude. But you should have the first three modes nearly degenerate in frequencies so my guess would be that something went probably wrong.

I have two more advices for DFPT users :
1. At the end of the DFPT run, for each qpoint, you have the fourteen components of the derivative of the total energy (just grep 'Fourteen'). At the end of this section, you have the value of the 2nd derivative of the total energy "2DETotal = XXXX Ha" then, two lines below, you have the same 2nd derivative of the total energy calculated in a non-variational way "non-var 2DETotal = XXXX Ha". These two values must be close!! Usually in my calculations, the difference between the two values is lower than 10^-5 Ha. Sometimes, in very complex systems (actinides, no symmetries, defects, etc.), I have a difference of ~10^-2 Ha. If the difference is larger than that, then your calculation is wrong. If it happens, you may first want to check the convergence of your ground state wave functions (tolwfr at least 10^-20). Usually it solves the problem.

2. The other important point is that the kpoint grid needs to satisfy the space group symmetry. However, when it does not, abinit doesn't say anything, the calculation runs smoothly but the results are wrong. The indication for this "bug" is that the acoustic sum rule is absolutely not satisfied (which I believe is the main subject of this topic ). In order to solve this issue there is a little bypass that I found and that you should use while we are correcting the issue: just use rfstrs=3 in your calculation at the Gamma q point. If the kpoint grid does not fully respect the space group symmetry, then in this case abinit will stop and will display an error message saying that the kpoint grid must be changed. Just change your kpoint grid so that it satisfies this rule. And always use rfstrs=3 in order to be sure that you have the correct k point grid.

If your calculations follow these two points, you should have better phonon spectra

As for the GGA implementation, the development will start very soon!

[roginovicci]

Dear Dr. Boris Dorado!

Thank you very much for the detailed answer! But concerning your suggestion I've found all my input data is correct. Namely:
1. I've checked in output file each 'Fourteen' section and found the values is pretty close:

Code: Select all

displacement of atom   1   along direction   1
 2DEtotal=    0.4149788857E+02 Ha. Also 2DEtotal=    0.112921497530E+04 eV
    (2DErelax=   -9.2443763861E+01 Ha. 2DEnonrelax=    1.3394165243E+02 Ha)
    (  non-var. 2DEtotal :    4.1497888575E+01 Ha)

displacement of atom   1   along direction   3
 2DEtotal=    0.8913498291E+02 Ha. Also 2DEtotal=    0.242548623530E+04 eV
    (2DErelax=   -1.9045119793E+02 Ha. 2DEnonrelax=    2.7958618084E+02 Ha)
    (  non-var. 2DEtotal :    8.9134982909E+01 Ha)

displacement of atom   5   along direction   1
 2DEtotal=    0.6236169058E+01 Ha. Also 2DEtotal=    0.169694790043E+03 eV
    (2DErelax=   -1.7093290100E+03 Ha. 2DEnonrelax=    1.7155651791E+03 Ha)
    (  non-var. 2DEtotal :    6.2361690572E+00 Ha)

 displacement of atom   5   along direction   2
2DEtotal=    0.3851608987E+02 Ha. Also 2DEtotal=    0.104807610628E+04 eV
    (2DErelax=   -1.7302963337E+03 Ha. 2DEnonrelax=    1.7688124236E+03 Ha)
    (  non-var. 2DEtotal :    3.8516089872E+01 Ha)
 displacement of atom   5   along direction   3
2DEtotal=    0.4594562599E+02 Ha. Also 2DEtotal=    0.125024406541E+04 eV
    (2DErelax=   -3.3547325007E+03 Ha. 2DEnonrelax=    3.4006781266E+03 Ha)
    (  non-var. 2DEtotal :    4.5945625986E+01 Ha)

But that section also point on overlaping
16 Contribution from 1st-order change of wavefunctions overlap

2. Using variable rfstrs=3 is forbidden for PAW, but I've checked the input with NCPP, all went smooth, more over the k-point grid I use is the one achieved with abinit input file containing
kptopt 1
prtkpt 1
kptrlen 60

both variants I use persists in output:

Code: Select all

    
   9     6   0   0      5.0000E-01      5.2200E+01      12     1
          0   6   0      5.0000E-01
          0   0   4      5.0000E-01
 
  26     6   0   0      5.0000E-01      5.2200E+01      12     2
          0   6   0      5.0000E-01
          0   0   4      5.0000E-01
 
   27     6   0   0      5.0000E-01      5.6229E+01      18     2
          0   6   0      5.0000E-01
          0   0   6      5.0000E-01

By the way what the difference between set #9 and #26 what 'iset' variable is?

We still have accoustic at gamma is big:

Code: Select all

 Phonon frequencies in cm-1    :
- -4.695097E+01  2.110835E+01  2.110835E+01  7.149604E+01  1.067714E+02
-  1.405374E+02  1.405374E+02  2.220122E+02  2.665549E+02  2.665549E+02 ...

But at least as for version 10.5 ASR works quite well and negative values vanishes.

[roginovicci]

Some additional info.

I'm using JTH PAW datasets, ABINIT complies on high values in Psp strength Dij matrix

Code: Select all

 pawio_print_ij: WARNING - 
  The matrix seems to have high value(s) !
  (  3 components have a value greater than   50.0). 
  It can cause instabilities during SCF convergence.
  Action: you should check your atomic dataset (psp file) 
          and look for "high" projector functions...  

Truly saying have no idea how to fix it.

The other thing is concerning to other PAW dataset taking from gpaw project. Using it provides acoustic phonon frequency to be same order:

Code: Select all

 Phonon frequencies in cm-1    :
-  2.426490E+01  2.463850E+01  2.463850E+01  5.949732E+01  1.212202E+02
-  1.547501E+02  1.547501E+02  2.303079E+02  2.600356E+02  2.600356E+02

Should we conclude here, that gpaw "potentials" works better in selected structure?

[Boris]

Hello Eugene

Your result with the GPAW PAW dataset is very interesting. The frequencies at Gamma are all degenerated and have values of around 0,7 THz, which seems way better.

We will look into this, thanks!

Boris

[roginovicci]

The paw dataset I've used in calculations above was generated with gpaw project the following way:
For oxigen:
gpaw-setup -l -f LDA -r 1.1,1.05,1.05 -c "[He]" -w -p --compensation-charge-radius=1. -e "0.78;0.97;0.0" O
for silicium:
gpaw-setup -l -f LDA -r 1.8,1.9,1.8 -c "[Ne]" -w -p --compensation-charge-radius=1.6 -e "0.55;0.99;0.0" Si

These "potentials" gives no overlapping. But require rather high ecut = 50Ha!

[Boris]

The high cutoff is probably because of the O dataset. The one I generated with atompaw, as well as the JTH dataset, both require a 50 Ha cutoff.

[roginovicci]

That's sound strange as my result (for JTH datasets) gives energy convergence of 0.1mHa at ecut=40Ha (k-grid is ngkpt )

Code: Select all

           
             ecut1     2.50000000E+01 Hartree
             ecut2     3.00000000E+01 Hartree
             ecut3     3.50000000E+01 Hartree
             ecut4     4.00000000E+01 Hartree
             ecut5     4.50000000E+01 Hartree
           etotal1    -1.4590089049E+02
           etotal2    -1.4590617116E+02
           etotal3    -1.4590888493E+02
           etotal4    -1.4590931988E+02
           etotal5    -1.4590935345E+02

In case of gpaw "potential" the paw radius is approximately 1a.u. that's for sure the reason why ecut is high. In case of JTH "potential" the one for Oxigen is far softer, Rpaw=1.4, so high ecut is mostly because system is very compact the minimal distance is about 3a.u, that's very small distance between O and Si atoms.

Anyway, I will try to perform geometry optimization and response calculations with higher ecut and let you know.

[roginovicci]

OK, I've performed calculations with higher ecut=50Ha. Now all acoustic modes is negative and of an order of tenth cm-1 (-45cm-1 max)

Code: Select all

 Phonon frequencies in cm-1    :                                                                                                                                          
- -4.589082E+01 -1.724093E+01 -1.724093E+01  5.334798E+01  1.029513E+02  

After applying ASR the negative frequencies wanishes and the difference with lower ecut in range of 1cm-1. Thus thread should be marked as SOLVED?

Dear Dr. Boris Dorado, may I ask a question? You was saying you able to generate atompaw "potentials". I've found JTH paw datasets produces warning wile scf cycles:

Code: Select all

 pawio_print_ij: WARNING -                                                                                                                                                
  The matrix seems to have high value(s) !                                                                                                                                
  (  3 components have a value greater than   50.0).                                                                                                                      
  It can cause instabilities during SCF convergence.                                                                                                                      
  Action: you should check your atomic dataset (psp file)                                                                                                                 
          and look for "high" projector functions...  

I suggest that's because scheme used to generate PS partial-waves for Si atom is "polynom2". I've change "custom" option to "VANDERBILT" and generated new "potential" and WARNING I've mention has left, and convergence seems to go faster (more testing needed). But I did not found how to set up E_loc value in section governing the scheme used to get VPS(r) (local) pseudopotential from all-electron effective potential Veff(r) using "ultrasoft" pseudization scheme (without norm conservation constraint). I've found the E_loc value acts only on "Vloc energy" in output file. The high values also produces ghost states some time. Should I use E_loc as high to reduce "Vloc energy" or is there some other rule to choose E_loc.

I do understand the questions is out of topic, but I'll be appreciate a lot for any hint.

Eugene.