# Step 1: Setup of the Silicon model system

The first task in this tutorial is to build the model system we will use for both the classical and QM/MM simulations. We will use an approximately 2D model system in the thin strip geometry

## 1.1 Building the bulk unit cell (30 minutes)

### Import the relevant modules and functions

We start by importing all the functions we will need. Create a new
script named `make_crack.py`

and add the following lines:

```
from ase.lattice import bulk
from ase.lattice.cubic import Diamond
from ase.constraints import FixAtoms
import ase.units as units
from quippy import set_fortran_indexing
from quippy.potential import Potential, Minim
from quippy.elasticity import youngs_modulus, poisson_ratio
from quippy.io import write
from quippy.crack import (print_crack_system,
G_to_strain,
thin_strip_displacement_y,
find_crack_tip_stress_field)
```

Note that some routines come from ASE and others from quippy. We will use ASE for basic atomic manipulations, and quippy to provide the interatomic potentials plus some special purpose functionality.

Note

For interactive use, it is convenient to import everything from the
entire quippy package with `from qlab import *`

as described
in the Practical considerations section. We chose not to do that in these scripts to
make it clear where each function we are using is defined, and to make it easier
to look them up in the online documentation.

### Definition of the simulation parameters

Let’s first define the parameters needed to construct our model system. There are three possible crack systems. For now, we will use the first (uncommented) one, \((111)[01\bar{1}]\), which means a crack propagating on the \((111)\) cleavage plane (the lowest surface energy of all silicon surfaces) with the crack front along the \([01\bar{1}]\) direction:

```
# System 1. (111)[0-11]
crack_direction = (-2, 1, 1) # Miller index of x-axis
cleavage_plane = (1, 1, 1) # Miller index of y-axis
crack_front = (0, 1, -1) # Miller index of z-axis
# # System 2. (110)[001]
# crack_direction = (1,-1,0)
# cleavage_plane = (1,1,0)
# crack_front = (0,0,1)
# # System 3. (110)[1-10]
# crack_direction = (0,0,-1)
# cleavage_plane = (1,1,0)
# crack_front = (1,-1,0)
```

If you have time later, you can come back to this point and change to one of the other fracture systems. Next we need various geometric parameters:

```
width = 200.0*units.Ang # Width of crack slab
height = 100.0*units.Ang # Height of crack slab
vacuum = 100.0*units.Ang # Amount of vacuum around slab
crack_seed_length = 40.0*units.Ang # Length of seed crack
strain_ramp_length = 30.0*units.Ang # Distance over which strain is ramped up
initial_G = 5.0*(units.J/units.m**2) # Initial energy flow to crack tip
```

Note the explicit unit conversion: some of this is unnecessary as we
are using the ase.units module where ```
Ang = eV =
1
```

. The energy release rate initial_G is given in the
widely used units of J/m^{2}.

Next we define some parameters related to the classical interatomic potential:

```
relax_fmax = 0.025*units.eV/units.Ang # Maximum force criteria
param_file = 'params.xml' # XML file containing
# interatomic potential parameters
mm_init_args = 'IP SW' # Initialisation arguments
# for the classical potential
```

And finally the output file:

```
output_file = 'crack.xyz' # File to which structure will be written
```

You should download the `params.xml`

file, which contains
the parameters for the SW potential (and also for DFTB, needed for
Step 3: LOTF hybrid MD simulation of fracture in Si)

### Finding the equilibrium lattice constant for Si

To find the Si equilibrium lattice constant a0 with the SW potential,
let’s first build the 8-atom diamond cubic cell for silicon, with an initial
guess at lattice constant of 5.44 A. This can be done using the
`bulk()`

function from the `ase.structure`

module:

```
si_bulk = ... # Build the 8-atom diamond cubic cell for Si
```

The variable si_bulk is an `Atoms`

object. It
has various attributes and methods that will be introduced as necessary
during this tutorial.

Once you have created your si_bulk object, run the `make_crack.py`

script from within ipython with the `run`

command. Providing you
have imported everything from the `qlab`

module, will then be
able to interactively visualise the Si unit cell with the
`view()`

function from the `qlab`

module, which you
should type in at the ipython prompt:

```
In [5]: view(si_bulk)
```

This will pop up an AtomEye [Li2003] window showing the 8-atom
silicon cell, with the unit cell boundary drawn with a thick black
line. You can rotate the system with the left mouse button, translate
by holding Control and tracking, or translate within the periodic
boundaries by holding Shift and dragging. Zoom in and out by
dragging with the right mouse button (or scroll wheel, if you have
one). Press b to toggle the display of bonds. For more help on
AtomEye see its web page or the documentation
for the `qlab`

and `atomeye`

modules.

