Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
@qzhu2017
qzhu2017 committed Jan 17, 2024
1 parent 425d883 commit c14fbc9
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 56 deletions.
92 changes: 46 additions & 46 deletions pyxtal/elasticity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,13 +102,13 @@ def Voigt_6x6_to_full_3x3x3x3(C):
----------
C : array_like
6x6 stiffness matrix (Voigt notation).
Returns
-------
C : array_like
3x3x3x3 stiffness matrix.
"""

C = np.asarray(C)
C_out = np.zeros((3,3,3,3), dtype=float)
for i, j, k, l in itertools.product(range(3), range(3), range(3), range(3)):
Expand Down Expand Up @@ -251,7 +251,7 @@ def invariants(s, syy=None, szz=None, syz=None, sxz=None, sxy=None,

def rotate_cubic_elastic_constants(C11, C12, C44, A, tol=1e-6):
"""
Return rotated elastic moduli for a cubic crystal given the elastic
Return rotated elastic moduli for a cubic crystal given the elastic
constant in standard C11, C12, C44 notation.
Parameters
Expand All @@ -270,7 +270,7 @@ def rotate_cubic_elastic_constants(C11, C12, C44, A, tol=1e-6):
A = np.asarray(A)

# Is this a rotation matrix?
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
np.eye(3, dtype=float)) > tol):
raise RuntimeError('Matrix *A* does not describe a rotation.')

Expand Down Expand Up @@ -301,7 +301,7 @@ def rotate_cubic_elastic_constants(C11, C12, C44, A, tol=1e-6):

def rotate_elastic_constants(C, A, tol=1e-6):
"""
Return rotated elastic moduli for a general crystal given the elastic
Return rotated elastic moduli for a general crystal given the elastic
constant in Voigt notation.
Parameters
Expand All @@ -320,7 +320,7 @@ def rotate_elastic_constants(C, A, tol=1e-6):
A = np.asarray(A)

# Is this a rotation matrix?
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
np.eye(3, dtype=float)) > tol):
raise RuntimeError('Matrix *A* does not describe a rotation.')

Expand Down Expand Up @@ -361,7 +361,7 @@ def rotate(self, A):
A = np.asarray(A)

# Is this a rotation matrix?
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
np.eye(3, dtype=float)) > self.tol):
raise RuntimeError('Matrix *A* does not describe a rotation.')

Expand Down Expand Up @@ -392,7 +392,7 @@ def _rotate_explicit(self, A):
A = np.asarray(A)

# Is this a rotation matrix?
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
if np.sometrue(np.abs(np.dot(np.array(A), np.transpose(np.array(A))) -
np.eye(3, dtype=float) ) > self.tol):
raise RuntimeError('Matrix *A* does not describe a rotation.')

Expand Down Expand Up @@ -438,7 +438,7 @@ def compliance(self):

###

def measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None,
def measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None,
logfile=None, **kwargs):
"""
Brute-force measurement of elastic constants for a triclinic (general)
Expand Down Expand Up @@ -473,7 +473,7 @@ def measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None,
for j in range(3):
a.set_cell(cell, scale_atoms=True)
a.set_positions(r0)

D = np.eye(3)
D[i, j] += 0.5*delta
D[j, i] += 0.5*delta
Expand Down Expand Up @@ -548,7 +548,7 @@ def measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None,
[ 0, 0, 0, 4, 0, 20],
[10, 14, 17, 0, 5, 0],
[ 0, 0, 0, 20, 0, 6]]),

