Calculation error of elastic tensors

#1

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

#2

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?

#3

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

#4

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.

#5

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)

#6

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

#7

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

#8

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.

#9

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

#10

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.

#11

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.

#12

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.