diff --git a/pymatgen/analysis/elasticity/elastic.py b/pymatgen/analysis/elasticity/elastic.py index 21867936bfc..73a243aac64 100644 --- a/pymatgen/analysis/elasticity/elastic.py +++ b/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) diff --git a/tests/analysis/elasticity/test_elastic.py b/tests/analysis/elasticity/test_elastic.py index bc5d9813538..82f92f9f1c0 100644 --- a/tests/analysis/elasticity/test_elastic.py +++ b/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):