Now, we initialise the Stillinger-Weber (SW) classical interatomic
potential using quippy’s `Potential`

class

```
mm_pot = Potential('IP SW', param_filename='params.xml')
```

The equilibrium lattice constant a0 can now be found by minimising the
cell degrees of freedom with respect to the virial tensor calculated by the
SW potential. First, we need to attach a calculator (i.e. the SW
potential, mm_pot we just created) to the si_bulk object,
using the method `set_calculator()`

:

```
si_bulk. ... # Attach the SW potential to si_bulk
```

This means that subsequent requests to calculate energy or forces of si_bulk will be performed using our SW potential.

The minimisation can now be carried out by making a
`Minim`

class from the si_bulk Atoms,
requesting that both atomic positions and cell degrees of freedom
should be relaxed. Then run the minimisation until the maximum force
is below `fmax=1e-2`

, using the `run()`

method

```
minim = ... # Initialise the minimiser from si_bulk
print('Minimising bulk unit cell')
minim. ... # Run the minimisation
```

The lattice constant a0 can be easily obtained from the relaxed
lattice vectors using the `cell()`

attribute of
the si_bulk object, which returns a \(3 \times 3\) matrix
containing the lattice vectors as rows in Cartesian coordinates,
i.e. `si_bulk.cell[0,0]`

is the x coordinate of the first lattice
vector.

```
a0 = ... # Get the lattice constant
print('Lattice constant %.3f A\n' % a0)
```

As a check, you should find a value for a0 of around 5.431 A.

Once you have obtained a0, you should replace the si_bulk object with a new bulk cell using this lattice constant, so that the off-diagonal components of the lattice are exactly zero:

```
si_bulk = ... # Make a new 8-atom bulk cell with correct a0
si_bulk. ... # re-attach the SW potential as a calculator
```

### Milestone 1.1

At this point your script should look something like `make_crack_1.py`

.

## 1.2 Calculation of elastic and surface properties of silicon (30 minutes)

### Calculation of the Young’s modulus and the Poisson ratio

Following the discussion above section, we need to
calculate some elastic properties of our model silicon. To calculate the Young’s
modulus E along the direction perpendicular to the cleavage plane, and the
Poisson ratio \(\nu\) in the \(xy\) plane, we need the \(6 \times
6\) matrix of the elastic constants \(C_{ij}\). This matrix c can be
calculated using the `get_elastic_constants()`

method of the mm_pot Potential object.

```
c = mm_pot. ... # Get the 6x6 C_ij matrix
print('Elastic constants (GPa):')
print((c / units.GPa).round(0))
print('')
```

Here, the `GPa`

constant from the ase.units module module is used to
convert from pressure units of eV/A^{3} into GPa.

The Young’s modulus E and the Poisson ratio nu can now be calculated,
given c, the cleavage_plane and the crack_direction (defined in the
parameters section above), using the functions
`youngs_modulus()`

and
`poisson_ratio()`

from the
`quippy.elasticity`

module.

```
E = ... # Get E
print('Young\'s modulus %.1f GPa' % (E / units.GPa))
nu = ... # Get nu
print('Poisson ratio % .3f\n' % nu)
```

As a check, for the \((111)[01\bar{1}]\) crack system, you should get a Young’s modulus of 142.8 GPa and a Poisson ratio of 0.265.

### Calculation of the surface energy of the cleavage plane

To calculate the surface energy per unit area gamma of the cleavage_plane, we build a Si slab unit cell aligned with the requested crystallographic orientation. The orientation of the crack system can be printed using the following command:

```
print_crack_system(crack_direction, cleavage_plane, crack_front)
```

The new unit slab can be obtained using the
`ase.lattice.cubic.Diamond`

from the `ase.lattice`

module, which is used as follows:

```
unit_slab = Diamond(directions=[crack_direction,
cleavage_plane,
crack_front],
size=(1, 1, 1),
symbol='Si',
pbc=True,
latticeconstant=a0)
print('Unit slab with %d atoms per unit cell:' % len(unit_slab))
print(unit_slab.cell)
print('')
```

You can visualise the new cell with `view(unit_slab)`

(type this at
the ipython prompt after running the script as it is so far, don’t
add it to the script file):

We now shift the unit_slab vertically so that we will open up a surface along a \((111)\) glide plane, cutting vertically aligned bonds (see e.g. this image). This choice gives the lowest energy surface. We then map the positions back into the unit cell:

```
unit_slab.positions[:, 1] += (unit_slab.positions[1, 1] -
unit_slab.positions[0, 1]) / 2.0
unit_slab.set_scaled_positions(unit_slab.get_scaled_positions())
```

The `positions`

is a (N,3) array containing
the Cartesian coordinates of the atoms, and
`set_scaled_positions()`

