From c14fbc93c4342fac353dd3665ae2437470cdc284 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Wed, 17 Jan 2024 10:24:58 -0500 Subject: [PATCH] format --- pyxtal/elasticity.py | 92 ++++++++++++++++++++++---------------------- pyxtal/lattice.py | 16 ++++---- pyxtal/supergroup.py | 2 +- requirements.txt | 1 - 4 files changed, 55 insertions(+), 56 deletions(-) diff --git a/pyxtal/elasticity.py b/pyxtal/elasticity.py index d957bbdd..8f5e8eb3 100644 --- a/pyxtal/elasticity.py +++ b/pyxtal/elasticity.py @@ -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)): @@ -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 @@ -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.') @@ -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 @@ -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.') @@ -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.') @@ -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.') @@ -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) @@ -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 @@ -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], @@ -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): """ @@ -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. @@ -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). @@ -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 ----- @@ -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 @@ -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] + @@ -1038,21 +1038,21 @@ 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 @@ -1060,14 +1060,14 @@ def poisson_ratio(C, l, m): 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 @@ -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 @@ -1091,13 +1091,13 @@ 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 @@ -1105,18 +1105,18 @@ def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6): 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: @@ -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) @@ -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 diff --git a/pyxtal/lattice.py b/pyxtal/lattice.py index f6885755..fd018786 100644 --- a/pyxtal/lattice.py +++ b/pyxtal/lattice.py @@ -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']: @@ -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) @@ -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') @@ -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 @@ -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]) diff --git a/pyxtal/supergroup.py b/pyxtal/supergroup.py index d31c3110..c81bdfa1 100644 --- a/pyxtal/supergroup.py +++ b/pyxtal/supergroup.py @@ -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: diff --git a/requirements.txt b/requirements.txt index 988d77f2..952858f6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,3 @@ networkx>=2.3 scipy>=1.7.3 importlib_metadata>=1.4 ase>=3.18.0 -pyshtools>=4.10.3