How to calculate elastic constants with SIESTA?

The answer to this question can of course be found in SIESTA mail list archives. Let me repeat my own answer which has been posted there recently:

“SIESTA calculates stress tensor and you know stress/strain relationship from the textbooks. You have to relax your system properly to have your stress tensor be close to zero. Then manually deform lattice vectors in selected directions and obtain the value of stress. For example, if you want to calculate C11, you have to deform lattice vector a (I suppose it is (a,0,0)) to a(1+e), where e is your value of strain. Then you make the calculation with relaxing atomic positions but NOT lattice vectors (VariableCell false) and obtain the value of stress tensor Sigma_ij and find C11:

C11 = Sigma11/e”

This gives you the main idea. Here I want to describe this procedure more properly and focus on how to overcome some common problems.

Problem 1. Which value of strain to choose?
It depends on the material. You have to be sure that you're in the linear regime and Hook's law works, so the smaller the strain the better chance that you'll not be mistaken. BUT! Small values of strain mean less accuracy (I'll explain why), so the value of strain shouldn't be too small. I recommend the value of 0.5% or 1%. You may also (if you want) perform convergence test of elastic constant vs value of strain. To my experience, values in the range of 0.3%--1% give good results.

Problem 2. How to choose the components of strain tensor?
That's simple. If you want to calculate C11, the strain tensor eij should be

e 0 0
0 0 0
0 0 0

Where e is the value of strain (0.01 typically) Such tensor can also be used for computing C12 and C13. For C22 it should be

0 0 0
0 e 0
0 0 0

And the same for C33. As for shear elastic constants (C44, C55 and C66) it has nondiagonal elements. For example for C66 it will be:

0 e 0
0 0 0
0 0 0

Then to calculate elastic constant you have to deform your lattice tensor (containing three lattice vectors) like this:

aij = a0ij(1+eij)

And then run the calculation with relaxation of atomic positions (VariableCell false). SIESTA gives the values of stress tensor. For example, you want to calculate C66. Only e12 is nonzero, so the expression for Sigma12 will be:

Sigma12 = C1212e12 = C66e
then C66 = Sigma12/e

This is in the ideal case assuming that the initial lattice was fully relaxed and the stress tensor was completely zero for relaxed cell. But it is not the case, yes?

Problem 2. I've relaxed the cell but the stress tensor is not zero. Does it influence the results?

Yes, but not much. You may improve the results if you provide “zero stress correction”. For example, for C66 elastic constant:

C66 = (Sigma12-Sigma012)/e,

Where Sigma0ij is the stress tensor of the relaxed cell. It should be zero but it doesn't. This correction will work only if Sigma012 << Sigma12! If not, I advice to increase e, so that it will work. Note also, if you use some geometry constraints, you will not be able to compute elastic constants or the error will be big!

Problem 3. When I perform the deformation of +1% and -1%, I obtain different values of stress tensor.
Yes, this happen. See Problem 2 for solution - introducing zero stress correction will improve the results (maybe:). I advice that you always do two calculations for e>0 and e<0. In ideal case the results will be the same, but in real world you'll get the difference. Taking the average value of two elastic constants and introducing zero stress correction will improve the accuracy.

The last but not least: even if you have a cubic system with C11, C12 and C44 independent elastic constants, compute all of them! Calculate C11, C22 etc. to reduce the error - if you have time to do so.

Hope this short information helps. Feel free to correct me if I'm wrong or ask questions in the comments (you don't need to register!)