Elastic constant calculation

Question

I want to calculate elastic constants with VASP "by hand" (not by, e.g., DFPT) as I have a big structure (64 atoms). For this, I am making tests on $\ce{Si}$ at the moment.

So, I have done it for $C_{11}$ for which I have obtained $C_{11}=\pu{1543.75 kbar}$, as compared with $\pu{1543.48 kbar}$ with DFPT.

Now, I want to calculate $C_{44}$. I have used, e.g., the following strain matrix:

0.0  0.0  0.0
0.0  0.0  0.0015
0.0  0.0015 0.0

which leads to the cell vector matrix:

     5.4664034019782948    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.4664034019782948    0.0150000000000000
     0.0000000000000000    0.0150000000000000    5.4664034019782948

I have run the SCF calculations for which I have the following results (energies in eV):

-0.00150 -43.37379723
-0.00075 -43.37495130
 0.00000 -43.37531894
 0.00075 -43.37495130
 0.00150 -43.37379723

When I fit the curve with the polynomial $f(x)=a+bx^2$, to find $a$ and $b$, and calculate $C_{44}$ as $C_{44}=\frac{2b}{V}$ (as I do for $C_{11}$), I find:

$b=\pu{679.563 eV}$
$V=\pu{163.34 Å^3}$
$C_{44}=\pu{13326 kbar}$
which is obviously wrong as $C_{44}=\pu{1001.53 kbar}$ from DFPT, if one does not account for the ionic relaxation (=$\pu{-251.1746 kbar}$).

Does anyone know what's wrong in what I am doing?

Thank you,
Pascal

Answer 1

The deformation parameter ($q$) in your strain matrix ($q$=0.0015) seems to be 10 times smaller than the increment ($q$=0.015) in the lattice-vector matrix.

To correct for this, you can fit to $E=a+b(10 \times x)^2$ and to get $a = -43.3753$, $b = 6.7956$, and $C_{44}=$133.3108 kilobar.

To find C$_{44}$, the strain matrix used actually is

0.0   0.0   0.0  
0.0   0.0   q/2  
0.0   q/2   0.0  

To account for the denominator 2 in the strain vectors, You can fit to $E=a+b(2\times 10 \times x)^2$. Now, we get $b = 1.6989$, and $C_{44}=$33.3277 kilobar ($a$ is same as before).

If the fitting is done using a dimensionless deformation parameter (in units of the lattice constant), by fitting to $E=a+b(2\times 10 \times x/5.4664034019782948)^2$, we get $b=50.7660$ and $C_{44}=$995.8838 kilobar which is close to the DFPT value. This last step makes sure correctly that the dimension of $b$ is eV and the dimension of (2$b$/$V$) to be eV/$Å^2$.

Sometimes, the strain matrix used for $C_{44}$ includes a small diagonal element and has the form

q**2/(4-q**2) 0.0   0.0  
0.0           0.0   q/2  
0.0           q/2   0.0  

See for example this link. If you try with this strain matrix, you may get results in better agreement with DFPT.