Hello! I am wondering about the calculation of the piezoelectric modulus. For example, in RbSnF₃ (mp-998193), it is listed as 3.51. I wanted to reproduce this value using the piezoelectric tensor, so I wrote the following Python code:
import numpy as np
e = np.array([
[-3.188, 0.4841, -0.2276, 0.2722, -0.2063, 0.3177],
[0.3411, -2.762, -0.4208, -0.7191, 0.1107, 0.0267],
[-0.5021, -0.5808, -0.0881, -0.0464, 0.0141, 0.059]
])
def longitudinal_modulus(n, tensor):
"""
Compute longitudinal piezoelectric modulus
tensor: 3x6 Voigt
"""
n = n / np.linalg.norm(n)
sigma = np.array([
n[0]**2,
n[1]**2,
n[2]**2,
2*n[1]*n[2],
2*n[0]*n[2],
2*n[0]*n[1]
])
d_long = n[0]*np.dot(tensor[0,:], sigma) + \
n[1]*np.dot(tensor[1,:], sigma) + \
n[2]*np.dot(tensor[2,:], sigma)
return d_long
n_theta = 50 # divisions in theta
n_phi = 100 # divisions in phi
max_d = -np.inf
best_n = None
for i in range(n_theta):
theta = np.pi * i / (n_theta - 1) # polar angle
for j in range(n_phi):
phi = 2 * np.pi * j / n_phi # azimuthal angle
# convert to Cartesian unit vector
nx = np.sin(theta) * np.cos(phi)
ny = np.sin(theta) * np.sin(phi)
nz = np.cos(theta)
n_vec = np.array([nx, ny, nz])
d_val = longitudinal_modulus(n_vec, e)
if d_val > max_d:
max_d = d_val
best_n = n_vec
print("Maximum longitudinal modulus ||e_ij||_max =", max_d, "C/m²")
print("Direction vector for maximum:", best_n)
However, it gives the result of 3.27. Close, but not quite… Is there any error in the formula my code implements? Does the data used to calculate the piezoelectric modulus on MP differ from the data the piezoelectric tensor is based on? Thank you! 