minor change.

%% Cell type:code id: tags:

``` python

from math import cos, sin, asin, acos, gcd, atan2

import numpy as np

from numpy import array, dot, degrees, cross

from numpy.linalg import inv, det, solve, norm

import pandas as pd

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# import the python scripts:

import gb_code.csl_generator as csl

import gb_code.gb_generator as gbc

%matplotlib notebook

```

%% Cell type:markdown id: tags:

## Produce Lists of CSL boundaries for any given rotation axis (hkl) :

%% Cell type:code id: tags:

``` python

# for example: [1, 0, 0], [1, 1, 0] or [1, 1, 1]

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

# list Sigma boundaries < 50

csl.print_list(axis, 50)

```

%% Output

    Sigma:     1  Theta:   0.00

    Sigma:     3  Theta:  60.00

    Sigma:     7  Theta:  38.21

    Sigma:    13  Theta:  27.80

    Sigma:    19  Theta:  46.83

    Sigma:    21  Theta:  21.79

    Sigma:    31  Theta:  17.90

    Sigma:    37  Theta:  50.57

    Sigma:    39  Theta:  32.20

    Sigma:    43  Theta:  15.18

    Sigma:    49  Theta:  43.57

%% Cell type:markdown id: tags:

## Select a sigma and get the characteristics of the GB:

%% Cell type:code id: tags:

``` python

# pick a sigma for this axis, ex: 19.

sigma = 19

theta, m, n = csl.get_theta_m_n_list(axis, sigma)[0]

R = csl.rot(axis, theta)

# Minimal CSL cells. The plane orientations and the orthogonal cells

# will be produced from these original cells.

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

print('Angle:', degrees(theta), '\n', 'Sigma:', sigma, '\n',

      'Minimal cells:', '\n', M1, '\n', M2, '\n')

```

%% Output

    Angle: 46.82644889274108

     Sigma: 19

     Minimal cells:

     [[ 0  3  1]

     [-3 -2  1]

     [ 2  0  1]]

     [[ 2  3  1]

     [-3  0  1]

     [ 0 -2  1]]

%% Cell type:markdown id: tags:

## Produce Lists of GB planes for the chosen boundary :

%% Cell type:code id: tags:

``` python

# the higher the limit the higher the indices of GB planes produced.

lim = 3

V1, V2, M, Gb = csl.Create_Possible_GB_Plane_List(axis, m, n, lim)

# the following data frame shows the created list of GB planes and their corresponding types

df = pd.DataFrame(

    {'GB1': list(V1),

     'GB2': list(V2),

     'Type': Gb

    })

df.head()

```

%% Output

                GB1           GB2   Type

    0    [-3, 2, 0]    [-3, 0, 2]  Mixed

    1    [0, -3, 2]    [2, -3, 0]  Mixed

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

    3    [3, -2, 0]    [3, 0, -2]  Mixed

    4  [-1, -1, -1]  [-1, -1, -1]  Twist

%% Cell type:markdown id: tags:

## Criteria for finding the GB plane of interest*:

###  1- Based on the type of GB plane:

%% Cell type:markdown id: tags:

#### _*The following criteria searches the generated data frame. To extend the search you can increase the limit (lim) in the above cell._

%% Cell type:code id: tags:

``` python

# Gb types: Symmetric Tilt, Tilt, Twist, Mixed

df[df['Type'] == 'Tilt'].head()

```

%% Output

                 GB1          GB2  Type

    104   [-8, 7, 1]   [-8, 1, 7]  Tilt

    113   [1, -8, 7]   [7, -8, 1]  Tilt

    116  [-1, 8, -7]  [-7, 8, -1]  Tilt

    119  [8, -7, -1]  [8, -1, -7]  Tilt

    277   [7, 1, -8]   [1, 7, -8]  Tilt

%% Cell type:markdown id: tags:

### 2 - Based on the minimum number of atoms in the orthogonal cell:

%% Cell type:markdown id: tags:

####  _This can be of interest for DFT calculations that require smaller cells. The search may take a few minutes if your original limit is large as it must calculate all the orthogonal cells to know the number of atoms._

%% Cell type:code id: tags:

``` python

basis = 'fcc'

Number = np.zeros(len(V1))

Number_atoms = []

for i in range((len(V1))):

    Number_atoms.append(csl.Find_Orthogonal_cell(basis,axis,m,n,V1[i])[2])

```

%% Cell type:code id: tags:

``` python

# show me Gb planes that have orthogonal cells with less than max_num_atoms, here: 500.

df['Number'] = Number_atoms

max_num_atoms = 300

df[df['Number'] < max_num_atoms]

```

%% Output

                  GB1           GB2            Type  Number

    2       [1, 1, 1]     [1, 1, 1]           Twist     228

    4    [-1, -1, -1]  [-1, -1, -1]           Twist     228

    6      [3, -5, 2]   [5, -3, -2]  Symmetric Tilt     228

    12    [-3, 5, -2]    [-5, 3, 2]  Symmetric Tilt     228

    104    [-8, 7, 1]    [-8, 1, 7]            Tilt     228

    113    [1, -8, 7]    [7, -8, 1]            Tilt     228

    116   [-1, 8, -7]   [-7, 8, -1]            Tilt     228

    119   [8, -7, -1]   [8, -1, -7]            Tilt     228

    277    [7, 1, -8]    [1, 7, -8]            Tilt     228

    288   [-7, -1, 8]   [-1, -7, 8]            Tilt     228

