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
```

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