ElasticTensor jobs issues

[Leimeroth]

Summary

ElasticTensor suffers from several problems that were discussed and partly fixed in #246 already. Basically the conclusion was that the tensor would need a secondary symmetrization step that is currently not implemented. Since I recently calculated some elastic tensors again and the chosen parameters seemed to overly influence the results I was wondering whether there is another issue with the way the ElasticTensor is calculated. Additionally while the lacking symmetrization step will lead to slight differences a simple relations like C66 = (C11-C12)/2 for HCP should not give completely wrong answers, but still it does even for a simple calculation with an interatomic potential where stresses should be accurate, no k-mesh changes etc. So I compared results to those of pymatgen and found that instead of an error of 20 GPa pymatgen has one 2.7, which seems far more reasonable for fitted values. Also it allows symmetrization, which reduces the error to somethin like 1e-10.

pyiron Version and Platform

Newest github version

Expected Behavior

At least somewhat close results between pyiron ElasticTensor and pymatgen for a simple test case.

Actual Behavior

Pymatgen and pyiron give very different answers when calculating an elastic tensor. Since pymatgen seems to follow simple symmetry laws much more closely I suspect some issue with the pyiron protocol.

Steps to Reproduce

Unpack and run attached notebook.

Further Information, Files, and Links

TestET.zip

[jan-janssen]

I would recommend to use ElasticMatrixJob instead of ElasticTensor: https://github.com/pyiron/pyiron_gpl/blob/main/pyiron_gpl/elastic/elastic.py

Simple example:

import pyiron_gpl
from pyiron_atomistics import Project

pr = Project("elastic")
pr.remove_jobs(silently=True, recursive=True)

potential = "2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1"
structure = pr.create_ase_bulk("Mg")

lmp_mini = pr.create_job(pr.job_type.Lammps, "lmp_mini")
lmp_mini.structure = structure
lmp_mini.potential = potential
lmp_mini.calc_minimize(pressure=0.0)
lmp_mini.run()

lmp_elastic = pr.create_job(pr.job_type.Lammps, "lmp_mg_elastic")
lmp_elastic.structure = lmp_mini.get_structure()
lmp_elastic.potential = potential
elastic_tensor = lmp_elastic.create_job(pr.job_type.ElasticTensor, "elastic_mg")
elastic_tensor.run()

lmp_elastic = pr.create_job(pr.job_type.Lammps, "lmp_mg_elastic2")
lmp_elastic.structure = lmp_mini.get_structure()
lmp_elastic.potential = potential
elastic_matrix = lmp_elastic.create_job(pr.job_type.ElasticMatrixJob, "elastic_mg2")
elastic_matrix.run()

print(elastic_tensor['output/elastic_tensor'], elastic_matrix["output/elasticmatrix"]["C"])

[jan-janssen]

array([[ 7.93560420e+01,  2.52026820e+01,  2.56915107e+01, 2.82521551e-01,  1.46283347e-01, -1.30715207e+01],
       [ 2.58065280e+01,  7.78505073e+01,  2.66521625e+01, -3.20806650e-01,  7.95525250e-01, -1.08234795e+01],
       [ 1.61339533e+01,  1.57041745e+01,  7.49109879e+01, -1.60893194e-01,  3.75390940e-01, -7.60173857e+00],
       [-1.17196275e-02, -1.49410145e-01, -3.73451781e-02, 1.34870536e+01,  4.30497647e-01,  1.74548664e-01],
       [-1.02014412e-01, -9.21805568e-02,  1.25646772e-01, -1.85145297e-01,  1.38863244e+01,  1.32139892e-01],
       [ 1.24940814e+01,  1.10353797e+01,  1.22774143e+01, 1.67962736e-01,  4.56404352e-01,  2.17401300e+01]])

vs.

array([[68.71081125, 26.25107837, 16.48307692,  0.        ,  0.        , 0.        ],
       [26.25107837, 68.71081125, 16.48307692,  0.        ,  0.        , 0.        ],
       [16.48307692, 16.48307692, 70.47791543,  0.        ,  0.        , 0.        ],
       [ 0.        ,  0.        ,  0.        , 13.47390771,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , 13.47390771,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        , 21.22986644]])

[Leimeroth]

Symmetrized pymatgen values: C11: 147.5651423032322, C12: 68.1513539891288, C13: 74.05067707520763,C33: 167.72366499865117, C44: 44.82160029471396, C66: 39.70689415682077

ElasticMatrix: C11: 147.1155101588847, C12: 81.26022428115812, C13: 77.56639448136447,C33: 168.55305184061808, C44: 44.202016915046904, C66: 32.92764293886329

