diff --git a/gb_code/csl_generator.py b/gb_code/csl_generator.py index 71b5ca6..509ec41 100644 --- a/gb_code/csl_generator.py +++ b/gb_code/csl_generator.py @@ -31,7 +31,8 @@ def get_cubic_sigma(uvw, m, n=1): """ - CSL analytical formula based on the book: 'Interfaces in crystalline materials', + CSL analytical formula based on the book: + 'Interfaces in crystalline materials', Sutton and Balluffi, clarendon press, 1996. generates possible sigma values. arguments: @@ -81,6 +82,7 @@ def get_theta_m_n_list(uvw, sigma): thetas.sort(key=lambda x: x[0]) return thetas + def print_list(uvw, limit): """ prints a list of smallest sigmas/angles for a given axis(uvw). @@ -88,10 +90,11 @@ def print_list(uvw, limit): for i in range(limit): tt = get_theta_m_n_list(uvw, i) if len(tt) > 0: - theta, m, n = tt[0] + theta, _, _ = tt[0] print("Sigma: {0:3d} Theta: {1:5.2f} " .format(i, degrees(theta))) + def rot(a, Theta): """ produces a rotation matrix. @@ -112,7 +115,7 @@ def rot(a, Theta): # Helpful Functions: -#-------------------# +# -------------------# def integer_array(A, tol=1e-7): """ @@ -120,6 +123,7 @@ def integer_array(A, tol=1e-7): """ return np.all(abs(np.round(A) - A) < tol) + def angv(a, b): """ returns the angle between two vectors. @@ -127,6 +131,7 @@ def angv(a, b): ang = acos(np.round(dot(a, b)/norm(a)/norm(b), 8)) return round(degrees(ang), 7) + def ang(a, b): """ returns the cos(angle) between two vectors. @@ -134,6 +139,7 @@ def ang(a, b): ang = np.round(dot(a, b)/norm(a)/norm(b), 7) return abs(ang) + def CommonDivisor(a): """ returns the common factor of vector a and the reduced vector. @@ -146,6 +152,7 @@ def CommonDivisor(a): CommFac.append(i) return(a.astype(int), (abs(np.prod(CommFac)))) + def SmallestInteger(a): """ returns the smallest multiple integer to make an integer array. @@ -157,12 +164,13 @@ def SmallestInteger(a): break return (testV, i) if integer_array(testV) else None + def integerMatrix(a): """ returns an integer matrix from row vectors. """ Found = True - b = np.zeros((3,3)) + b = np.zeros((3, 3)) a = np.array(a) for i in range(3): for j in range(1, 2000): @@ -175,6 +183,7 @@ def integerMatrix(a): print("Can not make integer matrix!") return (b) if Found else None + def SymmEquivalent(arr): """ returns cubic symmetric eqivalents of the given 2 dimensional vector. @@ -213,6 +222,7 @@ def SymmEquivalent(arr): Result = np.array(Result) return np.unique(Result, axis=0) + def Tilt_Twist_comp(v1, uvw, m, n): """ returns the tilt and twist components of a given GB plane. @@ -270,7 +280,7 @@ def Create_Possible_GB_Plane_List(uvw, m=5, n=1, lim=5): [1, 1, 0], [0, 1, 1], [1, 0, 1], - ], dtype='float') + ], dtype='float') # Find GB plane coordinates: for i in range(len(indice_0)): if CommonDivisor(indice_0[i])[1] <= 1: @@ -302,6 +312,7 @@ def Create_Possible_GB_Plane_List(uvw, m=5, n=1, lim=5): return (V1, V2, MeanPlanes, GBtype) + def Create_minimal_cell_Method_1(sigma, uvw, R): """ finds Minimal cell by means of a numerical search. @@ -376,6 +387,7 @@ def MiniCell_search(indices, MiniCell_1, R, sigma): else: return Found + def Basis(basis): """ defines the basis. @@ -407,6 +419,7 @@ def Basis(basis): return basis + def Find_Orthogonal_cell(basis, uvw, m, n, GB1): """ finds Orthogonal cells from the CSL minimal cells. @@ -473,7 +486,7 @@ def Find_Orthogonal_cell(basis, uvw, m, n, GB1): OrthoCell_1 = OrthoCell_1.astype(float) OrthoCell_2 = OrthoCell_2.astype(float) - if basis == 'sc' or basis == 'diamond' : + if basis == 'sc' or basis == 'diamond': return ((OrthoCell_1.astype(float), OrthoCell_2.astype(float), Num.astype(int))) @@ -487,17 +500,18 @@ def Find_Orthogonal_cell(basis, uvw, m, n, GB1): else: return None + def print_list_GB_Planes(uvw, basis, m, n, lim=3): """ prints lists of GB planes given an axis, basis, m and n. """ uvw = np.array(uvw) - V1, V2, Mean, Type = Create_Possible_GB_Plane_List(uvw, m, n, lim) + V1, V2, _, Type = Create_Possible_GB_Plane_List(uvw, m, n, lim) for i in range(len(V1)): Or = Find_Orthogonal_cell(basis, uvw, m, n, V1[i]) if Or: - print ("{0:<20s} {1:<20s} {2:<20s} {3:<10s}" - .format(str(V1[i]), str(V2[i]), Type[i], str(Or[2]))) + print("{0:<20s} {1:<20s} {2:<20s} {3:<10s}" + .format(str(V1[i]), str(V2[i]), Type[i], str(Or[2]))) # ___CSL/DSC vector construction___# @@ -608,6 +622,7 @@ def face_centering(a): return z_f if det(z_f) == 0.25 else None + def DSC_vec(basis, sigma, minicell): """ a discrete shift complete (DSC) @@ -626,6 +641,7 @@ def DSC_vec(basis, sigma, minicell): D = dot(D_sc, face_centering(D_sc)) return D + def CSL_vec(basis, minicell): """ CSL minimal cell for sc, fcc and bcc. @@ -642,6 +658,7 @@ def CSL_vec(basis, minicell): C = dot(C_sc, face_centering(C_sc)) return C + def DSC_on_plane(D, Pnormal): """ projects the given DSC network on a given plane. @@ -668,8 +685,10 @@ def CSL_density(basis, minicell, plane): density = 1/(hnorm * det(C)) return (abs(density), 1/hnorm) -# An auxilary function to help reduce the size of the small orthogonal cell that -# is decided otherwise based on Sc for fcc and bcc +# An auxilary function to help reduce the size of the small orthogonal cell +# that is decided otherwise based on Sc for fcc and bcc + + def Ortho_fcc_bcc(basis, O1, O2): ortho1 = np.zeros((3, 3)) ortho2 = np.zeros((3, 3)) @@ -685,15 +704,17 @@ def Ortho_fcc_bcc(basis, O1, O2): ortho1[:, i] = O1[:, i] / Min_d ortho2[:, i] = O2[:, i] / Min_d return (ortho1, ortho2) - # Writing to a yaml file that will be read by gb_generator + + def Write_to_io(axis, m, n, basis): """ an input file for gb_generator.py that can be customized. It also contains the output from the usage of csl_generator.py. """ - my_dict = {'GB_plane': str([axis[0], axis[1], axis[2]]), 'lattice_parameter': '4', + my_dict = {'GB_plane': str([axis[0], axis[1], axis[2]]), + 'lattice_parameter': '4', 'overlap_distance': '0.0', 'which_g': 'g1', 'rigid_trans': 'no', 'a': '10', 'b': '5', 'dimensions': '[1,1,1]', @@ -743,8 +764,8 @@ def Write_to_io(axis, m, n, basis): f.write('# File type, either VASP or LAMMPS input \n') f.write(list(my_dict.keys())[8] + ': ' + list(my_dict.values())[8] + '\n\n\n') - f.write('# The following is your csl_generator output. YOU DO NOT NEED ' - 'TO CHANGE THEM! \n\n') + f.write('# The following is your csl_generator output.' + ' YOU DO NOT NEED TO CHANGE THEM! \n\n') f.write('axis'+': ' + str([axis[0], axis[1], axis[2]]) + '\n') f.write('m' + ': ' + str(m) + '\n') f.write('n' + ': ' + str(n) + '\n') @@ -763,7 +784,6 @@ def main(): uvw = np.array([int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3])]) uvw = CommonDivisor(uvw)[0] - if len(sys.argv) == 4: limit = 100 print(" List of possible CSLs for {} axis sorted by Sigma " @@ -792,7 +812,7 @@ def main(): sigma = int(sys.argv[5]) try: - t, m, n = get_theta_m_n_list(uvw, sigma)[0] + _, m, n = get_theta_m_n_list(uvw, sigma)[0] Write_to_io(uvw, m, n, basis) print("----------List of possible CSL planes for Sigma {}---------" @@ -814,7 +834,7 @@ def main(): try: - t, m, n = get_theta_m_n_list(uvw, sigma)[0] + _, m, n = get_theta_m_n_list(uvw, sigma)[0] lim = int(sys.argv[6]) if lim > 10: @@ -839,10 +859,3 @@ def main(): if __name__ == '__main__': main() - - - - - - - diff --git a/gb_code/gb_generator.py b/gb_code/gb_generator.py index b578937..74cc195 100644 --- a/gb_code/gb_generator.py +++ b/gb_code/gb_generator.py @@ -1,5 +1,5 @@ -#!/usr/bin/env python +# !/usr/bin/env python """ This module produces GB structures. You need to run csl_generator first to get the info necessary for your grain boundary of interest. @@ -7,8 +7,9 @@ The code runs in one mode only and takes all the necessary options to write a final GB structure from the input io_file, which is written after you run csl_generator.py. You must customize the io_file that comes with some default - values. For ex.: input the GB_plane of interest from running the CSl_generator - in the second mode. Once you have completed customizing the io_file, run: + values. For ex.: input the GB_plane of interest from running the + CSl_generator in the second mode. Once you have completed customizing the + io_file, run: 'gb_generator.py io_file' """ @@ -17,7 +18,7 @@ import numpy as np from numpy import dot, cross from numpy.linalg import det, norm -import gb_code.csl_generator as cslgen +import csl_generator as cslgen import warnings @@ -71,8 +72,12 @@ def ParseGB(self, axis, basis, LatP, m, n, gb): self.gbplane = np.array(gb) try: - self.ortho1, self.ortho2, self.Num = cslgen.Find_Orthogonal_cell( - self.basis, self.axis, self.m, self.n, self.gbplane) + self.ortho1, self.ortho2, self.Num = \ + cslgen.Find_Orthogonal_cell(self.basis, + self.axis, + self.m, + self.n, + self.gbplane) except: print(""" @@ -111,28 +116,23 @@ def WriteGB(self, overlap=0.0, rigid=False, except: print('Make sure the a and b integers are there!') sys.exit() - - xdel, ydel, x_indice, y_indice = self.Find_overlapping_Atoms() - print ("<<------ {} atoms are being removed! ------>>" - .format(len(xdel))) + self.Expand_Super_cell() + xdel, _, x_indice, y_indice = self.Find_overlapping_Atoms() + print("<<------ {} atoms are being removed! ------>>" + .format(len(xdel))) if self.whichG == "G1" or self.whichG == "g1": self.atoms1 = np.delete(self.atoms1, x_indice, axis=0) - xdel[:, 0] = xdel[:, 0] + norm(self.ortho1[:, 0]) - self.atoms1 = np.vstack((self.atoms1, xdel)) elif self.whichG == "G2" or self.whichG == "g2": self.atoms2 = np.delete(self.atoms2, y_indice, axis=0) - ydel[:, 0] = ydel[:, 0] - norm(self.ortho1[:, 0]) - self.atoms2 = np.vstack((self.atoms2, ydel)) else: print("You must choose either 'g1', 'g2' ") sys.exit() - self.Expand_Super_cell() if not self.trans: count = 0 - print ("<<------ 1 GB structure is being created! ------>>") + print("<<------ 1 GB structure is being created! ------>>") if self.File == "LAMMPS": self.Write_to_Lammps(count) elif self.File == "VASP": @@ -150,14 +150,14 @@ def WriteGB(self, overlap=0.0, rigid=False, except: print('Make sure the a and b integers are there!') sys.exit() - print ("<<------ 0 atoms are being removed! ------>>") + print("<<------ 0 atoms are being removed! ------>>") self.Expand_Super_cell() self.Translate(a, b) else: self.Expand_Super_cell() count = 0 - print ("<<------ 1 GB structure is being created! ------>>") + print("<<------ 1 GB structure is being created! ------>>") if self.File == "LAMMPS": self.Write_to_Lammps(count) elif self.File == "VASP": @@ -204,7 +204,8 @@ def CSL_Ortho_unitcell_atom_generator(self): Atoms = [] tol = 0.001 if V > 5e6: - print("Warning! It may take a very long time to produce this cell!") + print("Warning! It may take a very long time" + "to produce this cell!") # produce Atoms: for i in range(V): @@ -250,7 +251,8 @@ def CSL_Bicrystal_Atom_generator(self): self.atoms1 = dot(self.rot1, self.atoms1.T).T self.atoms2 = dot(self.rot2, self.atoms2.T).T - self.atoms2[:, 0] = self.atoms2[:, 0] - norm(Or_2[0, :]) + # tol = 0.01 + self.atoms2[:, 0] = self.atoms2[:, 0] - norm(Or_2[0, :]) # - tol # print(self.atoms2, norm(Or_2[0, :]) ) return @@ -288,10 +290,19 @@ def Find_overlapping_Atoms(self): """ finds the overlapping atoms. """ - IndX = np.where([self.atoms1[:, 0] < 1])[1] + periodic_length = norm(self.ortho1[:, 0]) * self.dim[0] + periodic_image = self.atoms2 + [periodic_length * 2, 0, 0] + # select atoms contained in a smaller window around the GB and its + # periodic image + IndX = np.where([(self.atoms1[:, 0] < 1) | + (self.atoms1[:, 0] > (periodic_length - 1))])[1] IndY = np.where([self.atoms2[:, 0] > -1])[1] - X_new = self.atoms1[self.atoms1[:, 0] < 1] - Y_new = self.atoms2[self.atoms2[:, 0] > -1] + IndY_image = np.where([periodic_image[:, 0] < + (periodic_length + 1)])[1] + X_new = self.atoms1[IndX] + Y_new = np.concatenate((self.atoms2[IndY], periodic_image[IndY_image])) + IndY_new = np.concatenate((IndY, IndY_image)) + # create a meshgrid search routine x = np.arange(0, len(X_new), 1) y = np.arange(0, len(Y_new), 1) indice = (np.stack(np.meshgrid(x, y)).T).reshape(len(x) * len(y), 2) @@ -300,7 +311,12 @@ def Find_overlapping_Atoms(self): indice_y = indice[norms < self.overD][:, 1] X_del = X_new[indice_x] Y_del = Y_new[indice_y] - return (X_del, Y_del, IndX[indice_x], IndY[indice_y]) + + if (len(X_del) != len(Y_del)): + print("Warning! the number of deleted atoms" + "in the two grains are not equal!") + # print(type(IndX), len(IndY), len(IndY_image)) + return (X_del, Y_del, IndX[indice_x], IndY_new[indice_y]) def Translate(self, a, b): @@ -311,7 +327,7 @@ def Translate(self, a, b): tol = 0.001 if (1 - cslgen.ang(self.gbplane, self.axis) < tol): - M1, M2 = cslgen.Create_minimal_cell_Method_1( + M1, _ = cslgen.Create_minimal_cell_Method_1( self.sigma, self.axis, self.R) D = (1 / self.sigma * cslgen.DSC_vec(self.basis, self.sigma, M1)) Dvecs = cslgen.DSC_on_plane(D, self.gbplane) @@ -320,8 +336,8 @@ def Translate(self, a, b): shift2 = TransDvecs[:, 1] / 2 a = b = 3 else: - #a = 10 - #b = 5 + # a = 10 + # b = 5 if norm(self.ortho1[:, 1]) > norm(self.ortho1[:, 2]): shift1 = (1 / a) * (norm(self.ortho1[:, 1]) * @@ -379,23 +395,23 @@ def Write_to_Vasp(self, trans): xlo = -1 * np.round(norm(self.ortho1[:, 0]) * dimx * self.LatP, 8) xhi = np.round(norm(self.ortho1[:, 0]) * dimx * self.LatP, 8) - LenX= xhi - xlo + LenX = xhi - xlo ylo = 0.0 yhi = np.round(norm(self.ortho1[:, 1]) * dimy * self.LatP, 8) - LenY= yhi - ylo + LenY = yhi - ylo zlo = 0.0 zhi = np.round(norm(self.ortho1[:, 2]) * dimz * self.LatP, 8) LenZ = zhi - zlo Wf = np.concatenate((X_new, Y_new)) - with open(name + plane + '_' + overD + '_' +Trans, 'w') as f: + with open(name + plane + '_' + overD + '_' + Trans, 'w') as f: f.write('#POSCAR written by GB_code \n') f.write('1 \n') f.write('{0:.8f} 0.0 0.0 \n'.format(LenX)) f.write('0.0 {0:.8f} 0.0 \n'.format(LenY)) f.write('0.0 0.0 {0:.8f} \n'.format(LenZ)) - f.write('{} {} \n'.format(len(X),len(Y))) + f.write('{} {} \n'.format(len(X), len(Y))) f.write('Cartesian\n') np.savetxt(f, Wf, fmt='%.8f %.8f %.8f') f.close() @@ -439,7 +455,7 @@ def Write_to_Lammps(self, trans): Wf = np.concatenate((W1, W2)) FinalMat = np.concatenate((Counter.T, Wf), axis=1) - with open(name + plane + '_' + overD + '_' +Trans, 'w') as f: + with open(name + plane + '_' + overD + '_' + Trans, 'w') as f: f.write('#Header \n \n') f.write('{} atoms \n \n'.format(NumberAt)) f.write('2 atom types \n \n') @@ -489,24 +505,24 @@ def main(): if overlap > 0 and rigid: gbI.WriteGB( - overlap = overlap, whichG = whichG, rigid = rigid, a = a, - b = b, dim1 = dim1, dim2 = dim2, dim3 = dim3, file = file + overlap=overlap, whichG=whichG, rigid=rigid, a=a, + b=b, dim1=dim1, dim2=dim2, dim3=dim3, file=file ) elif overlap > 0 and not rigid: gbI.WriteGB( - overlap = overlap, whichG = whichG, rigid = rigid, - dim1 = dim1, dim2 = dim2, dim3 = dim3, file = file + overlap=overlap, whichG=whichG, rigid=rigid, + dim1=dim1, dim2=dim2, dim3=dim3, file=file ) elif overlap == 0 and rigid: gbI.WriteGB( - overlap = overlap, rigid = rigid, a = a, - b = b, dim1 = dim1, dim2 = dim2, dim3 = dim3, - file = file + overlap=overlap, rigid=rigid, a=a, + b=b, dim1=dim1, dim2=dim2, dim3=dim3, + file=file ) elif overlap == 0 and not rigid: gbI.WriteGB( - overlap = overlap, rigid = rigid, - dim1 = dim1, dim2 = dim2, dim3 = dim3, file = file + overlap=overlap, rigid=rigid, + dim1=dim1, dim2=dim2, dim3=dim3, file=file ) else: print(__doc__)