Add optional Voigt order output, test written output

Author: harveydevereux

janus_core/calculations/elasticity.py

@@ -74,6 +74,8 @@ class Elasticity(BaseCalculation):
     write_kwargs
         Keyword arguments to pass to ase.io.write to save generated structures.
         Default is {}.
+    write_voigt
+        Whether to write out in Voigt notation, Default is True.
 
     Attributes
     ----------
@@ -102,6 +104,7 @@ def __init__(
         write_results: bool = True,
         write_structures: bool = False,
         write_kwargs: OutputKwargs | None = None,
+        write_voigt: bool = True,
     ) -> None:
         """
         Initialise class.
@@ -149,6 +152,8 @@ def __init__(
         write_kwargs
             Keyword arguments to pass to ase.io.write to save generated structures.
             Default is {}.
+        write_voigt
+            Whether to write out in Voigt notation, Default is True.
         """
         read_kwargs, minimize_kwargs, write_kwargs = none_to_dict(
             read_kwargs, minimize_kwargs, write_kwargs
@@ -160,6 +165,7 @@ def __init__(
         self.write_results = write_results
         self.write_structures = write_structures
         self.write_kwargs = write_kwargs
+        self.write_voigt = write_voigt
 
         if (
             (self.minimize or self.minimize_all)
@@ -286,7 +292,13 @@ def run(self) -> ElasticTensor:
                         "homogeneous_poisson",
                     )
                 ]
-                for cijkl in self.elastic_tensor.flatten():
+
+                vals = (
+                    self.elastic_tensor.voigt
+                    if self.write_voigt
+                    else self.elastic_tensor
+                )
+                for cijkl in vals.flatten():
                     values.append(cijkl)
                 print(" ".join(map(str, values)), file=out)
         return self.elastic_tensor
@@ -311,21 +323,18 @@ def _calculate_elasticity(self) -> None:
         ]
         for i, deformed_structure in enumerate(self.deformed_structures):
             deformed_structure.calc = copy(self.struct.calc)
+            progress = f"{i} / {len(self.deformed_structures)}"
             if self.minimize_all:
                 if self.logger:
                     self.logger.info(
-                        "Minimising deformed structure %u / %u",
-                        i,
-                        len(self.deformed_structures),
+                        "Calculating stress for deformed structure " + progress
                     )
                 optimizer = GeomOpt(deformed_structure, **self.minimize_kwargs)
                 optimizer.run()
 
             if self.logger:
                 self.logger.info(
-                    "Calculating stress for deformed structure %u / %u",
-                    i,
-                    len(self.deformed_structures),
+                    "Calculating stress for deformed structure " + progress,
                 )
 
             # Always append first original structure

tests/test_elasticity.py

@@ -5,21 +5,110 @@
 from pathlib import Path
 
 from ase.build import bulk
+import numpy as np
 from pytest import approx
 
 from janus_core.calculations.elasticity import Elasticity
 from janus_core.helpers.mlip_calculators import choose_calculator
+from tests.utils import assert_log_contains
 
 DATA_PATH = Path(__file__).parent / "data"
 MODEL_PATH = Path(__file__).parent / "models" / "mace_mp_small.model"
 
 
 def test_calc_elasticity(tmp_path):
     """Test calculating elasticity for Aluminium."""
+    elasticity_path = tmp_path / "elasticity-elastic_tensor.dat"
+    log_file = tmp_path / "elasticity.log"
+
     struct = bulk("Al", crystalstructure="fcc")
     struct.calc = choose_calculator(arch="mace_mp", model=MODEL_PATH)
 
-    elasticity = Elasticity(struct, file_prefix=tmp_path / "elasticity")
+    elasticity = Elasticity(
+        struct, file_prefix=tmp_path / "elasticity", log_kwargs={"filename": log_file}
+    )
     elastic_tensor = elasticity.run()
 
     assert elastic_tensor.k_reuss == approx(80.34637735135381, rel=1e-3)
+
+    # Check geometry optimization run by default
+    assert_log_contains(
+        log_file,
+        includes=["Using filter", "Using optimizer", "Starting geometry optimization"],
+    )
+
+    assert elasticity_path.exists()
+    written_elasticity = np.loadtxt(elasticity_path)
+    # Defaulting to 6x6 Voigt notation
+    assert len(written_elasticity) == 45
+
+    for i, prop in enumerate(
+        (
+            "k_reuss",
+            "k_voigt",
+            "k_vrh",
+            "g_reuss",
+            "g_voigt",
+            "g_vrh",
+            "y_mod",
+            "universal_anisotropy",
+            "homogeneous_poisson",
+        )
+    ):
+        assert elastic_tensor.property_dict[prop] == approx(
+            written_elasticity[i], rel=1e-3
+        )
+
+    assert written_elasticity[9:] == approx(elastic_tensor.voigt.flatten(), rel=1e-3)
+
+
+def test_no_optimize_no_write_voigt(tmp_path):
+    """Test calculating elasticity for Aluminium without optimization."""
+    elasticity_path = tmp_path / "elasticity-elastic_tensor.dat"
+    log_file = tmp_path / "elasticity.log"
+
+    struct = bulk("Al", crystalstructure="fcc")
+    struct.calc = choose_calculator(arch="mace_mp", model=MODEL_PATH)
+
+    elasticity = Elasticity(
+        struct,
+        file_prefix=tmp_path / "elasticity",
+        log_kwargs={"filename": log_file},
+        minimize=False,
+        write_voigt=False,
+    )
+    elastic_tensor = elasticity.run()
+
+    assert elastic_tensor.k_reuss == approx(80.34637735135381, rel=1e-3)
+
+    # Check geometry optimization run by default
+    assert_log_contains(
+        log_file,
+        excludes=["Using filter", "Using optimizer", "Starting geometry optimization"],
+    )
+
+    assert elasticity_path.exists()
+    written_elasticity = np.loadtxt(elasticity_path)
+    # Selected to full tensor 3x3x3x3
+    assert len(written_elasticity) == 90
+
+    assert written_elasticity[0] == approx(elastic_tensor.k_reuss, rel=1e-3)
+
+    for i, prop in enumerate(
+        (
+            "k_reuss",
+            "k_voigt",
+            "k_vrh",
+            "g_reuss",
+            "g_voigt",
+            "g_vrh",
+            "y_mod",
+            "universal_anisotropy",
+            "homogeneous_poisson",
+        )
+    ):
+        assert elastic_tensor.property_dict[prop] == approx(
+            written_elasticity[i], rel=1e-3
+        )
+
+    assert written_elasticity[9:] == approx(elastic_tensor.flatten(), rel=1e-3)