ElasticTensor: C11: 161.22549494586735, C12: 69.33523917507553, C13: 95.33237480466062, C33: 167.051460432368, C44: 52.830961255120506, C66: 26.68943522375893

So basically 3 implementations 3 answers. Pymatgen and ElasticMatrix agree well besides C12 (and C66, but C12 and C66 are not independent) Pymatgen offers the added benefits of rotating the tensor to IEEE convention very easily. In the case of ElasticTensor I find it irritating that even C11 differs so much.

[samwaseda]

Can you post the code for pymatgen? Maybe we're gonna simply use their code.

[Leimeroth]

Basically I did:

import pymatgen.analysis.elasticity as elastic
import pymatgen.core.tensors as pmg_tensors
pmgS = s.to_pymatgen()
deformed_pmgStructures = elastic.DeformedStructureSet(pmgS, symmetry=True)

counter = 0
for deformed_pmgS in defomred_pmgStructures.deformed_structures:
    deformed_s = pr.create.structure.from_pymatgen(deformed_pmgS)
    job = pr.create.job.Lammps(f"Deformed_{counter}")
    job.structure = deformed_s
    job.potential = relax.potential
    job.run()
    counter +=1

cauchy_stresses = []
# PK2_stresses are recommended in
# https://pymatgen.org/pymatgen.analysis.elasticity.elastic.html diff fit function
# Probably because strains are finite
PK2_stresses = []
strains = []
for i in range(counter):
    job = pr.load(f"Deformed_{i}")
    cauchy_stress = -job["output/generic/pressures"][0]
    cauchy_stresses.append(cauchy_stress)
    stress = elastic.Stress(cauchy_stress)
    deformation = defomred_pmgStructures.deformations[i]
    strains.append(deformation.green_lagrange_strain)
    PK2_stresses.append(stress.piola_kirchoff_2(deformation))

tensor = elastic.diff_fit(strains, PK_stresses)[0]
# Symmetrization step
tensor_sym = tensor.fit_to_structure(pmgS)
# Voigt notation
tv_pmg = tensor_sym.voigt

``

[Leimeroth]

Probably necessary to run tensor_populated = tensor.populate(pmgS, prec=0.1) had some weird results for cubic structures when using symmtery reduced datasets without it.

[jyang2009]

Dear all,

I'm trying to do a elastic tensor calculation with lammps and was reading this issue. My question is: is it general that one needs to perform symmetrization for the elastic tensor of a lammps calculation? Why lammps by itself cannot produce the correctly symmetrized tensor? I'm also reading the example of calculating the elastic tensor given by lammps (here), there they don't seem to suggest that additional symmetry operation is needed. Thanks a lot in advance.

Best, Jing

[Leimeroth]

is it general that one needs to perform symmetrization for the elastic tensor of a lammps calculation?

In many codes calculating elastic tensors no symmetrization happens, because deformations are applied in a way that only calculated independent elastic constants once. If the elastic tensor itself is calculated without respect to symmetry some small errors between the values that should be dependent due to the crystal symmetry could occur for various reasons. For example in DFT calculations the k-point mesh could cause trouble. Numeric errors can also occur in the process of fitting the stress or energy. As noted in my original question, the current pyiron protocol seems to have a more serious issue though. That is why I brought up the hcp example, for which some small deviations can probably be expected due to the numerical issues, but the actual results massively violate the symmetry constraints.

Why lammps by itself cannot produce the correctly symmetrized tensor?

Lammps itself only calculates the stresses or energies. The fitting of elastic constants is done in some other code or within the script they provide. I never checked the results of the Lammps script, but feel free to try it with something like a hcp structure, manually check whether symmetry constraints are violated and report back. I would be very interested to see the results. They could be added to my funny table showing that three implementations actually give three different results. Maybe a fourth one gives yet another answer.

Also, I just noticed that lammps also does some symmetrization for cubic systems, see these lines in their script

# Average moduli for cubic crystals

variable C11cubic equal (${C11all}+${C22all}+${C33all})/3.0
variable C12cubic equal (${C12all}+${C13all}+${C23all})/3.0
variable C44cubic equal (${C44all}+${C55all}+${C66all})/3.0

Pymatgen just does it in a more general fashion for arbitrary crystal systems by analyzing the symmetry beforehand.

[jyang2009]

Dear Niklas,

Thanks a lot for your insights. This is very helpful. The lammps script indeed gives results that violate the symmetry constraints. Specifically for the hcp structure I got a large non-zero C16, which is why I started looking into this. I also found this link which discusses how to use pymatgen to symmetrize the elastic tensor. I copy it here in case someone is interested.

Best, Jing