Add elastic properties calculator

matcalc/elastic.py

@@ -0,0 +1,87 @@
+"""Calculator for phonon properties."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from pymatgen.analysis.elasticity import ElasticTensor, Strain, DeformedStructureSet
+from pymatgen.io.ase import AseAtomsAdaptor
+
+from .base import PropCalc
+from .relaxation import RelaxCalc
+
+if TYPE_CHECKING:
+    from ase.calculators.calculator import Calculator
+
+
+class ElasticCalc(PropCalc):
+    """Calculator for elastic properties."""
+
+    def __init__(
+        self,
+        calculator: Calculator,
+        norm_strains=0.01,
+        shear_strains=0.01,
+        relax_structure=True,
+        fmax=0.1,
+    ):
+        """
+        Args:
+            calculator: ASE Calculator to use.
+            norm_strains: strain value to apply to each normal mode.
+            shear_strains: strain value to apply to each shear mode.
+            relax_structure: whether to relax the provided structure with the given calculator.
+            fmax: maximum force in the relaxed structure (if relax_structure).
+        """
+        self.calculator = calculator
+        self.norm_strains = norm_strains
+        self.shear_strains = shear_strains
+        self.relax_structure = relax_structure
+        self.fmax = fmax
+
+    def calc(self, structure) -> dict:
+        """
+        Calculates elastic properties of Pymatgen structure with units determined by the calculator.
+
+        Args:
+            structure: Pymatgen structure.
+
+        Returns: {
+            "elastic_tensor": Elastic tensor as a pymatgen ElasticTensor object,
+            "g_vrh": Voigt-Reuss-Hill shear modulus based on elastic tensor,
+            "k_vrh": Voigt-Reuss-Hill bulk modulus based on elastic tensor,
+            "y_mod": Young's modulus based on elastic tensor,
+        }
+        """
+        if self.relax_structure:
+            rcalc = RelaxCalc(self.calculator, fmax=self.fmax)
+            structure = rcalc.calc(structure)["final_structure"]
+
+        adaptor = AseAtomsAdaptor()
+        deformed_structure_set = DeformedStructureSet(
+            structure,
+            [
+                self.norm_strains,
+            ],
+            [
+                self.shear_strains,
+            ],
+        )
+        stresses = []
+        for deformed_structure in deformed_structure_set:
+            atoms = adaptor.get_atoms(deformed_structure)
+            atoms.calc = self.calculator
+            stresses.append(atoms.get_stress(voigt=False))
+
+        strains = [Strain.from_deformation(deformation) for deformation in deformed_structure_set.deformations]
+        atoms = adaptor.get_atoms(structure)
+        atoms.calc = self.calculator
+        elastic_tensor = ElasticTensor.from_independent_strains(
+            strains, stresses, eq_stress=atoms.get_stress(voigt=False)
+        )
+        return {
+            "elastic_tensor": elastic_tensor,
+            "g_vrh": elastic_tensor.g_vrh,
+            "k_vrh": elastic_tensor.k_vrh,
+            "y_mod": elastic_tensor.y_mod,
+        }

tests/test_elastic.py

@@ -0,0 +1,20 @@
+"""Tests for ElasticCalc class"""
+from __future__ import annotations
+
+import pytest
+
+from matcalc.elastic import ElasticCalc
+
+
+def test_ElasticCalc(LiFePO4, M3GNetUPCalc):
+    """Tests for ElasticCalc class"""
+    calculator = M3GNetUPCalc
+    ecalc = ElasticCalc(calculator, norm_strains=0.02, shear_strains=0.04, fmax=0.01)
+
+    # Test LiFePO4 with relaxation
+    results = ecalc.calc(LiFePO4)
+    assert results["elastic_tensor"].shape == (3, 3, 3, 3)
+    assert results["elastic_tensor"][0][1][1][0] == pytest.approx(0.6441543434291928, rel=0.0001)
+    assert results["k_vrh"] == pytest.approx(1.109278785217532, rel=0.0001)
+    assert results["g_vrh"] == pytest.approx(0.5946891263210372, rel=0.0001)
+    assert results["y_mod"] == pytest.approx(1513587180.4865916, rel=0.0001)