'triclinic': np.array([[ 1, 7, 8, 9, 10, 11],
[ 7, 2, 12, 13, 14, 15],
[ 8, 12, 3, 16, 17, 18],
Expand Down Expand Up @@ -718,7 +718,7 @@ def generate_strained_configs(at0, symmetry='triclinic', N_steps=5, delta=1e-2):
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2,
def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2,
optimizer=None, verbose=True, GPa=True, tag='test',
graphics=False, logfile=None, **kwargs):
"""
Expand All @@ -737,7 +737,7 @@ def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2,
N_steps : int
Number of atomic configurations to generate for each strain pattern.
Default is 5. Absolute strain values range from -delta*N_steps/2
to +delta*N_steps/2.
to +delta*N_steps/2.
delta : float
Strain increment for analytical derivatives of stresses.
Default is 1e-2.
Expand All @@ -753,7 +753,7 @@ def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2,
If True, convert to GPa
graphics : bool
If True, use :mod:`matplotlib.pyplot` to plot the stress vs. strain
curve for each C_ij component fitted. Default True.
curve for each C_ij component fitted. Default True.
logfile : bool
Log file to write optimizer output to. Default None (i.e. suppress
output).
Expand All @@ -768,7 +768,7 @@ def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2,
If scipy.stats module is available then error estimates for each C_ij
component are obtained from the accuracy of the linear regression.
Otherwise an array of np.zeros((6,6)) is returned.
Notes
-----
Expand Down Expand Up @@ -909,7 +909,7 @@ def do_fit(index1, index2, stress, strain, patt):
#print("Cell\n", at.get_cell())
strs += "Strain {:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}\n".format(*strain_info)
strs += "Stress (GPa) {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}\n".format(*(stress_info/units.GPa))
f.write(strs)
f.write(strs)
print(strs)

# Do the linear regression
Expand Down Expand Up @@ -1014,14 +1014,14 @@ def youngs_modulus(C, l):
Notes
-----
Formula is from W. Brantley, Calculated elastic constants for stress problems associated
Formula is from W. Brantley, Calculated elastic constants for stress problems associated
with semiconductor devices. J. Appl. Phys., 44, 534 (1973).
"""
"""

S = inv(C) # Compliance matrix
lhat = l/norm(l) # Normalise directions

# Youngs modulus in direction l, ratio of stress sigma_l
# Youngs modulus in direction l, ratio of stress sigma_l
# to strain response epsilon_l
E = 1.0/(S[0,0] - 2.0*(S[0,0]-S[0,1]-0.5*S[3,3])*(lhat[0]*lhat[0]*lhat[1]*lhat[1] +
lhat[1]*lhat[1]*lhat[2]*lhat[2] +
Expand All @@ -1038,36 +1038,36 @@ def poisson_ratio(C, l, m):
Notes
-----
Formula is from W. Brantley, Calculated elastic constants for stress problems associated
Formula is from W. Brantley, Calculated elastic constants for stress problems associated
with semiconductor devices. J. Appl. Phys., 44, 534 (1973).
"""

S = inv(C) # Compliance matrix
lhat = l/norm(l) # Normalise directions
mhat = m/norm(m)

# Poisson ratio v_lm: response in m direction to strain in
# Poisson ratio v_lm: response in m direction to strain in
# l direction, v_lm = - epsilon_m/epsilon_l
v = -((S[0,1] + (S[0,0]-S[0,1]-0.5*S[3,3])*(lhat[0]*lhat[0]*mhat[0]*mhat[0] +
lhat[1]*lhat[1]*mhat[1]*mhat[1] +
lhat[2]*lhat[2]*mhat[2]*mhat[2])) /
(S[0,0] - 2.0*(S[0,0]-S[0,1]-0.5*S[3,3])*(lhat[0]*lhat[0]*lhat[1]*lhat[1] +
lhat[1]*lhat[1]*lhat[2]*lhat[2] +
lhat[2]*lhat[2]*mhat[2]*mhat[2])) /
(S[0,0] - 2.0*(S[0,0]-S[0,1]-0.5*S[3,3])*(lhat[0]*lhat[0]*lhat[1]*lhat[1] +
lhat[1]*lhat[1]*lhat[2]*lhat[2] +
lhat[0]*lhat[0]*lhat[2]*lhat[2])))
return v


def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6):
"""
Calculate elastic moduli from 6x6 elastic constant matrix C_{ij}.
The elastic moduli calculated are: Young's muduli, Poisson's ratios,
shear moduli, bulk mudulus and bulk mudulus tensor.
If a direction l is specified, the system is rotated to have it as its
x direction (see Notes for details). If R is specified the system is
rotated according to it.
Parameters
----------
C : array_like
Expand All @@ -1077,7 +1077,7 @@ def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6):
of the original system)
R : array_like, optional
3x3 rotation matrix.
Returns
-------
E : array_like
Expand All @@ -1091,32 +1091,32 @@ def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6):
Bulk modulus.
K : array_like
3x3 matrix with bulk modulus tensor.
Other Parameters
----------------
tol : float, optional
tolerance for checking validity of rotation and comparison
of vectors.
Notes
---
It works by rotating the elastic constant tensor to the desired
direction, so it should be valid for arbitrary crystal structures.
If only l is specified there is an infinite number of possible
rotations. The chosen one is a rotation along the axis orthogonal
to the plane defined by the vectors (1, 0, 0) and l.
Bulk modulus tensor as defined in
O. Rand and V. Rovenski, "Analytical Methods in Anisotropic
Elasticity", Birkh\"auser (2005), pp. 71.
"""

if R is not None:
R = np.asarray(R)