and
`get_scaled_positions()`

are necessary to ensure
all the atoms are mapped back inside the unit cell before we open
up a surface. This is the result of applying the shift (do another
`view(unit_slab)`

to update your AtomEye viewer).

Note how the top and bottom layers now correspond to \((111)\) glide planes, so that the cell boundary now corresponds to a shuffle plane as required.

We now make a copy of the unit_slab and create a surface unit cell
with surfaces parallel to the cleavage_plane. We can use the
`ase.atoms.Atoms.center()`

method which, besides centring the
atoms in the unit cell, allows some vacuum to be added on both sides
of the slab along a specified axis (use `axis=0`

for the x-axis,
or `axis=1`

for the y-axis, etc.). The amount of vacuum you add is
not critical, but could be taken from the vacuum parameter in the
parameters section above:

```
surface = unit_slab.copy()
surface. ... # Add vacuum along y axis
```

You should get a surface unit cell which looks something like this:

Here, the atoms have been coloured by coordination by pressing the k key. The green atoms on the surfaces are three-fold coordinated.

Now that we have both the bulk unit slab and the surface unit cell,
the surface energy gamma for the cleavage plane can be calculated
using the SW potential. Once a calculator (e.g. mm_pot) is attached
to an `Atoms`

object, the potential energy of the
atomic system can be calculated with
`get_potential_energy()`

. It is useful to know
that the number of atoms in an Atoms object can be obtained by the
list-method len (e.g. len(si_bulk) gives the number of atoms in
si_bulk), and that the volume of a cell can be calculated with
`get_volume()`

:

```
surface. ... # Attach SW potential to surface atoms
E_surf = ... # Get potential energy of surface system
E_per_atom_bulk = ... # Get potential energy per atom for bulk slab
area = ... # Calculate surface area using volume and cell
gamma = ... # Calculate surface energy
print('Surface energy of %s surface %.4f J/m^2\n' %
(cleavage_plane, gamma / (units.J / units.m ** 2)))
```

As a check, you should obtain \(\gamma_{(111)}\) = 1.36 J/m^{2}. You may want to verify that this result is converged
with respect to the number of layers in the system (note the cutoff
distance of the SW potential, which you can obtain with
`mm_pot.cutoff()`

, is about 3.93 A, just beyond the second neighbour
distance).

### Milestone 1.2

At this point your script should look something like `make_crack_2.py`

## 1.3 Setup of the crack slab supercell (30 minutes)

### Replicating the unit cell to form a slab supercell

Now, we have all the ingredients needed to build the full crack slab system and to apply the requested strain field.

We start by building the full slab system. First, we need to find the number
of unit_slab cells along x and y that approximately match width and
height (see parameters section).
Note that the python function `int()`

can be used to
convert a floating point number into an integer, truncating towards zero:

```
nx = ... # Find number of unit_slab cells along x
ny = ... # Find number of unit_slab cells along y
```

To make sure that the slab is centered on a bond along the y direction, the number of units cell in this direction, ny, must be even:

```
if ny % 2 == 1:
ny += 1
```

The crack supercell is now simply obtained by replicating unit_slab \(nx \times ny \times 1\) times along the three axes:

```
crack_slab = unit_slab * (nx, ny, 1)
```

As we did before for the surface system, vacuum has to be introduced along
the x and y axes (*Hint:* use the `center()`

method twice)

```
crack_slab. ... # Add vacuum along x
crack_slab. ... # Add vacuum along y
```

The crack_slab is now centered on the origin in the xy plane to make it simpler to apply strain:

```
crack_slab.positions[:, 0] -= crack_slab.positions[:, 0].mean()
crack_slab.positions[:, 1] -= crack_slab.positions[:, 1].mean()
```

and its original width and height values are saved, and will later be used to measure the strain:

```
orig_width = (crack_slab.positions[:, 0].max() -
crack_slab.positions[:, 0].min())
orig_height = (crack_slab.positions[:, 1].max() -
crack_slab.positions[:, 1].min())
print(('Made slab with %d atoms, original width and height: %.1f x %.1f A^2' %
(len(crack_slab), orig_width, orig_height)))
```

The original y coordinates of the top and bottom of the slab and the original x coordinates of the left and right surfaces are also saved:

```
top = crack_slab.positions[:, 1].max()
bottom = crack_slab.positions[:, 1].min()
left = crack_slab.positions[:, 0].min()
right = crack_slab.positions[:, 0].max()
```

At this point, your crack_slab should look something like this:

### Setting constraints to fix the edge atoms

