documentation, changes to notebooks. Fixed an error in gb_generator.

   "metadata": {},

   "outputs": [],

   "source": [

    "# the higher the limit the higher the indices of GB plnaes produced.\n",

    "# the higher the limit the higher the indices of GB planes produced.\n",

    "\n",

    "lim = 5\n",

    "\n",

#!/usr/bin/env python

import sys

import random

from math import degrees, atan, sqrt, pi, ceil, cos, acos, sin, gcd, radians

import numpy as np

from numpy import array, dot, cross

from numpy.linalg import det, norm, inv

Usage = """

 The code runs in two modes.

"""

This module is a collection of functions that produce CSL properties.

When run from the terminal, the code runs in two modes.

 First mode:

  "python CSLgenerator.py axis(u v w) [limit]" ----->  Where the u v w are the

  'python CSLgenerator.py axis(u v w) [limit]' ----->  Where the u v w are the

  indices of the rotation axis such as 1 0 0, 1 1 1, 1 1 0 and so on. The limit

  is the maximum Sigma of interest.

  (the limit by default: 100)

 Second mode:

 "python CSLgenerator.py axis(u v w) basis sigma [limit]" -----> Where basis is

 'python CSLgenerator.py axis(u v w) basis sigma [limit]' -----> Where basis is

  either fcc, bcc, diamond or sc. You read the sigma of interest from the first

  mode run. The limit here refers to CSL GB inclinations. The bigger the limit,

  the higher the indices of CSL planes.

  (the limit by default: 2)

  your chosen axis, basis and sigma will be written to an io_file which will

  later be read by gb_geberator.py

  later be read by gb_geberator.py.

"""

# CSL analytical formula based on

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

import sys

import random

from math import degrees, atan, sqrt, pi, ceil, cos, acos, sin, gcd, radians

import numpy as np

from numpy import dot, cross

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',

     Sutton and Balluffi, clarendon press, 1996.

     generates possible sigma values given an axis and two int m and n.

    """

    u, v, w = uvw

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

    sigma = m*m + n*n * sqsum

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

    """

    generates possible theta values given an axis and two int m and n.

    """

    u, v, w = uvw

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

    if m > 0:

    else:

        return pi

# produce the smallest angle for a given sigma

def get_theta_m_n_list(uvw, sigma):

    """

    returnes a list of sigma boundaries for a given sigma threshold.

    the smallest angles have been taken into account.

    """

    if sigma == 1:

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

    thetas = []

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

    return thetas

# produce a list of sigma for a given axis

def print_list(uvw, limit):

    """

    prints a list of sigmas/angles for a given axis.

    """

    for i in range(limit):

        tt = get_theta_m_n_list(uvw, i)

        if len(tt) > 0:

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

                  .format(i, degrees(theta)))

# The rotation matrix

def rot(a, Theta):

    """

    produces a rotation matrix from the angle and axis.

    """

    c = cos(Theta)

    s = sin(Theta)

    a = a / norm(a)

    ax, ay, az = a

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

    return np.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),

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

# Helpful Functions:

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

# make sure you have an integer array

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

    """

    returns an integer array.

    """

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

# angle between two vectors

def angv(a, b):

    """

    returns the angle between two vectors.

    """

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

    return round(degrees(ang), 7)

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

    return abs(ang)

# returns reduced vector and the common factor

def CommonDivisor(a):

    """

    returns the reduced vector and the common factor of vector a.

    """

    CommFac = []

    a = np.array(a)

    for i in range(2, 100):

            CommFac.append(i)

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

# returns the smallest multiple integer to make an integer array

# (Look at the 200 limit!)

def SmallestInteger(a):

    """

    returns the smallest multiple integer to make an integer array.

    """

    a = np.array(a)

    for i in range(1, 200):

        testV = i * a

            break

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

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

# of each vector

def SymmEquivalent(arr):

    """

    returns cubic symmetric eqivalents of the given 2 dimensional vector.

    """

    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]]

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

    """

    returns the tilt and twist components of a given GB plane.

    """

    theta = get_cubic_theta(uvw, m, n)

    R = rot(uvw, theta)

    v2 = np.round(dot(R, v1), 6).astype(int)

                     .format(tilt, twist))

# Generate GB planes and specify the character #

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

    """

    generate GB planes and specifies the character.

    """

    uvw = np.array(uvw)

    Theta = get_cubic_theta(uvw, m, n)

    Sigma = get_cubic_sigma(uvw, m, n)

    return (V1, V2, MeanPlanes, GBtype)

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

