DFPT (IBRION = 8)

[hamza_bouafia_LGP]

Dear all,
I have two questions concerning DFPT calculations:
1) after launching the calculation, the first calculation converges according to EDIFF but after linear response progess dE (convergence) does not exceed 10E-4, how to increase this precision? or is this value (default) very sufficient?
2) the calculation is long, especially for large cells,if there is a calculation stop (due to a problem), how to restart the calculation so as not to redo all the calculation steps.
Thank you in advance
Best,
H. BOUAFIA

[merzuk.kaltak]

Dear Hamza,

For IBRION=8 you may set the POTIM tag to influence the magnitude of the displacement.
This might improve the precision slightly, but can also cause numerical instabilities.
If you mean something else, then please provide some input and output files and point us to the quantities you are looking at.

Concerning your second question:
Unfortunately, continuing an IBRION=8 calculation from an arbitrary displacement step is not possible yet. Here, I refer you to the phonopy package, which has a good vasp interface and produces the symmetrically distinct displacement for a given POSCAR. The corresponding vasp calculations for each displacement can then be done independently from each other.

[hamza_bouafia_LGP]

Thank you very much for your reply:
here is part of the output:

Code: Select all

 entering main loop
           N           E                                   dE                       d eps       ncg            rms             rms(c)
DAV:   1     0.281026179159E+05    0.28103E+05   -0.37631E+05 32784   0.295E+03
DAV:   2     0.837302049800E+04   -0.19730E+05   -0.17590E+05 26896   0.946E+02
DAV:   3     0.102418767560E+04   -0.73488E+04   -0.61959E+04 23808   0.383E+02
DAV:   4    -0.717248788510E+03   -0.17414E+04   -0.16681E+04 25872   0.185E+02
DAV:   5    -0.911021076819E+03   -0.19377E+03   -0.19144E+03 35408   0.686E+01    0.164E+02
DAV:   6    -0.948130097367E+03   -0.37109E+02   -0.63311E+02 36656   0.438E+01    0.728E+01
DAV:   7    -0.969693918365E+03   -0.21564E+02   -0.12688E+02 36512   0.419E+01    0.405E+01
DAV:   8    -0.980827370209E+03   -0.11133E+02   -0.37069E+01 43776   0.144E+01    0.194E+01
DAV:   9    -0.983231304405E+03   -0.24039E+01   -0.85767E+00 44016   0.522E+00    0.611E+00
DAV:  10    -0.983690197254E+03   -0.45889E+00   -0.18458E+00 44032   0.382E+00    0.242E+00
DAV:  11    -0.983775678052E+03   -0.85481E-01   -0.52368E-01 44000   0.152E+00    0.127E+00
DAV:  12    -0.983781442814E+03   -0.57648E-02   -0.92092E-02 43792   0.105E+00    0.321E-01
DAV:  13    -0.983784603301E+03   -0.31605E-02   -0.18386E-02 44016   0.374E-01    0.225E-01
DAV:  14    -0.983785571721E+03   -0.96842E-03   -0.36702E-03 43872   0.215E-01    0.903E-02
DAV:  15    -0.983785843149E+03   -0.27143E-03   -0.86400E-04 43312   0.840E-02    0.295E-02
DAV:  16    -0.983785884155E+03   -0.41005E-04   -0.17134E-04 43392   0.448E-02    0.155E-02
DAV:  17    -0.983785899272E+03   -0.15117E-04   -0.40384E-05 43904   0.194E-02    0.612E-03
DAV:  18    -0.983785903736E+03   -0.44642E-05   -0.10901E-05 43264   0.108E-02    0.307E-03
DAV:  19    -0.983785904228E+03   -0.49199E-06   -0.24478E-06 43648   0.477E-03    0.196E-03
DAV:  20    -0.983785904202E+03    0.25626E-07   -0.63153E-07 43920   0.256E-03    0.178E-03
DAV:  21    -0.983785904046E+03    0.15643E-06   -0.18088E-07 42816   0.129E-03    0.154E-03
DAV:  22    -0.983785904005E+03    0.41124E-07   -0.37380E-08 42608   0.657E-04    0.151E-03
DAV:  23    -0.983785903973E+03    0.31287E-07   -0.11403E-09 32896   0.337E-04    0.151E-03
DAV:  24    -0.983785903980E+03   -0.63592E-08    0.48165E-09 18832   0.182E-04
   1 F= -.98378590E+03 E0= -.98378590E+03  d E =-.533462E-23  mag=     0.0000
 Linear response reoptimize wavefunctions to high precision
DAV:   1    -0.983785903987E+03   -0.69267E-08    0.50556E-09 16384   0.140E-04
DAV:   2    -0.983785903987E+03    0.14552E-10    0.51074E-09 16224   0.135E-04
DAV:   3    -0.983785903986E+03    0.11642E-09    0.51554E-09 15648   0.106E-04
 Linear response DOF=          12
 Linear response progress:
  Degree of freedom:   1/ 12
            N       E                                         dE                 d eps             ncg        rms             rms(c)
