Skip to content

Commit

Permalink
MAINT: Fixed bug with incorrect Voigt notation and cleaned code
Browse files Browse the repository at this point in the history
  • Loading branch information
@Fraser-Birks
Fraser-Birks committed Nov 15, 2023
1 parent 6d6463c commit 9b4b00e
Show file tree
Hide file tree
Showing 2 changed files with 70 additions and 126 deletions.
104 changes: 32 additions & 72 deletions matscipy/cauchy_born.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
from itertools import permutations

from scipy.stats.qmc import LatinHypercube, scale

from matscipy.elasticity import full_3x3_to_Voigt_6_strain, Voigt_6_to_full_3x3_strain
import ase


Expand Down Expand Up @@ -68,12 +68,12 @@ def predict_gradient(self, grad_phi):
return grad_phi @ self.model

def save(self):
"""Saves the regression model to file: CB_model.npy"""
np.save('CB_model.npy', self.model)
"""Saves the regression model to file: CB_model.txt"""
np.savetxt('CB_model.txt', self.model)

def load(self):
"""Loads the regression model from file: CB_model.npy"""
self.model = np.load('CB_model.npy')
"""Loads the regression model from file: CB_model.txt"""
self.model = np.loadtxt('CB_model.txt')


class CubicCauchyBorn:
Expand Down Expand Up @@ -140,6 +140,8 @@ def set_sublattices(self, atoms, A, read_from_atoms=False):

# generate U (right hand stretch tensor) to apply
if read_from_atoms:
if ('lattice1mask' not in atoms.arrays) or ('lattice2mask' not in atoms.arrays):
raise KeyError('Lattice masks not found in atoms object')
lattice1mask = atoms.arrays['lattice1mask']
lattice2mask = atoms.arrays['lattice2mask']
else:
Expand Down Expand Up @@ -183,6 +185,8 @@ def set_sublattices(self, atoms, A, read_from_atoms=False):
if all(element == lattice1mask[0] for element in lattice1mask):
lattice1mask = force_diff_lattice[:, 2] > 0
lattice2mask = force_diff_lattice[:, 2] < 0
if all(element == lattice1mask[0] for element in lattice1mask):
raise RuntimeError('No forces detected in any direction when shift applied - cannot assign sublattices')

# set new array in atoms
try: # make new arrays if they didn't already exist
Expand All @@ -209,10 +213,11 @@ def switch_sublattices(self, atoms):
----------
atoms : ASE atoms object with defined sublattices
"""
tmpl1 = self.lattice1mask.copy()
tmpl2 = self.lattice2mask.copy()
self.lattice1mask = tmpl2.copy()
self.lattice2mask = tmpl1.copy()
self.lattice1mask, self.lattice2mask = self.lattice2mask.copy(), self.lattice1mask.copy()
#tmpl1 = self.lattice1mask.copy()
#tmpl2 = self.lattice2mask.copy()
#self.lattice1mask = tmpl2.copy()
#self.lattice2mask = tmpl1.copy()
lattice1maskatoms = atoms.arrays['lattice1mask']
lattice2maskatoms = atoms.arrays['lattice2mask']
lattice1maskatoms[:] = self.lattice1mask
Expand Down Expand Up @@ -256,13 +261,7 @@ def f_gl(E_vec, calc, unitcell):
"""
# green lagrange version of function
# turn E into matrix
E = np.zeros([3, 3])
E[0, 0] = E_vec[0]
E[1, 1] = E_vec[1]
E[2, 2] = E_vec[2]
E[2, 1], E[1, 2] = E_vec[3], E_vec[3]
E[2, 0], E[0, 2] = E_vec[4], E_vec[4]
E[0, 1], E[1, 0] = E_vec[5], E_vec[5]
E = Voigt_6_to_full_3x3_strain(E_vec)
# get U^2
Usqr = 2 * E + np.eye(3)
# square root matrix to get U
Expand All @@ -286,8 +285,7 @@ def f_gl(E_vec, calc, unitcell):
# relax cell
opt = LBFGS(unitcell_copy, logfile=None)
opt.run(fmax=1e-10)
relaxed_shift[:] = unitcell_copy.get_positions(
)[1] - unitcell_copy.get_positions()[0]
relaxed_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]

# get shift
shift_diff = relaxed_shift - initial_shift
Expand Down Expand Up @@ -954,13 +952,7 @@ def eval_eps(eps_vec, grad_f, hess_f):
return predictions

# turn E into voigt vector
E_voigt = np.zeros([np.shape(E)[0], 6])
E_voigt[:, 0] = E[:, 0, 0]
E_voigt[:, 1] = E[:, 1, 1]
E_voigt[:, 2] = E[:, 2, 2]
E_voigt[:, 3] = E[:, 1, 2]
E_voigt[:, 4] = E[:, 0, 2]
E_voigt[:, 5] = E[:, 0, 1]
E_voigt = full_3x3_to_Voigt_6_strain(E)

