import json
import math
import numpy as np
from scipy import optimize
__author__ = "Romain Gaillac and François-Xavier Coudert"
__version__ = "2025.08.26"
__license__ = "MIT"
[docs]
def dirVec(theta, phi):
"""Return a unit vector associated with angles theta and phi"""
return [math.sin(theta) * math.cos(phi), math.sin(theta) * math.sin(phi), math.cos(theta)]
[docs]
def dirVec1(theta, phi, chi):
"""Return the first unit vector associated with angles (theta, phi, chi)"""
return [math.sin(theta) * math.cos(phi), math.sin(theta) * math.sin(phi), math.cos(theta)]
[docs]
def dirVec2(theta, phi, chi):
"""Return the second unit vector associated with angles (theta, phi, chi)"""
return [
math.cos(theta) * math.cos(phi) * math.cos(chi) - math.sin(phi) * math.sin(chi),
math.cos(theta) * math.sin(phi) * math.cos(chi) + math.cos(phi) * math.sin(chi),
-math.sin(theta) * math.cos(chi),
]
# Functions to minimize/maximize
[docs]
def minimize(func, dim):
"""Find the global minimum of a function with 2 (theta, phi) or 3 (theta, phi, chi) parameters"""
if dim == 2:
# Mechanical properties are centrosymmetric so we reduce the range for phi
r = ((0, np.pi), (0, np.pi))
n = 25
elif dim == 3:
# Mechanical properties are centrosymmetric so we reduce the range for phi and chi
r = ((0, np.pi), (0, np.pi), (0, np.pi))
n = 10
# With good heuristics for our number of points, brute-force is not too bad
# because we are in low dimension.
return optimize.brute(func, r, Ns=n, full_output=True, finish=optimize.fmin)[0:2]
[docs]
def maximize(func, dim):
"""Find the global maximum of a function with 2 (theta, phi) or 3 (theta, phi, chi) parameters"""
res = minimize(lambda x: -func(x), dim)
return (res[0], -res[1])
[docs]
class Elastic:
"""
An elastic tensor, along with methods to access it.
It will raise a ValueError when the input is invalid, and a TypeError
if the input is valid but corresponds to a 2D material.
"""
def __init__(self, s):
"""Initialize the elastic tensor from a string"""
# Argument can be a 6-line string, a list of list, or a string representation of the list of list
# If the argument is JSON, decode it
if isinstance(s, str):
try:
s = json.loads(s)
except:
pass
if isinstance(s, np.ndarray) or isinstance(s, list):
# If we already have a list or numpy array, fine
mat = s
elif isinstance(s, str):
# Remove braces and pipes
s = s.replace("|", " ").replace("(", " ").replace(")", " ")
# Remove empty lines
lines = [line for line in s.split("\n") if line.strip()]
if len(lines) == 3:
raise TypeError("is a 2D material")
if len(lines) != 6:
raise ValueError("should have three or six rows")
# Convert to list of list of floats
try:
mat = [list(map(float, line.split())) for line in lines]
except:
raise ValueError("not all entries are numbers")
else:
raise ValueError("invalid argument as matrix")
# Make it into a square matrix
try:
mat = np.array(mat)
except:
# Is it upper triangular?
if list(map(len, mat)) == [6, 5, 4, 3, 2, 1]:
mat = [[0] * i + mat[i] for i in range(6)]
mat = np.array(mat)
if list(map(len, mat)) == [3, 2, 1]:
mat = [[0] * i + mat[i] for i in range(3)]
mat = np.array(mat)
# Is it lower triangular?
if list(map(len, mat)) == [1, 2, 3, 4, 5, 6]:
mat = [mat[i] + [0] * (5 - i) for i in range(6)]
mat = np.array(mat)
if list(map(len, mat)) == [1, 2, 3]:
mat = [mat[i] + [0] * (2 - i) for i in range(3)]
mat = np.array(mat)
if not isinstance(mat, np.ndarray):
raise ValueError("should be a square or triangular matrix")
if mat.shape == (3, 3):
raise TypeError("is a 2D material")
if mat.shape != (6, 6):
raise ValueError("should be a square or triangular matrix")
# Check that is is symmetric, or make it symmetric
if np.linalg.norm(np.tril(mat, -1)) == 0:
mat = mat + np.triu(mat, 1).transpose()
if np.linalg.norm(np.triu(mat, 1)) == 0:
mat = mat + np.tril(mat, -1).transpose()
if np.linalg.norm(mat - mat.transpose()) > 1e-3:
raise ValueError("should be symmetric, or triangular")
if np.linalg.norm(mat - mat.transpose()) > 0:
# It was almost symmetric: symmetrize it completely
mat = 0.5 * (mat + mat.transpose())
# Store it
self.CVoigt = mat
# Put it in a more useful representation
try:
self.SVoigt = np.linalg.inv(self.CVoigt)
except:
raise ValueError("matrix is singular")
VoigtMat = [[0, 5, 4], [5, 1, 3], [4, 3, 2]]
def SVoigtCoeff(p, q):
return 1.0 / ((1 + p // 3) * (1 + q // 3))
self.Smat = [
[
[
[
SVoigtCoeff(VoigtMat[i][j], VoigtMat[k][l]) * self.SVoigt[VoigtMat[i][j]][VoigtMat[k][l]]
for i in range(3)
]
for j in range(3)
]
for k in range(3)
]
for l in range(3)
]
[docs]
def is2D(self):
"""Return True if the material is 2D."""
return False
[docs]
def isOrthorhombic(self):
"""Return True if the elastic tensor is orthorhombic."""
def iszero(x):
return abs(x) < 1.0e-3
return (
iszero(self.CVoigt[0][3])
and iszero(self.CVoigt[0][4])
and iszero(self.CVoigt[0][5])
and iszero(self.CVoigt[1][3])
and iszero(self.CVoigt[1][4])
and iszero(self.CVoigt[1][5])
and iszero(self.CVoigt[2][3])
and iszero(self.CVoigt[2][4])
and iszero(self.CVoigt[2][5])
and iszero(self.CVoigt[3][4])
and iszero(self.CVoigt[3][5])
and iszero(self.CVoigt[4][5])
)
[docs]
def isCubic(self):
"""Return True if the elastic tensor is cubic."""
def iszero(x):
return abs(x) < 1.0e-3
return (
iszero(self.CVoigt[0][3])
and iszero(self.CVoigt[0][4])
and iszero(self.CVoigt[0][5])
and iszero(self.CVoigt[1][3])
and iszero(self.CVoigt[1][4])
and iszero(self.CVoigt[1][5])
and iszero(self.CVoigt[2][3])
and iszero(self.CVoigt[2][4])
and iszero(self.CVoigt[2][5])
and iszero(self.CVoigt[3][4])
and iszero(self.CVoigt[3][5])
and iszero(self.CVoigt[4][5])
and iszero(self.CVoigt[0][0] - self.CVoigt[1][1])
and iszero(self.CVoigt[0][0] - self.CVoigt[2][2])
and iszero(self.CVoigt[0][0] - self.CVoigt[1][1])
and iszero(self.CVoigt[0][0] - self.CVoigt[2][2])
and iszero(self.CVoigt[3][3] - self.CVoigt[4][4])
and iszero(self.CVoigt[3][3] - self.CVoigt[5][5])
and iszero(self.CVoigt[0][1] - self.CVoigt[0][2])
and iszero(self.CVoigt[0][1] - self.CVoigt[1][2])
)
[docs]
def Young(self, x):
"""Return the Young's modulus in the direction x=(theta, phi)."""
a = dirVec(x[0], x[1])
r = sum(
[
a[i] * a[j] * a[k] * a[l] * self.Smat[i][j][k][l]
for i in range(3)
for j in range(3)
for k in range(3)
for l in range(3)
]
)
return 1 / r
[docs]
def Young_2(self, x, y):
"""Return the Young's modulus in the direction x, y."""
a = dirVec(x, y)
r = sum(
[
a[i] * a[j] * a[k] * a[l] * self.Smat[i][j][k][l]
for i in range(3)
for j in range(3)
for k in range(3)
for l in range(3)
]
)
return 1 / r
[docs]
def LC(self, x):
"""Return the linear compressibility in the direction x=(theta, phi)."""
a = dirVec(x[0], x[1])
r = sum([a[i] * a[j] * self.Smat[i][j][k][k] for i in range(3) for j in range(3) for k in range(3)])
return 1000 * r
[docs]
def LC_2(self, x, y):
"""Return the linear compressibility in the direction x, y."""
a = dirVec(x, y)
r = sum([a[i] * a[j] * self.Smat[i][j][k][k] for i in range(3) for j in range(3) for k in range(3)])
return 1000 * r
[docs]
def shear(self, x):
"""Return the shear modulus in the direction x=(theta, phi, chi)."""
a = dirVec(x[0], x[1])
b = dirVec2(x[0], x[1], x[2])
r = sum(
[
a[i] * b[j] * a[k] * b[l] * self.Smat[i][j][k][l]
for i in range(3)
for j in range(3)
for k in range(3)
for l in range(3)
]
)
return 1 / (4 * r)
[docs]
def Poisson(self, x):
"""Return the Poisson's ratio in the direction x=(theta, phi, chi)."""
a = dirVec(x[0], x[1])
b = dirVec2(x[0], x[1], x[2])
r1 = sum(
[
a[i] * a[j] * b[k] * b[l] * self.Smat[i][j][k][l]
for i in range(3)
for j in range(3)
for k in range(3)
for l in range(3)
]
)
r2 = sum(
[
a[i] * a[j] * a[k] * a[l] * self.Smat[i][j][k][l]
for i in range(3)
for j in range(3)
for k in range(3)
for l in range(3)
]
)
return -r1 / r2
[docs]
def averages(self):
"""Return the Voigt, Reuss, and Hill averages for the elastic tensor."""
A = (self.CVoigt[0][0] + self.CVoigt[1][1] + self.CVoigt[2][2]) / 3
B = (self.CVoigt[1][2] + self.CVoigt[0][2] + self.CVoigt[0][1]) / 3
C = (self.CVoigt[3][3] + self.CVoigt[4][4] + self.CVoigt[5][5]) / 3
a = (self.SVoigt[0][0] + self.SVoigt[1][1] + self.SVoigt[2][2]) / 3
b = (self.SVoigt[1][2] + self.SVoigt[0][2] + self.SVoigt[0][1]) / 3
c = (self.SVoigt[3][3] + self.SVoigt[4][4] + self.SVoigt[5][5]) / 3
KV = (A + 2 * B) / 3
GV = (A - B + 3 * C) / 5
KR = 1 / (3 * a + 6 * b)
GR = 5 / (4 * a - 4 * b + 3 * c)
KH = (KV + KR) / 2
GH = (GV + GR) / 2
return [
[KV, 1 / (1 / (3 * GV) + 1 / (9 * KV)), GV, (1 - 3 * GV / (3 * KV + GV)) / 2],
[KR, 1 / (1 / (3 * GR) + 1 / (9 * KR)), GR, (1 - 3 * GR / (3 * KR + GR)) / 2],
[KH, 1 / (1 / (3 * GH) + 1 / (9 * KH)), GH, (1 - 3 * GH / (3 * KH + GH)) / 2],
]
[docs]
def eigenvalues(self):
"""Return the eigenvalues of the stiffness matrix."""
return np.sort(np.linalg.eig(self.CVoigt)[0])
[docs]
def shear2D(self, x):
"""Find the minimum and maximum shear modulus for a given direction x=(theta, phi)."""
ftol = 0.001
xtol = 0.01
def func1(z):
return self.shear([x[0], x[1], z[0]])
r1 = optimize.minimize(func1, np.pi / 2.0, args=(), method="Powell", options={"xtol": xtol, "ftol": ftol})
def func2(z):
return -self.shear([x[0], x[1], z[0]])
r2 = optimize.minimize(func2, np.pi / 2.0, args=(), method="Powell", options={"xtol": xtol, "ftol": ftol})
return (float(r1.fun), -float(r2.fun))
[docs]
def shear3D(self, x, y, guess1=np.pi / 2.0, guess2=np.pi / 2.0):
"""Find the minimum and maximum shear modulus for a given direction (x, y)."""
tol = 0.0005
def func1(z):
return self.shear([x, y, z[0]])
r1 = optimize.minimize(func1, guess1, args=(), method="L-BFGS-B")
def func2(z):
return -self.shear([x, y, z[0]])
r2 = optimize.minimize(func2, guess2, args=(), method="L-BFGS-B")
return (float(r1.fun), -float(r2.fun), float(r1.x[0]), float(r2.x[0]))
[docs]
def Poisson2D(self, x):
"""Find the minimum and maximum Poisson's ratio for a given direction x=(theta, phi)."""
ftol = 0.001
xtol = 0.01
def func1(z):
return self.Poisson([x[0], x[1], z[0]])
r1 = optimize.minimize(func1, np.pi / 2.0, args=(), method="Powell", options={"xtol": xtol, "ftol": ftol})
def func2(z):
return -self.Poisson([x[0], x[1], z[0]])
r2 = optimize.minimize(func2, np.pi / 2.0, args=(), method="Powell", options={"xtol": xtol, "ftol": ftol})
return (min(0, float(r1.fun)), max(0, float(r1.fun)), -float(r2.fun))
[docs]
def Poisson3D(self, x, y, guess1=np.pi / 2.0, guess2=np.pi / 2.0):
"""Find the minimum and maximum Poisson's ratio for a given direction (x, y)."""
tol = 0.005
def func1(z):
return self.Poisson([x, y, z[0]])
r1 = optimize.minimize(func1, guess1, args=(), method="L-BFGS-B")
def func2(z):
return -self.Poisson([x, y, z[0]])
r2 = optimize.minimize(func2, guess2, args=(), method="L-BFGS-B")
return (min(0, float(r1.fun)), max(0, float(r1.fun)), -float(r2.fun), float(r1.x[0]), float(r2.x[0]))
[docs]
class ElasticOrtho(Elastic):
"""An elastic tensor, for the specific case of an orthorhombic system"""
def __init__(self, arg):
"""Initialize from a matrix, or from an Elastic object"""
if isinstance(arg, str):
Elastic.__init__(self, arg)
elif isinstance(arg, Elastic):
self.CVoigt = arg.CVoigt
self.SVoigt = arg.SVoigt
self.Smat = arg.Smat
else:
raise TypeError("ElasticOrtho constructor argument should be string or Elastic object")
[docs]
def Young(self, x):
"""Return the Young's modulus in the direction x=(theta, phi)."""
ct2 = math.cos(x[0]) ** 2
st2 = 1 - ct2
cf2 = math.cos(x[1]) ** 2
sf2 = 1 - cf2
s11 = self.Smat[0][0][0][0]
s22 = self.Smat[1][1][1][1]
s33 = self.Smat[2][2][2][2]
s44 = 4 * self.Smat[1][2][1][2]
s55 = 4 * self.Smat[0][2][0][2]
s66 = 4 * self.Smat[0][1][0][1]
s12 = self.Smat[0][0][1][1]
s13 = self.Smat[0][0][2][2]
s23 = self.Smat[1][1][2][2]
return 1 / (
ct2**2 * s33
+ 2 * cf2 * ct2 * s13 * st2
+ cf2 * ct2 * s55 * st2
+ 2 * ct2 * s23 * sf2 * st2
+ ct2 * s44 * sf2 * st2
+ cf2**2 * s11 * st2**2
+ 2 * cf2 * s12 * sf2 * st2**2
+ cf2 * s66 * sf2 * st2**2
+ s22 * sf2**2 * st2**2
)
[docs]
def LC(self, x):
"""Return the linear compressibility in the direction x=(theta, phi)."""
ct2 = math.cos(x[0]) ** 2
cf2 = math.cos(x[1]) ** 2
s11 = self.Smat[0][0][0][0]
s22 = self.Smat[1][1][1][1]
s33 = self.Smat[2][2][2][2]
s12 = self.Smat[0][0][1][1]
s13 = self.Smat[0][0][2][2]
s23 = self.Smat[1][1][2][2]
return 1000 * (ct2 * (s13 + s23 + s33) + (cf2 * (s11 + s12 + s13) + (s12 + s22 + s23) * (1 - cf2)) * (1 - ct2))
[docs]
def shear(self, x):
"""Return the shear modulus in the direction x=(theta, phi, chi)."""
ct = math.cos(x[0])
ct2 = ct * ct
st2 = 1 - ct2
cf = math.cos(x[1])
sf = math.sin(x[1])
sf2 = sf * sf
cx = math.cos(x[2])
cx2 = cx * cx
sx = math.sin(x[2])
sx2 = 1 - cx2
s11 = self.Smat[0][0][0][0]
s22 = self.Smat[1][1][1][1]
s33 = self.Smat[2][2][2][2]
s44 = 4 * self.Smat[1][2][1][2]
s55 = 4 * self.Smat[0][2][0][2]
s66 = 4 * self.Smat[0][1][0][1]
s12 = self.Smat[0][0][1][1]
s13 = self.Smat[0][0][2][2]
s23 = self.Smat[1][1][2][2]
r = (
ct2 * ct2 * cx2 * s44 * sf2
+ cx2 * s44 * sf2 * st2 * st2
+ 4 * cf**3 * ct * cx * (-2 * s11 + 2 * s12 + s66) * sf * st2 * sx
+ 2
* cf
* ct
* cx
* sf
* (
ct2 * (s44 - s55)
+ (4 * s13 - 4 * s23 - s44 + s55 - 4 * s12 * sf2 + 4 * s22 * sf2 - 2 * s66 * sf2) * st2
)
* sx
+ s66 * sf2 * sf2 * st2 * sx2
+ cf**4 * st2 * (4 * ct2 * cx2 * s11 + s66 * sx2)
+ ct2 * (2 * cx2 * (2 * s33 + sf2 * (-4 * s23 - s44 + 2 * s22 * sf2)) * st2 + s55 * sf2 * sx2)
+ cf**2
* (
ct2 * ct2 * cx2 * s55
+ ct2 * (-2 * cx2 * (4 * s13 + s55 - 2 * (2 * s12 + s66) * sf2) * st2 + s44 * sx2)
+ st2 * (cx2 * s55 * st2 + 2 * (2 * s11 - 4 * s12 + 2 * s22 - s66) * sf2 * sx2)
)
)
return 1 / r
[docs]
def Poisson(self, x):
"""Return the Poisson's ratio in the direction x=(theta, phi, chi)."""
ct = math.cos(x[0])
ct2 = ct * ct
st2 = 1 - ct2
cf = math.cos(x[1])
sf = math.sin(x[1])
cx = math.cos(x[2])
sx = math.sin(x[2])
s11 = self.Smat[0][0][0][0]
s22 = self.Smat[1][1][1][1]
s33 = self.Smat[2][2][2][2]
s44 = 4 * self.Smat[1][2][1][2]
s55 = 4 * self.Smat[0][2][0][2]
s66 = 4 * self.Smat[0][1][0][1]
s12 = self.Smat[0][0][1][1]
s13 = self.Smat[0][0][2][2]
s23 = self.Smat[1][1][2][2]
return (
-(ct**2 * cx**2 * s33 * st2)
- cf**2 * cx**2 * s13 * st2 * st2
- cx**2 * s23 * sf**2 * st2 * st2
+ ct * cx * s44 * sf * st2 * (ct * cx * sf + cf * sx)
- ct**2 * s23 * (ct * cx * sf + cf * sx) ** 2
- cf**2 * s12 * st2 * (ct * cx * sf + cf * sx) ** 2
- s22 * sf**2 * st2 * (ct * cx * sf + cf * sx) ** 2
+ cf * ct * cx * s55 * st2 * (cf * ct * cx - sf * sx)
- cf * s66 * sf * st2 * (ct * cx * sf + cf * sx) * (cf * ct * cx - sf * sx)
- ct**2 * s13 * (cf * ct * cx - sf * sx) ** 2
- cf**2 * s11 * st2 * (cf * ct * cx - sf * sx) ** 2
- s12 * sf**2 * st2 * (cf * ct * cx - sf * sx) ** 2
) / (
ct**4 * s33
+ 2 * cf**2 * ct**2 * s13 * st2
+ cf**2 * ct**2 * s55 * st2
+ 2 * ct**2 * s23 * sf**2 * st2
+ ct**2 * s44 * sf**2 * st2
+ cf**4 * s11 * st2 * st2
+ 2 * cf**2 * s12 * sf**2 * st2 * st2
+ cf**2 * s66 * sf**2 * st2 * st2
+ s22 * sf**4 * st2 * st2
)
[docs]
class Elastic2D:
"""Elastic tensor for 2D material, along with methods to access it"""
def __init__(self, s):
"""Initialize the elastic tensor from a string"""
# Argument can be a 3-line string, a list of list, or a string representation of the list of list
try:
if isinstance(json.loads(s), list):
s = json.loads(s)
except:
pass
if isinstance(s, str):
# Remove braces and pipes
s = s.replace("|", " ").replace("(", " ").replace(")", " ")
# Remove empty lines
lines = [line for line in s.split("\n") if line.strip()]
if len(lines) != 3:
raise ValueError("should have three rows")
# Convert to float
try:
mat = [list(map(float, line.split())) for line in lines]
except:
raise ValueError("not all entries are numbers")
elif isinstance(s, list):
# If we already have a list, simply use it
mat = s
else:
raise ValueError("invalid argument as matrix")
# Make it into a square matrix
try:
mat = np.array(mat)
except:
# Is it upper triangular?
if list(map(len, mat)) == [3, 2, 1]:
mat = [[0] * i + mat[i] for i in range(3)]
mat = np.array(mat)
# Is it lower triangular?
if list(map(len, mat)) == [1, 2, 3]:
mat = [mat[i] + [0] * (2 - i) for i in range(3)]
mat = np.array(mat)
if not isinstance(mat, np.ndarray):
raise ValueError("should be a square or triangular matrix")
if mat.shape != (3, 3):
raise ValueError("should be a square or triangular matrix")
# Check that is is symmetric, or make it symmetric
if np.linalg.norm(np.tril(mat, -1)) == 0:
mat = mat + np.triu(mat, 1).transpose()
if np.linalg.norm(np.triu(mat, 1)) == 0:
mat = mat + np.tril(mat, -1).transpose()
if np.linalg.norm(mat - mat.transpose()) > 1e-3:
raise ValueError("should be symmetric, or triangular")
if np.linalg.norm(mat - mat.transpose()) > 0:
mat = 0.5 * (mat + mat.transpose())
# Store it
self.CVoigt = mat
# Put it in a more useful representation
try:
self.SVoigt = np.linalg.inv(self.CVoigt)
except:
raise ValueError("matrix is singular")
# Store the important components
self.s11 = self.SVoigt[0][0]
self.s12 = self.SVoigt[0][1]
self.s16 = self.SVoigt[0][2]
self.s22 = self.SVoigt[1][1]
self.s26 = self.SVoigt[1][2]
self.s66 = self.SVoigt[2][2]
[docs]
def is2D(self):
"""Return True if the material is 2D."""
return True
[docs]
def Young(self, theta):
"""Return the Young's modulus in the direction theta."""
ct = math.cos(theta)
st = math.sin(theta)
return 1 / (
self.s11 * ct**4
+ self.s22 * st**4
+ 2 * self.s16 * ct**3 * st
+ 2 * self.s26 * ct * st**3
+ (2 * self.s12 + self.s66) * ct**2 * st**2
)
[docs]
def shear(self, theta):
"""Return the shear modulus in the direction theta."""
ct = math.cos(theta)
st = math.sin(theta)
calc = (
(self.s11 + self.s22 - 2 * self.s12) * ct**2 * st**2
+ self.s66 / 4 * (ct**4 + st**4 - 2 * st**2 * ct**2)
+ self.s16 * (st**3 * ct - ct**3 * st)
+ self.s26 * (ct**3 * st - st**3 * ct)
)
return 1 / (4 * calc)
[docs]
def Poisson(self, theta):
"""Return the Poisson's ratio in the direction theta."""
ct = math.cos(theta)
st = math.sin(theta)
num = (
(self.s11 + self.s22 - self.s66) * ct**2 * st**2
+ self.s12 * (ct**4 + st**4)
+ self.s16 * (st**3 * ct - ct**3 * st)
+ self.s26 * (ct**3 * st - st**3 * ct)
)
denom = (
(2 * self.s12 + self.s66) * ct**2 * st**2
+ self.s11 * ct**4
+ self.s22 * st**4
+ 2 * self.s16 * st**3 * ct
+ 2 * self.s26 * ct**3 * st
)
return -num / denom
[docs]
def eigenvalues(self):
"""Return the eigenvalues of the stiffness matrix."""
return np.sort(np.linalg.eig(self.CVoigt)[0])