Changes to EOS and elasticity calculators

ElasticityCalc
- Accept floats as strains
- Fix small bug with list of length 1 and residuals
- Set default use_equilibrium to True (which is default when floats are used for strains)
- Modify testing values, add test for float strains.

EOS
- Minor linting, modify test values for new M3GNet version.

Other
- matgl to 0.9.1 in requirements

matcalc/elasticity.py

@@ -24,31 +24,35 @@ class ElasticityCalc(PropCalc):
     def __init__(
         self,
         calculator: Calculator,
-        norm_strains: tuple[float, ...] = (0.001, 0.003, 0.005, 0.01),
-        shear_strains: tuple[float, ...] = (0.001, 0.003, 0.005, 0.01),
+        norm_strains: tuple[float, ...] | float = (-0.01, -0.005, 0.005, 0.01),
+        shear_strains: tuple[float, ...] | float = (-0.06, -0.03, 0.03, 0.06),
         fmax: float = 0.1,
         relax_structure: bool = True,
-        use_equilibrium: bool = False,
+        use_equilibrium: bool = True,
     ) -> None:
         """
         Args:
             calculator: ASE Calculator to use.
-            norm_strains: strain values to apply to each normal mode.
-                Defaults to (0.001, 0.003, 0.005, 0.01).
-            shear_strains: strain values to apply to each shear mode.
-                Defaults to (0.001, 0.003, 0.005, 0.01).
+            norm_strains: single or multiple strain values to apply to each normal mode.
+                Defaults to (-0.01, -0.005, 0.005, 0.01).
+            shear_strains: single or multiple strain values to apply to each shear mode.
+                Defaults to (-0.06, -0.03, 0.03, 0.06).
             fmax: maximum force in the relaxed structure (if relax_structure). Defaults to 0.1.
             relax_structure: whether to relax the provided structure with the given calculator.
                 Defaults to True.
-            use_equilibrium: whether to use the equilibrium stress and strain.
-                Defaults to False.
+            use_equilibrium: whether to use the equilibrium stress and strain. Ignored and set
+                to True if either norm_strains or shear_strains has length 1 or is a float.
+                Defaults to True.
         """
         self.calculator = calculator
-        self.norm_strains = norm_strains
-        self.shear_strains = shear_strains
+        self.norm_strains = tuple(np.array([1]) * np.array(norm_strains))
+        self.shear_strains = tuple(np.array([1]) * np.array(shear_strains))
         self.relax_structure = relax_structure
         self.fmax = fmax
