From 1e7f93afc987c2798c1241f3fd1a0ac098037a11 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Thu, 16 Nov 2023 12:12:18 -0800 Subject: [PATCH] add units to ElasticTensor methods fix typos like lower case voigt --- pymatgen/analysis/elasticity/elastic.py | 91 +++++++++++++------------ pymatgen/analysis/elasticity/strain.py | 29 ++++---- pymatgen/analysis/nmr.py | 4 +- pymatgen/analysis/piezo.py | 2 +- pymatgen/core/tensors.py | 12 ++-- 5 files changed, 72 insertions(+), 66 deletions(-) diff --git a/pymatgen/analysis/elasticity/elastic.py b/pymatgen/analysis/elasticity/elastic.py index 537879e0253..d61423018d8 100644 --- a/pymatgen/analysis/elasticity/elastic.py +++ b/pymatgen/analysis/elasticity/elastic.py @@ -9,7 +9,7 @@ import itertools import math import warnings -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal import numpy as np import sympy as sp @@ -23,6 +23,8 @@ from pymatgen.core.units import Unit if TYPE_CHECKING: + from collections.abc import Sequence + from pymatgen.core import Structure @@ -75,7 +77,7 @@ def calculate_stress(self, strain): strain = np.array(strain) if strain.shape == (6,): strain = Strain.from_voigt(strain) - assert strain.shape == (3, 3), "Strain must be 3x3 or voigt-notation" + assert strain.shape == (3, 3), "Strain must be 3x3 or Voigt notation" stress_matrix = self.einsum_sequence([strain] * (self.order - 1)) / factorial(self.order - 1) return Stress(stress_matrix) @@ -121,24 +123,22 @@ def wrapper(self, *args, **kwargs): class ElasticTensor(NthOrderElasticTensor): """ - This class extends Tensor to describe the 3x3x3x3 - second-order elastic tensor, C_{ijkl}, with various - methods for estimating other properties derived from - the second order elastic tensor. + This class extends Tensor to describe the 3x3x3x3 second-order elastic tensor, + C_{ijkl}, with various methods for estimating other properties derived from the second + order elastic tensor (e. g. bulk modulus, shear modulus, Young's modulus, Poisson's ratio) + in units of eV/A^3. """ def __new__(cls, input_array, tol: float = 1e-4): """ - Create an ElasticTensor object. The constructor throws an error if - the shape of the input_matrix argument is not 3x3x3x3, i. e. in true - tensor notation. Issues a warning if the input_matrix argument does - not satisfy standard symmetries. Note that the constructor uses - __new__ rather than __init__ according to the standard method of - subclassing numpy ndarrays. + Create an ElasticTensor object. The constructor throws an error if the shape of + the input_matrix argument is not 3x3x3x3, i. e. in true tensor notation. Issues a + warning if the input_matrix argument does not satisfy standard symmetries. Note + that the constructor uses __new__ rather than __init__ according to the standard + method of subclassing numpy ndarrays. Args: - input_array (3x3x3x3 array-like): the 3x3x3x3 array-like - representing the elastic tensor + input_array (3x3x3x3 array-like): the 3x3x3x3 array-like representing the elastic tensor tol (float): tolerance for initial symmetry test of tensor """ @@ -148,33 +148,32 @@ def __new__(cls, input_array, tol: float = 1e-4): @property def compliance_tensor(self): """ - Returns the Voigt-notation compliance tensor, - which is the matrix inverse of the - Voigt-notation elastic tensor. + Returns the Voigt notation compliance tensor, which is the matrix + inverse of the Voigt notation elastic tensor. """ s_voigt = np.linalg.inv(self.voigt) return ComplianceTensor.from_voigt(s_voigt) @property def k_voigt(self) -> float: - """Returns the K_v bulk modulus.""" + """Returns the K_v bulk modulus (in eV/A^3).""" return self.voigt[:3, :3].mean() @property def g_voigt(self) -> float: - """Returns the G_v shear modulus.""" + """Returns the G_v shear modulus (in eV/A^3).""" return ( 2 * self.voigt[:3, :3].trace() - np.triu(self.voigt[:3, :3]).sum() + 3 * self.voigt[3:, 3:].trace() ) / 15.0 @property def k_reuss(self) -> float: - """Returns the K_r bulk modulus.""" + """Returns the K_r bulk modulus (in eV/A^3).""" return 1 / self.compliance_tensor.voigt[:3, :3].sum() @property def g_reuss(self) -> float: - """Returns the G_r shear modulus.""" + """Returns the G_r shear modulus (in eV/A^3).""" return 15 / ( 8 * self.compliance_tensor.voigt[:3, :3].trace() - 4 * np.triu(self.compliance_tensor.voigt[:3, :3]).sum() @@ -183,12 +182,12 @@ def g_reuss(self) -> float: @property def k_vrh(self) -> float: - """Returns the K_vrh (Voigt-Reuss-Hill) average bulk modulus.""" + """Returns the K_vrh (Voigt-Reuss-Hill) average bulk modulus (in eV/A^3).""" return 0.5 * (self.k_voigt + self.k_reuss) @property def g_vrh(self) -> float: - """Returns the G_vrh (Voigt-Reuss-Hill) average shear modulus.""" + """Returns the G_vrh (Voigt-Reuss-Hill) average shear modulus (in eV/A^3).""" return 0.5 * (self.g_voigt + self.g_reuss) @property @@ -452,7 +451,7 @@ def from_pseudoinverse(cls, strains, stresses): stresses (Nx3x3 array-like): list or array of stresses strains (Nx3x3 array-like): list or array of strains """ - # convert the stress/strain to Nx6 arrays of voigt-notation + # convert the stress/strain to Nx6 arrays of voigt notation warnings.warn( "Pseudo-inverse fitting of Strain/Stress lists may yield " "questionable results from vasp data, use with caution." @@ -486,11 +485,11 @@ def from_independent_strains(cls, strains, stresses, eq_stress=None, vasp=False, if len(set(ss_dict) - set(strain_states)) > 0: warnings.warn("Extra strain states in strain-stress pairs are neglected in independent strain fitting") c_ij = np.zeros((6, 6)) - for i in range(6): - istrains = ss_dict[strain_states[i]]["strains"] - istresses = ss_dict[strain_states[i]]["stresses"] - for j in range(6): - c_ij[i, j] = np.polyfit(istrains[:, i], istresses[:, j], 1)[0] + for ii in range(6): + strains = ss_dict[strain_states[ii]]["strains"] + stresses = ss_dict[strain_states[ii]]["stresses"] + for jj in range(6): + c_ij[ii, jj] = np.polyfit(strains[:, ii], stresses[:, jj], 1)[0] if vasp: c_ij *= -0.1 # Convert units/sign convention of vasp stress tensor c = cls.from_voigt(c_ij) @@ -525,16 +524,15 @@ class ElasticTensorExpansion(TensorCollection): (e. g. symmetrization, voigt conversion, etc.). """ - def __init__(self, c_list): + def __init__(self, c_list: Sequence) -> None: """ Initialization method for ElasticTensorExpansion. Args: - c_list (list or tuple): sequence of Tensor inputs - or tensors from which the elastic tensor - expansion is constructed. + c_list (list or tuple): sequence of Tensor inputs or tensors from which + the elastic tensor expansion is constructed. """ - c_list = [NthOrderElasticTensor(c, check_rank=4 + i * 2) for i, c in enumerate(c_list)] + c_list = [NthOrderElasticTensor(c, check_rank=4 + idx * 2) for idx, c in enumerate(c_list)] super().__init__(c_list) @classmethod @@ -547,14 +545,14 @@ def from_diff_fit(cls, strains, stresses, eq_stress=None, tol: float = 1e-10, or return cls(c_list) @property - def order(self): + def order(self) -> int: """ Order of the elastic tensor expansion, i. e. the order of the highest included set of elastic constants. """ return self[-1].order - def calculate_stress(self, strain): + def calculate_stress(self, strain) -> float: """ Calculate's a given elastic tensor's contribution to the stress using Einstein summation. @@ -562,7 +560,7 @@ def calculate_stress(self, strain): return sum(c.calculate_stress(strain) for c in self) def energy_density(self, strain, convert_GPa_to_eV=True): - """Calculates the elastic energy density due to a strain.""" + """Calculates the elastic energy density due to a strain in eV/A^3 or GPa.""" return sum(c.energy_density(strain, convert_GPa_to_eV) for c in self) def get_ggt(self, n, u): @@ -579,7 +577,7 @@ def get_ggt(self, n, u): 2 * gk ) - def get_tgt(self, temperature=None, structure=None, quad=None): + def get_tgt(self, temperature: float | None = None, structure: Structure = None, quad=None): """ Gets the thermodynamic Gruneisen tensor (TGT) by via an integration of the GGT weighted by the directional heat @@ -592,7 +590,7 @@ def get_tgt(self, temperature=None, structure=None, quad=None): Args: temperature (float): Temperature in kelvin, if not specified will return non-cv-normalized value - structure (float): Structure to be used in directional heat + structure (Structure): Structure to be used in directional heat capacity determination, only necessary if temperature is specified quad (dict): quadrature for integration, should be @@ -601,6 +599,7 @@ def get_tgt(self, temperature=None, structure=None, quad=None): """ if temperature and not structure: raise ValueError("If using temperature input, you must also include structure") + assert structure # for mypy quad = quad or DEFAULT_QUAD points = quad["points"] @@ -677,7 +676,9 @@ def omega(self, structure: Structure, n, u): vel = (1e9 * self[0].einsum_sequence([n, u, n, u]) / (weight / vol)) ** 0.5 return vel / l0 - def thermal_expansion_coeff(self, structure: Structure, temperature, mode="debye"): + def thermal_expansion_coeff( + self, structure: Structure, temperature: float, mode: Literal["dulong - petit", "debye"] = "debye" + ): """ Gets thermal expansion coefficient from third-order constants. @@ -915,7 +916,7 @@ def find_eq_stress(strains, stresses, tol: float = 1e-10): def get_strain_state_dict(strains, stresses, eq_stress=None, tol: float = 1e-10, add_eq=True, sort=True): """ - Creates a dictionary of voigt-notation stress-strain sets + Creates a dictionary of voigt notation stress-strain sets keyed by "strain state", i. e. a tuple corresponding to the non-zero entries in ratios to the lowest nonzero value, e.g. [0, 0.1, 0, 0.2, 0, 0] -> (0,1,0,2,0,0) @@ -973,9 +974,9 @@ def generate_pseudo(strain_states, order=3): """Generates the pseudo-inverse for a given set of strains. Args: - strain_states (6xN array like): a list of voigt-notation - "strain-states", i. e. perturbed indices of the strain - as a function of the smallest strain e. g. (0, 1, 0, 0, 1, 0) + strain_states (6xN array like): a list of Voigt-notation strain-states, + i. e. perturbed indices of the strain as a function of the smallest + strain e. g. (0, 1, 0, 0, 1, 0) order (int): order of pseudo-inverse to calculate Returns: @@ -1009,7 +1010,7 @@ def generate_pseudo(strain_states, order=3): def get_symbol_list(rank, dim=6): """ - Returns a symbolic representation of the voigt-notation + Returns a symbolic representation of the Voigt-notation tensor that places identical symbols for entries related by index transposition, i. e. C_1121 = C_1211 etc. diff --git a/pymatgen/analysis/elasticity/strain.py b/pymatgen/analysis/elasticity/strain.py index 7d4c3531bd8..00e428a1be0 100644 --- a/pymatgen/analysis/elasticity/strain.py +++ b/pymatgen/analysis/elasticity/strain.py @@ -17,6 +17,8 @@ from pymatgen.core.tensors import SquareTensor, symmetry_reduce if TYPE_CHECKING: + from collections.abc import Sequence + from numpy.typing import ArrayLike from pymatgen.core.structure import Structure @@ -79,20 +81,20 @@ def apply_to_structure(self, structure: Structure): return def_struct @classmethod - def from_index_amount(cls, matrixpos, amt): + def from_index_amount(cls, matrix_pos, amt): """ Factory method for constructing a Deformation object from a matrix position and amount. Args: - matrixpos (tuple): tuple corresponding the matrix position to + matrix_pos (tuple): tuple corresponding the matrix position to have a perturbation added amt (float): amount to add to the identity matrix at position - matrixpos + matrix_pos """ - f = np.identity(3) - f[matrixpos] += amt - return cls(f) + ident = np.identity(3) + ident[matrix_pos] += amt + return cls(ident) class DeformedStructureSet(collections.abc.Sequence): @@ -101,7 +103,13 @@ class that generates a set of independently deformed structures that can be used to calculate linear stress-strain response. """ - def __init__(self, structure: Structure, norm_strains=None, shear_strains=None, symmetry=False): + def __init__( + self, + structure: Structure, + norm_strains: Sequence[float] = (-0.01, -0.005, 0.005, 0.01), + shear_strains: Sequence[float] = (-0.06, -0.03, 0.03, 0.06), + symmetry=False, + ) -> None: """ Construct the deformed geometries of a structure. Generates m + n deformed structures according to the supplied parameters. @@ -109,14 +117,11 @@ def __init__(self, structure: Structure, norm_strains=None, shear_strains=None, Args: structure (Structure): structure to undergo deformation norm_strains (list of floats): strain values to apply - to each normal mode. + to each normal mode. Defaults to (-0.01, -0.005, 0.005, 0.01). shear_strains (list of floats): strain values to apply - to each shear mode. + to each shear mode. Defaults to (-0.06, -0.03, 0.03, 0.06). symmetry (bool): whether or not to use symmetry reduction. """ - norm_strains = norm_strains or [-0.01, -0.005, 0.005, 0.01] - shear_strains = shear_strains or [-0.06, -0.03, 0.03, 0.06] - self.undeformed_structure = structure self.deformations: list[Deformation] = [] self.def_structs: list[Structure] = [] diff --git a/pymatgen/analysis/nmr.py b/pymatgen/analysis/nmr.py index 118bbc8a685..b5a0c8cb31a 100644 --- a/pymatgen/analysis/nmr.py +++ b/pymatgen/analysis/nmr.py @@ -49,7 +49,7 @@ def __new__(cls, cs_matrix, vscale=None): or a 1x3 array of the primary sigma values corresponding to the principal axis system vscale (6x1 array-like): 6x1 array-like scaling the - Voigt-notation vector with the tensor entries + Voigt notation vector with the tensor entries """ t_array = np.array(cs_matrix) @@ -139,7 +139,7 @@ def __new__(cls, efg_matrix, vscale=None): or a 1x3 array of the primary values corresponding to the principal axis system vscale (6x1 array-like): 6x1 array-like scaling the - voigt-notation vector with the tensor entries + Voigt notation vector with the tensor entries """ t_array = np.array(efg_matrix) diff --git a/pymatgen/analysis/piezo.py b/pymatgen/analysis/piezo.py index 8929d7a83f0..a29d393abf4 100644 --- a/pymatgen/analysis/piezo.py +++ b/pymatgen/analysis/piezo.py @@ -18,7 +18,7 @@ class PiezoTensor(Tensor): - """This class describes the 3x6 piezo tensor in Voigt-notation.""" + """This class describes the 3x6 piezo tensor in Voigt notation.""" def __new__(cls, input_array, tol: float = 1e-3): """ diff --git a/pymatgen/core/tensors.py b/pymatgen/core/tensors.py index f54d5ac89c6..8ae564102f8 100644 --- a/pymatgen/core/tensors.py +++ b/pymatgen/core/tensors.py @@ -53,7 +53,7 @@ def __new__(cls, input_array, vscale=None, check_rank=None): input_array: (array-like with shape 3^N): array-like representing a tensor quantity in standard (i. e. non-voigt) notation vscale: (N x M array-like): a matrix corresponding - to the coefficients of the voigt-notation tensor + to the coefficients of the Voigt-notation tensor check_rank: (int): If not None, checks that input_array's rank == check_rank. Defaults to None. """ @@ -270,7 +270,7 @@ def symmetrized(self): @property def voigt_symmetrized(self): - """Returns a "voigt"-symmetrized tensor, i. e. a voigt-notation + """Returns a "voigt"-symmetrized tensor, i. e. a Voigt-notation tensor such that it is invariant w.r.t. permutation of indices. """ if not (self.rank % 2 == 0 and self.rank >= 2): @@ -635,7 +635,7 @@ def as_dict(self, voigt: bool = False) -> dict: Args: voigt (bool): flag for whether to store entries in - voigt-notation. Defaults to false, as information + Voigt notation. Defaults to false, as information may be lost in conversion. Returns (dict): @@ -668,11 +668,11 @@ class TensorCollection(collections.abc.Sequence, MSONable): or for having a tensor expansion. """ - def __init__(self, tensor_list, base_class=Tensor): + def __init__(self, tensor_list: Sequence, base_class=Tensor) -> None: """:param tensor_list: List of tensors. :param base_class: Class to be used. """ - self.tensors = [base_class(t) if not isinstance(t, base_class) else t for t in tensor_list] + self.tensors = [tensor if isinstance(tensor, base_class) else base_class(tensor) for tensor in tensor_list] def __len__(self): return len(self.tensors) @@ -848,7 +848,7 @@ def __new__(cls, input_array, vscale=None): input_array (3x3 array-like): the 3x3 array-like representing the content of the tensor vscale (6x1 array-like): 6x1 array-like scaling the - voigt-notation vector with the tensor entries + Voigt-notation vector with the tensor entries """ obj = super().__new__(cls, input_array, vscale, check_rank=2) return obj.view(cls)