%% Cell type:markdown id: tags:

### 3- Based on proximity to a particular type boundary, for example the Symmetric tilts:

%% Cell type:markdown id: tags:

####  _This is how I created various steps on grain boundaries with vicinal orientations in the following work: __(https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.2.043601)__ See also here:_

__(https://www.mpie.de/2955247/GrainBoundaryDynamics)__

%% Cell type:code id: tags:

``` python

SymmTiltGbs = []

for i in range(len(V1)):

    if str(Gb[i]) == 'Symmetric Tilt':

        SymmTiltGbs.append(V1[i])

# Find GBs less than Delta (here 6) degrees from any of the symmetric tilt boundaries in this system

Delta = 6

Min_angles = []

for i in range(len(V1)):

    angles = []

    for j in range(len(SymmTiltGbs)):

        angles.append(csl.angv(V1[i],SymmTiltGbs[j]))

    Min_angles.append(min(angles))

```

%% Cell type:code id: tags:

``` python

df['Angles'] = Min_angles

df[df['Angles'] < Delta]

```

%% Output

                   GB1            GB2            Type  Number    Angles

    6       [3, -5, 2]    [5, -3, -2]  Symmetric Tilt     228  0.000000

    12     [-3, 5, -2]     [-5, 3, 2]  Symmetric Tilt     228  0.000000

    54      [2, 3, -5]    [-2, 5, -3]  Symmetric Tilt     988  0.000000

    57     [5, -2, -3]     [3, 2, -5]  Symmetric Tilt     988  0.000000

    70      [-5, 2, 3]    [-3, -2, 5]  Symmetric Tilt     988  0.000000

    73     [-2, -3, 5]     [2, -5, 3]  Symmetric Tilt     988  0.000000

    250  [-10, 14, -7]    [-16, 8, 5]           Mixed   26220  5.350633

    254   [10, -14, 7]   [16, -8, -5]           Mixed   26220  5.350633

    259    [8, -16, 5]  [14, -10, -7]           Mixed   26220  5.350633

    262   [-8, 16, -5]   [-14, 10, 7]           Mixed   26220  5.350633

%% Cell type:markdown id: tags:

## Select a GB plane and go on:

### You only need to pick the GB1 plane, from any of the three criteria in the cells above

%% Cell type:code id: tags:

``` python

GB_plane = [-2, -3, 5]

# lattice parameter

LatP = 4

basis = 'fcc'

```

%% Cell type:code id: tags:

``` python

# just a piece of info, how much of this mixed boundary is tilt or twist?

csl.Tilt_Twist_comp(GB_plane, axis, m, n)

```

%% Output

    Pure tilt boundary with a tilt component:  46.83

%% Cell type:code id: tags:

``` python

# instantiate a GB:

my_gb = gbc.GB_character()

# give all the characteristics

my_gb.ParseGB (axis, basis, LatP, m, n, GB_plane)

my_gb.ParseGB(axis, basis, LatP, m, n, GB_plane)

# Create the bicrystal

my_gb.CSL_Bicrystal_Atom_generator()

# Write a GB :

# by default overlap = 0.0, rigid = False,

#           dim1, dim2, dim3 = [1, 1, 1]

#           file = 'LAMMPS'

# read the io_file and/or README for more info, briefly: when you give overlap > 0

# you need to decide 'whichG': atoms to be removed from either G1 or G2.

# if rigid = True, then you need two integers a and b to make a mesh on the GB plane

# a and b can be put to 10 and 5 (just a suggestion) respectively for any GB except

# TWIST for the TWISTS the code will handle it internally regardless of your input a and b

# ex:

#my_gb.WriteGB(overlap=0.3, whichG='g1', rigid= True, a=3, b=2, dim2=5 , file='VASP')

my_gb.WriteGB(overlap=0.3, whichG='g1', dim1=3, dim2=3, dim3=2, file='VASP')

```

%% Output

    <<------ 13 atoms are being removed! ------>>

    <<------ 1 GB structure is being created! ------>>

%% Cell type:code id: tags:

``` python

# extract atom positions of the two grains:

X = my_gb.atoms1

Y = my_gb.atoms2

```

%% Cell type:code id: tags:

``` python

# 3d plot of the gb, that can be shown in any cartesian view direction:

def plot_gb(X,Y, view_dir = [0,1,0]):

    fig = plt.figure()

    ax = fig.add_subplot(111, projection='3d')

    ax.scatter(X[:,0],X[:,1],X[:,2],'o', s = 20, facecolor = 'y',edgecolor='none', alpha=0.2 )

    ax.scatter(Y[:,0],Y[:,1],Y[:,2],'o',s = 20,facecolor = 'b',edgecolor='none', alpha=0.2)

    # Show [0, 0, 0] as a red point

    ax.scatter(0,0,0,'s',s = 200,facecolor = 'r',edgecolor='none')

    ax.set_proj_type('ortho')

    ax.grid(False)

    az = degrees(atan2(view_dir[1],view_dir[0]))

    el = degrees(asin(view_dir[2]/norm(view_dir)))

    ax.view_init(azim = az, elev = el)

    ax.set_xlabel('X axis')

    ax.set_ylabel('Y axis')

    ax.set_zlabel('Z axis')

    return

```

%% Cell type:code id: tags:

``` python

%matplotlib notebook

plot_gb(X,Y,[0,1,0])

```

%% Output