During the MD simulations, the positions of the top and bottom rows of atoms will be kept fixed. More precisely, these rows of atoms will only be moved rigidly when the strain is applied and will not move in response to forces from the interatomic potential (see the discussion of the thin strip geometry above). To do this, we initialise a fixed_mask array that is True for each atom whose position needs to be fixed, and False otherwise:

```
fixed_mask = ((abs(crack_slab.positions[:, 1] - top) < 1.0) |
(abs(crack_slab.positions[:, 1] - bottom) < 1.0))
```

Note that the `|`

operator is shorthand for a logical ‘or’
operation. After re-running the latest version of your script and
executing `view(crack_slab)`

, you can colour the atoms by
fixed_mask using the `aux_property_coloring()`

function

```
aux_property_coloring(fixed_mask)
```

which colours the atoms where fixed_mask is True in red and those where it is False in blue, like this:

Now we can use the `FixAtoms`

class to
fix the positions of the atoms according to the mask fixed_mask, and
then attach the constraint to crack_slab using
`set_constraint()`

:

```
const = ... # Initialise the constraint
crack_slab. ... # Attach the constraint to crack_slab
print('Fixed %d atoms\n' % fixed_mask.sum())
```

To create the crack seed, we now apply the initial strain ramp. First,
we need to convert the chosen energy release rate initial_G into a
strain. This can be done using the `G_to_strain()`

function which implements the thin strip equation described
above. The strain is then used to displace
the y coordinate of the atomic positions according to the strain
ramp produced by the `thin_strip_displacement_y()`

function. Here, the crack_seed_length and the strain_ramp_length
parameters should be used. The objective is that atoms to the left of
`left + crack_seed_length`

should be rigidly shifted vertically, and
those to the right of ```
left + crack_seed_length +
strain_ramp_length
```

should be uniformly strained, with a transition
region in between.

```
strain = ... # Convert G into strain
crack_slab.positions[:, 1] += ... # update the atoms positions along y
print('Applied initial load: strain=%.4f, G=%.2f J/m^2' %
(strain, initial_G / (units.J / units.m**2)))
```

This is the resulting crack slab, for the \((111)\) case:

### Relaxation of the crack slab

To obtain a good starting point for the MD, we need to perform an
approximate geometry optimisation of the slab, keeping the top and
bottom rows of atoms fixed. Once again, our mm_pot needs to be
attached to crack_slab and the minimiser
`Minim`

initialised (note that here it does
not make sense to relax the cell since we have vacuum in two
directions). We can then perform the minimisation until the maximum
force is below the relax_fmax defined in the parameters
section:

```
print('Relaxing slab')
crack_slab. ... # Attach the calculator to crack_slab
minim = ... # Initialise the minimiser
minim. ... # Run the minimisation until forces are relax_fmax
```

Here’s what your minimised crack slab should look like:

### Locating the crack tip

Before starting the next steps, it is useful to find the initial
position of the crack tip. This is provided by the
`find_crack_tip_stress_field()`

function:

```
crack_pos = find_crack_tip_stress_field(crack_slab, calc=mm_pot)
print 'Found crack tip at position %s' % crack_pos
```

This function works by fitting the components of the Irwin crack stress field to the per-atom stresses calculated by the classical SW potential, allowing the origin of the analytical stress field to move during the fit. Then, we simply set this point to be the current crack position.

### Saving the output file

Finally, we can save all the calculated materials properties inside the
crack_slab `Atoms`

object, before writing it to disk:

```
crack_slab.info['nneightol'] = 1.30 # set nearest neighbour tolerance
crack_slab.info['LatticeConstant'] = a0
crack_slab.info['C11'] = c[0, 0]
crack_slab.info['C12'] = c[0, 1]
crack_slab.info['C44'] = c[3, 3]
crack_slab.info['YoungsModulus'] = E
crack_slab.info['PoissonRatio_yx'] = nu
crack_slab.info['SurfaceEnergy'] = gamma
crack_slab.info['OrigWidth'] = orig_width
crack_slab.info['OrigHeight'] = orig_height
crack_slab.info['CrackDirection'] = crack_direction
crack_slab.info['CleavagePlane'] = cleavage_plane
crack_slab.info['CrackFront'] = crack_front
crack_slab.info['strain'] = strain
crack_slab.info['G'] = initial_G
crack_slab.info['CrackPos'] = crack_pos
crack_slab.info['is_cracked'] = False
```

We can save our results, including all the extra properties and information, in extendedxyz in the output_file, whose name is defined in the parameters section:

```
print('Writing crack slab to file %s' % output_file)
write(output_file, crack_slab)
```

### Milestone 1.3

At this point your final script should look something like
Step 1 solution — make_crack.py, and your XYZ file like `crack.xyz`

.

When you are ready, proceed to Step 2: Classical MD simulation of fracture in Si.