re-structuring, changing README accordingly

 "cells": [

  {

   "cell_type": "code",

   "execution_count": 1,

   "execution_count": 2,

   "metadata": {},

   "outputs": [],

   "source": [

    "from mpl_toolkits.mplot3d import Axes3D\n",

    "\n",

    "# import the python scripts:\n",

    "import csl_generator as csl\n",

    "import gb_code.csl_generator as csl\n",

    "\n",

    "%matplotlib notebook"

   ]

This python package helps you create orthogonal grain boundary supercells for atomistic calculations. The code is based on the coincident site lattice (CSL) formulations for cubic materials (sc, bcc, fcc, diamond). I intend to extend it to hcp structures soon. The code produces a final structure to be read in [LAMMPS](https://lammps.sandia.gov/) or [VASP](https://www.vasp.at/).

# Overview

There are two main scripts: [_csl_generator.py_](./csl_generator.py) and [_gb_generator.py_](./csl_generator.py) which you need to use in this order to produce the final grain boundary (GB) structure.

There are two main scripts: [_csl_generator.py_](./csl_generator.py) and [_gb_generator.py_](./csl_generator.py) which you need to use in this order to produce the final grain boundary (GB) structure. These scripts are both modules (a collection of functions/classes) and can be executed

from the command line.

In this README I will explain the steps to use the code in the Linux Terminal and I have also attached two _jupyter notebooks_ ([Usage_of_GB_code.ipynb](./Usage_of_GB_code.ipynb), [Dichromatic_pattern_CSL.ipynb](./Dichromatic_pattern_CSL.ipynb)) which describe how the code can be accessed and used in the notebooks by various examples. These notebooks have extra functionality. The former is for the general usage of the code with some tips to locate GBs of interest, the latter depicts how CSL construction can be used for different purposes.

You can use [this link](https://mybinder.org/v2/gh/oekosheri/GB_code/master) for an interactive Jupyter notebook environment provided by Binder. [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/oekosheri/GB_code/master)

To use it locally, you will need python3 and numpy1.14 for the main scripts and additionally matplotlib and pandas to use the auxilliary Jupyter notebooks. For installation simply clone or download the code in your terminal.

To use it locally, you will need python3 and numpy1.14 for the main scripts and additionally matplotlib and pandas to use the auxilliary Jupyter notebooks.

# Installation guide

For installation simply clone or download the code in your terminal and in the main directory of the package type:

```

> pip install .

```

This will copy the modules to your active python site-packages, thereby making them importable in any python script and will put the scrpits in the python bin, thereby making them executable in the shell.

# Usage

To pick a grain boundary 5 degrees of freedom need to be fixed: rotation axis, rotation angle and GB plane orientation.

  _ex:_

```

> python csl_generator.py 1 1 1 50

> csl_generator.py 1 1 1 50

   List of possible CSLs for [1 1 1] axis sorted by Sigma

Sigma:     1  Theta:   0.00

  _ex:_

```

> python csl_generator.py 1 1 1 diamond 13

> csl_generator.py 1 1 1 diamond 13

----------List of possible CSL planes for Sigma 13---------

 GB1-------------------GB2-------------------Type----------Number of Atoms

As an example, we change the default gb_plane to [2,  1, -2], overlap_distance to 0.3 and rigid_trans to 'yes' in the io_file.

To produce the GB of interest we go on with: [_gb_generator.py_](./gb_generator.py)

```

> python gb_generator.py io_file

> gb_generator.py io_file

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

<<------ 50 GB structures are being created! ------>>

```

%% 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 csl_generator as csl

import gb_generator as gbc

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)

# 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

#!/usr/bin/env python

# This files marks the containing directory as Python package.

 No newline at end of file