Calculation error of elastic tensors

I downloaded VASP files of some materials, and used them to calculate the elastic tensors by VASP.

But the tensors I got are different with which I got from Materials Project. Could you please tell me why?

the INCAR I use:

SYSTEM = diamond
PREC = Accurate
ENCUT = 600.0
EDIFF = 1e-6
ISMEAR = 0
SIGMA = 0.05
POTIM = 0.100
LCHARG = FALSE
LWAVE = FALSE
IBRION=6
ISIF=3
NFREE=4
EDIFFG = 1e-5
EDIFFG = -0.001

In general, elastic tensor calculations are very sensitive to changes in parameters. Your smearing type (ISMEAR=0) is different than the one we use, which might explain the difference. Also, using the primitive vs conventional cell can make a difference. How large is the change between your and our result?

Just for reference, all of our elastic tensor methods are contained in the atomate package.

Yes, just to clarify further, we would need a specific example of the material, your specific input crystal (poscar) and your output tensor to be able to comment further. There are many many variables here, include a difference in setting, differences in deformations applied, limits of numerical accuracy, potentially difference of pseudopotentials etc.

I’d be surprised if tensors agreed exactly so some small variation is to be expected. And as Joey says, these are very sensitive to convergence also.

The material ID is mp-1017566, whose stiffness tensor is (GPa)
[97 -20 -20 0 0 0
-20 97 -20 0 0 0
-20 -20 97 0 0 0
0 0 0 102 0 0
0 0 0 0 102 0
0 0 0 0 0 102]

The POSCAR file was dowloaded at https://materialsproject.org/materials/mp-1017566/

which is
Ge1 Pb1 O3
1.0
3.888437 0.000000 0.000000
0.000000 3.888437 0.000000
0.000000 0.000000 3.888437
Ge Pb O
1 1 3
direct
0.500000 0.500000 0.500000 Ge
0.000000 0.000000 0.000000 Pb
0.000000 0.500000 0.500000 O
0.500000 0.000000 0.500000 O
0.500000 0.500000 0.000000 O

And I have also dowloaded the INCAR file from https://materialsproject.org/materials/mp-1017566/ , in which I can see your smearing type (ISMEAR=-5)

Here are two results (GPa) I calculated by VASP:
1).
[243.8 109.7 109.7 0 0 0
109.7 243.8 109.7 0 0 0
109.7 109.7 243.8 0 0 0
0 0 0 137.9 0 0
0 0 0 0 137.9 0
0 0 0 0 0 137.9]
(INCAR file can be seen in the topic title)
2).
[249.8 114.4 114.4 0 0 0
114.4 249.8 114.4 0 0 0
114.4 114.4 249.8 0 0 0
0 0 0 136.7 0 0
0 0 0 0 136.7 0
0 0 0 0 0 136.7]
(INCAR file was only changed from ISMEAR=0 to ISMEAR=-5)

Thanks for your reply! I have already replied the details.

Thanks for your reply! I have already replied the POSCAR and INCAR files.

There could be several causes of this difference since elastic calculations are very sensitive. The smearing, the ENCUT, the way VASP does finite-difference using the POTIM, the pseudopotentials you use all affect this calculation. On top of this, I don’t believe VASP will relax your structure and run an elastic calculation in one run. These need to be two separate calculations.

I would recommend starting by relaxing the structure. Note, VASP can only relax on energy or forces, not both as you appear to be trying to do. Then let VASP perform an elastic calculation. That should get you much closer to our value, although it’s likely to differ still since we don’t use VASP’s built-in algorithm.

It’s also possible the MP tensor is incorrect here … we should check.

Ok, @shyamd pointed out this publication https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3344236/figure/f4-ijms-13-04632/ which seems to suggest the experimental bulk modulus (~40 GPA for a glass) might be closer to the MP value that the value from the tensor you pasted above, which gives a VRH bulk modulus of ~ 100 GPa, though this is probably not sufficient to say one way or the other for sure.