RMM:   1    -0.732166311022E+00   -0.73217E+00    0.14873E+00 42737   0.173E+00
RMM:   2    -0.690532602964E+00    0.41634E-01   -0.10088E-01 62770   0.305E-01    0.425E-01
RMM:   3    -0.692969825360E+00   -0.24372E-02   -0.19137E-02 71897   0.109E-01    0.260E-01
RMM:   4    -0.697959721298E+00   -0.49899E-02   -0.12460E-02 75492   0.821E-02    0.419E-02
RMM:   5    -0.698789057266E+00   -0.82934E-03    0.25689E-03 80995   0.542E-02    0.191E-02
RMM:   6    -0.700919025866E+00   -0.21300E-02   -0.63299E-03 81167   0.504E-02    0.129E-02
RMM:   7    -0.700322859036E+00    0.59617E-03    0.12000E-02 81910   0.475E-02    0.652E-03
RMM:   8    -0.701943060480E+00   -0.16202E-02    0.52733E-03 82110   0.465E-02    0.320E-03
RMM:   9    -0.701084380203E+00    0.85868E-03    0.18297E-02 82244   0.458E-02    0.172E-03
RMM:  10    -0.701402452670E+00   -0.31807E-03    0.18299E-02 82431   0.452E-02    0.199E-03
 first order change of Fermi-energy  2.918196390653272E-002
 force on displaced ion       1 direction 1 : -59.704   0.000  -0.000
 

############################
so as you notice
for the first step (scf) we see that the convergence given in INCAR has been confirmed (1.E-8) but for the linear response, it does not exceed 1.E-5

Code: Select all

############################ INCAR
ISMEAR =  0           
SIGMA  =  0.05      
IBRION =  8       
EDIFF  =  1E-08       
PREC   =  Accurate    
ENCUT  =  520        
IALGO  =  38       
LREAL=.FALSE.      
ISPIN  =  2
MAGMOM = 16*0 96*0 4*4 8*-4 4*4
NELM   =  300        
NELMIN =  2            
 GGA  =  PS  
LDAU      = .TRUE.
LDAUTYPE  = 2
LDAUL     = -1 -1  3
LDAUU     = 0.00 0.00 5.67
LDAUJ     = 0.00 0.00 0.00
LMAXMIX   = 6
NSW = 1
EDIFFG =  1E-08

[merzuk.kaltak]

For a small test system I was not able to reproduce this behavior with vasp-6.4.2.
Please upload (KPOINTS, POTCAR, POSCAR, OUTCAR) and let us know which vasp version you use.

[hamza_bouafia_LGP]

Ia m sorry for the delay, because I redid the calculations.
Because of the calculation time, I reduced the cell to 1*1*1 instead of 2*1*3 and used lreal = auto instead of .FALSE.
As you see in the out file, some cycles of linear response dE do not exceed E-5.

[merzuk.kaltak]

Dear Hamza,

I have looked into the relevant file src/ilinear_response.F. The corresponding lines of code that are responsible for the break criterion in the ionic displacement steps are the following

Code: Select all

 844 !======================== end of loop ENDLSC ===========================
 845 ! end of the selfconsistent calculation loop
 846 ! and write commands
 847 !=======================================================================
 848        INFO%LABORT=.FALSE.
 849 
 850        IF(ABS(DESUM)*SCALE**2<INFO%EDIFF.AND.ABS(DE)*SCALE**2<INFO%EDIFF) INFO%LABORT=.TRUE.
 851        ! rms stable for four steps, ok stop, no way to improve
 852        IF(N>RESOLVE_DEG_NTIMES+4) THEN
 853           IF (ABS((RMS(N)-RMS(N-1))/RMS(N))<0.1_q .AND. ABS((RMS(N)-RMS(N-2))/RMS(N))<0.1_q .AND. ABS((RMS(N)-RMS(N-3))/RMS(N))<0.1_q) INFO%LABORT=.TRUE.
 854        ENDIF
 855        !-----do not stop before minimum number of iterations is reached
 856        IF (N < ABS(INFO%NELMIN)) INFO%LABORT=.FALSE.

As you can see vasp aborts the electronic minimization step in the ionic displacement based on two criteria:

• Sum of eigenvalues and change in energy times SCALE**2 < EDIFF. Note SCALE=0.02 at this stage.

• the residual vector rms does not change for more than 10% for 4 consecutive iterations

The first criterion can be influenced indirectly by setting EDIFF in the INCAR.
In contrast, the second criterion can not be changed from the INCAR.
The reason why these lines were introduced was the fact that we typically observed stagnation of rms after several steps. Running the minimizer for some additional iterations does not change the results significantly (if at all).
For instance, you can change the code to check also for stagnation in rms(c) (similar to src/elinear_response.F) like so

Code: Select all

        IF(N>RESOLVE_DEG_NTIMES+4) THEN 
          IF ( ABS((RMS(N)-RMS(N-1))/RMS(N))<0.1_q .AND. ABS((RMS(N)-RMS(N-2))/RMS(N))<0.1_q .AND.  ABS((RMS(N)-RMS(N-3))/RMS(N))<0.1_q )THEN
              !require also the same for charge residuum
              IF (ABS((RMST(N)-RMST(N-1))/RMST(N))<1E-1 .AND. ABS((RMST(N)-RMST(N-2))/RMST(N))<1E-1 .AND. ABS((RMST(N)-RMST(N-3))/RMST(N))<0.1_q) THEN 
                 INFO%LABORT=.TRUE.
              ENDIF
          ENDIF
       ENDIF

Note, this requires to change all occurrences of RMST in the subroutine LR_MAIN to RMST(N) and to define the variable as

Code: Select all

    REAL(q) :: RMST(INFO%NELM)            ! magnitude of residual vector in charge

If you run this code for your system you will notice no significant improvement, except that the code runs for more iterations. This demonstrates the choice for rms as the relevant convergence criterion.

N.B.: The abortion of the minimization based on stagnation of rms is implemented only in the linear response code. A conventional single point calculation does not respect this convergence criterion.

[hamza_bouafia_LGP]

Thank you very much for the help and for all this explanation.