# Is this a rotation matrix?
if np.sometrue(np.abs(np.dot(np.array(R), np.transpose(np.array(R))) -
if np.sometrue(np.abs(np.dot(np.array(R), np.transpose(np.array(R))) -
np.eye(3, dtype=float)) > tol):
raise RuntimeError('Matrix *R* does not describe a rotation.')
else:
Expand Down Expand Up @@ -1163,7 +1163,7 @@ def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6):

# Bulk modulus
B = 1/np.sum(S[0:3, 0:3])

# Bulk modulus tensor
Crt = Voigt_6x6_to_full_3x3x3x3(Cr)
K = np.einsum('ijkk', Crt)
Expand All @@ -1177,18 +1177,18 @@ def elastic_properties(C):
Kv = C[:3,:3].mean()
Gv = (C[0,0]+C[1,1]+C[2,2] - (C[0,1]+C[1,2]+C[2,0]) + 3*(C[3,3]+C[4,4]+C[5,5]))/15
Ev = 1/((1/(3*Gv))+(1/(9*Kv)))
vv = 0.5*(1-((3*Gv)/(3*Kv+Gv)));
vv = 0.5*(1-((3*Gv)/(3*Kv+Gv)));

S = np.linalg.inv(C)
Kr = 1/((S[0,0]+S[1,1]+S[2,2])+2*(S[0,1]+S[1,2]+S[2,0]))
Gr = 15/(4*(S[0,0]+S[1,1]+S[2,2])-4*(S[0,1]+S[1,2]+S[2,0])+3*(S[3,3]+S[4,4]+S[5,5]))
Er = 1/((1/(3*Gr))+(1/(9*Kr)))
vr = 0.5*(1-((3*Gr)/(3*Kr+Gr)))

Kh = (Kv+Kr)/2
Gh = (Gv+Gr)/2
Eh = (Ev+Er)/2
vh = (vv+vr)/2
Kr = 1/((S[0,0]+S[1,1]+S[2,2])+2*(S[0,1]+S[1,2]+S[2,0]))
Gr = 15/(4*(S[0,0]+S[1,1]+S[2,2])-4*(S[0,1]+S[1,2]+S[2,0])+3*(S[3,3]+S[4,4]+S[5,5]))
Er = 1/((1/(3*Gr))+(1/(9*Kr)))
vr = 0.5*(1-((3*Gr)/(3*Kr+Gr)))

Kh = (Kv+Kr)/2
Gh = (Gv+Gr)/2
Eh = (Ev+Er)/2
vh = (vv+vr)/2
return Kv, Gv, Ev, vv, Kr, Gr, Er, vr, Kh, Gh, Eh, vh


16 changes: 8 additions & 8 deletions pyxtal/lattice.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ def get_dofs(self, ltype):
get the number of degree of freedom
"""
if ltype in ["triclinic"]:
dofs = [3, 3]
dofs = [3, 3]
elif ltype in ["monoclinic"]:
dofs = [3, 1]
elif ltype in ['orthorhombic']:
Expand Down Expand Up @@ -984,7 +984,7 @@ def find_transition_to_orthoslab(self, c=(0,0,1), a=(1,0,0), m=5):

tol = 1e-3
direction = np.array(c)

# find the simplest a-direction
if np.dot(np.array(a), direction) < tol:
a_hkl = np.array(a)
Expand All @@ -1001,7 +1001,7 @@ def find_transition_to_orthoslab(self, c=(0,0,1), a=(1,0,0), m=5):
a_hkl = a_hkls[np.argmin(np.abs(a_hkls).sum(axis=1))]
a_vector = np.dot(a_hkl, self.matrix)
#print('a_hkl', a_hkl)

# find the simplest b-direction
b_hkl = None
min_angle_ab = float('inf')
Expand All @@ -1010,20 +1010,20 @@ def find_transition_to_orthoslab(self, c=(0,0,1), a=(1,0,0), m=5):
for l in range(-m, m+1):
hkl = np.array([h, k, l])
if [h, k, l] != [0, 0, 0] and not has_reduction(hkl):
if abs(np.dot(hkl, direction)) < tol:
if abs(np.dot(hkl, direction)) < tol:
vector = np.dot(hkl, self.matrix)
angle1 = angle(vector, a_vector, radians=False)
if abs(90-angle1) < min_angle_ab:
min_angle_ab = abs(90-angle1)
b_hkl = hkl
b_vector = vector

#print('b_hkl', b_hkl, min_angle_ab)
# change the sign
if abs(angle(np.cross(a_hkl, b_hkl), direction))>tol:
if abs(angle(np.cross(a_hkl, b_hkl), direction))>tol:
b_hkl *= -1
b_vector *= -1

## update the c_direction
ab_plane = np.cross(a_vector, b_vector)#; print('ab_plane', ab_plane)
c_hkl = None
Expand All @@ -1040,7 +1040,7 @@ def find_transition_to_orthoslab(self, c=(0,0,1), a=(1,0,0), m=5):
min_angle_c = angle1
c_hkl = hkl
c_vector = vector

#print(a_hkl, b_hkl, c_hkl)
return np.vstack([a_hkl, b_hkl, c_hkl])

Expand Down
2 changes: 1 addition & 1 deletion pyxtal/supergroup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -457,7 +457,7 @@ def calc_disps(self, split_id, solution, d_tol):
if mask is None or len(mask)<3:
def fun(translation, mapping, splitter, mask):
return self.symmetrize_dist(splitter, mapping, mask, translation)[0]

res = minimize(fun, translations[id], args=(mappings[id], splitter, mask),
method='Nelder-Mead', options={'maxiter': 10})
if res.fun < max_disp:
Expand Down
1 change: 0 additions & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,3 @@ networkx>=2.3
scipy>=1.7.3
importlib_metadata>=1.4
ase>=3.18.0
pyshtools>=4.10.3

0 comments on commit c14fbc9

Please sign in to comment.