I have run several times after changing the smearing, the ENCUT and POTIM according to the INCAR I downloaded at MP. But the results hardly changed.

Meanwhile, I calculated the elastic tensor for MP-48 and got a result similar to that MP gives. So I think maybe the INCAR file I used has no problem.

Also, I tried to strech MP-1017566 by VASP and found the material would shrink in the transverse direction, which means it has a positive Poisson’s ratio at this direction. But according to the stiffness tensor from MP:
[97 -20 -20 0 0 0
-20 97 -20 0 0 0
-20 -20 97 0 0 0
0 0 0 102 0 0
0 0 0 0 102 0
0 0 0 0 0 102],
this material is not supposed to have a positive Poisson’s ratio at any direction.

We don’t use VASP’s elastic constant routine. By using VASP’s built-in elastic routine, you’re doing something different than what we perform and asking why you’re getting a different answer. As Joey mentioned, our workflow is in atomate, and you’re welcome to run it directly and if you’re still getting different results we could help. You can also try manually running our procedure as follows:

1.) Relax the structure to an energy convergence of ~ 1E-6 using a K-point density of 7,000 per reciprocal Angstrom^3 and turn off symmetry so VASP can’t mess with the K-mesh.
This has to be done in a separate calculation. Are you doing this?
2.) Generate 4 strained structures corresponding to -1%, -0.5%, 0.5%, 1% for each of the 6 unique deformation directions. I believe VASP only uses 1 per direction.
3.) Run static calculations for each of these new deformed structures using an ENCUT of 700, k-point density of 7,000 per reciprocal angstrom^3, turning off symmetry, using ADDGRID=True to get an even more precise stress tensor. In addition, the K-point mesh for these calculations is for the symmetry of the deformed structure. VASP uses the k-point mesh of the original structure, as far as I can tell.
4.) We then fit these back to the elastic tensor using an over-determined matrix inverse, which lets us separate out any potential non-linearities.

Our procedure is excessive when a material is hard and insulating, such as mp-48 where you are likely to get a good elastic tensor with a simple procedure as what is built into VASP. As it gets softer and more semiconducting, the effect of all these differences gets magnified. Especially for something that is potentially soft. For the case of PbGeO3 in particular, I would be worried about the strains that VASP uses for generating the deformations, whether you have a high enough ENCUT + ADDGRID to get accurate stress tensors and whether you’re relaxing sufficiently to compute an elastic tensor in the first place over the other INCAR parameters.

I’m sorry for replying after so long.

I have calculated the Poisson’s ratios and elastic tensors of several materials reapetedly with the INCAR and KPOINTS changed according to your suggests. But the reasults are still different with that from MP.

Another problem is that, we find the maximum Poisson’s ratio of MP-1017566 is negative by rotating the stiffness tensor you provide. But after stretching the material by VASP, we got a positive Poisson’s ratio, which is impossible.

So, could you please give us your input files or stretch this material at one direction and calculate one of its Poisson’s ratio?

Thank you!

Hi @MilMou,

I agree it’s important to figure out what’s going on here. If it helps, here are the stress/strains we used to fit this specific elastic tensor:

cauchy_stresses = [ 
    [ [ 0.485835569, 0.0, -0.0 ], [ -0.0, -0.100418809, -0.0 ], [ -0.0, -0.0, -0.100418809 ] ], 
    [ [ -0.177069249, -0.0, -0.0 ], [ -0.0, -0.166530955, 1.024832653 ], [ -0.0, 1.024832653, -0.18702881500000001  ] ], 
    [ [  -0.10041880899999997, 0.0, 0.0 ], [ 0.0, 0.4858355689999999, 0.0 ], [ 0.0, 0.0, -0.10041880899999997  ] ], 
    [ [  -0.10041880899999997, 0.0, 0.0 ], [ 0.0, -0.10041880899999997, 0.0 ], [ 0.0, 0.0, 0.4858355689999999  ] ], 
    [ [  -0.16653095499999998, 0.0, -1.0248326529999998 ], [ 0.0, -0.17706924899999998, 0.0 ], [ -1.0248326529999998, 0.0, -0.187028815  ] ], 
    [ [  -0.17706924899999998, 0.0, 0.0 ], [ 0.0, -0.16653095499999998, -1.0248326529999998 ], [ 0.0, -1.0248326529999998, -0.187028815  ] ], 
    [ [  -0.16653095499999998, 0.0, 1.0248326529999998 ], [ 0.0, -0.17706924899999998, 0.0 ], [ 1.0248326529999998, 0.0, -0.187028815  ] ], 
    [ [  -0.187028815, -1.0248326529999998, 0.0 ], [ -1.0248326529999998, -0.16653095499999998, 0.0 ], [ 0.0, 0.0, -0.17706924899999998  ] ], 
    [ [  -0.187028815, 1.0248326529999998, 0.0 ], [ 1.0248326529999998, -0.16653095499999998, 0.0 ], [ 0.0, 0.0, -0.17706924899999998 ] ]
]

strains = [
    [ [ 0.004999903894075897, 0.0, 0.0 ], [ 0.0, -1.1102230246251565E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251565E-16  ] ], 
    [ [  -1.1102230246251565E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251565E-16, 0.005000051366872475 ], [ 0.0, 0.005000051366872475, -8.041918347911903E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, 0.004999903894075896, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, 0.004999903894075896  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, -0.005000051366872474 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ -0.005000051366872474, 0.0, -8.041918347911901E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251563E-16, -0.005000051366872474 ], [ 0.0, -0.005000051366872474, -8.041918347911901E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.005000051366872474 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ 0.005000051366872474, 0.0, -8.041918347911901E-8  ] ], 
    [ [  -8.041918347911901E-8, -0.005000051366872474, 0.0 ], [ -0.005000051366872474, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16  ] ], 
    [ [  -8.041918347911901E-8, 0.005000051366872474, 0.0 ], [ 0.005000051366872474, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16 ] ]
]

And this is for structure:

Full Formula (Ge1 Pb1 O3)
Reduced Formula: GePbO3
abc   :   3.893560   3.893560   3.893560
angles:  90.000000  90.000000  90.000000
Sites (5)
  #  SP       a     b    c    magmom
---  ----  ----  ----  ---  --------
  0  Ge     0.5   0.5  0.5        -0
  1  Pb     0     0    0          -0
  2  O     -0     0.5  0.5         0
  3  O      0.5   0.5  0          -0
  4  O      0.5  -0    0.5        -0

Best,

Matt

Thanks for your reply!

We know that stress tensor = stiffness tensor * strain tensor
But when I used the tensors you provided, the nine results stiffness tensors are not the same.
Are the nine sets of strain and stress tensors one-to-one?

They are independent deformations, since each strain tensor has really only 1 independent non-zero entry. You have to solve for the elastic tensor using all 9 sets at once.

Sorry for my mistake!

Actually, I also used the stiffness tensor your website provids to check the 9 sets of data. I think every stress tensor equals the stiffness tensor multiplied by the strain tensor. But the result is that the 1st set of strain and stress keep this equaltion hold and the 2nd set of data do not

Are you saying the first stress/strain tensors don’t work with the MP Elastic tensor? They all look fine to me. Here is the code I used to verify that the deviation is less than 0.2 GPa.

import numpy as np
from pymatgen import MPRester
from pymatgen.analysis.elasticity import ElasticTensor, Stress, Strain

task_id = "mp-1017566"