def Create_minimal_cell_Method_1(sigma, uvw, R):

    """

    finds Minimal cell by means of a numerical search.

    (An alternative analytical method can be used too).

    """

    uvw = np.array(uvw)

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

    MiniCell_1[:, 2] = uvw

    indice = (np.stack(np.meshgrid(x, y, z)).T).reshape(len(x) ** 3, 3)

    # remove 0 vectors and uvw from the list

    indice_0 = indice[np.where(np.sum(abs(indice), axis=1) != 0)]

    condition1 = ((abs(dot(indice_0, uvw) / norm(indice_0, axis=1) /

                  norm(uvw))).round(7))

    else:

        return (Found)

# Defining the basis

def Basis(basis):

    """

    defines the basis.

    """

    # Cubic basis

    if str(basis) == 'fcc':

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

    return (basis)

# Find Orthogonal cells from the minimal cells:

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

    """

    finds Orthogonal cells from the CSL minimal cells.

    """

    # Changeable limit

    lim = 15

    uvw = np.array(uvw)

    else:

        return None

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

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)

    for i in range(len(V1)):

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

# A DSC cell

def DSC_vec(basis, sigma, minicell):

    """

    a DSC network for the given sigma and minimal cell.

    """

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

    if basis == 'sc':

        D = D_sc.copy()

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

    """

    your CSL cell for sc, fcc and bcc.

    """

    C_sc = minicell.copy()

    if basis == 'sc':

        C = C_sc.copy()

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

    """

    projects the given DSC network on the given plane.

    """

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

    Pnormal = np.array(Pnormal)

    for i in range(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 "

    """

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

    """

    plane = np.array(plane)

    C = CSL_vec(basis, minicell)

    h = dot(C.T,plane)

    if (len(sys.argv) != 4 and len(sys.argv) != 5 and len(sys.argv) != 6 and

            len(sys.argv) != 7):

        print(Usage)

        print(__doc__)

    else:

        uvw = np.array([int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3])])

        except:

            print("""

                Careful! limit must be an integer ...

                  """, Usage)

                  """, __doc__)

    if len(sys.argv) == 6:

#!/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.

 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:

 'python gb_generator.py io_file'

 """

import sys

import numpy as np

from numpy.linalg import det, norm

import csl_generator as cslgen

Usage = """

 The code runs in one mode only and takes all necessary options to

 write a final GB structure from the input io_file, which is written

 after you run csl_generator.py. Add the GB_plane of interest from the

 CS planes list (GB1) after running csl_generator.py in second mode.

 "python gb_generator.py io_file"

 """

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

class GB_character:

    """

    a grain boundary class, encompassing all the characteristics of a GB.

    """

    def __init__(self):

        self.axis = np.array([1, 0, 0])

        self.sigma = 1

        self.trans = False

    def ParseGB(self, axis, basis, LatP, m, n, gb):

        """

        parses the GB input.

        """

        self.axis = np.array(axis)

        self.m = int(m)

        self.n = int(n)

            sys.exit()

    def WriteGB(self, *args):

        """

        parses the arguments and writes the final structure to the file.

        """

        self.overD = float(args[0])

        if self.overD > 0:

            self.whichG = str(args[1])

        else:

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

            sys.exit()

    # The fastest way of building CSL unitcells:

    def CSL_Ortho_unitcell_atom_generator(self):

        """

        populates a unitcell from the orthogonal vectors.

        """

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

        """

        builds the unitcells for both grains g1 and g2.

        """

        Or_1 = self.ortho1.T

        Or_2 = self.ortho2.T

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

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

        return

    # Expand the smallest CSL unitcell:

    def Expand_Super_cell(self):

        """

        expands the smallest CSL unitcell to the given dimensions.

        """

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

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

        c = norm(self.ortho1[:, 2])

        self.atoms2 = np.array(Y_new)

        return

    # find overlapping atoms:

    def Find_overlapping_Atoms(self):

        """

        finds the overlapping atoms.

        """

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

        """

        translates the GB on a mesh created by a, b integers and write to LAMMPS.

        """

        tol = 0.001

        if (1 - cslgen.ang(self.gbplane, self.axis) < tol):

                self.atoms1 = atoms1_new

                self.Write_to_Lammps(count)

    # Write a single GB without translations to LAMMPS:

    def Write_to_Lammps(self, trans):

        """

        write a single GB without translations to LAMMPS.

        """

        name = 'input_G'

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

        if self.overD > 0:

            gbI.WriteGB(overlap, rigid, dim1, dim2, dim3)

    else:

        print(Usage)

        print(__doc__)

    return