The jupyter notebooks were modified, WriteGB was modified, PEP 8 was modified

 "cells": [

  {

   "cell_type": "code",

   "execution_count": 11,

   "execution_count": 3,

   "metadata": {},

   "outputs": [],

   "source": [

  },

  {

   "cell_type": "code",

   "execution_count": 3,

   "execution_count": 4,

   "metadata": {},

   "outputs": [

    {

  },

  {

   "cell_type": "code",

   "execution_count": 4,

   "execution_count": 18,

   "metadata": {},

   "outputs": [

    {

    "\n",

    "R = csl.rot(axis, theta)\n",

    "\n",

    "# Minimal CSL cells. The plane orientations and the orthogonal cells\n",

    "# will be produced from these original cells.\n",

    "M1, M2 = csl.Create_minimal_cell_Method_1(sigma,axis,R)\n",

    "\n",

    "print('Angle:',degrees(theta),'\\n','Sigma:',sigma,'\\n', 'Minimal cells:','\\n', M1,'\\n',M2,'\\n', )"

  },

  {

   "cell_type": "code",

   "execution_count": 5,

   "execution_count": 9,

   "metadata": {},

   "outputs": [],

   "source": [

  },

  {

   "cell_type": "code",

   "execution_count": 6,

   "execution_count": 10,

   "metadata": {},

   "outputs": [

    {

       "4    [2, -1, 0]    [2, 0, -1]  Mixed"

      ]

     },

     "execution_count": 6,

     "execution_count": 10,

     "metadata": {},

     "output_type": "execute_result"

    }

  },

  {

   "cell_type": "code",

   "execution_count": 7,

   "execution_count": 11,

   "metadata": {},

   "outputs": [

    {

       "2     [1, 1, 1]     [1, 1, 1]  Twist"

      ]

     },

     "execution_count": 7,

     "execution_count": 11,

     "metadata": {},

     "output_type": "execute_result"

    }

  },

  {

   "cell_type": "code",

   "execution_count": 17,

   "execution_count": 15,

   "metadata": {},

   "outputs": [],

   "source": [

  },

  {

   "cell_type": "code",

   "execution_count": 9,

   "execution_count": 16,

   "metadata": {},

   "outputs": [],

   "source": [

  },

  {

   "cell_type": "code",

   "execution_count": 18,

   "execution_count": 19,

   "metadata": {},

   "outputs": [

    {

    {

     "data": {

      "text/plain": [

       "[<mpl_toolkits.mplot3d.art3d.Line3D at 0x125393860>]"

       "[<mpl_toolkits.mplot3d.art3d.Line3D at 0x11f0d1860>]"

      ]

     },

     "execution_count": 18,

     "execution_count": 19,

     "metadata": {},

     "output_type": "execute_result"

    }

   ],

   "source": [

    "%matplotlib notebook\n",

    "ax = PlotPlane(v1, lim=4)\n",

    "\n",

    "# function to plot a line \n",

  },

  {

   "cell_type": "code",

   "execution_count": 19,

   "execution_count": 20,

   "metadata": {},

   "outputs": [

    {

       "268   [7, 0, -7]   [3, 5, -8]  Tilt"

      ]

     },

     "execution_count": 19,

     "execution_count": 20,

     "metadata": {},

     "output_type": "execute_result"

    }

  },

  {

   "cell_type": "code",

   "execution_count": 14,

   "execution_count": 21,

   "metadata": {},

   "outputs": [],

   "source": [

  },

  {

   "cell_type": "code",

   "execution_count": 15,

   "execution_count": 24,

   "metadata": {},

   "outputs": [

    {

       "(-20, 3)"

      ]

     },

     "execution_count": 15,

     "execution_count": 24,

     "metadata": {},

     "output_type": "execute_result"

    }

   "source": [

    "ax = PlotPlane(v1, lim=27)\n",

    "\n",

    "# I intended to decompose the mixed grain boundary [7, 16, -29] into two lower energy facets.\n",

    "# 1 unit vector of  = 7 units of {1 2 3} plane + 2 units of {0 1 4}. \n",

    "# I intended to decompose the mixed grain boundary (7, 16, -29) into two lower energy facets.\n",

    "# 1 unit vector of (7, 16, -29) = 7 units of (1 2 3) plane + 2 units of (0 1 4). \n",

    "\n",

    "PlotLine(dot(R,[7, 16, -29]), color='r')\n",

    "PlotLine(dot(R,[1, 2, -3]), end=7, color='r')\n",

    "PlotLine(dot(R,[0, 1, -4]), origin=dot(R,7*np.array([1, 2, -3])), end=2, color='r') \n",

    "PlotLine(dot(R,[1, 2, -3]), length=7, color='r')\n",

    "PlotLine(dot(R,[0, 1, -4]), origin=dot(R,7*np.array([1, 2, -3])), length=2, color='r') \n",

    "ax.set_xlim(-5, 20)\n",

    "ax.set_ylim(-5, 20)\n",

    "ax.set_zlim(-20, 3)\n",

  later be read by gb_geberator.py

 """

# CSL analytical formula based on

# 'Interfaces in crystalline materials', Sutton and Balluffi, clarendon press, 1996.

def get_cubic_sigma(uvw, m, n=1):

    u, v, w = uvw

    sqsum = u*u + v*v + w*w

    else:

        return pi

# produce the smallest angle for a given sigma

def get_theta_m_n_list(uvw, sigma):

    if sigma == 1:

        return [(0., 0., 0.)]

                thetas.sort(key=lambda x: x[0])

    return thetas

# produce a list of sigma for a given axis

def print_list(uvw, limit):

    for i in range(limit):

            print("Sigma:   {0:3d}  Theta:  {1:5.2f} "

                  .format(i, degrees(theta)))

# The rotation matrix

def rot(a, Theta):

    c = cos(Theta)

    s = sin(Theta)

    a = a / norm(a)

    ax = a[0]

    ay = a[1]

    az = a[2]

    ax, ay, az = a

    return array([[c + ax * ax * (1 - c), ax * ay * (1 - c) - az * s,

                   ax * az * (1 - c) + ay * s],

                  [ay * ax * (1 - c) + az * s, c + ay * ay * (1 - c),

# Helpful Functions:

#-------------------#

# make sure you have an integer array

def integer_array(A, tol=1e-7):

    return np.all(abs(np.round(A) - A) < tol)

# angle between two vectors

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):

    ang = np.round(dot(a, b)/norm(a)/norm(b), 7)

    return abs(ang)

# returns reduced vector and the common factor

def CommonDivisor(a):

    CommFac = []

    a = np.array(a)

            CommFac.append(i)

    return(a.astype(int), (abs(np.prod(CommFac))))

# returns smallest multiple integer to make an integer array

# returns the smallest multiple integer to make an integer array

# (Look at the 200 limit!)

def SmallestInteger(a):

    a = np.array(a)

    for i in range(1, 200):

            break

    return (testV, i) if integer_array(testV) else None

# array is a 2 dimensional array of vectors, returns cubic symmetric eqivalents

# arr is a 2 dimensional array of vectors, returns cubic symmetric eqivalents

# of each vector

def SymmEquivalent(array):

def SymmEquivalent(arr):

    Sym = np.zeros([24, 3, 3])

    Sym[0, :] = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

    Sym[1, :] = [[1, 0, 0], [0, -1, 0], [0, 0, -1]]

    Sym[22, :] = [[0, 0, 1], [0, -1, 0], [1, 0, 0]]

    Sym[23, :] = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]

    array = np.atleast_2d(array)

    arr = np.atleast_2d(arr)

    Result = []

    for i in range(len(Sym)):

        for j in range(len(array)):

            Result.append(dot(Sym[i, :], array[j]))

        for j in range(len(arr)):

            Result.append(dot(Sym[i, :], arr[j]))

    Result = np.array(Result)

    return (np.unique(Result, axis=0))

# returns the tilt and twist component of a given GB plane:

def Tilt_Twist_comp(v1, uvw, m, n):

    theta = get_cubic_theta(uvw, m, n)

    R = rot(uvw, theta)

# Generate GB planes and specify the character #

def Create_Possible_GB_Plane_List(uvw, m=5, n=1, lim=5):

    uvw = np.array(uvw)

    Theta = get_cubic_theta(uvw, m, n)

    R1 = rot(uvw, Theta)

    # List and character of possible GB planes:

    x = np.arange(-lim, lim + 1, 1)

    y = x

    z = x

                   [1, 0, 1],

                   ], dtype='float')

    # Find GB plane coordinates:

    for i in range(len(indice_0)):

        if (CommonDivisor(indice_0[i])[1] <= 1):

            V1[i, :] = (indice_0[i, 0] * Min_1[:, 0] +

    MeanPlanes = (V1 + V2) / 2

    # Check the type of GB plane: Symmetric tilt, tilt, twist

    for i in range(len(V1)):

        if ang(V1[i], uvw) < tol:

    return (V1, V2, MeanPlanes, GBtype)

# ------------------------Find Minimal cell by my method-----------------:

# Find Minimal cell by my method (An alternative analytical method can be used too):

def Create_minimal_cell_Method_1(sigma, uvw, R):

    uvw = np.array(uvw)

    MiniCell_1 = np.zeros([3, 3])

    MiniCell_1[:, 2] = uvw

    Found = False

    count = 0

    while (not Found) and count < len(TestVecs) - 1:

        if (1 - ang(TestVecs[count], MiniCell_1[:, 2]) > tol):

            # and  (ang(TestVecs[i],uvw) > tol):

            MiniCell_1[:, 1] = (TestVecs[count])

            count += 1

            for j in range(len(TestVecs)):

                if (1 - ang(TestVecs[j], MiniCell_1[:, 2]) > tol) and (

                        1 - ang(TestVecs[j], MiniCell_1[:, 1]) > tol):

                    if (ang(TestVecs[j],

    else:

        return (Found)

# Defining the basis

def Basis(basis):

    # Cubic basis

    if str(basis) == 'fcc':

        basis = np.array([[0, 0, 0],

                         [0.5, 0, 0.5],

    return (basis)

# Find Orthogonal cells:

# Find Orthogonal cells from the minimal cells:

def Find_Orthogonal_cell(basis, uvw, m, n, GB1):

    # Changeable limit

    lim = 15

    uvw = np.array(uvw)

    # OrthoCell_2 = OrthoCell_2

    # Test

    # OrthoCell_3 = (dot(R, OrthoCell_1))

    Volume_1 = (round(det(OrthoCell_1), 5))

    Volume_2 = (round(det(OrthoCell_2), 5))

    Num = Volume_1 * len(Basis(basis)) * 2

    else:

        return None

# print lists of GB planes:

# print lists of GB planes given an axis, basis, m and n :

def print_list_GB_Planes(uvw, basis, m, n, lim=3):

    uvw = np.array(uvw)

# ___CSL/DSC vector construction___#

# According to Grimmer, DSC and CSL lattices for bcc and fcc were made from the

# Sc lattice via body_centering and face_centering:

# According to Grimmer et al. (Acta Cryst. (1974). A30, 197-207) ,

#  DSC and CSL lattices for bcc and fcc were made from the Sc lattice

# via body_centering and face_centering:

def odd_even(M1):

    d_e = np.array([['a', 'a', 'a'],

                    ['a', 'a', 'a'],

    return (z_f if det(z_f) == 0.25 else None)

# A DSC cell

def DSC_vec(basis, sigma, minicell):

    D_sc = np.round(sigma * (inv(minicell).T), 6).astype(int)

    if basis == 'sc':

        D = dot(D_sc, face_centering(D_sc))

    return (D)

# Your CSL cell for sc, fcc and bcc. For diamond use fcc.

def CSL_vec(basis, minicell):

    C_sc = minicell.copy()

    if basis == 'sc':

        C = dot(C_sc, face_centering(C_sc))

    return (C)

# The projected DSC cell on the plane of interest.

def DSC_on_plane(D, Pnormal):

    D_proj = np.zeros((3, 3))

                        Pnormal/norm(Pnormal))

    return (D_proj)

# The density of CSL points on a given plane

def CSL_density(basis, minicell, plane):

    " returns CSL density of the given plane and its d_spacing "

    plane = np.array(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

def Ortho_fcc_bcc(basis, O1, O2):

    ortho1 = np.zeros((3, 3))

    ortho2 = np.zeros((3, 3))

        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):

    my_dict = {'GB_plane': str([axis[0], axis[1], axis[2]]), 'lattice_parameter': '4',

                'csl_generator as GB1 \n \n')

        f.write(list(my_dict.keys())[0] + ': ' + list(my_dict.values())[0] +

                '\n')

        f.write('# lattic parameter in Angstrom \n \n')

        f.write('# lattice parameter in Angstrom \n \n')

        f.write(list(my_dict.keys())[1] + ': ' + list(my_dict.values())[1] +

                '\n')

        f.write('# atoms that are closer than this fraction of the lattice '

                'on the GB_plane or not (yes or no)\n')

        f.write('# When yes, for any GB aside from twist GBs, the two inplane \n'

        '# CSL vecs will be divided by integers a and b to produce a*b initial \n'

        '# CSL vectors will be divided by integers a and b to produce a*b initial \n'

        '# configurations. The default values produce 50 initial structures \n'

        '# if you choose no for rigid_trans, you do not need to care about a and b. \n'

        '# twist boundaries are handled internally \n\n')

            print("----------List of possible CSL planes for Sigma {}---------"

                  .format(sigma))

            print(" GB1---------------GB2-------------Type----------"

            print(" GB1-------------------GB2-------------------Type----------"

                  "Number of Atoms ")

            print_list_GB_Planes(uvw, basis, m, n, lim)

            print("----------List of possible CSL planes for Sigma {}---------"

                  .format(sigma))

            print(" GB1---------------GB2-------------Type----------"

            print(" GB1-------------------GB2-------------------Type----------"

                  "Number of Atoms ")

            print_list_GB_Planes(uvw, basis, m, n, lim)

 "python gb_generator.py io_file"

 """

# A grain boundary class, encompassing all the characteristics of a GB

class GB_character:

    def __init__(self):

            sys.exit()

    def WriteGB(self, *args):

        self.overD = float(args[0])

        if self.overD > 0:

            self.whichG = str(args[1])

                    print('Make sure the input arguments are right!')

                    sys.exit()

                self.dim = np.array([int(args[2]), int(args[3]), int(args[4])])

                self.Expand_Super_cell()

                count = 0

                print ("<<------ 1 GB structure is being created! ------>>")

                self.Write_to_Lammps(count)

        else:

            print('Overlap distance is not inputted correctly!')

            print('Overlap distance is not inputted incorrectly!')

            sys.exit()

    # The fastest way of building CSL unitcells:

    def CSL_Ortho_unitcell_atom_generator(self):

        # The fastest way of building unitcells:

        Or = self.ortho.T

        LoopBound = np.zeros((3, 2), dtype=float)

        transformed = []

            self.Atoms = None

        return

    # Build the unitcells for both grains g1 and g2:

    def CSL_Bicrystal_Atom_generator(self):

        Or_1 = self.ortho1.T

        Or_2 = self.ortho2.T

        self.rot1 = np.array([Or_1[0, :] / norm(Or_1[0, :]),

                             Or_1[1, :] / norm(Or_1[1, :]),

                             Or_1[2, :] / norm(Or_1[2, :])])

        # print(self.atoms2, norm(Or_2[0, :]) )

        return

    # Expand the smallest CSL unitcell:

    def Expand_Super_cell(self):

        a = norm(self.ortho1[:, 0])

        b = norm(self.ortho1[:, 1])

        self.atoms2 = np.array(Y_new)

        return

    # find overlapping atoms:

    def Find_overlapping_Atoms(self):

        IndX = np.where([self.atoms1[:, 0] < 1])[1]

        IndY = np.where([self.atoms2[:, 0] > -1])[1]

        X_new = self.atoms1[self.atoms1[:, 0] < 1]

        Y_del = Y_new[indice_y]

        return (X_del, Y_del, IndX[indice_x], IndY[indice_y])

    # translate the GB on a mesh created by a, b integers and write to LAMMPS:

    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(

                self.atoms1 = atoms1_new

                self.Write_to_Lammps(count)

    # Write a single GB without translations to LAMMPS:

    def Write_to_Lammps(self, trans):

        name = 'input_G'

        plane = str(self.gbplane[0])+str(self.gbplane[1])+str(self.gbplane[2])

        if self.overD > 0:

def main():

    import yaml

    if len(sys.argv) == 2:

        io_file = sys.argv[1]

        file = open(io_file, 'r')