Elasticity

Compute the full elastic tensor and derived moduli from finite-strain stress calculations.

The ELASTICITY task computes the second-order elastic tensor of a crystal by:

1. Fully relaxing the structure with OPT.
2. Applying a set of normal and shear strains to generate deformed structures via pymatgen.analysis.elasticity.DeformedStructureSet.
3. Computing the stress tensor of each deformed structure.
4. Fitting a linear stress-strain relationship to extract the 6×6 Voigt elastic tensor.

The implementation is adapted from MatCalc.

Function signature

from mlip_arena.tasks import ELASTICITY

result = ELASTICITY(
    atoms=atoms,
    calculator=calc,
    optimizer="BFGSLineSearch",
    optimizer_kwargs=None,
    filter="FrechetCell",
    filter_kwargs=None,
    criterion=None,
    normal_strains=np.linspace(-0.01, 0.01, 4),
    shear_strains=np.linspace(-0.06, 0.06, 4),
    persist_opt=True,
    cache_opt=False,
)

Parameters

body.atoms · ase.Atoms · required

Input structure. A copy is made internally before relaxation.

body.calculator · ase.calculators.calculator.BaseCalculator · required

Calculator for energy, force, and stress evaluations.

body.optimizer · Optimizer | str · default: "BFGSLineSearch"

Optimizer passed to the initial OPT task. See structure optimization for valid values.

body.optimizer_kwargs · dict | None · default: "None"

Extra keyword arguments forwarded to the optimizer constructor.

body.filter · Filter | str | None · default: "FrechetCell"

Cell filter for the initial full relaxation. Defaults to "FrechetCell" to allow both positions and cell to relax.

body.filter_kwargs · dict | None · default: "None"

Extra keyword arguments forwarded to the filter constructor.

body.criterion · dict | None · default: "None"

Convergence criterion dict forwarded to the OPT task (e.g. {"fmax": 0.001}).

body.normal_strains · list[float] | ndarray · default: "np.linspace(-0.01, 0.01, 4)"

Normal strain magnitudes applied along each of the three Cartesian directions. Default covers ±1% in 4 steps. More points improve the linear fit accuracy.

body.shear_strains · list[float] | ndarray · default: "np.linspace(-0.06, 0.06, 4)"

Shear strain magnitudes applied for the off-diagonal components of the strain tensor. Default covers ±6% in 4 steps. Larger range needed because MLIPs often show weaker shear stiffness.

body.persist_opt · boolean · default: "true"

If True, the OPT result is persisted to Prefect’s result store.

body.cache_opt · boolean · default: "false"

If True, the OPT result is cached and reused if the same structure and calculator are provided again.

Return value

Returns a dict on success, or a Prefect State if the initial relaxation fails:

Key	Type	Description
elastic_tensor	pymatgen.analysis.elasticity.ElasticTensor	Full 3×3×3×3 (Voigt 6×6) elastic tensor in GPa, zeroed below numerical tolerance
residuals_sum	float	Sum of least-squares residuals from all stress-strain fits — a measure of fit quality

The ElasticTensor object provides derived moduli:

et = result["elastic_tensor"]

print("K_Voigt (GPa):", et.k_voigt)
print("K_Reuss (GPa):", et.k_reuss)
print("G_Voigt (GPa):", et.g_voigt)
print("G_Reuss (GPa):", et.g_reuss)
print("Young's modulus (GPa):", et.y_mod / 1e9)   # convert Pa → GPa
print("Poisson ratio:", et.universal_anisotropy)

Example

1. Import and set up

   import numpy as np
   from ase.build import bulk
   from mlip_arena.tasks import ELASTICITY
   from mlip_arena.tasks.utils import get_calculator
   from mlip_arena.models import MLIPEnum
   
   atoms = bulk("Cu", "fcc", a=3.6)
   calc = get_calculator(MLIPEnum.MACE_MP)
   
2. Run ELASTICITY

   result = ELASTICITY(
       atoms=atoms,
       calculator=calc,
       filter="FrechetCell",
       criterion={"fmax": 0.001},
       normal_strains=np.linspace(-0.01, 0.01, 5),
       shear_strains=np.linspace(-0.06, 0.06, 5),
   )
   
3. Extract elastic moduli

   et = result["elastic_tensor"]
   
   print("Elastic tensor (GPa):")
   print(et.voigt)
   
   print(f"Bulk modulus (Voigt):  {et.k_voigt:.1f} GPa")
   print(f"Bulk modulus (Reuss):  {et.k_reuss:.1f} GPa")
   print(f"Shear modulus (Voigt): {et.g_voigt:.1f} GPa")
   print(f"Shear modulus (Reuss): {et.g_reuss:.1f} GPa")
   print(f"Fit residuals sum: {result['residuals_sum']:.4f}")
   

For accurate elastic constants, use a tight convergence criterion (e.g. fmax=0.001) for the initial relaxation. Poorly relaxed structures produce inconsistent stress-strain data.