"""
Piezo sensitivity analysis module.
"""
from pymatgen.core.tensors import Tensor
import pymatgen.io.phonopy
import numpy as np
import warnings
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer as sga
from monty.dev import requires
try:
from phonopy import Phonopy
from phonopy.harmonic import dynmat_to_fc as dyntofc
except ImportError:
Phonopy = None
__author__ = "Handong Ling"
__copyright__ = "Copyright 2019, The Materials Project"
__version__ = "1.0"
__maintainer__ = "Handong Ling"
__email__ = "hling@lbl.gov"
__status__ = "Development"
__date__ = "Feb, 2019"
[docs]class BornEffectiveCharge:
"""
This class describes the Nx3x3 born effective charge tensor
"""
def __init__(self, structure, bec, pointops, tol=1e-3):
"""
Create an BornEffectiveChargeTensor object defined by a
structure, point operations of the structure's atomic sites.
Note that the constructor uses __new__ rather than __init__
according to the standard method ofsubclassing numpy ndarrays.
Args:
input_matrix (Nx3x3 array-like): the Nx3x3 array-like
representing the born effective charge tensor
"""
self.structure = structure
self.bec = bec
self.pointops = pointops
self.BEC_operations = None
if not np.sum(self.bec) < tol:
warnings.warn(
"Input born effective charge tensor does "
"not satisfy charge neutrality"
)
[docs] def get_BEC_operations(self, eigtol=1e-05, opstol=1e-03):
"""
Returns the symmetry operations which maps the tensors
belonging to equivalent sites onto each other in the form
[site index 1, site index 2, [Symmops mapping from site
index 1 to site index 2]]
Args:
eigtol (float): tolerance for determining if two sites are
related by symmetry
opstol (float): tolerance for determining if a symmetry
operation relates two sites
Return:
list of symmetry operations mapping equivalent sites and
the indexes of those sites.
"""
bec = self.bec
struc = self.structure
ops = sga(struc).get_symmetry_operations(cartesian=True)
uniquepointops = []
for op in ops:
uniquepointops.append(op)
for atom in range(len(self.pointops)):
for op in self.pointops[atom]:
if op not in uniquepointops:
uniquepointops.append(op)
passed = []
relations = []
for site in range(len(bec)):
unique = 1
eig1, vecs1 = np.linalg.eig(bec[site])
index = np.argsort(eig1)
neweig = np.real([eig1[index[0]], eig1[index[1]], eig1[index[2]]])
for index in range(len(passed)):
if np.allclose(neweig, passed[index][1], atol=eigtol):
relations.append([site, index])
unique = 0
passed.append([site, passed[index][0], neweig])
break
if unique == 1:
relations.append([site, site])
passed.append([site, neweig])
BEC_operations = []
for atom in range(len(relations)):
BEC_operations.append(relations[atom])
BEC_operations[atom].append([])
for op in uniquepointops:
new = op.transform_tensor(self.bec[relations[atom][1]])
# Check the matrix it references
if np.allclose(new, self.bec[relations[atom][0]], atol=opstol):
BEC_operations[atom][2].append(op)
self.BEC_operations = BEC_operations
[docs] def get_rand_BEC(self, max_charge=1):
"""
Generate a random born effective charge tensor which obeys a structure's
symmetry and the acoustic sum rule
Args:
max_charge (float): maximum born effective charge value
Return:
np.array Born effective charge tensor
"""
struc = self.structure
symstruc = sga(struc)
symstruc = symstruc.get_symmetrized_structure()
l = len(struc)
BEC = np.zeros((l, 3, 3))
for atom in range(len(self.BEC_operations)):
if self.BEC_operations[atom][0] == self.BEC_operations[atom][1]:
temp_tensor = Tensor(np.random.rand(3, 3) - 0.5)
temp_tensor = sum(
[temp_tensor.transform(symm_op) for symm_op in self.pointops[atom]]
) / len(self.pointops[atom])
BEC[atom] = temp_tensor
else:
tempfcm = np.zeros([3, 3])
for op in self.BEC_operations[atom][2]:
tempfcm += op.transform_tensor(BEC[self.BEC_operations[atom][1]])
BEC[self.BEC_operations[atom][0]] = tempfcm
if len(self.BEC_operations[atom][2]) != 0:
BEC[self.BEC_operations[atom][0]] = BEC[
self.BEC_operations[atom][0]
] / len(self.BEC_operations[atom][2])
# Enforce Acoustic Sum
disp_charge = np.einsum("ijk->jk", BEC) / l
add = np.zeros([l, 3, 3])
for atom in range(len(self.BEC_operations)):
if self.BEC_operations[atom][0] == self.BEC_operations[atom][1]:
temp_tensor = Tensor(disp_charge)
temp_tensor = sum(
[temp_tensor.transform(symm_op) for symm_op in self.pointops[atom]]
) / len(self.pointops[atom])
add[self.BEC_operations[atom][0]] = temp_tensor
else:
temp_tensor = np.zeros([3, 3])
for op in self.BEC_operations[atom][2]:
temp_tensor += op.transform_tensor(
add[self.BEC_operations[atom][1]]
)
add[self.BEC_operations[atom][0]] = temp_tensor
if len(self.BEC_operations[atom]) != 0:
add[self.BEC_operations[atom][0]] = add[
self.BEC_operations[atom][0]
] / len(self.BEC_operations[atom][2])
BEC = BEC - add
return BEC * max_charge
[docs]class InternalStrainTensor:
"""
This class describes the Nx3x3x3 internal tensor defined by a
structure, point operations of the structure's atomic sites.
"""
def __init__(self, structure, ist, pointops, tol=1e-3):
"""
Create an InternalStrainTensor object.
Args:
input_matrix (Nx3x3x3 array-like): the Nx3x3x3 array-like
representing the internal strain tensor
"""
self.structure = structure
self.ist = ist
self.pointops = pointops
self.IST_operations = None
obj = self.ist
if not (obj - np.transpose(obj, (0, 1, 3, 2)) < tol).all():
warnings.warn(
"Input internal strain tensor does " "not satisfy standard symmetries"
)
[docs] def get_IST_operations(self, opstol=1e-03):
"""
Returns the symmetry operations which maps the tensors
belonging to equivalent sites onto each other in the form
[site index 1, site index 2, [Symmops mapping from site
index 1 to site index 2]]
Args:
opstol (float): tolerance for determining if a symmetry
operation relates two sites
Return:
list of symmetry operations mapping equivalent sites and
the indexes of those sites.
"""
struc = self.structure
ops = sga(struc).get_symmetry_operations(cartesian=True)
uniquepointops = []
for op in ops:
uniquepointops.append(op)
for atom in range(len(self.pointops)):
for op in self.pointops[atom]:
if op not in uniquepointops:
uniquepointops.append(op)
IST_operations = []
for atom in range(len(self.ist)):
IST_operations.append([])
for j in range(0, atom):
for op in uniquepointops:
new = op.transform_tensor(self.ist[j])
# Check the matrix it references
if np.allclose(new, self.ist[atom], atol=opstol):
IST_operations[atom].append([j, op])
self.IST_operations = IST_operations
[docs] def get_rand_IST(self, max_force=1):
"""
Generate a random internal strain tensor which obeys a structure's
symmetry and the acoustic sum rule
Args:
max_force(float): maximum born effective charge value
Return:
InternalStrainTensor object
"""
l = len(self.structure)
IST = np.zeros((l, 3, 3, 3))
for atom in range(len(self.IST_operations)):
temp_tensor = np.zeros([3, 3, 3])
for op in self.IST_operations[atom]:
temp_tensor += op[1].transform_tensor(IST[op[0]])
if len(self.IST_operations[atom]) == 0:
temp_tensor = Tensor(np.random.rand(3, 3, 3) - 0.5)
for dim in range(3):
temp_tensor[dim] = (temp_tensor[dim] + temp_tensor[dim].T) / 2
temp_tensor = sum(
[temp_tensor.transform(symm_op) for symm_op in self.pointops[atom]]
) / len(self.pointops[atom])
IST[atom] = temp_tensor
if len(self.IST_operations[atom]) != 0:
IST[atom] = IST[atom] / len(self.IST_operations[atom])
return IST * max_force
[docs]class ForceConstantMatrix:
"""
This class describes the NxNx3x3 force constant matrix defined by a
structure, point operations of the structure's atomic sites, and the
shared symmetry operations between pairs of atomic sites.
"""
def __init__(self, structure, fcm, pointops, sharedops, tol=1e-3):
"""
Create an ForceConstantMatrix object.
Args:
input_matrix (NxNx3x3 array-like): the NxNx3x3 array-like
representing the force constant matrix
"""
self.structure = structure
self.fcm = fcm
self.pointops = pointops
self.sharedops = sharedops
self.FCM_operations = None
[docs] def get_FCM_operations(self, eigtol=1e-05, opstol=1e-05):
"""
Returns the symmetry operations which maps the tensors
belonging to equivalent sites onto each other in the form
[site index 1a, site index 1b, site index 2a, site index 2b,
[Symmops mapping from site index 1a, 1b to site index 2a, 2b]]
Args:
eigtol (float): tolerance for determining if two sites are
related by symmetry
opstol (float): tolerance for determining if a symmetry
operation relates two sites
Return:
list of symmetry operations mapping equivalent sites and
the indexes of those sites.
"""
struc = self.structure
ops = sga(struc).get_symmetry_operations(cartesian=True)
uniquepointops = []
for op in ops:
uniquepointops.append(op)
for atom in range(len(self.pointops)):
for op in self.pointops[atom]:
if op not in uniquepointops:
uniquepointops.append(op)
passed = []
relations = []
for atom1 in range(len(self.fcm)):
for atom2 in range(atom1, len(self.fcm)):
unique = 1
eig1, vecs1 = np.linalg.eig(self.fcm[atom1][atom2])
index = np.argsort(eig1)
neweig = np.real([eig1[index[0]], eig1[index[1]], eig1[index[2]]])
for entry in range(len(passed)):
if np.allclose(neweig, passed[entry][2], atol=eigtol):
relations.append(
[atom1, atom2, passed[entry][0], passed[entry][1]]
)
unique = 0
break
if unique == 1:
relations.append([atom1, atom2, atom2, atom1])
passed.append([atom1, atom2, np.real(neweig)])
FCM_operations = []
for entry in range(len(relations)):
good = 0
FCM_operations.append(relations[entry])
FCM_operations[entry].append([])
good = 0
for op in uniquepointops:
new = op.transform_tensor(
self.fcm[relations[entry][2]][relations[entry][3]]
)
if np.allclose(
new, self.fcm[relations[entry][0]][relations[entry][1]], atol=opstol
):
FCM_operations[entry][4].append(op)
good = 1
if (
relations[entry][0] == relations[entry][3]
and relations[entry][1] == relations[entry][2]
):
good = 1
if (
relations[entry][0] == relations[entry][2]
and relations[entry][1] == relations[entry][3]
):
good = 1
if good == 0:
FCM_operations[entry] = [
relations[entry][0],
relations[entry][1],
relations[entry][3],
relations[entry][2],
]
FCM_operations[entry].append([])
for op in uniquepointops:
new = op.transform_tensor(
self.fcm[relations[entry][2]][relations[entry][3]]
)
if np.allclose(
new.T,
self.fcm[relations[entry][0]][relations[entry][1]],
atol=opstol,
):
FCM_operations[entry][4].append(op)
self.FCM_operations = FCM_operations
return FCM_operations
[docs] def get_unstable_FCM(self, max_force=1):
"""
Generate an unsymmeterized force constant matrix
Args:
max_charge (float): maximum born effective charge value
Return:
numpy array representing the force constant matrix
"""
struc = self.structure
operations = self.FCM_operations
# set max force in reciprocal space
numsites = len(struc.sites)
D = (1 / max_force) * 2 * (np.ones([numsites * 3, numsites * 3]))
for op in operations:
same = 0
transpose = 0
if op[0] == op[1] and op[0] == op[2] and op[0] == op[3]:
same = 1
if op[0] == op[3] and op[1] == op[2]:
transpose = 1
if transpose == 0 and same == 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = np.zeros(
[3, 3]
)
D[3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3] = np.zeros(
[3, 3]
)
for symop in op[4]:
tempfcm = D[3 * op[2]:3 * op[2] + 3, 3 * op[3]:3 * op[3] + 3]
tempfcm = symop.transform_tensor(tempfcm)
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] += tempfcm
if len(op[4]) != 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] / len(op[4])
D[3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
].T
continue
else:
temp_tensor = Tensor(np.random.rand(3, 3) - 0.5) * max_force
temp_tensor_sum = sum(
[
temp_tensor.transform(symm_op)
for symm_op in self.sharedops[op[0]][op[1]]
]
)
temp_tensor_sum = temp_tensor_sum / (len(self.sharedops[op[0]][op[1]]))
if op[0] != op[1]:
for pair in range(len(op[4])):
temp_tensor2 = temp_tensor_sum.T
temp_tensor2 = op[4][pair].transform_tensor(temp_tensor2)
temp_tensor_sum = (temp_tensor_sum + temp_tensor2) / 2
else:
temp_tensor_sum = (temp_tensor_sum + temp_tensor_sum.T) / 2
D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] = temp_tensor_sum
D[
3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3
] = temp_tensor_sum.T
return D
[docs] def get_symmetrized_FCM(self, unsymmetrized_fcm, max_force=1):
"""
Generate a symmeterized force constant matrix from an unsymmeterized matrix
Args:
unsymmetrized_fcm (numpy array): unsymmeterized force constant matrix
max_charge (float): maximum born effective charge value
Return:
3Nx3N numpy array representing the force constant matrix
"""
operations = self.FCM_operations
D = unsymmetrized_fcm
for op in operations:
same = 0
transpose = 0
if op[0] == op[1] and op[0] == operations[2] and op[0] == op[3]:
same = 1
if op[0] == op[3] and op[1] == op[2]:
transpose = 1
if transpose == 0 and same == 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = np.zeros(
[3, 3]
)
for symop in op[4]:
tempfcm = D[3 * op[2]:3 * op[2] + 3, 3 * op[3]:3 * op[3] + 3]
tempfcm = symop.transform_tensor(tempfcm)
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] += tempfcm
if len(op[4]) != 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] / len(op[4])
D[3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
].T
continue
else:
temp_tensor = Tensor(
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3]
)
temp_tensor_sum = sum(
[
temp_tensor.transform(symm_op)
for symm_op in self.sharedops[op[0]][op[1]]
]
)
if len(self.sharedops[op[0]][op[1]]) != 0:
temp_tensor_sum = temp_tensor_sum / (
len(self.sharedops[op[0]][op[1]])
)
# Apply the proper transformation if there is an equivalent already
if op[0] != op[1]:
for pair in range(len(op[4])):
temp_tensor2 = temp_tensor_sum.T
temp_tensor2 = op[4][pair].transform_tensor(temp_tensor2)
temp_tensor_sum = (temp_tensor_sum + temp_tensor2) / 2
else:
temp_tensor_sum = (temp_tensor_sum + temp_tensor_sum.T) / 2
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = temp_tensor_sum
D[3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3] = temp_tensor_sum.T
return D
[docs] def get_stable_FCM(self, fcm, fcmasum=10):
"""
Generate a symmeterized force constant matrix that obeys the objects symmetry
constraints, has no unstable modes and also obeys the acoustic sum rule through an
iterative procedure
Args:
fcm (numpy array): unsymmeterized force constant matrix
fcmasum (int): number of iterations to attempt to obey the acoustic sum
rule
Return:
3Nx3N numpy array representing the force constant matrix
"""
check = 0
count = 0
while check == 0:
# if resymmetrizing brings back unstable modes 20 times, the method breaks
if count > 20:
check = 1
break
eigs, vecs = np.linalg.eig(fcm)
maxeig = np.max(-1 * eigs)
eigsort = np.argsort(np.abs(eigs))
for i in range(3, len(eigs)):
if eigs[eigsort[i]] > 1e-06:
eigs[eigsort[i]] = -1 * maxeig * np.random.rand()
diag = np.real(np.eye(len(fcm)) * eigs)
fcm = np.real(np.matmul(np.matmul(vecs, diag), vecs.T))
fcm = self.get_symmetrized_FCM(fcm)
fcm = self.get_asum_FCM(fcm)
eigs, vecs = np.linalg.eig(fcm)
unstable_modes = 0
eigsort = np.argsort(np.abs(eigs))
for i in range(3, len(eigs)):
if eigs[eigsort[i]] > 1e-06:
unstable_modes = 1
if unstable_modes == 1:
count = count + 1
continue
else:
check = 1
return fcm
# acoustic sum
[docs] def get_asum_FCM(self, fcm, numiter=15):
"""
Generate a symmeterized force constant matrix that obeys the objects symmetry
constraints and obeys the acoustic sum rule through an iterative procedure
Args:
fcm (numpy array): 3Nx3N unsymmeterized force constant matrix
numiter (int): number of iterations to attempt to obey the acoustic sum
rule
Return:
numpy array representing the force constant matrix
"""
# set max force in reciprocal space
operations = self.FCM_operations
numsites = len(self.structure)
D = np.ones([numsites * 3, numsites * 3])
for num in range(numiter):
X = np.real(fcm)
# symmetry operations
pastrow = 0
total = np.zeros([3, 3])
for col in range(numsites):
total = total + X[0:3, col * 3:col * 3 + 3]
total = total / (numsites)
for op in operations:
same = 0
transpose = 0
if op[0] == op[1] and op[0] == op[2] and op[0] == op[3]:
same = 1
if op[0] == op[3] and op[1] == op[2]:
transpose = 1
if transpose == 0 and same == 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = np.zeros(
[3, 3]
)
for symop in op[4]:
tempfcm = D[
3 * op[2]:3 * op[2] + 3, 3 * op[3]:3 * op[3] + 3
]
tempfcm = symop.transform_tensor(tempfcm)
D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] += tempfcm
if len(op[4]) != 0:
D[3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] / len(op[4])
D[3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3] = D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
].T
continue
else:
# Get the difference in the sum up to this point
currrow = op[0]
if currrow != pastrow:
total = np.zeros([3, 3])
for col in range(numsites):
total = (
total
+ X[
currrow * 3:currrow * 3 + 3, col * 3:col * 3 + 3
]
)
for col in range(currrow):
total = (
total
- D[
currrow * 3:currrow * 3 + 3, col * 3:col * 3 + 3
]
)
total = total / (numsites - currrow)
pastrow = currrow
# Apply the point symmetry operations of the site
temp_tensor = Tensor(total)
temp_tensor_sum = sum(
[
temp_tensor.transform(symm_op)
for symm_op in self.sharedops[op[0]][op[1]]
]
)
if len(self.sharedops[op[0]][op[1]]) != 0:
temp_tensor_sum = temp_tensor_sum / (
len(self.sharedops[op[0]][op[1]])
)
# Apply the proper transformation if there is an equivalent already
if op[0] != op[1]:
for pair in range(len(op[4])):
temp_tensor2 = temp_tensor_sum.T
temp_tensor2 = op[4][pair].transform_tensor(temp_tensor2)
temp_tensor_sum = (temp_tensor_sum + temp_tensor2) / 2
else:
temp_tensor_sum = (temp_tensor_sum + temp_tensor_sum.T) / 2
D[
3 * op[0]:3 * op[0] + 3, 3 * op[1]:3 * op[1] + 3
] = temp_tensor_sum
D[
3 * op[1]:3 * op[1] + 3, 3 * op[0]:3 * op[0] + 3
] = temp_tensor_sum.T
fcm = fcm - D
return fcm
[docs] @requires(Phonopy, "phonopy not installed!")
def get_rand_FCM(self, asum=15, force=10):
"""
Generate a symmeterized force constant matrix from an unsymmeterized matrix
that has no unstable modes and also obeys the acoustic sum rule through an
iterative procedure
Args:
force (float): maximum force constant
asum (int): number of iterations to attempt to obey the acoustic sum
rule
Return:
NxNx3x3 np.array representing the force constant matrix
"""
numsites = len(self.structure.sites)
structure = pymatgen.io.phonopy.get_phonopy_structure(self.structure)
pnstruc = Phonopy(structure, np.eye(3), np.eye(3))
dyn = self.get_unstable_FCM(force)
dyn = self.get_stable_FCM(dyn)
dyn = np.reshape(dyn, (numsites, 3, numsites, 3)).swapaxes(1, 2)
dynmass = np.zeros([len(self.structure), len(self.structure), 3, 3])
masses = []
for j in range(numsites):
masses.append(self.structure.sites[j].specie.atomic_mass)
dynmass = np.zeros([numsites, numsites, 3, 3])
for m in range(numsites):
for n in range(numsites):
dynmass[m][n] = dyn[m][n] * np.sqrt(masses[m]) * np.sqrt(masses[n])
supercell = pnstruc.get_supercell()
primitive = pnstruc.get_primitive()
converter = dyntofc.DynmatToForceConstants(primitive, supercell)
dyn = np.reshape(np.swapaxes(dynmass, 1, 2), (numsites * 3, numsites * 3))
converter.set_dynamical_matrices(dynmat=[dyn])
converter.run()
fc = converter.get_force_constants()
return fc
[docs]def get_piezo(BEC, IST, FCM, rcond=0.0001):
"""
Generate a random piezoelectric tensor based on a structure and corresponding
symmetry
Args:
BEC (numpy array): Nx3x3 array representing the born effective charge tensor
IST (numpy array): Nx3x3x3 array representing the internal strain tensor
FCM (numpy array): NxNx3x3 array representing the born effective charge tensor
rcondy (float): condition for excluding eigenvalues in the pseudoinverse
Return:
3x3x3 calculated Piezo tensor
"""
numsites = len(BEC)
temp_fcm = np.reshape(np.swapaxes(FCM, 1, 2), (numsites * 3, numsites * 3))
eigs, vecs = np.linalg.eig(temp_fcm)
K = np.linalg.pinv(
-temp_fcm,
rcond=np.abs(eigs[np.argsort(np.abs(eigs))[2]])
/ np.abs(eigs[np.argsort(np.abs(eigs))[-1]])
+ rcond,
)
K = np.reshape(K, (numsites, 3, numsites, 3)).swapaxes(1, 2)
return np.einsum("ikl,ijlm,jmno->kno", BEC, K, IST) * 16.0216559424
[docs]@requires(Phonopy, "phonopy not installed!")
def rand_piezo(struc, pointops, sharedops, BEC, IST, FCM, anumiter=10):
"""
Generate a random piezoelectric tensor based on a structure and corresponding
symmetry
Args:
struc (pymatgen structure): structure whose symmetry operations the piezo tensor must obey
pointops: list of point operations obeyed by a single atomic site
sharedops: list of point operations shared by a pair of atomic sites
BEC (numpy array): Nx3x3 array representing the born effective charge tensor
IST (numpy array): Nx3x3x3 array representing the internal strain tensor
FCM (numpy array): NxNx3x3 array representing the born effective charge tensor
anumiter (int): number of iterations for acoustic sum rule convergence
Return:
list in the form of [Nx3x3 random born effective charge tenosr,
Nx3x3x3 random internal strain tensor, NxNx3x3 random force constant matrix, 3x3x3 piezo tensor]
"""
bec = BornEffectiveCharge(struc, BEC, pointops)
bec.get_BEC_operations()
rand_BEC = bec.get_rand_BEC()
ist = InternalStrainTensor(struc, IST, pointops)
ist.get_IST_operations()
rand_IST = ist.get_rand_IST()
fcm = ForceConstantMatrix(struc, FCM, pointops, sharedops)
fcm.get_FCM_operations()
rand_FCM = fcm.get_rand_FCM()
P = get_piezo(rand_BEC, rand_IST, rand_FCM) * 16.0216559424 / struc.volume
return (rand_BEC, rand_IST, rand_FCM, P)