Commit b1fca9b3 authored by Sherri Hadian's avatar Sherri Hadian
Browse files

fixing bug in Find_overlapping_atoms + some pep8 modifications

parent 4a41879a
......@@ -31,7 +31,8 @@ from numpy.linalg import det, norm, inv
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()
#!/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 @@ to get the info necessary for your grain boundary of interest.
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 sys
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
"""
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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 @@ class GB_character:
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__)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment