Add diffusive thermal conductivity model proposed by Agne et al. (#3546)

* add agne model

* add test and source in doc string

* add test and source in doc string

* add types

* add due.dcite decorator for Doi("10.1039/C7EE03256K")

* refactor test

---------

Co-authored-by: anaik <[email protected]>
Co-authored-by: Janosh Riebesell <[email protected]>

Author: naik-aakash

pymatgen/analysis/elasticity/elastic.py

@@ -21,10 +21,13 @@
 from pymatgen.analysis.elasticity.stress import Stress
 from pymatgen.core.tensors import DEFAULT_QUAD, SquareTensor, Tensor, TensorCollection, get_uvec
 from pymatgen.core.units import Unit
+from pymatgen.util.due import Doi, due
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
 
+    from numpy.typing import ArrayLike
+
     from pymatgen.core import Structure
 
 
@@ -198,7 +201,7 @@ def y_mod(self) -> float:
         """
         return 9.0e9 * self.k_vrh * self.g_vrh / (3 * self.k_vrh + self.g_vrh)
 
-    def directional_poisson_ratio(self, n, m, tol: float = 1e-8):
+    def directional_poisson_ratio(self, n: ArrayLike, m: ArrayLike, tol: float = 1e-8) -> float:
         """
         Calculates the poisson ratio for a specific direction
         relative to a second, orthogonal direction.
@@ -215,15 +218,15 @@ def directional_poisson_ratio(self, n, m, tol: float = 1e-8):
         v *= -1 / self.compliance_tensor.einsum_sequence([n] * 4)
         return v
 
-    def directional_elastic_mod(self, n):
+    def directional_elastic_mod(self, n) -> float:
         """Calculates directional elastic modulus for a specific vector."""
         n = get_uvec(n)
         return self.einsum_sequence([n] * 4)
 
     @raise_if_unphysical
-    def trans_v(self, structure: Structure):
+    def trans_v(self, structure: Structure) -> float:
         """
-        Calculates transverse sound velocity (in SI units) using the
+        Calculates transverse sound velocity using the
         Voigt-Reuss-Hill average bulk modulus.
 
         Args:
@@ -233,19 +236,17 @@ def trans_v(self, structure: Structure):
             float: transverse sound velocity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
         n_atoms = structure.composition.num_atoms
         weight = float(structure.composition.weight)
-        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * volume)
+        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * structure.volume)
         if self.g_vrh < 0:
             raise ValueError("k_vrh or g_vrh is negative, sound velocity is undefined")
         return (1e9 * self.g_vrh / mass_density) ** 0.5
 
     @raise_if_unphysical
-    def long_v(self, structure: Structure):
+    def long_v(self, structure: Structure) -> float:
         """
-        Calculates longitudinal sound velocity (in SI units)
-        using the Voigt-Reuss-Hill average bulk modulus.
+        Calculates longitudinal sound velocity using the Voigt-Reuss-Hill average bulk modulus.
 
         Args:
             structure: pymatgen structure object
@@ -254,18 +255,16 @@ def long_v(self, structure: Structure):
             float: longitudinal sound velocity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
         n_atoms = structure.composition.num_atoms
         weight = float(structure.composition.weight)
-        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * volume)
+        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * structure.volume)
         if self.g_vrh < 0:
             raise ValueError("k_vrh or g_vrh is negative, sound velocity is undefined")
         return (1e9 * (self.k_vrh + 4 / 3 * self.g_vrh) / mass_density) ** 0.5
 
     @raise_if_unphysical
-    def snyder_ac(self, structure: Structure):
-        """
-        Calculates Snyder's acoustic sound velocity (in SI units).
+    def snyder_ac(self, structure: Structure) -> float:
+        """Calculates Snyder's acoustic sound velocity.
 
         Args:
             structure: pymatgen structure object
@@ -274,20 +273,19 @@ def snyder_ac(self, structure: Structure):
             float: Snyder's acoustic sound velocity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
         n_atoms = structure.composition.num_atoms
-        num_density = 1e30 * n_sites / volume
+        site_density = 1e30 * n_sites / structure.volume
         tot_mass = sum(e.atomic_mass for e in structure.species)
         avg_mass = 1.6605e-27 * tot_mass / n_atoms
         return (
             0.38483
             * avg_mass
             * ((self.long_v(structure) + 2 * self.trans_v(structure)) / 3) ** 3.0
-            / (300 * num_density ** (-2 / 3) * n_sites ** (1 / 3))
+            / (300 * site_density ** (-2 / 3) * n_sites ** (1 / 3))
         )
 
     @raise_if_unphysical
