Problems of elastic constant calculations, IBRION=6

[byin]

Dear admin,

Here are problems in the calculation of elastic constants with IBRION=6.


1.
In the documentary, there is "For NFREE=1, only a single displacement is applied (it is strongly recommended to avoid NFREE=1)". Is this suggestion applied only to the calculation of Hessian matrix, or both (Hessian and elastic constants)? And why shall we avoid NFREE=1?


2.
In the source code finite_diff.F, the elastic constants are calculated using standard differentiation formulas

Code: Select all

!
! contract calculated forces and stress using standard differentiation formulas
!
       ALLOCATE(SUM_FORCES(DOF,3,NIONS),SUM_SIF(DOF,3,3))
       SELECT CASE(NDISPL)
       CASE(1)
          DO N=1,DOF
             SUM_FORCES(N,:,:)=(DISPL_FORCES(1,N,:,:)-INITIAL_FORCE)/STEP
             SUM_SIF(N,:,:)   =-(DISPL_SIF(1,N,:,:)-INITIAL_SIF)/STEP
          END DO
       CASE(2)
          SUM_FORCES = (1._q/(2._q*STEP))*(DISPL_FORCES(1,:,:,:)-DISPL_FORCES(2,:,:,:))
          SUM_SIF    =-(1._q/(2._q*STEP))*(DISPL_SIF(1,:,:,:)   -DISPL_SIF(2,:,:,:))
       CASE(4)
          SUM_FORCES = (1._q/(12._q*STEP))* &
               (8._q*DISPL_FORCES(2,:,:,:)-8._q*DISPL_FORCES(3,:,:,:) &
               -DISPL_FORCES(1,:,:,:)+     DISPL_FORCES(4,:,:,:))
          SUM_SIF =-(1._q/(12._q*STEP))* &
               (8._q*DISPL_SIF(2,:,:,:)-8._q*DISPL_SIF(3,:,:,:) &
               -DISPL_SIF(1,:,:,:)+     DISPL_SIF(4,:,:,:))
       END SELECT

But in the documentary, it is said that the elastic constants are derived from the strain-stress relationship, ref.[4]. In ref. [4], the Cij tensor is calculated by the least-square fitting of 27 parameters. Please have a check.


3.
I tried NFREE=1, POTIM=0.001 in the calculation of fcc Pt, one-atom primitive cell. Here are the stress results, 1+6+3 in total:
grep "in kB" OUTCAR
in kB 0.00925 0.00925 0.00917 -0.00011 -0.00011 -0.00010
in kB -3.07460 -2.19066 -2.19074 -0.00010 -0.00010 -0.00010
in kB -2.18765 -3.08326 -2.19296 -0.00010 -0.00011 -0.00010
in kB -2.19291 -2.18815 -3.08274 -0.00010 -0.00011 -0.00011
in kB -0.00129 -0.00219 -0.00441 -0.67219 -0.00011 -0.00010
in kB -0.00067 0.00829 0.00737 0.00891 -0.68360 -0.00011
in kB 0.00750 -0.00061 0.00828 -0.00013 0.00809 -0.68277
in kB 0.00388 0.00391 0.00381 -0.00010 -0.00016 0.00714
in kB 0.00433 0.00432 0.00425 -0.00013 -0.00010 0.00015
in kB 0.00436 0.00435 0.00427 -0.00010 -0.00012 -0.00010

According to the code above, C11=-(-3.07460-0.00925)/0.001=3083.85.
But in the OUTCAR, the results is

ELASTIC MODULI (kBar)
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX 3086.9237 2202.0964 2202.0984 -0.0015 -0.0030 -0.0003
YY 2199.0845 3095.5888 2204.3237 -0.0011 0.0040 -0.0063
ZZ 2204.3508 2199.5853 3094.9934 -0.0052 -0.0001 0.0050
XY 10.5403 11.4316 13.5783 672.0889 0.0002 0.0001
YZ 9.9187 0.9527 1.7942 -9.0204 683.4927 0.0093
ZX 1.7496 9.8541 0.8925 0.0271 -8.1934 682.6604
--------------------------------------------------------------------------------

Why there is a difference here, 3083.85 vs 3086.9237? Although small, why this happens?


4.
Here, what we need is the change of stress tensor caused by the applied strain. As it is well known, due to the insufficient of cutoff energy, there may exist Pulay stress. We can express s_DFT+s_Pulay=s_true. In the calculation of the change of stress, will the Pulay stress part cancel out, leading to a reliable result? What is your opinion?


Looking forward to your reply.

Best,
Binglun

[support_vasp]

Hi,

We're sorry that we didn’t answer your question. This does not live up to the quality of support that we aim to provide. The team has since expanded. If we can still help with your problem, please ask again in a new post, linking to this one, and we will answer as quickly as possible.

Best wishes,

VASP