-        self.use_equilibrium = use_equilibrium
+        if len(self.norm_strains) > 1 and len(self.shear_strains) > 1:
+            self.use_equilibrium = use_equilibrium
+        else:
+            self.use_equilibrium = True
 
     def calc(self, structure: Structure) -> dict[str, float | ElasticTensor | Structure]:
         """
@@ -140,6 +144,6 @@ def _elastic_tensor_from_strains(
             for j in range(6):
                 fit = np.polyfit(strain[:, i], stress[:, j], 1, full=True)
                 c_ij[i, j] = fit[0][0]
-                residuals_sum += fit[1][0]
+                residuals_sum += fit[1][0] if len(fit[1]) > 0 else 0.0
         c = ElasticTensor.from_voigt(c_ij)
         return c.zeroed(tol), residuals_sum

matcalc/eos.py

@@ -5,8 +5,8 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from sklearn.metrics import r2_score
 from pymatgen.analysis.eos import BirchMurnaghan
+from sklearn.metrics import r2_score
 
 from .base import PropCalc
 from .relaxation import RelaxCalc
@@ -98,7 +98,7 @@ def calc(self, structure: Structure) -> dict:
         bm = BirchMurnaghan(volumes=volumes, energies=energies)
         bm.fit()
 
-        volumes, energies = zip(*sorted(list(zip(volumes, energies)), key=lambda i: i[0]))
+        volumes, energies = zip(*sorted(zip(volumes, energies), key=lambda i: i[0]))
 
         return {
             "eos": {"volumes": volumes, "energies": energies},

requirements-ci.txt

@@ -5,5 +5,5 @@ coveralls
 mypy
 ruff
 black
-matgl==0.9.0
+matgl==0.9.1
 chgnet==0.2.0

tests/test_elasticity.py

@@ -21,12 +21,12 @@ def test_elastic_calc(Li2O, M3GNetCalc):
     # Test Li2O with equilibrium structure
     results = ecalc.calc(Li2O)
     assert results["elastic_tensor"].shape == (3, 3, 3, 3)
-    assert results["elastic_tensor"][0][1][1][0] == pytest.approx(0.48107847906067014, rel=1e-3)
-    assert results["bulk_modulus_vrh"] == pytest.approx(0.6463790972016057, rel=1e-3)
-    assert results["shear_modulus_vrh"] == pytest.approx(0.4089633387159433, rel=1e-3)
-    assert results["youngs_modulus"] == pytest.approx(1013205376.4173204, rel=1e-3)
-    assert results["residuals_sum"] == pytest.approx(3.511205625736905e-08, rel=1e-2)
-    assert results["structure"].lattice.a == pytest.approx(3.3089642445687586, rel=1e-4)
+    assert results["elastic_tensor"][0][1][1][0] == pytest.approx(0.5014895636122672, rel=1e-3)
+    assert results["bulk_modulus_vrh"] == pytest.approx(0.6737897607182401, rel=1e-3)
+    assert results["shear_modulus_vrh"] == pytest.approx(0.4179219576918434, rel=1e-3)
+    assert results["youngs_modulus"] == pytest.approx(1038959096.5809333, rel=1e-3)
+    assert results["residuals_sum"] == pytest.approx(3.8487476828544434e-08, rel=1e-2)
+    assert results["structure"].lattice.a == pytest.approx(3.2885851104196875, rel=1e-4)
 
     # Test Li2O without the equilibrium structure
     ecalc = ElasticityCalc(
@@ -38,4 +38,17 @@ def test_elastic_calc(Li2O, M3GNetCalc):
     )
 
     results = ecalc.calc(Li2O)
-    assert results["residuals_sum"] == pytest.approx(2.673485844932801e-08, rel=1e-2)
+    assert results["residuals_sum"] == pytest.approx(2.9257237571340992e-08, rel=1e-2)
+
+    # Test Li2O with float
+    ecalc = ElasticityCalc(
+        calculator,
+        fmax=0.1,
+        norm_strains=0.004,
+        shear_strains=0.004,
+        use_equilibrium=True,
+    )
+
+    results = ecalc.calc(Li2O)
+    assert results["residuals_sum"] == 0.0
+    assert results["bulk_modulus_vrh"] == pytest.approx(0.6631894154825593, rel=1e-3)

tests/test_eos.py

@@ -14,17 +14,17 @@ def test_eos_calc(Li2O, LiFePO4, M3GNetCalc):
     results = pcalc.calc(Li2O)
 
     assert {*results} == {"eos", "r2_score_bm", "bulk_modulus_bm"}
-    assert results["bulk_modulus_bm"] == pytest.approx(69.868, rel=1e-2)
+    assert results["bulk_modulus_bm"] == pytest.approx(65.57980045603279, rel=1e-2)
     assert {*results["eos"]} == {"volumes", "energies"}
     assert results["eos"]["volumes"] == pytest.approx(
-        [18.70, 19.97, 21.30, 22.69, 24.14, 25.65, 27.22, 28.85, 30.55, 32.31, 34.14],
+        [18.38, 19.63, 20.94, 22.3, 23.73, 25.21, 26.75, 28.36, 30.02, 31.76, 33.55],
         rel=1e-3,
     )
     assert results["eos"]["energies"] == pytest.approx(
-        [-13.51, -13.78, -13.98, -14.11, -14.17, -14.19, -14.17, -14.12, -14.04, -13.94, -13.81],
+        [-13.52, -13.77, -13.94, -14.08, -14.15, -14.18, -14.16, -14.11, -14.03, -13.94, -13.83],
         rel=1e-3,
     )
-
+    pcalc = EOSCalc(calculator, relax_structure=False)
     results = list(pcalc.calc_many([Li2O, LiFePO4]))
     assert len(results) == 2
-    assert results[1]["bulk_modulus_bm"] == pytest.approx(58.716, rel=1e-2)
+    assert results[1]["bulk_modulus_bm"] == pytest.approx(54.5953851822073, rel=1e-2)