Skip to content

Commit

Permalink
add units to ElasticTensor methods
Browse files Browse the repository at this point in the history
fix typos like lower case voigt
  • Loading branch information
janosh committed Nov 16, 2023
1 parent bc5edc5 commit 1e7f93a
Show file tree
Hide file tree
Showing 5 changed files with 72 additions and 66 deletions.
91 changes: 46 additions & 45 deletions pymatgen/analysis/elasticity/elastic.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -23,6 +23,8 @@
from pymatgen.core.units import Unit

if TYPE_CHECKING:
from collections.abc import Sequence

from pymatgen.core import Structure


Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
"""
Expand All @@ -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()
Expand All @@ -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
Expand Down Expand Up @@ -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."
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -547,22 +545,22 @@ 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.
"""
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):
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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"]
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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.
Expand Down
29 changes: 17 additions & 12 deletions pymatgen/analysis/elasticity/strain.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand All @@ -101,22 +103,25 @@ 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.
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] = []
Expand Down
4 changes: 2 additions & 2 deletions pymatgen/analysis/nmr.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down
2 changes: 1 addition & 1 deletion pymatgen/analysis/piezo.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down
Loading

0 comments on commit 1e7f93a

Please sign in to comment.