# Elastic constant calculation

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

• The "Dear all" at the beginning of my question has vanished somehow... Apr 2, 2022 at 17:05
• I don't have an answer for you, but you might try using VaspKit (vaspkit.com), which will create input files to get elastic properties. Apr 2, 2022 at 17:11
• Thanks, that could be an option. Pascal Apr 2, 2022 at 18:39
• @Pascal StackExchange has certain filters in place to remove "conversational" elements from posts. Unlike a lot of forums, SE is tried to limit posts to just questions and answers. Further clarification and conversation is usually reserved for comments or chat.
– Tyberius
Apr 2, 2022 at 18:40
• May be you can also try using pymatgen, the analysis module (analysis.elasticity.elastic module). The good thing about this package is you can see from the source code how they calculate the elastic constant. Apr 3, 2022 at 2:47

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.

• Thank you Raghunathan for your help. That help me a lot. For sure, I used the first matrix you are talking about. I probably mixed up 0.0015 and 0.0150... Best, Apr 3, 2022 at 8:34