cauchy_stresses = [ 
    [ [ 0.485835569, 0.0, -0.0 ], [ -0.0, -0.100418809, -0.0 ], [ -0.0, -0.0, -0.100418809 ] ], 
    [ [ -0.177069249, -0.0, -0.0 ], [ -0.0, -0.166530955, 1.024832653 ], [ -0.0, 1.024832653, -0.18702881500000001  ] ], 
    [ [  -0.10041880899999997, 0.0, 0.0 ], [ 0.0, 0.4858355689999999, 0.0 ], [ 0.0, 0.0, -0.10041880899999997  ] ], 
    [ [  -0.10041880899999997, 0.0, 0.0 ], [ 0.0, -0.10041880899999997, 0.0 ], [ 0.0, 0.0, 0.4858355689999999  ] ], 
    [ [  -0.16653095499999998, 0.0, -1.0248326529999998 ], [ 0.0, -0.17706924899999998, 0.0 ], [ -1.0248326529999998, 0.0, -0.187028815  ] ], 
    [ [  -0.17706924899999998, 0.0, 0.0 ], [ 0.0, -0.16653095499999998, -1.0248326529999998 ], [ 0.0, -1.0248326529999998, -0.187028815  ] ], 
    [ [  -0.16653095499999998, 0.0, 1.0248326529999998 ], [ 0.0, -0.17706924899999998, 0.0 ], [ 1.0248326529999998, 0.0, -0.187028815  ] ], 
    [ [  -0.187028815, -1.0248326529999998, 0.0 ], [ -1.0248326529999998, -0.16653095499999998, 0.0 ], [ 0.0, 0.0, -0.17706924899999998  ] ], 
    [ [  -0.187028815, 1.0248326529999998, 0.0 ], [ 1.0248326529999998, -0.16653095499999998, 0.0 ], [ 0.0, 0.0, -0.17706924899999998 ] ]
]

strains = [
    [ [ 0.004999903894075897, 0.0, 0.0 ], [ 0.0, -1.1102230246251565E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251565E-16  ] ], 
    [ [  -1.1102230246251565E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251565E-16, 0.005000051366872475 ], [ 0.0, 0.005000051366872475, -8.041918347911903E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, 0.004999903894075896, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, 0.004999903894075896  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, -0.005000051366872474 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ -0.005000051366872474, 0.0, -8.041918347911901E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.0 ], [ 0.0, -1.1102230246251563E-16, -0.005000051366872474 ], [ 0.0, -0.005000051366872474, -8.041918347911901E-8  ] ], 
    [ [  -1.1102230246251563E-16, 0.0, 0.005000051366872474 ], [ 0.0, -1.1102230246251563E-16, 0.0 ], [ 0.005000051366872474, 0.0, -8.041918347911901E-8  ] ], 
    [ [  -8.041918347911901E-8, -0.005000051366872474, 0.0 ], [ -0.005000051366872474, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16  ] ], 
    [ [  -8.041918347911901E-8, 0.005000051366872474, 0.0 ], [ 0.005000051366872474, -1.1102230246251563E-16, 0.0 ], [ 0.0, 0.0, -1.1102230246251563E-16 ] ]
]



mpr = MPRester()

elastic_tensor = mpr.query({"task_id":task_id},["elasticity"])[0]["elasticity"]["elastic_tensor"]
elastic_tensor = ElasticTensor.from_voigt(elastic_tensor)
cauchy_stresses = [Stress(c) for c in cauchy_stresses]
strains = [Strain(s) for s in strains]

test_stresses = [elastic_tensor.calculate_stress(s) for s in strains]

diff = [c-s for c,s in zip(cauchy_stresses,test_stresses)]
max_diff = np.max(np.abs(diff))
print(f"Max deviation is: {max_diff} GPa")

No, I mean the first stress/strain tensors do work with the MP Elastic tensor but the second stress/strain tensors do not.

I see your point. It looks like all of the shear components show up in the normal modes. That being said, this shouldn’t be allowed according to the symmetry rules, so they get set to 0. The stresses do show good linearity since the + and - deformations (see the last two for instance) mirror each other well. These shear stresses won’t contribute much to the Poisson ratio though. The most likely culprit is the fact that the e11,e22,e33 are all less than e44,e55,e66.