Molecular Dynamics Simulation of Fracture in Quartz

This tutorial was prepared for use at a hands-on session at the Advanced Oxide Interfaces Workshop, ICTP, Trieste, May 2011.

Setting up your working environment

In this tutorial we are going to be using the QUIP (short for QUantum mechanics and Interatomic Potentials) code, which is a molecular dynamics code written in Fortran 95 and equipped with a Python interface named quippy.

As well as working with the standard unix command line we’ll also be using ipython, an interactive Python shell. In the listings below, shell commands are prefixed by a $ to identify them; everything else is a Python command which should be typed into the ipython shell. If you’re not familiar with Python don’t worry, it’s quite easy to pickup and the syntax is very similar to Fortran.

This tutorial assumes you have installed QUIP and quippy; see the Installation of QUIP and quippy section if you need some help on doing this.

Starting a fracture simulation

To maximise the amount of progress that we can make in the hands-on session, the first task is to get a fracture simulation running.

We’re going to be modelling a crack in \(\alpha\)-quartz opening on the basal plane \((0001)\), with the crack front parallel to the \([\bar{1}10]\) direction. The input structure is an XYZ file contains a \(150~\AA \times 50 \AA\) crack system with around 4000 atoms. The y axis is aligned with the \((0001)\) plane, and the z axis with the \([\bar{1}10]\) direction. The system contains a seed crack under uniaxial tension in the y direction with a strain of 23%.


Download the following input files:

  1. - the input structure, pictured above
  2. quartz_crack.xml - parameters for potential and molecular dynamics
  3. - structure of the quartz primitive cell

Before we begin, we need to make sure that the QUIP crack program is compiled by running the following commands:

$ export QUIP_ARCH=linux_x86_64_gfortran # or whichever arch is relevant
$ make Programs/crack

Then press enter to accept default settings.


The value you should assign to the variable QUIP_ARCH depends on your platform (see Installation of QUIP and quippy for more details). QUIP_ROOT refers to the directory where the the QUIP source code is located.

To make the command crack availale, copy the executable ${QUIP_ROOT}/build/${QUIP_ARCH} to a directory on your PATH, e.g. ~/bin.

Similarly, to compile eval run:

$ make QUIP_Programs/eval

It is highly recommended to change the name when copying the eval prgram to a directory on PATH to avoid a clash with the builtin eval command.

Start the simulation by running the crack program, which takes the “stem” of the input filenames as its only argument:

$ crack quartz_crack > quartz_crack.out &

Note that we’re redirecting the output to quartz_crack.out and running crack in the background. As well as outputting status information to quartz_crack.out, the simulation will periodically write the MD trajectory to a file named

You can monitor the progress of the simulation by opening the trajectory file with my modified version of AtomEye:

$ A

In AtomEye, you can move the system within the periodic boundaries to centre it in the middle of the cell by holding Shift and dragging with your left mouse button, or by pressing Shift+z. Zoom in and out by dragging with the right mouse button. Press b to toggle the display of bonds. You can use the Insert and Delete keys to move forwards or backwards through the trajectory - note the frame number in the title bar of the window. You can quit AtomEye by pressing q.

For more help on AtomEye see its web page or my notes on this modified version. In particular, you might find this startup file useful (copy to ~/.A to use).

Have a look at the output in quartz_crack.out. To see how the simulation is progressing, you can search this for lines starting with “D”:

$ grep "^D " quartz_crack.out

The first number gives the time in femtoseconds (the integration timestep of 0.5 fs is chosen to accurately sample the highest frequency phonon mode in the system), and the second and third are the instantaneous and time-averaged temperatures in Kelvin. The average temperature is computed using an exponential moving average.

After a few hundred femtoseconds, you should see the crack start to propagate in a fairly steady way, similar to the snapshot shown below (as a rough guide you can expect the simulation to progress at a rate of about 2 picoseconds per hour).


The simulation is run in the NVT ensemble, where the number of atoms N, volume V and temperature T are kept constant. The system is relatively small so a thermostat is used to regulate the temperature. The accelerating crack gives out a lot of energy, and since the thermostat used here is relatively weak, to avoid unduly influencing the dynamics, the temperature will initially rise before settling down once the crack reaches its equilbirum velocity.

Keep your crack simulation running as you proceed with the next parts of the tutorial (avoid closing the terminal in which you started the job!).

Quartz Primitive Unit Cell

While the fracture simulation is running, we will try to predict the expected behaviour of a crack in quartz according to continuum elasticity. We’ll start by performing some calculations on the quartz primitive cell using quippy.

Start ipython and import all the quippy routines:

$ ipython
In [1]: from quippy import *

Construct an \(\alpha\)-quartz primitive unit cell and save it to an XYZ file:

aq = alpha_quartz(a=4.84038097073, c=5.3285240037, u=0.464175616171,
  x=0.411742710542, y=0.278727453998, z=0.109736032769)