# return the shift vectors
if method == 'taylor':
Expand Down Expand Up @@ -1219,13 +1211,7 @@ def check_for_refit(
# print('Es',E)
# we get back the strain vectors in the lattice frame
# turn E into voigt vectors
E_voigt = np.zeros([np.shape(E)[0], 6])
E_voigt[:, 0] = E[:, 0, 0]
E_voigt[:, 1] = E[:, 1, 1]
E_voigt[:, 2] = E[:, 2, 2]
E_voigt[:, 3] = E[:, 1, 2]
E_voigt[:, 4] = E[:, 0, 2]
E_voigt[:, 5] = E[:, 0, 1]
E_voigt = full_3x3_to_Voigt_6_strain(E)

# pass set of Es to evaluate to refit_regression
self.refit_regression(atoms, E_voigt)
Expand Down Expand Up @@ -1336,13 +1322,7 @@ def eval_shift(self, E_vec, unitcell):
"""
# green lagrange version of function
# turn E into matrix
E = np.zeros([3, 3])
E[0, 0] = E_vec[0]
E[1, 1] = E_vec[1]
E[2, 2] = E_vec[2]
E[2, 1], E[1, 2] = E_vec[3], E_vec[3]
E[2, 0], E[0, 2] = E_vec[4], E_vec[4]
E[0, 1], E[1, 0] = E_vec[5], E_vec[5]
E = Voigt_6_to_full_3x3_strain(E_vec)
# get U^2
Usqr = 2 * E + np.eye(3)
# square root matrix
Expand Down Expand Up @@ -1462,39 +1442,19 @@ def evaluate_shift_gradient_regression(self, E, dE):
"""

# turn E and dE into voigt vectors
E_voigt = np.zeros([np.shape(E)[0], 6])
E_voigt[:, 0] = E[:, 0, 0]
E_voigt[:, 1] = E[:, 1, 1]
E_voigt[:, 2] = E[:, 2, 2]
E_voigt[:, 3] = E[:, 1, 2]
E_voigt[:, 4] = E[:, 0, 2]
E_voigt[:, 5] = E[:, 0, 1]

dE_voigt = np.zeros([np.shape(dE)[0], 6])
dE_voigt[:, 0] = dE[:, 0, 0]
dE_voigt[:, 1] = dE[:, 1, 1]
dE_voigt[:, 2] = dE[:, 2, 2]
dE_voigt[:, 3] = dE[:, 1, 2]
dE_voigt[:, 4] = dE[:, 0, 2]
dE_voigt[:, 5] = dE[:, 0, 1]
E_voigt = full_3x3_to_Voigt_6_strain(E)
dE_voigt = full_3x3_to_Voigt_6_strain(dE)

# Get the rotated versions of the strain and strain deriv tensors
epsx = E_voigt
epsy = self.permutation(E_voigt, 1)
epsz = self.permutation(E_voigt, 2)

depsx = dE_voigt
depsy = self.permutation(dE_voigt, 1)
depsz = self.permutation(dE_voigt, 2)
predictions = np.zeros([np.shape(epsx)[0], 3])

# Predict the shift gradients using the fitted regression model
predictions[:, 0] = self.RM.predict_gradient(
self.grad_basis_function_evaluation(epsx, depsx))
predictions[:, 1] = self.RM.predict_gradient(
self.grad_basis_function_evaluation(epsy, depsy))
predictions[:, 2] = self.RM.predict_gradient(
self.grad_basis_function_evaluation(epsz, depsz))
eps = [E_voigt, self.permutation(E_voigt, 1), self.permutation(E_voigt, 2)]

deps = [dE_voigt, self.permutation(dE_voigt, 1), self.permutation(dE_voigt, 2)]
predictions = np.zeros([np.shape(eps[0])[0], 3])

for i in range(3):
# Predict the shift gradients using the fitted regression model
predictions[:, i] = self.RM.predict_gradient(
self.grad_basis_function_evaluation(eps[i], deps[i]))
return predictions

def get_shift_gradients(
Expand Down Expand Up @@ -1590,11 +1550,11 @@ def get_shift_gradients(
return dshifts

def save_regression_model(self):
"""Saves the regression model to file: CB_model.npy"""
"""Saves the regression model to file: CB_model.txt"""
self.RM.save()

def load_regression_model(self):
"""Loads the regression model from file: CB_model.npy"""
"""Loads the regression model from file: CB_model.txt"""
self.RM = RegressionModel()
self.RM.load()

Expand Down
92 changes: 38 additions & 54 deletions tests/test_cauchy_born_corrector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
from scipy.linalg import sqrtm
import ase.io
import matscipytest
from matscipy.elasticity import Voigt_6_to_full_3x3_strain


class TestPredictCauchyBornShifts(matscipytest.MatSciPyTestCase):
Expand Down Expand Up @@ -44,64 +45,49 @@ def test_set_sublattice(self):

def test_fit_taylor_model(self):
self.cb.fit_taylor()
# print(self.cb.grad_f)
# print(self.cb.hess_f)
grad_f = np.array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -
1.09995082e+00, -2.22044605e-12, -3.33066907e-12])
hess_f = np.array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.03958794e-01, 0.00000000e+00, 1.11022302e-08],
print(self.cb.grad_f)
print(self.cb.hess_f)
grad_f = np.array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-5.49975395e-01, 0.00000000e+00, -2.22044605e-12])
hess_f = np.array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
1.01979464e-01, 1.66533454e-08, 1.66533454e-08],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
9.12311482e-01, -5.55111512e-09, 0.00000000e+00],
4.56155597e-01, -1.66533454e-08, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
9.12311465e-01, 0.00000000e+00, 0.00000000e+00],
[2.03958794e-01, 9.12311493e-01, 9.12311465e-01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[0.00000000e+00, -5.55111512e-09, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[1.11022302e-08, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
4.56155608e-01, 0.00000000e+00, 0.00000000e+00],
[1.01979464e-01, 4.56155597e-01, 4.56155608e-01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[1.66533454e-08, -1.66533454e-08, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[1.66533454e-08, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
assert np.allclose(self.cb.grad_f, grad_f, atol=1e-6)
assert np.allclose(self.cb.hess_f, hess_f, atol=1e-6)

def E_cart3D(self, x, y, z, eps=None):
E = np.zeros([len(x), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps = np.reshape(eps,[1,6])
eps_vec = np.repeat(eps,len(x),axis=0)
E = Voigt_6_to_full_3x3_strain(eps_vec)
return E

def E_cylind3D(self, r, theta, z, eps=None):
E = np.zeros([len(r), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps = np.reshape(eps,[1,6])
eps_vec = np.repeat(eps,len(r),axis=0)
E = Voigt_6_to_full_3x3_strain(eps_vec)
return E

def F_cart3D(self, x, y, z, eps=None):
E = np.zeros([len(x), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps = np.reshape(eps,[1,6])
eps_vec = np.repeat(eps,len(x),axis=0)
E = Voigt_6_to_full_3x3_strain(eps_vec)
U = np.zeros_like(E)
for i in range(np.shape(E)[0]):
U[i, :, :] = sqrtm((2*E[i, :, :])+np.eye(3))
return U # with no rigid rotation, U is F

def F_cylind3D(self, r, theta, z, eps=None):
E = np.zeros([len(r), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps = np.reshape(eps,[1,6])
eps_vec = np.repeat(eps,len(r),axis=0)
E = Voigt_6_to_full_3x3_strain(eps_vec)
U = np.zeros_like(E)
for i in range(np.shape(E)[0]):
U[i, :, :] = sqrtm((2*E[i, :, :])+np.eye(3))
Expand Down Expand Up @@ -234,13 +220,13 @@ def test_regression_model_refit(self):
# print(shift_err_before,shift_err_after)

def E_cart3D_with_de(self, x, y, z, eps=None, de=0):
E = np.zeros([len(x), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]+de
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps_arr = eps.copy()
eps_arr[1] += de
eps_arr = np.reshape(eps_arr,[1,6])
eps_vec = np.repeat(eps_arr,len(x),axis=0)
print('eps_vec', eps_vec[0,:])
E = Voigt_6_to_full_3x3_strain(eps_vec)
print('E out', E[0,:,:])
return E

def test_regression_model_gradient_E(self):
Expand Down Expand Up @@ -282,13 +268,11 @@ def test_regression_model_gradient_E(self):
assert np.allclose(nu_grad_fd, nu_grad, 1e-8)

def F_cart3D_with_de(self, x, y, z, eps=None, de=0):
E = np.zeros([len(x), 3, 3])
E[:, 0, 0] = eps[0]
E[:, 1, 1] = eps[1]+de
E[:, 2, 2] = eps[2]
E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
eps_arr = eps.copy()
eps_arr[1] += de
eps_arr = np.reshape(eps_arr,[1,6])
eps_vec = np.repeat(eps_arr,len(x),axis=0)
E = Voigt_6_to_full_3x3_strain(eps_vec)
U = np.zeros_like(E)
for i in range(np.shape(E)[0]):
U[i, :, :] = sqrtm((2*E[i, :, :])+np.eye(3))
Expand Down

0 comments on commit 9b4b00e

Please sign in to comment.