LCALCEPS result different in 6.2 vs 5.4.4

[aohara]

Dear support / VASP community,

While working with a student who is learning to calculate dielectric tensors (static, not frequency dependent) to use in eventual defect calculations, I found an odd behavior in VASP version 6.2 vs 5.4.4 and was able to reproduce it on two different computing centers.

We are calculating the dielectric tensor of the primitive cell of diamond with the PBE potential on an 8x8x8 Gamma centered k-point grid. The plane-wave cutoff is 450 eV and ISMEAR = 0 with SIGMA = 0.005.

In VASP 5.4.4, if I use LCALCEPS I get 5.711869 and if I use LEPSILON I get 5.886591 which is pretty good agreement (especially given that these two methods can have different convergence criterion in terms of k-points, etc).

However, if I use VASP 6.2 then for LCALCEPS I get 3.29291 and with LEPSILON I get 5.886592.

So with the two different version, I get essentially no difference between the LEPSILON calculations, but I get a rather significant different with LCALCEPS.

Did something change in the implementation of LCALCEPS between VASP 5.4.4 and VASP 6.x? Is there something else that I'm missing that would cause this?

Thanks for any insight that can be provided.

Below, I provide the text of my INCAR, POSCAR, and KPOINTS files:

-------------------------------------------------------------------------------------------------
INCAR:
-------------------------------------------------------------------------------------------------

Code: Select all

SYSTEM = Diamond

! Initialization flags
ISTART = 0       ! 0=new, 1=continue
ICHARG = 2       ! 1=file, 2=atomic, 11=fixed from file (nscf)
INIWAV = 1       ! Needed for fresh starts, random initial wave function

! Electronic Loop Controls
ENCUT  = 450     ! Plane-wave cutoff
EDIFF  = 1E-5    ! stopping-criterion
ALGO   = Normal  ! Controls diagonalization algorithm
NELMIN = 2       ! Minimum number of steps
NELM   = 80      ! Maximum number of steps
LREAL = .FALSE.  ! Real space project yes/no
ISPIN = 1        ! 1 = non-spin-polarized, 2 = spin-polarized
LCALCEPS = .TRUE.

! Hybrid functional
!LHFCALC = .TRUE.
!LFOCKACE = .TRUE.
!HFSCREEN = 0.2
!AEXX = 0.25
!PRECFOCK = Normal

! Ionic Relaxation Controls
!NSW    = 0       ! Number of ionic steps
IBRION = 6      ! Controls ionic algorithm: -1 = off, 2 = conj. grad., 1 = quasinewton
!ISIF   = 2       ! Controls how relaxation occurs: 2 = ions only, 3 = ions plus lattice vectors
!EDIFFG = -1E-2   ! Ionic convergence criterion. Positive is energy (eV), negative is force (eV/A)

! Density of States and Integration Flags
ISMEAR = 0      ! Method for determining partial occupancies/integration
                 ! -5 is the tetrahedron method, 0 is Gaussian (we'll mostly use these
SIGMA = 0.005    ! Broadening ("smearing") width for Gaussian
!LORBIT = 11      ! DOSCAR and lm decomposed PROCAR file
! EMIN/EMAX defualt to sensible values (min/max value for DOS)
NEDOS  = 5000    ! Number of grid points in DOSCAR (density of states)
                 ! For certain plots, make these densor

! Parallelization flags: important on large calculations
NCORE = 1
KPAR  = 6

! Output related flags
LWAVE  = .TRUE.  ! write WAVECAR
LCHARG = .TRUE.  ! write CHGCAR, CHG
LVTOT  = .FALSE.  ! write LOCPOT (full potential)
LVHAR  = .FALSE.  ! write LOCPOT (Hartree only)

-------------------------------------------------------------------------------------------------
POSCAR:
-------------------------------------------------------------------------------------------------

Code: Select all

Diamond
3.574131
        0.0000000000         0.5000000000         0.5000000000
        0.5000000000         0.0000000000         0.5000000000
        0.5000000000         0.5000000000         0.0000000000
   C
    2
Direct
     0.000000000         0.000000000         0.000000000
     0.250000000         0.250000000         0.250000000

-------------------------------------------------------------------------------------------------
KPOINTS:
-------------------------------------------------------------------------------------------------

Code: Select all

Automatic mesh
0
Gamma
8 8 8
0 0 0

[henrique_miranda]

Thank you for the bug report!
I was able to reproduce this issue and I am currently investigating the cause for it.
I will update you as soon as I find something.

EDIT1: decreasing EDIFF (to 1e-8 for example) allows you to recover a result closer to the one from vasp.5.4.4.
I am trying to find why it used to converge to the right result even with such a high-tolerance (1e-5) as the one you are using.

[henrique_miranda]

Ok, I will try to resume here the differences:
From VASP 5.4.4 we have found an issue in the PROJXYZ routine (leading to numeric inaccuracies) which has been fixed in VASP 6.1.0.
In your case, the difference between the two results is very drastic because you are using a too low tolerance for energy convergence.
Here are some results that I obtained:

vasp.5.4.4

Code: Select all

EDIFF=1e-5  "epsilon_scf" = 5.71189506
EDIFF=1e-8  "epsilon_scf" = 5.80725066
EDIFF=1e-10 "epsilon_scf" = 5.80971898

vasp.6.2.0

Code: Select all

EDIFF=1e-5  "epsilon_scf" = 3.29291019
EDIFF=1e-8  "epsilon_scf" = 5.79111341
EDIFF=1e-10 "epsilon_scf" = 5.80729297

Also worth pointing out that In vasp.6.2.0 the convergence is reached with a much smaller number of iterations.

Hope this clarifies the issue

[aohara]

Henrique,

Thank you for such a detailed follow up analysis. I'm glad that it isn't a bug in the new VASP 6.2 and relates to additional criterion threshold, since I've been encouraging the other members of our group to start migrating to our VASP 6 license.

On a technical side, to improve my understanding, what is the reason for needing such a tight convergence threshold on the electronic degrees of freedom for doing the dielectric calculation? Especially since these convergence criterion (~10-8) are much tighter than are usually needed for routine calculations and relaxations and in particular for the response to E-fields vs perturbation theory (since the LEPSILON seems to be much less sensitive), what drives the sensitivity? Is it primarily due to the fact that the small fields affect the wave functions and their energies by such a small amount that you need to be below the threshold of noise?

Thanks again,
Andy

[henrique_miranda]

Dear Andy,

The bug that was fixed from 5.4.4 to 6.2.0 is related to the phase factors of the projectors and this lead to some numerical noise in the calculations. In this particular example that you have sent, the convergence tolerance criteria for the energy is quite large (especially for a system with 2 atoms). My impression is that because of the bug in vasp 5.4.4 it took more iterative steps to reach the convergence criteria, which somehow lead the result to be closer to the correct value.

Notice that the convergence criteria are on the energy and not on the wavefunctions or the final quantity you are interested in (dielectric function). The energy converges much quicker. Nonetheless, it is important that you always check that your convergence criterium is strict enough to not affect your final quantity.

In the case of a DFPT calculation (LEPSILON) you might notice the following lone in the output:

Code: Select all

Linear response reoptimize wavefunctions to high precision

The wavefunctions are re-optimized to high precision to ensure that final results are accurate even if a large EDIFF was specified in the input.