Here a and c are the lattice parameters (in \(\AA\)), u is the internal displacement of the silicon atoms and x, y, and z the internal coordinates of the oxygens (see image below).


The variable aq is referred to as Atoms object. It contains various attributes representing the system of interest which you can inspect from the ipython command line, e.g. to print out the lattice vectors, which are stored as columns in the lattice matrix

print aq.lattice

Or to print the atomic number and position of the first atom

print aq.z[1], aq.pos[1]

If you want to find out the syntax for a quippy function, type its name following by a question mark, e.g.

alpha_quartz ?

will print out the function signature and list of parameter types, as well as a brief description of what alpha_quartz() does. If you are interested in learning more about using QUIP and quippy, you could take a look at the online Introductory Tutorial.

Keep your ipython session open as we’ll need it later. Open the file you created with AtomEye. You should see something like this:


You can also download for comparison.


The main advantage of AtomEye is that it scales well to very larger systems, and does a pretty good job of understanding periodic boundary conditions. However, it can be a bit confusing for small cells like this. Since the quartz primitive cell contains only 9 atoms, AtomEye doubles the unit cell along the short lattice vectors (a and b) resulting in 36 atoms.

Calculating Elastic Properties

In order to predict the onset and velocity of fracture, we need to determine the values of the Young’s modulus \(E\), the Poisson ratio \(\nu\), surface energy density \(\gamma\) and the Rayleigh Wave speed \(c_R\). These are all properties of the classical potential which in this case is a short-ranged version of the TS polarisable potential (see this paper in J. Chem. Phys. for more about this potential).

Apart from \(\gamma\), all of these quantities can be obtained from the \(6 \times 6\) matrix of elastic constants \(C_{ij}\), so we will start by calculating this.

\(C\) is defined by \(\bm\sigma = C\bm\epsilon\) where \(\bm\sigma\) and \(\bm\epsilon\) are six component stress and strain vectors following the Voigt convention:

\[ \begin{align}\begin{aligned}\bm\sigma = \left( \sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{yz}, \sigma_{xz}, \sigma_{xy} \right)\\\bm\epsilon = \left( \varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, 2\,\varepsilon_{yz}, 2\,\varepsilon_{xz}, 2\,\varepsilon_{xy} \right)\end{aligned}\end{align} \]

The simplest way to calculate \(C\) with QUIP is to use the command line QUIP eval program. You will need the file you made earlier, as well as an XML file containing the parameters of the classical potential:

$ eval init_args="IP TS" param_file=TS_params.xml cij

Here init_args describes the kind of potential to use, at_file is the file containing the unit cell and param_file is the potential parameters. cij tells the eval program that we want it to calculate the elastic constants; this is done using the virial stress tensor. Have a look inside TS_params.xml to see the values which give the potential parameters controlling short range repsulsion, Yukawa-screened Coulomb interaction and dipole polarisability. For example you can see that oxygen (atomic number 8) is polarisable in this model and silicon (atomic number 14) is not.

Make a file called cij.dat containing the matrix output by eval (remove the “CIJ” at the beginning of each line). You can load this matrix into quippy using the command

C = loadtxt("cij.dat")

Extension: For trigonal crystals like \(\alpha\)-quartz, there are only actually 6 independent values in the \(C\) matrix:

\[\begin{split}C = \left( \begin{array}{cccccc} C_{11} & C_{12} & C_{13} & C_{14} & 0 & 0\\ C_{12} & C_{11} & C_{13} & -C_{14} & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ C_{14} & -C_{14} & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & C_{14} \\ 0 & 0 & 0 & 0 & C_{14} & C_{66} \end{array} \right)\end{split}\]

where \(C_{66}\) is given by \(\frac{1}{2}(C_{11} - C_{12})\). We can exploit this symmetry to get all these values using only two strain patterns: \(\epsilon_{xx}\) and \(\epsilon_{yy}+\epsilon_{zz}\).

If you like, you could try using the following quippy code to evaluate \(C\) taking the cystal symmetry into account to see how the results differ from those obtained with eval.

p = Potential('IP TS', param_filename='TS_params.xml')
C = elastic_constants(p, aq, graphics=False, sym='trigonal_low')

Which components are most different? Why do you think this is? Think about the effect of internal relaxation: compare with the values of the \(C^0_{ij}\) tensor obtained if internal relaxation is not allowed (use the c0ij option to the eval program). Why do you think some components are particularily sensitive to internal relaxation?

Young’s Modulus and Poisson Ratio

To calculate the effective Young’s modulus and Poisson ratio for our quartz fracture simulation, we need to rotate the elastic constant matrix to align it with the axes of our \((0001)[\bar{1}10]\) fracture geometry. You can create the required rotation matrix using:

