Problems of elastic constant calculations, IBRION=6
Posted: Fri Sep 15, 2017 11:17 am
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
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
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