-    def snyder_opt(self, structure: Structure):
+    def snyder_opt(self, structure: Structure) -> float:
         """
         Calculates Snyder's optical sound velocity (in SI units).
 
@@ -298,18 +296,17 @@ def snyder_opt(self, structure: Structure):
             float: Snyder's optical sound velocity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
-        num_density = 1e30 * n_sites / volume
+        site_density = 1e30 * n_sites / structure.volume
         return (
             1.66914e-23
             * (self.long_v(structure) + 2 * self.trans_v(structure))
             / 3.0
-            / num_density ** (-2 / 3)
+            / site_density ** (-2 / 3)
             * (1 - n_sites ** (-1 / 3))
         )
 
     @raise_if_unphysical
-    def snyder_total(self, structure: Structure):
+    def snyder_total(self, structure: Structure) -> float:
         """
         Calculates Snyder's total sound velocity (in SI units).
 
@@ -322,7 +319,7 @@ def snyder_total(self, structure: Structure):
         return self.snyder_ac(structure) + self.snyder_opt(structure)
 
     @raise_if_unphysical
-    def clarke_thermalcond(self, structure: Structure):
+    def clarke_thermalcond(self, structure: Structure) -> float:
         """
         Calculates Clarke's thermal conductivity (in SI units).
 
@@ -333,16 +330,15 @@ def clarke_thermalcond(self, structure: Structure):
             float: Clarke's thermal conductivity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
         tot_mass = sum(e.atomic_mass for e in structure.species)
         n_atoms = structure.composition.num_atoms
         weight = float(structure.composition.weight)
         avg_mass = 1.6605e-27 * tot_mass / n_atoms
-        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * volume)
+        mass_density = 1.6605e3 * n_sites * weight / (n_atoms * structure.volume)
         return 0.87 * 1.3806e-23 * avg_mass ** (-2 / 3) * mass_density ** (1 / 6) * self.y_mod**0.5
 
     @raise_if_unphysical
-    def cahill_thermalcond(self, structure: Structure):
+    def cahill_thermalcond(self, structure: Structure) -> float:
         """
         Calculates Cahill's thermal conductivity (in SI units).
 
@@ -353,21 +349,47 @@ def cahill_thermalcond(self, structure: Structure):
             float: Cahill's thermal conductivity (in SI units)
         """
         n_sites = len(structure)
-        volume = structure.volume
-        num_density = 1e30 * n_sites / volume
-        return 1.3806e-23 / 2.48 * num_density ** (2 / 3) * (self.long_v(structure) + 2 * self.trans_v(structure))
+        site_density = 1e30 * n_sites / structure.volume
+        return 1.3806e-23 / 2.48 * site_density ** (2 / 3) * (self.long_v(structure) + 2 * self.trans_v(structure))
+
+    @due.dcite(
+        Doi("10.1039/C7EE03256K"),
+        description="Minimum thermal conductivity in the context of diffuson-mediated thermal transport",
+    )
+    @raise_if_unphysical
+    def agne_diffusive_thermalcond(self, structure: Structure) -> float:
+        """
+        Calculates Agne's diffusive thermal conductivity (in SI units).
+
+        Please cite the original authors if using this method
+        M. T. Agne, R. Hanus, G. J. Snyder, Energy Environ. Sci. 2018, 11, 609-616.
+        DOI: https://doi.org/10.1039/C7EE03256K
+
+        Args:
+            structure: pymatgen structure object
+
+        Returns:
+            float: Agne's diffusive thermal conductivity (in SI units)
+        """
+        n_sites = len(structure)
+        site_density = 1e30 * n_sites / structure.volume
+        return (
+            0.76
+            * (site_density ** (2 / 3))
+            * 1.3806e-23
+            * ((1 / 3) * (2 * self.trans_v(structure) + self.long_v(structure)))
+        )
 
     @raise_if_unphysical
-    def debye_temperature(self, structure: Structure):
+    def debye_temperature(self, structure: Structure) -> float:
         """
-        Estimates the Debye temperature from longitudinal and
-        transverse sound velocities.
+        Estimates the Debye temperature from longitudinal and transverse sound velocities.
 
         Args:
             structure: pymatgen structure object
 
         Returns:
-            float: debye temperature (in SI units)
+            float: Debye temperature (in SI units)
         """
         v0 = structure.volume * 1e-30 / len(structure)
         vl, vt = self.long_v(structure), self.trans_v(structure)

tests/analysis/elasticity/test_elastic.py

@@ -117,52 +117,59 @@ def test_directional_poisson_ratio(self):
 
     def test_structure_based_methods(self):
         # trans_velocity
-        assert self.elastic_tensor_1.trans_v(self.structure) == approx(1996.35019877)
+        struct = self.structure
+        assert self.elastic_tensor_1.trans_v(struct) == approx(1996.35019877)
         # long_velocity
-        assert self.elastic_tensor_1.long_v(self.structure) == approx(3534.68123832)
+        assert self.elastic_tensor_1.long_v(struct) == approx(3534.68123832)
         # Snyder properties
-        assert self.elastic_tensor_1.snyder_ac(self.structure) == approx(18.06127074)
-        assert self.elastic_tensor_1.snyder_opt(self.structure) == approx(0.18937465)
-        assert self.elastic_tensor_1.snyder_total(self.structure) == approx(18.25064540)
+        assert self.elastic_tensor_1.snyder_ac(struct) == approx(18.06127074)
+        assert self.elastic_tensor_1.snyder_opt(struct) == approx(0.18937465)
+        assert self.elastic_tensor_1.snyder_total(struct) == approx(18.25064540)
         # Clarke
-        assert self.elastic_tensor_1.clarke_thermalcond(self.structure) == approx(0.3450307)
+        assert self.elastic_tensor_1.clarke_thermalcond(struct) == approx(0.3450307)
         # Cahill
-        assert self.elastic_tensor_1.cahill_thermalcond(self.structure) == approx(0.37896275)
+        cahill_thermal_cond = self.elastic_tensor_1.cahill_thermalcond(struct)
+        assert cahill_thermal_cond == approx(0.37896275)
+        # Agne
+        agne_thermal_cond = self.elastic_tensor_1.agne_diffusive_thermalcond(struct)
+        assert agne_thermal_cond == approx(0.23808966)
+        # Test Agne / Cahill factor
+        assert agne_thermal_cond / cahill_thermal_cond == approx(0.6282666)
         # Debye
-        assert self.elastic_tensor_1.debye_temperature(self.structure) == approx(198.8037985019)
+        assert self.elastic_tensor_1.debye_temperature(struct) == approx(198.8037985019)
 
         # structure-property dict
-        sprop_dict = self.elastic_tensor_1.get_structure_property_dict(self.structure)
-        assert sprop_dict["long_v"] == approx(3534.68123832)
-        for val in sprop_dict.values():
+        struct_prop_dict = self.elastic_tensor_1.get_structure_property_dict(struct)
+        assert struct_prop_dict["long_v"] == approx(3534.68123832)
+        for val in struct_prop_dict.values():
             assert not isinstance(val, FloatWithUnit)
-        for k, v in sprop_dict.items():
-            if k == "structure":
-                assert v == self.structure
+        for key, val in struct_prop_dict.items():
+            if key == "structure":
+                assert val == struct
             else:
-                f = getattr(self.elastic_tensor_1, k)
-                if callable(f):
-                    assert getattr(self.elastic_tensor_1, k)(self.structure) == approx(v)
+                attr = getattr(self.elastic_tensor_1, key)
+                if callable(attr):
+                    assert getattr(self.elastic_tensor_1, key)(struct) == approx(val)
                 else:
-                    assert getattr(self.elastic_tensor_1, k) == approx(v)
+                    assert getattr(self.elastic_tensor_1, key) == approx(val)
 
         # Test other sprop dict modes
-        sprop_dict = self.elastic_tensor_1.get_structure_property_dict(self.structure, include_base_props=False)
-        assert "k_vrh" not in sprop_dict
+        struct_prop_dict = self.elastic_tensor_1.get_structure_property_dict(struct, include_base_props=False)
+        assert "k_vrh" not in struct_prop_dict
 
         # Test ValueError being raised for structure properties
         test_et = deepcopy(self.elastic_tensor_1)
         test_et[0][0][0][0] = -100000
         prop_dict = test_et.property_dict
-        for attr_name in sprop_dict:
+        for attr_name in struct_prop_dict:
             if attr_name not in ([*prop_dict, "structure"]):
                 with pytest.raises(
                     ValueError, match="Bulk or shear modulus is negative, property cannot be determined"
                 ):
-                    getattr(test_et, attr_name)(self.structure)
+                    getattr(test_et, attr_name)(struct)
         with pytest.raises(ValueError, match="Bulk or shear modulus is negative, property cannot be determined"):
-            test_et.get_structure_property_dict(self.structure)
-        noval_sprop_dict = test_et.get_structure_property_dict(self.structure, ignore_errors=True)
+            test_et.get_structure_property_dict(struct)
+        noval_sprop_dict = test_et.get_structure_property_dict(struct, ignore_errors=True)
         assert noval_sprop_dict["snyder_ac"] is None
 
     def test_new(self):