R = rotation_matrix(aq, (0,0,0,1), (-1,1,0))
print R

Next we transform \(C\) using the rotation matrix, and calculate the compliance matrix \(S = C^{-1}\):

C_eff = transform_elasticity(C, R)
S_eff = numpy.linalg.inv(C_eff)

C transforms as a rank-4 tensor since the full relation between stress and strain is \(\sigma_{ij} = c_{ijkl} \epsilon_{kl}\), meaning that the transformed tensor is

\[c'_{ijkl} = \sum_{p,q,r,s = 1}^{3} R_{pi}R_{qj}R_{rk}R_{sl} c_{pqrs}\]

Finally we can work out the values of the effective Young’s modulus and Poisson ratio in the standard way from the components of the compliance matrix:

E = 1/S_eff[2,2]
print "Young's modulus", E, "GPa"

nu = -S_eff[1,2]/S_eff[2,2]
print "Poisson ratio", nu

Extension: Try rotating the quartz primitive cell directly:

aqr = transform(aq, R)

Open in AtomEye and confirm that it is oriented so that the \((0001)\) surface is aligned with the vertical (y) axis like the fracture system or, pictured below. Use eval to directly compute the elastic constant matrix of the rotated cell. How well does this matrix compare to C_eff?

Surface Energy

The file contains a 54 atom unit cell for the (0001) surface of \(\alpha\)-quartz (shown below). Use either the QUIP eval program or quippy to calculate the surface energy density \(\gamma\) predicted by our classical potential for this surface.


If you use eval, use the E command line option to get the program to calculate the potential energy of the input cell. If you use quippy, you should construct a Potential object using the XML parameters from the file TS_params.xml as shown above. You can then calculate the energy using the calc() function - see the tutorial section on Potentials and DynamicalSystems for more details.


The energy unit of QUIP are electron volt (eV), and the distance units are Angstrom (\(\AA\)).

Hint: you can use the expression

\[\gamma = \frac{E_{surf} - E_{bulk}}{2A}\]

where \(E_{surf}\) and \(E_{bulk}\) are the energies of surface and bulk configurations containing the same number of atoms, and \(A\) is the area of the open surface (with a factor of two because there are two open surfaces in this unit cell).

You should get a value for \(\gamma\) of around 3.5 J/m2.

Extension: what effect does relaxing the atomic positions have on the surface energy? (use the relax argument to the eval program, or the minim() function in quippy). What happens if you anneal the surface using the md program?

$ md pot_init_args="IP TS" params_in_file=TS_params.xml \ dt=0.5 N_steps=1000

Which surface energy do you think is more relevant for predicting the onset of fracture, relaxed or unrelaxed?

Energy Release Rate for Static Fracture

According to continuum elasticity, the strain energy release rate of an advancing crack is defined by

\[G = - \frac{\partial U_E}{\partial c}\]

where \(U_E\) is the total strain energy and \(c\) is the crack length. The well-known Griffith criteria uses an energy-balance argument to equate the critical value of \(G\) at which fracture can first occur to the energy required to create two new surfaces.

According to Griffith, we should expect crack propagation to become favourable for \(G > 2\gamma\), where \(\gamma\) is the surface energy density.

Calculating G for a thin strip

We are using the thin strip geometry illustrated below for our fracture simulations, with the top and bottom edges fixed.


The advantage of this setup is that the energy release rate G does not depend on the crack length, and can be found analytically by considering the energetics of an advancing crack.

The horizontal edges of the strip are given a uniform normal displacement \(\delta\), so the applied strain is \(\epsilon_0 = \delta / h\). Far ahead of the crack, the strip is in uniaxial tension: \(\epsilon_{yy} \to \epsilon_0\) as \(x \to \infty\).

The stress far ahead of the crack is given by \(\sigma_{0} = E' \epsilon_{0}\), and therefore the elastic energy per unit length and per unit thickness far ahead of the crack tip is

\[W_\infty = \frac{1}{2}E'\epsilon_0^2\cdot 2h = \frac{\delta^2 E'}{h}\]

where \(E'\) is the effective Young’s modulus.

Far behind the tip, the energy density is zero. Since no energy disappears through the clamped edges, if the crack is to advance by unit distance, a vertical strip of material with energy density \(W_\infty\) is effectively replaced by a strip with energy density zero.

The energy supplied to the crack tip is therefore equal to \(W_\infty\), so the energy release rate is simply

\[G = W_\infty = \frac{\delta^2 E'}{h}\]

In our simulations we will use periodic boundary conditions in the \(z\) direction, so we have plane strain loading (\(\epsilon_{zz} = 0\)), which means that the effective Young’s modulus \(E'\) is given by \(E/(1-\nu^2)\), where \(E\) is the Young’s modulus in the \(y\) relevant direction and \(\nu\) is the Poisson ratio, so finally we have

\[G = \frac{E \delta^2}{(1- \nu^2)h} = \frac{E \epsilon_0^2 h}{1 - \nu^2}\]

Use your values for the Young’s modulus, Poisson ratio and surface energy to calculate the value of \(G\) (in units of J/m2) and strain \(\epsilon\) at which a sample of \(\alpha\)-quartz with a height of 50 \(\AA\) is expected to fracture according to the continuum Griffith criterium. How does this compare to the initial strain we applied to our fracture specimen?

Changing the loading of the fracture system

Once your crack simulation should has run for a couple of picoseconds the crack should have reached its terminal velocity so you can stop the simulation (you can do this nicely creating an empty file named stop_run, or simply by killing the process).

We are going to take the current state of the simulation and rescale it homogeneously to change the applied load. We will then continue the simulation starting from the rescaled system. In this way we will be able to investigate the relationship between the loading G and the equilibrium crack velocity V.

You should find a file named in your job directory. This is a checkpoint file which contains a full snapshot of all the instantaneous positions and velocities of the atoms in the system. Load this file into quippy and rescale it as follows:

a = Atoms('')
params = CrackParams('quartz_crack.xml')

b = crack_rescale_homogeneous_xy(a, params, new_strain)

Replace new_strain with the target strain which should be between 0.15 and 0.30. If you inspect the new file you should see it’s identical to the orignal apart from the rescaling in the x and y directions (not along z since the system is periodic in that direction). Copy the bulk cell and XML parameters and start a new crack simulation:

$ cp
$ cp quartz_crack.xml quartz_crack_rescaled.xml
$ crack quartz_crack_rescaled > quartz_crack_rescaled.out &

Wait for the restarted simulation to reach a steady state and then estimate the crack velocity from looking at the XYZ file (right-click on an atom in AtomEye to print out its position; the interval between frames is 10 fs) or by plotting the crack position as a function of time by extracting lines starting with CRACK_TIP from the output file (you might find this script useful to do this; note however that the crack tip position is found using the atomic coordination numbers so if thre are defects in your cell it will not work correctly).

Energy Release Rate for Dynamic Fracture

Freund extended this approach to dynamic fracture by writing the total energy required to break bonds at the crack surface as the product of the static energy release rate \(G\) and a universal velocity-dependent function which he showed can be approximated as a linear function of the crack speed \(V\)

\[2\gamma \sim G \left( 1 - \frac{V}{c_R} \right)\]

The Rayleigh surface wave speed \(c_R\) sets the ultimate limit for the crack propagation speed. The expected velocity as a function of the loading \(G/2\gamma\) is then

\[\frac{V}{c_R} = 1 - \frac{2 \gamma}{G}\]

This is the Freund equation of motion. Calculating the Rayleigh wave speed - the speed with which elastic waves travel on a free surface - is fairly straightforward for isotropic materials. For anisotropic materials like quartz the calculation is more involved since the speed will be different for wave propagation in different crystallographic directions. For our case, the value turns out to be \(c_R \sim 9.3\) km/s using the elastic constants calculated for the short-range TS classical potential.

How does your measured value for the crack velocity compare to that predicted by the Freund equation of motion for the same G?

At the end of the session we will try to combine everybody’s results to produce a plot of \(V(G)\).

Extension tasks

  1. Look carefully at the MD trajectory from your fracture simulation in AtomEye. Do you notice anything about the bond-breaking events at the crack tip? How would you calculate the amount of charge which ends up on each of the cleaved surfaces?

  2. What is the effect of changing the simulation temperature? Trying changing the value of sim_temp in the XML parameter file to something much larger, e.g. 1000 K. Is the fracture qualitatively different? What happens if you turn off the thermostat (change the ensemble setting from NVT to NVE).

  3. The expression \(E_y = 1/S_{22}\) is an approximation and is not strictly valid for hexagonal materials. Check this approximation by performing a tensile test: take the rotated configuration aqr apply a series of small strains in the y direction, and compute the stress using the classical potential, using code based on the following:

    p = Potential('IP TS', param_filename='TS_params.xml')
    T = numpy.diag([1., 1+.eps_yy, 1.])            # matrix for strain in y
    aqr_t = transform(aqr, T)                      # apply transformation
    p.calc(aqr_t, virial=True)                     # compute virial tensor
    sig_yy = -v[2,2]/aqr_t.cell_volume()*GPA       # convert to stress

    Calculate an improved estimate for \(E_y\) by fitting a straight line to sig_yy as function of eps_yy. How good is the approximation we used above?

  4. It’s important to bear in mind that what we’ve done here is necessarily very approximate for a number of reasons:

    • This is a very small system
    • The applied loading is very large
    • The description of bond-breaking provided by classical potential is simplistic

    Which of these do you think is most important? How could these limitations be overcome in a more complete simulation?

Input files and Solutions