Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 26-refactoring-the-infrastructure
  • MPGCQM_workshop
  • app
  • backup-05-08
  • develop
  • kappa_sigma_notebook
  • master
  • staging
8 results

Target

Select target project
  • nomad-lab/analytics
1 result
Select Git revision
  • 26-refactoring-the-infrastructure
  • MPGCQM_workshop
  • app
  • backup-05-08
  • develop
  • kappa_sigma_notebook
  • master
  • staging
8 results
Show changes
View open merge request
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4236 deletions
3rdparty/quip/bin/module_name deleted 100755 → 0
View file @ 23bcf97f
#!/bin/bash
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2010.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
usage="process module name to appropriate module filename case and suffix
reads module names from stdin
Usage: $0 [ -lc | -uc ] [ -suffix .suf ]"
lc=0
uc=0
suffix=".mod"
uc_set=0
lc_set=0
while [[ $# > 0 ]]; do
if [[ $1 == "-uc" ]]; then
lc=0
uc=1
uc_set=1
else
if [[ $1 == "-lc" ]]; then
lc=1
uc=0
lc_set=1
else
if [[ $1 == "-suffix" ]]; then
shift
suffix=$1
else
echo "$usage" 1>&2
exit 1
fi
fi
fi
shift
done
if [[ $lc_set == 1 && $uc_set == 1 ]]; then
echo "-lc and -uc conflict" 1>&2
exit 1
fi
SED_ARGS="s/\$/$suffix/"
TR_ARGS="a a"
if [[ $lc == 1 ]]; then
TR_ARGS="[A-Z] [a-z]"
fi
if [[ $uc == 1 ]]; then
TR_ARGS="[a-z] [A-Z]"
fi
# echo "args $lc $uc $suffix" 1>&2
shift
for name in `cat`; do
echo $name | tr $TR_ARGS | sed "${SED_ARGS}"
done
3rdparty/quip/bin/protect_mod deleted 100755 → 0
View file @ 23bcf97f
#!/bin/bash
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2010.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
usage="make copies of module files defined in a source file
Usage: $0 [ -lc | -uc ] [ -suffix suffix ] -- filename [ filename2 ... ]"
mydir=`dirname $0`
if [[ $# -lt 1 ]]; then
echo "$usage" 1>&2
exit 1
fi
args=""
while [[ $# -gt 0 && $1 != "--" ]]; do
args="$args $1"
shift
done
# echo "args $args" 1>&2
if [[ $1 != "--" ]]; then
echo "$usage" 1>&2
exit 1
fi
shift
for file in $*; do
modules=`egrep -i '^[ ]*module[ ]+[^ ]+[ ]*$!?' $file | awk '{print $2}'`
for mod in `echo $modules | $mydir/module_name $args`; do
if [ -f $mod ]; then
mv $mod $mod.save
fi
done
done
3rdparty/quip/bin/unprotect_mod deleted 100755 → 0
View file @ 23bcf97f
#!/bin/bash
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2010.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
usage="reuse saved copies of module files created by a source file if they are identical to new module file
environment: OD_SKIP = number of bytes to skip in module file header
Usage: $0 [ -lc | -uc ] [ -suffix suffix ] -- filename [ filename2 ... ]" 1>&2
mydir=`dirname $0`
if [[ $# -lt 1 ]]; then
echo "$usage" 1>&2
exit 1
fi
args=""
while [[ $# -gt 0 && $1 != "--" ]]; do
args="$args $1"
shift
done
# echo "args $args" 1>&2
if [[ $1 != "--" ]]; then
echo "$usage" 1>&2
exit 1
fi
shift
for file in $*; do
modules=`egrep -i '^[ ]*module[ ]+[^ ]+[ ]*$!?' $file | awk '{print $2}'`
for mod in `echo $modules | $mydir/module_name $args`; do
if [ -f $mod.save ]; then
cp $mod $mod.txt
cp $mod.save $mod.save.txt
if [[ ! -z $OD_SKIP ]] ; then
od -t x -j $OD_SKIP < $mod > $mod.txt
od -t x -j $OD_SKIP < $mod.save > $mod.save.txt
fi
if [[ ! -z $RE_SKIP ]] ; then
egrep -v "$RE_SKIP" < $mod > $mod.txt
egrep -v "$RE_SKIP" < $mod.save > $mod.save.txt
fi
diff -q $mod.txt $mod.save.txt > /dev/null
out=$?
if [[ $out == 0 ]]; then
mv $mod.save $mod
else
rm -f $mod.save
fi
# rm -f $mod.txt $mod.save.txt
fi
done
done
3rdparty/quip/bin/update-authors deleted 100755 → 0
View file @ 23bcf97f
#!/usr/bin/env python
import sys, os
# obtain list of users with access to our repo
group = os.popen('ssh cvs.tcm.phy.cam.ac.uk group -l jrk33').read()
users = group.split(':')[1].split()
users = set(users) # remove duplicates
# download cvs-users mailing list from TCM
cvs_users = os.popen('ssh jrk33@pc41.tcm.phy.cam.ac.uk grep "^#%" -A1 ~lists/cvs-users').read().split('\n')
# Extract usernames and email addresses from cvs-users mailing list
emails = {}
for user, email in zip(cvs_users[::3], cvs_users[1::3]):
user = user[2:].strip()
email = email.strip()
if user in users:
emails[user] = email
# Fill in missing users with @tcm email address
for user in users:
if user not in emails:
emails[user] = '%s@tcm.phy.cam.ac.uk' % user
# Print two columns: user email
for user, email in emails.iteritems():
print '%-15s%s' % (user, email)
3rdparty/quip/bin/validate.py deleted 100755 → 0
View file @ 23bcf97f
#!/usr/bin/env python
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2010.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
from xml.parsers.xmlproc import xmlproc
from xml.parsers.xmlproc import xmlval
from xml.parsers.xmlproc import xmldtd
import sys
p = xmlval.XMLValidator()
p.parse_resource(sys.argv[1])
print 'No errors found.'
3rdparty/quip/doc/Makefile deleted 100644 → 0
View file @ 23bcf97f
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
SPHINXPROJ = quippy
SOURCEDIR = .
BUILDDIR = _build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
3rdparty/quip/doc/Tutorials/Introduction.ipynb deleted 100644 → 0
View file @ 23bcf97f
%% Cell type:markdown id: tags:
# Getting started with `Atoms`
%% Cell type:markdown id: tags:
This example will walk you through getting started with using `quippy` by playing with a small molecule. It complements the (older) [introductory tutorial](http://libatoms.github.io/QUIP/tutorial.html) which is longer and goes more in depth, though it is recommended that new users start here first. This example is especially recommended for those who want to get going and working with `Atoms` as quickly as possible.
%% Cell type:markdown id: tags:
First, make sure you have QUIP and `quippy` properly installed. If you can run the cell below, then congratulations! If not, then see the [Installation instructions](http://libatoms.github.io/QUIP/install.html).
%% Cell type:code id: tags:
``` python
import numpy as np
import quippy
from quippy import Atoms
```
%% Cell type:markdown id: tags:
Remember, this is an interactive tutorial. You are encouraged to run the cells yourself and even try editing the code to see what else you can do!
%% Cell type:markdown id: tags:
Now let's try getting our hands on an `Atoms` object; this object is required to do almost anything useful in `quippy`.
%% Cell type:code id: tags:
``` python
struct = quippy.Atoms('methane.xyz')
```
%% Cell type:markdown id: tags:
OK, now what do we do with it? Let's look at the Atoms source file and see what sort of information is encoded in it:
%% Cell type:code id: tags:
``` python
!cat methane.xyz
```
%% Output
5
cutoff=-1.00000000 nneightol=1.20000000 pbc="T T T" Lattice="10.00000000 0.00000000 0.00000000 0.00000000 10.00000000 0.00000000 0.00000000 0.00000000 10.00000000" Properties=species:S:1:pos:R:3:Z:I:1
C 0.00000000 0.00000000 0.00000000 6
H 0.00000000 0.00000000 1.07000000 1
H 0.62414000 0.79255000 -0.35666000 1
H -0.99844000 0.14425000 -0.35667000 1
H 0.37430000 -0.93680000 -0.35667000 1
%% Cell type:markdown id: tags:
We have positions, atomic symbols, atomic numbers, and a unit cell. There's also some other information in there that we don't need right now. But let's see how to access these properties from `quippy`:
%% Cell type:code id: tags:
``` python
print("Positions:")
print(struct.get_positions())
print("Unit cell:")
print(struct.get_cell())
print("Atomic numbers:")
print(struct.get_atomic_numbers())
```
%% Output
Positions:
[[ 0. 0. 0. ]
[ 0. 0. 1.07 ]
[ 0.62414 0.79255 -0.35666]
[-0.99844 0.14425 -0.35667]
[ 0.3743 -0.9368 -0.35667]]
Unit cell:
[[ 10. 0. 0.]
[ 0. 10. 0.]
[ 0. 0. 10.]]
Atomic numbers:
[6 1 1 1 1]
%% Cell type:markdown id: tags:
(Remember, units in quippy are always Ångstrøms and electronvolts: http://libatoms.github.io/QUIP/units.html)
%% Cell type:markdown id: tags:
You might want get an idea of what this structure actually looks like. The native `xyz` format is supported by many open-source molecular viewers, like [VMD](http://www.ks.uiuc.edu/Research/vmd/) and [Avogadro](https://avogadro.cc/). If you've gotten the `atomeye` module to work you can use that as well.
We're still looking for a simple way to view molecular structures in the notebook environment; once we (the developers) get that working with quippy we'll update this tutorial.
%% Cell type:markdown id: tags:
Let's try another way to get information about the structure: This is a molecule, so we can check the bonds and angles.
%% Cell type:code id: tags:
``` python
struct.get_distance(0, 1, mic=True)
```
%% Output
1.0700000000000001
%% Cell type:markdown id: tags:
OK, that's a pretty reasonable C-H bond length in a methane molecule. What about the other ones?
%% Cell type:code id: tags:
``` python
[struct.get_distance(0, i, mic=True) for i in range(1, 5)]
```
%% Output
[1.0700000000000001,
1.0699965409757173,
1.0700018621479124,
1.0700038406005841]
%% Cell type:markdown id: tags:
_(Oh, and in case you're wondering about the `mic=True` part, that just means to use the minimum image convention. It's only really relevant when we use periodic systems, but it's included here for completeness.)_\n\nAnd how about an H-C-H angle?
%% Cell type:code id: tags:
``` python
struct.get_angle(1, 0, 2)
```
%% Output
109.47090748192832
%% Cell type:markdown id: tags:
which is the correct "tetrahedral angle" between the hydrogens in a methane molecule (it's equal to $\arccos\left(-\frac{1}{3}\right)$).
%% Cell type:markdown id: tags:
Note that the atom indices in these functions are *zero_based*, meaning the first atom (here, the carbon) is referred to by the index 0.
%% Cell type:markdown id: tags:
## ASE and `quippy`
%% Cell type:markdown id: tags:
If you've read some of the documentation, you may be aware of the following alternative method for getting the distance between the two atoms:
%% Cell type:code id: tags:
``` python
struct.distance_min_image(0, 1)
```
%% Output
0.0
%% Cell type:markdown id: tags:
DANGER WILL ROBINSON! As you may have noticed, this isn't the correct bond distance. In some circumstances the code may even just crash. This is because the above function is derived from the underlying Fortran code, QUIP, rather than the `get_distance` function from before, which came from ASE. Fortran, unlike Python and C, starts counting from one instead of zero (cue the holy wars), so we just gave QUIP an atom index that doens't even exist -- and that can have unpredictable results. Now, the QUIP-derived functions are often useful, but for now it'll just be too confusing to keep this indexing distinction in our heads - let's stick with ASE until we have good reason to do otherwise.
%% Cell type:markdown id: tags:
So how do we tell whether a function is derived from QUIP or from ASE? Well, let's take a look at this function's help message:
%% Cell type:code id: tags:
``` python
help(struct.distance_min_image)
```
%% Output
Help on method <lambda> in module quippy._atoms:
<lambda> lambda self, *args, **kwargs method of quippy.atoms.Atoms instance
This interface calculates the distance between the nearest periodic images of two points (or atoms).
Return minimum image distance between two atoms or positions.
End points can be specified by any combination of atoms indices
'i' and 'j' and absolute coordinates 'u' and 'w'. If 'shift' is
present the periodic shift between the two atoms or points will
be returned in it.
Wrapper around Fortran interface ``distance_min_image`` containing multiple routines:
.. function :: distance_min_image(v,w,[shift])
:param v:
:type v: input rank-1 array('d') with bounds (3)
:param w:
:type w: input rank-1 array('d') with bounds (3)
:param shift:
:type shift: in/output rank-1 array('i') with bounds (3), optional
:returns: **ret_distance8_vec_vec** -- float
Routine is wrapper around Fortran routine ``distance8_vec_vec`` defined in file :git:`src/libAtoms/Atoms_types.f95`.
.. function :: distance_min_image(i,j,[shift])
:param i:
:type i: input int
:param j:
:type j: input int
:param shift:
:type shift: in/output rank-1 array('i') with bounds (3), optional
:returns: **ret_distance8_atom_atom** -- float
Routine is wrapper around Fortran routine ``distance8_atom_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.
.. function :: distance_min_image(i,v,[shift])
:param i:
:type i: input int
:param v:
:type v: input rank-1 array('d') with bounds (3)
:param shift:
:type shift: in/output rank-1 array('i') with bounds (3), optional
:returns: **ret_distance8_atom_vec** -- float
Routine is wrapper around Fortran routine ``distance8_atom_vec`` defined in file :git:`src/libAtoms/Atoms_types.f95`.
.. function :: distance_min_image(v,j,[shift])
:param v:
:type v: input rank-1 array('d') with bounds (3)
:param j:
:type j: input int
:param shift:
:type shift: in/output rank-1 array('i') with bounds (3), optional
:returns: **ret_distance8_vec_atom** -- float
Routine is wrapper around Fortran routine ``distance8_vec_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.
%% Cell type:markdown id: tags:
A bit confusing, yes, but take a look at the final line:
Routine is wrapper around Fortran routine ``distance8_vec_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.
It says this function is a _wrapper_ around something in the Fortran code, which means it's from QUIP and only works with _one-based_ indices (i.e. to this function, the first atom has the index 1).
%% Cell type:code id: tags:
``` python
struct.distance_min_image(1, 2)
```
%% Output
1.07
%% Cell type:markdown id: tags:
Whereas the ASE function has a different help string:
%% Cell type:code id: tags:
``` python
help(struct.get_distance)
```
%% Output
Help on method get_distance in module ase.atoms:
get_distance(self, a0, a1, mic=False, vector=False) method of quippy.atoms.Atoms instance
Return distance between two atoms.
Use mic=True to use the Minimum Image Convention.
vector=True gives the distance vector (from a0 to a1).
%% Cell type:markdown id: tags:
See that part, '`in module ase.atoms`'? That means this function comes from ASE and works with zero-based indices.
You should _always_ do this check before using any `Atoms` function you're not yet familiar with.
%% Cell type:markdown id: tags:
### Useful QUIP equivalents
%% Cell type:markdown id: tags:
Just in case you ever want to work with the Fortran QUIP routines, here are some equivalents of the queries we've seen above:
%% Cell type:code id: tags:
``` python
[struct.distance_min_image(1, i) for i in range(2, 6)]
```
%% Output
[1.07, 1.0699965409757173, 1.0700018621479124, 1.070003840600584]
%% Cell type:code id: tags:
``` python
print(struct.cosine(1, 2, 3))
print(np.arccos(struct.cosine(1, 2, 3)) * 180 / np.pi)
```
%% Output
-0.333328180365
109.470907482
%% Cell type:markdown id: tags:
Note that the _order_ of the arguments in `Atoms.cosine` is different as well: The central atom index is now the first argument.
%% Cell type:markdown id: tags:
The properties can also be accessed as one-based arrays (called `FortranArray`s). Note that these are also the _transpose_ of the ASE arrays, so the row and column indices are swapped.
%% Cell type:code id: tags:
``` python
struct.pos
```
%% Output
FortranArray([[ 0. , 0. , 0.62414, -0.99844, 0.3743 ],
[ 0. , 0. , 0.79255, 0.14425, -0.9368 ],
[ 0. , 1.07 , -0.35666, -0.35667, -0.35667]])
%% Cell type:code id: tags:
``` python
struct.lattice
```
%% Output
FortranArray([[ 10., 0., 0.],
[ 0., 10., 0.],
[ 0., 0., 10.]])
%% Cell type:code id: tags:
``` python
struct.Z
```
%% Output
FortranArray([6, 1, 1, 1, 1], dtype=int32)
%% Cell type:markdown id: tags:
## Changing `Atoms`
%% Cell type:markdown id: tags:
There's a more complete tutorial on manipulating Atoms objects in ASE at [the ASE website](https://wiki.fysik.dtu.dk/ase/tutorials/manipulating_atoms.html); this is a shorter example just to get you started.
%% Cell type:markdown id: tags:
First, let's try reorienting our methane molecule in space. The first hydrogen atom is already oriented along the Z axis; let's try to rotate the molecule so that the second one points along the X axis.
%% Cell type:code id: tags:
``` python
r_h2 = struct.get_distance(0, 2, vector=True)
r_h2[2] = 0.0
np.arccos(r_h2[0] / np.linalg.norm(r_h2))
```
%% Output
0.90371859179284519
%% Cell type:markdown id: tags:
So we need to rotate it by about -0.9 radians about the Z axis.
%% Cell type:markdown id: tags:
But first, we want to keep the old methane molecule for comparison. Let's copy it - note that this must be done _explicitly_; the Python assignment operator (`=`) only assigns a new _reference_ to the underlying object! This can be one of the hardest things to understand for programmers coming from a different language, so make sure to read up on Python references and do some of their tutorials if you find this distinction confusing.
%% Cell type:code id: tags:
``` python
struct_old = struct.copy()
```
%% Cell type:code id: tags:
``` python
struct.rotate([0.0, 0.0, 1.0], -1.0 * np.arccos(r_h2[0] / np.linalg.norm(r_h2)))
```
%% Cell type:code id: tags:
``` python
struct.get_positions()
```
%% Output
array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 1.07000000e+00],
[ 1.00880436e+00, -1.11022302e-16, -3.56660000e-01],
[ -5.04400083e-01, 8.73653852e-01, -3.56670000e-01],
[ -5.04404280e-01, -8.73653852e-01, -3.56670000e-01]])
%% Cell type:markdown id: tags:
Exercise: Now compare the positions of the rotated structure with those of the new structure; you may want to compute the RMS difference.
%% Cell type:markdown id: tags:
Now change the cell and make a supercell
%% Cell type:markdown id: tags:
## Generating structures from scratch
%% Cell type:markdown id: tags:
We could have also just used ASE functionality to pull a methane structure from its database:
%% Cell type:code id: tags:
``` python
import ase
from ase.build import molecule
```
%% Cell type:code id: tags:
``` python
me = molecule('CH4')
```
%% Cell type:code id: tags:
``` python
me.get_positions()
```
%% Output
array([[ 0. , 0. , 0. ],
[ 0.629118, 0.629118, 0.629118],
[-0.629118, -0.629118, 0.629118],
[ 0.629118, -0.629118, -0.629118],
[-0.629118, 0.629118, -0.629118]])
%% Cell type:markdown id: tags:
Like the introductory tutorial at http://libatoms.github.io/QUIP/tutorial.html - use the `diamond` generator, view a supercell, do fun things with the structure...
%% Cell type:markdown id: tags:
## Computing the energy
%% Cell type:markdown id: tags:
Introduce `Potential`. Use something like `IP SW` to get the energy of the silicon supercell. Direct user to the list of potentials and the next tutorial.
%% Cell type:code id: tags:
``` python
```
3rdparty/quip/doc/Tutorials/NRL_TB_to_xml/README deleted 100644 → 0
View file @ 23bcf97f
Set QUIP_ARCH environment variable
export QUIP_ARCH=your_quip_arch # for sh/bash
or
setenv QUIP_ARCH your_quip_arch # for csh/tcsh
Set QUIP_ROOT environment variable to the path to your top level QUIP
source directory (where the src/, subdirectory is)
Configure the build, TB must be enabled
make config
Create executables
make src/Programs/NRL_TB_to_xml
Fix TBMD format input file (let's call it nrl_tb_param_file), as downloaded
from http://cst-www.nrl.navy.mil/bind by adding after each atomic mass
the atomic number for that element
Create XML format NRL-TB parameter file for QUIP
${QUIP_ROOT}/build/${QUIP_ARCH}/NRL_TB_to_xml tightbind.parms.NRL_TB.element_name.xml < nrl_tb_param_file
3rdparty/quip/doc/Tutorials/TS_params.xml deleted 100644 → 0
View file @ 23bcf97f
<TS_params
betapol="0.75"
cutoff_coulomb="20.0"
cutoff_ms="18.0"
tolpol="0.01"
yuksmoothlength="10.0"
n_types="2"
pred_order="2"
maxipol="60"
tewald="F"
yukalpha="0.1" >
<per_type_data atomic_num="8" pol="14.131863" z="-1.4295594" type="1" />
<per_type_data atomic_num="14" pol="0.0" z="2.8591188" type="2" />
<per_pair_data C_pol="0.44302622" atnum_j="8" atnum_i="8" D_ms="0.00030700577" gamma_ms="12.165654" B_pol="1.1221903" R_ms="7.0252019" />
<per_pair_data C_pol="-1.5003213" atnum_j="8" atnum_i="14" D_ms="0.0020129372" gamma_ms="11.350477" B_pol="1.973181" R_ms="4.5780828" />
<per_pair_data C_pol="0.0" atnum_j="14" atnum_i="14" D_ms="0.33967532" gamma_ms="-0.17694797" B_pol="0.0" R_ms="-0.085202834" />
</TS_params>
3rdparty/quip/doc/Tutorials/adaptive-qmmm-references.rst deleted 100644 → 0
View file @ 23bcf97f
References
==========
.. [Bernstein2009] Bernstein, N., Kermode, J. R., & Csányi,
G. Hybrid atomistic simulation methods for materials
systems. `Rep. Prog. Phys.
72(2), 026501 (2009). <http://dx.doi.org/10.1088/0034-4885/72/2/026501>`_
.. [Buehler2006] Buehler, M., Van Duin, A., & Goddard, W. (2006). Multiparadigm
Modeling of Dynamical Crack Propagation in Silicon Using a Reactive Force
Field. `Phys. Rev. Lett., 96(9), 95505. (2006)
<http://dx.doi.org/10.1103/PhysRevLett.96.095505>`_
.. [Csanyi2004] Csányi, G., Albaret, T., Payne, M., & De Vita,
A. 'Learn on the Fly': A Hybrid Classical and Quantum-Mechanical
Molecular Dynamics Simulation. `Phys. Rev. Lett.,
93(17), 175503. (2004)
<http://prl.aps.org/abstract/PRL/v93/i17/e175503>`_
.. [Csanyi2005] Csányi, G., Albaret, T., Moras, G., Payne, M. C., & De Vita, A.
Multiscale hybrid simulation methods for material systems. `J. Phys.: Cond
Mat. 17 R691-R703 (2005).
<http://stacks.iop.org/0953-8984/17/i=27/a=R02?key=crossref.f572d9a616845900307c923f5f385044>`_
.. [DeVita1998] De Vita, A., & Car, R. A novel scheme for
accurate MD simulations of large systems. Mater. Res. Soc. Symp. Proc, 491, 473–480 (1998).
.. [Elstner1998] Elstner, M., Porezag, D., Jungnickel, G., Elsner, J.,
Haugk, M., Frauenheim, T., Suhai, S., et
al. Self-consistent-charge density-functional tight-binding
method for simulations of complex materials
properties. `Phys. Rev. B. 58 7260 (1998).
<http://prola.aps.org/abstract/PRB/v58/i11/p7260_1>`_
.. [Fineberg1991] Fineberg, J., Gross, S., Marder, M., & Swinney, H.
Instability in dynamic fracture. `Phys. Rev. Lett., 67(4),
457–460 (1991). <http://dx.doi.org/10.1103/PhysRevLett.67.457>`_
.. [FrenkelSmit2001] Frenkel, D., & Smit, B.
`Understanding Molecular Simulation <http://molsim.chem.uva.nl/frenkel_smit/>`_,
Academic. (New York, 2001).
.. [Freund1998] Freund, L.
`Dynamic fracture mechanics <http://ebooks.cambridge.org/ebook.jsf?bid=CBO9780511546761>`_,
Cambridge University Press (Cambridge, 1998)
.. [Griffith1921] Griffith, A. A.The Phenomena of Rupture and Flow in Solids.
`Philosophical Transactions of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 221(582-593), 163–198. (1921)
<http://www.dx.doi.org/10.1098/rsta.1921.0006>`_
.. [Irwin1948] Irwin, G. R. Fracturing of Metals. (1948)
.. [Kermode2008] Kermode, J. R., Albaret, T., Sherman, D., Bernstein,
N., Gumbsch, P., Payne, M. C., Csányi, G., and A. De Vita. Low-speed
fracture instabilities in a brittle crystal. `Nature, 455,
1224–1227 (2008). <http://www.nature.com/doifinder/10.1038/nature07297>`_
.. [Kermode2008a] Kermode, J. R. Multiscale Hybrid Simulation of Brittle
Fracture. `PhD Thesis, University of Cambridge (2008).
<http://www.jrkermode.co.uk/Publications>`_
.. [Lawn1993] Lawn, B. R. `Fracture of Brittle Solids --- Second Edition
<http://books.google.co.uk/books/about/Fracture_of_Brittle_Solids.html?id=CPxFY2GSdsoC>`_
(Cambridge Univ Pr, 1993)
.. [Li2003] Li, J. AtomEye: an efficient atomistic configuration
viewer. `Modell. Simul. Mater. Sci. Eng. (2003) <http://li.mit.edu/Archive/Papers/03/Li03a.pdf>`_.
See also websites for `original version <http://li.mit.edu/Archive/Graphics/A>`_, and
`modified version <http://jrkermode.co.uk/AtomEye>`_.
.. [Moras2010] Moras, G., Choudhury, R., Kermode, J. R., Csányi, G.,
Payne, M. C., & De Vita, A. Hybrid Quantum/Classical Modeling of
Material Systems: The Learn on the Fly Molecular Dynamics
Scheme. In T. Dumitrica (Ed.), `Trends in Computational
Nanomechanics: Transcending Length and Time Scales
(pp. 1–23). Springer (2010)
<http:///dx.doi.org/10.1007/978-1-4020-9785-0_1>`_
.. [Pizzagalli2013] Pizzagalli, L., Godet, J., Guénolé, J., Brochard, S., Holmstrom, E.,
Nordlund, K., & Albaret, T. A new parametrization of the Stillinger-Weber
potential for an improved description of defects and plasticity of silicon.
`J. Phys.: Cond. Mat. 25(5), 055801. (2013)
<http://www.dx.doi.org/10.1088/0953-8984/25/5/055801>`_
.. [Stillinger1985] Stillinger, F. H., & Weber, T. A. Computer simulation
of local order in condensed phases of silicon. `Phys. Rev. B.,
31(8), 5262–5271. (1985). <http://link.aps.org/doi/10.1103/PhysRevB.31.5262>`_
.. [Swadener2002] Swadener, J. G., Baskes, M. I., & Nastasi, M. Molecular
Dynamics Simulation of Brittle Fracture in Silicon. `Phys. Rev. Lett. 89
085503 (2002). <http://dx.doi.org/10.1103/PhysRevLett.89.085503>`_
.. [Vink2001] Vink, R. L. C., Barkema, G. T., Van der Weg, W. F., & Mousseau, N.
Fitting the Stillinger–Weber potential to amorphous silicon. `J. Non. Cryst. Sol., 282(2-3), 248–255. (2001).
<http://dx.doi.org/10.1016/S0022-3093(01)00342-8>`_
.. [Zimmerman2004] Zimmerman, J. A., Webb, E. B., Hoyt, J. J., Jones,
R. E., Klein, P. A., & Bammann, D. J. Calculation of stress in
atomistic simulation. `Modell. Simul. Mater. Sci. Eng. (2004).
<http://stacks.iop.org/0965-0393/12/S319>`_
3rdparty/quip/doc/Tutorials/adaptive-qmmm-solutions.rst deleted 100644 → 0
View file @ 23bcf97f
Solutions
=========
.. _make_crack:
Step 1 solution --- ``make_crack.py``
-------------------------------------
Download :download:`make_crack.py`
.. literalinclude:: make_crack.py
.. _run_crack_classical:
Step 2 solution --- ``run_crack_classical.py``
----------------------------------------------
Download :download:`run_crack_classical.py`
.. literalinclude:: run_crack_classical.py
.. _run_crack_lotf:
Step 3 solution --- ``run_crack_lotf.py``
-----------------------------------------
Download :download:`run_crack_lotf.py`
.. literalinclude:: run_crack_lotf.py
3rdparty/quip/doc/Tutorials/adaptive-qmmm-step0.rst deleted 100644 → 0
View file @ 23bcf97f
Introduction
************
Scope of the tutorial
=====================
In this tutorial, we will carry out classical and hybrid multiscale
QM/MM (quantum mechanics/molecular mechanics) molecular dynamics
simulations using the Learn on the Fly (LOTF) schmee for fracture in
silicon. For the classical simulations we will use the
Stillinger-Weber [Stillinger1985]_ interatomic potential, which
provides a generally good description of many properties of silicon,
but, not, however of brittle fracture as we will see. The focus of the
tutorial is on embedding techniques, so we will use an approximate
quantum mechanical method which demonstrates the approach while
remaining fast enough to carry out calculations during the time
available, namely Density Functional Tight Binding [Elstner1998]_.
The tutorial is divided into three sections. There are also some
extension tasks at the end which you can complete if you have time. In
the first section we will prepare the model fracture system for our
simulations. In the second part, we will carry out classical molecular
dynamics, and in the third part of the tutorial, the 'Learn on the
Fly' ([DeVita1998]_, [Csanyi2004]_) embedding scheme will be used to
carry out coupled classical and quantum calculations. One of the main
advantages of the LOTF approach is that the QM region can be moved
during the simulation (*adaptive* QM/MM). You can read more details
about multiscale embedding methods applied to materials systems in
[Bernstein2009]_, and more about brittle fracture in silicon
investigated with the 'Learn on the Fly' technique in [Kermode2008]_.
.. _practical:
Practical considerations
========================
In this tutorial we will carry out simulations using a combination of
two packages: `quippy <http://www.jrkermode.co.uk/quippy>`_, a Python
interface to the `libAtoms/QUIP <http://www.libatoms.org>`_ MD code
developed at King's College London, Cambridge University, the Naval
Research Lab in Washington and the Fraunhofer IWM in Freiburg, and the
Atomic Simulation Environment, `ASE <https://wiki.fysik.dtu.dk/ase>`_,
a Python framework developed at the Centre for Atomistic Materials
Design (`CAMd <http://www.camd.dtu.dk/>`_) at DTU, Copenhagen.
If you're not familiar with Python don't worry, it's quite easy to
pick up and the syntax is very similar to Fortran. We have provided
template Python code which you should copy and paste into your text
editor. Save your script with a ``.py`` file extension. Whenever you
see ``...`` it means there is something you need to add to the code in
order to complete it; usually there will be a comment to give you a
hint. The source code listings are all in boxes like this::
print 'Example code listing' # Print something out
print ... # You need to add something here
# e.g. replace the ... with a number
To run code, we will be using `ipython <http://ipython.org>`_, an
interactive Python shell. You can start `ipython` with a simple shell
command::
$ ipython
Start by importing everything from the `quippy` package with the
command::
In [1]: from qlab import *
.. note::
If you have not installed the `AtomEye` plugin then the command
above will give an :exc:`ImportError`. You can use::
from quippy import *
instead, which only imports `quippy` and not `AtomEye`, but you
will then not be able to visualise interactively using the
:func:`~qlab.view` function - you can always save configurations to
external files and visualise them with other tools such as `VMD
<http://www.ks.uiuc.edu/Research/vmd/>`_ or `OVITO
<http://www.ovito.org>`_ instead. See the :mod:`qlab` and
:mod:`atomeye` documentation for further details.
You should prepare your Python scripts in a text editor and then run
them from within `ipython` with the `%run
<http://ipython.org/ipython-doc/stable/interactive/tutorial.html#running-and-editing>`_
command. For example::
In [2]: run make_crack.py
If your script is taking a long time to run and you realise something
is wrong you can interrupt it with `Ctrl+C`.
The tutorial is structured so that you will build up your script as
you go along. You should keep the same `ipython` session open, and
simply run the script again each time you add something new. You will
find it easier to follow the tutorial if you use the same variable
names as we have used here.
You can also execute individual Python statements by typing them
directly into the `ipython` shell. This is very useful for debugging
and for interactive visualisation (described in more detail in the
:ref:`visualisation2` section). Finally, it's possible to copy and
paste code into `ipython` using the `%paste` or `%cpaste` commands.
You can follow the links in the tutorial to read more documentation for
functions and classes we will be using, and for `quippy` routines you
can also browse the source code by clicking on the `[source]` link
which you will find to the right of class and function definitions in
the online documentations. You can also explore functions
interactively from within `ipython` by following the name of an object
or function with a ``?``, e.g. ::
In [3]: Atoms ?
displays information about the :class:`~.Atoms` class, or ::
In [4]: write ?
displays information about the :func:`~quippy.io.write` function. See
the `ipython tutorial
<http://ipython.org/ipython-doc/stable/interactive/tutorial.html>`_
for more information.
Each subsection indicates the approximate amount of time you should
spend on it. At the end of each subsection there is a *Milestone*
where a script for that stage of the tutorial is provided. If you run
out of time, just skip ahead to the next milestone, download the
script and then continue with the next section.
Continue with :ref:`step1`.
3rdparty/quip/doc/Tutorials/adaptive-qmmm-step1.rst deleted 100644 → 0
View file @ 23bcf97f
.. _step1:
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 :ref:`thin strip geometry <thin_strip>`
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 :ref:`practical` 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.
.. _parameters:
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, :math:`(111)[01\bar{1}]`, which
means a crack propagating on the :math:`(111)` cleavage plane (the
lowest surface energy of all silicon surfaces) with the crack front
along the :math:`[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
<https://wiki.fysik.dtu.dk/ase/ase/units.html>`_ where ``Ang = eV =
1``. The energy release rate `initial_G` is given in the
widely used units of J/m\ :superscript:`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 :download:`params.xml` file, which contains
the parameters for the SW potential (and also for DFTB, needed for
:ref:`step3`)
.. _latticeconstant:
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
:func:`~ase.structure.bulk` function from the :mod:`ase.structure` module::
si_bulk = ... # Build the 8-atom diamond cubic cell for Si
The variable `si_bulk` is an :class:`~ase.atoms.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 :mod:`qlab` module, will then be
able to interactively visualise the Si unit cell with the
:func:`~qlab.view` function from the :mod:`qlab` module, which you
should type in at the `ipython` prompt::
In [5]: view(si_bulk)
.. image:: si_bulk.png
:align: center
:width: 300
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
<http://mt.seas.upenn.edu/Archive/Graphics/A>`_ or the documentation
for the :mod:`qlab` and :mod:`atomeye` modules.
Now, we initialise the Stillinger-Weber (SW) classical interatomic
potential using quippy's :class:`~quippy.potential.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 :meth:`~ase.atoms.Atoms.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
:class:`~quippy.potential.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 :meth:`~quippy.potential.Minim.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 :meth:`~ase.atoms.Atoms.cell` attribute of
the `si_bulk` object, which returns a :math:`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 :download:`make_crack_1.py`.
1.2 Calculation of elastic and surface properties of silicon (30 minutes)
-------------------------------------------------------------------------
.. _youngs_modulus_and_poisson_ratio:
Calculation of the Young's modulus and the Poisson ratio
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Following the discussion :ref:`above <thin_strip>` 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 :math:`\nu` in the :math:`xy` plane, we need the :math:`6 \times
6` matrix of the elastic constants :math:`C_{ij}`. This matrix `c` can be
calculated using the :meth:`~quippy.potential.Potential.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 :attr:`~ase.units.GPa` constant from the `ase.units module
<https://wiki.fysik.dtu.dk/ase/ase/units.html>`_ module is used to
convert from pressure units of eV/A\ :superscript:`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
:ref:`parameters section <parameters>` above), using the functions
:func:`~quippy.elasticity.youngs_modulus` and
:func:`~quippy.elasticity.poisson_ratio` from the
:mod:`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 :math:`(111)[01\bar{1}]` crack system, you
should get a Young's modulus of 142.8 GPa and a Poisson ratio of
0.265.
.. _surface_energy:
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
:func:`ase.lattice.cubic.Diamond <lattice.cubic.Diamond>`
from the :mod:`ase.lattice <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):
.. image:: unit_slab.png
:align: center
:width: 400
We now shift the `unit_slab` vertically so that we will open up a
surface along a :math:`(111)` glide plane, cutting vertically aligned
bonds (see e.g. `this image
<http://ej.iop.org/images/0295-5075/72/3/410/Full/img11.gif>`_). 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 :attr:`~ase.atoms.Atoms.positions` is a `(N,3)` array containing
the Cartesian coordinates of the atoms, and
:meth:`~ase.atoms.Atoms.set_scaled_positions` and
:meth:`~ase.atoms.Atoms.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).
.. image:: unit_slab_shifted.png
:align: center
:width: 400
Note how the top and bottom layers now correspond to :math:`(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
:meth:`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
:ref:`parameters section <parameters>` above::
surface = unit_slab.copy()
surface. ... # Add vacuum along y axis
You should get a surface unit cell which looks something like this:
.. image:: surface.png
:align: center
:width: 400
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 :class:`~ase.atoms.Atoms` object, the potential energy of the
atomic system can be calculated with
:meth:`~ase.atoms.Atoms.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
:meth:`~ase.atoms.Atoms.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 :math:`\gamma_{(111)}` = 1.36 J/m\
:superscript:`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 :download:`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 :ref:`parameters section <parameters>`).
Note that the python function :py:func:`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`
:math:`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 :meth:`~ase.atoms.Atoms.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:
.. image:: crack_slab_1.png
:align: center
:width: 600
.. (You might find it useful to press `Shift+z` to centre the AtomEye
.. view on a fractional lattice coordinate of `(.5, .5, .5)` rather
.. than the default of `(0., 0., 0.`).)
.. _crack_fixatoms:
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
:ref:`discussion of the thin strip geometry above<thin_strip>`). 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 :func:`~qlab.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:
.. image:: fixed_mask.png
:align: center
:width: 600
Now we can use the :class:`~constraints.FixAtoms` class to
fix the positions of the atoms according to the mask `fixed_mask`, and
then attach the constraint to `crack_slab` using
:meth:`~ase.atoms.Atoms.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 :func:`~quippy.crack.G_to_strain`
function which implements the :ref:`thin strip equation described
above <thin_strip_equation>`. The `strain` is then used to displace
the `y` coordinate of the atomic positions according to the strain
ramp produced by the :func:`~quippy.crack.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 :math:`(111)` case:
.. image:: crack_slab_2.png
:align: center
:width: 600
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
:class:`~quippy.potential.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 :ref:`parameters
section <parameters>`::
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:
.. image:: crack_slab_3.png
:align: center
:width: 600
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
:func:`~quippy.crack.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` :class:`~ase.atoms.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 :ref:`extendedxyz` in the `output_file`, whose name is
defined in the :ref:`parameters section <parameters>`::
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
:ref:`make_crack`, and your XYZ file like :download:`crack.xyz`.
When you are ready, proceed to :ref:`step2`.
3rdparty/quip/doc/Tutorials/adaptive-qmmm-step2.rst deleted 100644 → 0
View file @ 23bcf97f
.. _step2:
Step 2: Classical MD simulation of fracture in Si
=================================================
This part of the tutorial requires the initial crack structure ``crack.xyz``
produced in :ref:`step1`. If you had problems completing the first part, you
can :download:`download it here <crack.xyz>` instead. Start by creating a new
empty script, named ``run_crack_classical.py``.
2.1 Initialisation of the atomic system (20 minutes)
----------------------------------------------------
Import the relevant modules
^^^^^^^^^^^^^^^^^^^^^^^^^^^
As before, we import the necessary modules, classes and functions::
import numpy as np
from ase.constraints import FixAtoms
from ase.md.verlet import VelocityVerlet
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
import ase.units as units
from quippy.atoms import Atoms
from quippy.potential import Potential
from quippy.io import AtomsWriter
from quippy.crack import (get_strain,
get_energy_release_rate,
ConstantStrainRate,
find_crack_tip_stress_field)
Note that the molecular dynamics classes and functions come from
`ASE`, while the :class:`~.Atoms` object, potentials and
crack-specific functionality come from `quippy`.
Definition of the simulation parameters
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. _parameters2:
Let's define the simulation parameters. The meaning of each parameter is explained
in the comments on the right of each line::
input_file = 'crack.xyz' # File from which to read crack slab structure
sim_T = 300.0*units.kB # Simulation temperature
nsteps = 10000 # Total number of timesteps to run for
timestep = 1.0*units.fs # Timestep (NB: time base units are not fs!)
cutoff_skin = 2.0*units.Ang # Amount by which potential cutoff is increased
# for neighbour calculations
tip_move_tol = 10.0 # Distance tip has to move before crack
# is taken to be running
strain_rate = 1e-5*(1/units.fs) # Strain rate
traj_file = 'traj.nc' # Trajectory output file (NetCDF format)
print_interval = 10 # Number of time steps between
# writing output frames
param_file = 'params.xml' # Filename of XML file containing
# potential parameters
mm_init_args = 'IP SW' # Initialisation arguments for
# classical potential
Setup of the atomic structure
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
As a first step, we need to initialise the
:class:`~quippy.atoms.Atoms` object by loading the atomic structure created
in :ref:`step1` from the `input_file` ``crack.xyz``. Note that the global
`fortran_indexing` setting should be set to `False`. Otherwise `quippy` uses
atom indices in the range :math:`1 \ldots N`, which would not be consistent with
the python indexing used in ASE (:math:`0\ldots N-1`).
It is also necessary to read in the original width and height of the slab and
the original crack position, which were saved in ``crack.xyz`` at the end
of :ref:`Step 1 <step1>`::
print 'Loading atoms from file %s' % input_file
atoms = ... # Load atoms from file
orig_height = atoms.info['OrigHeight'] # Initialise original height
orig_crack_pos = atoms.info['CrackPos'].copy() # Initialise original crack position
Note that we make a copy of the ``CrackPos`` entry in the
:attr:`~quippy.atoms.info` dictionary, since otherwise
`orig_crack_pos` would continue to refer to the current crack position
as it is updated during the dynamical simulation.
Setup of the constraints
^^^^^^^^^^^^^^^^^^^^^^^^
Now we can set constraints on the atomic structure which will be
enforced during dynamics. In order to carry out the fracture MD
simulation, we need to fix the positions of the top and bottom atomic
rows (we call this constraint `fix_atoms`).
The `fix_atoms` constraint, which is exactly the same as the constraint used
for :ref:`relaxing the positions of the crack slab <crack_fixatoms>` above. In
order to do this, we need to find the `y` coordinate of the top, bottom atomic
rows. The `x` coordinates of the left and right edges of the slab might also
be useful later on. This can be easily done as before::
top = ... # Maximum y coordinate
bottom = ... # Minimum y coordinate
left = ... # Minimum x coordinate
right = ... # Maximum x coordinate
Now it is possible to define the `fixed_mask` array, which is `True`
for each atom whose position needs to be fixed, and `False` otherwise,
exactly as before, and to initialise the `fix_atoms` constraint, in
the same way we did it in `Step 1` (i.e., using the
:class:`~constraints.FixAtoms` class)::
fixed_mask = ... # Define the list of fixed atoms
fix_atoms = ... # Initialise the constraint
print('Fixed %d atoms\n' % fixed_mask.sum()) # Print the number of fixed atoms
This constraint needs to be attached to our `atoms` object
(see :meth:`~quippy.atoms.Atoms.set_constraint`)::
atoms. ... # Attach the constraints to atoms
To increase :math:`\epsilon_{yy}` of all atoms at a constant rate (see
the `strain_rate` and `timestep` :ref:`parameters <parameters2>`), we
use the :class:`~quippy.crack.ConstantStrainRate` class::
strain_atoms = ConstantStrainRate(orig_height, strain_rate*timestep)
You can look at the documentation for the
:class:`~quippy.crack.ConstantStrainRate` class to see how this works. The
:meth:`~constraints.adjust_positions` routine simply increases the strain of
all atoms. We will attach this to the dynamics in step 2.2 below.
Setup of the potential
^^^^^^^^^^^^^^^^^^^^^^
Before starting the MD simulation, the SW classical potential must be
initialised and attached to the `atoms` object. As in `Step 1`, we
use quippy's :class:`~quippy.potential.Potential` class, but now we
need to pass the `cutoff_skin` parameter, which is used to decide when
the neighbour list needs to be updated (see the attribute
:attr:`~quippy.potential.Potential.cutoff_skin`). Moreover, we request
the potential to compute per-atom stresses whenever we compute forces
using :meth:`~quippy.potential.Potential.set_default_quantities`, to
save time when locating the crack tip (discussed in more detail
:ref:`below <position_crack_tip>`). The
:meth:`~quippy.atoms.Atoms.set_calculator` method should then be used
to set the calculator to the SW potential::
mm_pot = ... # Initialise the SW potential with cutoff_skin
mm_pot.set_default_quantities(['stresses'])
atoms. ... # Set the calculator
Milestone 2.1
^^^^^^^^^^^^^
At this stage your script should look something like :download:`this <run_crack_classical_1.py>`.
2.2 Setup and run the classical MD (20 minutes)
-----------------------------------------------
Setting initial velocities and constructing the dynamics object
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
There are still a few things that need to be done before running the
MD fracture simulation. We will follow the standard `ASE molecular
dynamics methodology
<https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html>`_. We will set
the initial temperature of the system to `2*sim_T`: it will then
equilibrate to `sim_T`, by the Virial theorem::
MaxwellBoltzmannDistribution(atoms, 2.0*sim_T)
A MD simulation in the NVE ensemble, using the Velocity Verlet
algorithm, can be initialised with the ASE
:class:`~md.verlet.VelocityVerlet` class, which requires two
arguments: the atoms and the time step (which should come from the
`timestep` :ref:`parameter <parameters2>`::
dynamics = ... # Initialise the dynamics
Printing status information
^^^^^^^^^^^^^^^^^^^^^^^^^^^
Let's also define a function that prints the relevant information at
each time step of the MD simulation. The information can be saved
inside the :attr:`~quippy.atoms.Atoms.info` dictionary, so that it
also gets saved to the trajectory file `traj_file`.
The elapsed simulation time can be obtained with
``dynamics.get_time()`` (note that the time unit in ASE is
:math:`\mathrm{\AA}\sqrt{\mathrm{amu}/\mathrm{eV}}`, not `fs`). You
should use the :meth:`~ase.atoms.Atoms.get_kinetic_energy` method to
calculate the temperature (*Note*: you will need the :attr:`units.kB`
constant, which gives the value of the Boltzmann constant in eV/K),
and the functions :func:`~quippy.crack.get_strain` and
:func:`~quippy.crack.get_energy_release_rate` to return the current
strain energy release rate, respectively. ::
def printstatus():
if dynamics.nsteps == 1:
print """
State Time/fs Temp/K Strain G/(J/m^2) CrackPos/A D(CrackPos)/A
---------------------------------------------------------------------------------"""
log_format = ('%(label)-4s%(time)12.1f%(temperature)12.6f'+
'%(strain)12.5f%(G)12.4f%(crack_pos_x)12.2f (%(d_crack_pos_x)+5.2f)')
atoms.info['label'] = 'D' # Label for the status line
atoms.info['time'] = ... # Get simulation time
# and convert to fs
atoms.info['temperature'] = ... # Get temperature in K
atoms.info['strain'] = ... # Get strain
atoms.info['G'] = ... # Get energy release rate,
# and convert to J/m^2
crack_pos = ... # Find crack tip as in step 1
atoms.info['crack_pos_x'] = crack_pos[0]
atoms.info['d_crack_pos_x'] = crack_pos[0] - orig_crack_pos[0]
print log_format % atoms.info
This logger can be now attached to the `dynamics`, so that the information is
printed at every time step during the simulations::
dynamics.attach(printstatus)
Checking if the crack has advanced
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The same can be done to check during the simulation if the crack has advanced,
and to stop incrementing the strain if it has::
def check_if_cracked(atoms):
crack_pos = ... # Find crack tip position
# stop straining if crack has advanced more than tip_move_tol
if not atoms.info['is_cracked'] and (crack_pos[0] - orig_crack_pos[0]) > tip_move_tol:
atoms.info['is_cracked'] = True
del atoms.constraints[atoms.constraints.index(strain_atoms)]
The `check_if_cracked` function can now be attached to the dynamical
system, requesting an interval of 1 step (i.e. every time) and passing the
`atoms` object along to the function as an extra argument::
dynamics.attach(check_if_cracked, 1, atoms)
We also need to attach the `:meth:`quippy.crack.ConstrainStrainRate.apply_strain` method
of `strain_atoms` to the dynamics::
dynamics.attach(strain_atoms.apply_strain, 1, atoms)
Saving the trajectory
^^^^^^^^^^^^^^^^^^^^^
Finally, we need to initialise the trajectory file `traj_file` and to
save frames to the trajectory every `traj_interval` time steps. This
is done by creating a trajectory object with the
:func:`~quippy.io.AtomsWriter` function, and then attaching this
trajectory to the `dynamics`::
trajectory = ... # Initialise the trajectory
dynamics. ... # Attach the trajectory with an interval of
# traj_interval, passing atoms as an extra argument
We will save the trajectory in :ref:`netcdf` format. This is a binary
file format that is similar with the :ref:`extendedxyz` format we used
earlier, with the advantage of being more efficient for large files.
Running the dynamics
^^^^^^^^^^^^^^^^^^^^
After all this, a single command will run the MD for `nsteps` (see the `ASE
molecular dynamics methodology
<https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html>`_ for more information)::
dynamics.run(nsteps)
Milestone 2.2
^^^^^^^^^^^^^
If you have problems you can look at the complete version of the
:ref:`run_crack_classical` script. Leave your classical MD simulation
running and move on to the next section of the tutorial.
The first few lines produced by the ``run_crack_classical.py`` script should
look something like this::
Loading atoms from file crack.xyz
Fixed 240 atoms
State Time/fs Temp/K Strain G/(J/m^2) CrackPos/A D(CrackPos)/A
---------------------------------------------------------------------------------
D 1.0 560.097755 0.08427 5.0012 -30.61 (-0.00)
D 2.0 550.752265 0.08428 5.0024 -30.61 (-0.00)
D 3.0 535.568949 0.08429 5.0036 -30.61 (-0.00)
D 4.0 515.074874 0.08430 5.0047 -30.61 (-0.00)
D 5.0 489.977973 0.08431 5.0059 -30.61 (-0.00)
D 6.0 461.140488 0.08432 5.0071 -30.61 (-0.00)
D 7.0 429.546498 0.08433 5.0083 -30.61 (-0.00)
D 8.0 396.264666 0.08434 5.0095 -30.61 (-0.01)
D 9.0 362.407525 0.08435 5.0107 -30.61 (-0.01)
D 10.0 329.088872 0.08436 5.0119 -30.61 (-0.01)
Here we see the current time, temperature, strain, energy release rate `G`, the
`x` coordinate of the crack position, and the change in the crack position since
the beginning of the simulation. In the early stages of the calculation, the
strain and `G` are both increasing, and the temperature is rapidly falling
towards ``sim_T = 300`` as anticipated.
.. _visualisation2:
2.3 Visualisation and Analysis (as time permits)
------------------------------------------------
Start another `ipython` session is a new terminal with plotting
support enabled, using the shell command::
ipython --pylab
This will allow you to look at the progress of your classical fracture
simulation while it continues to run. All the example code given in
this section should be entered directly at the `ipython` prompt.
The first step is to import everything from `quippy` using the
:mod:`qlab` interactive module, then open your trajectory using the
:func:`~qlab.view` function::
from qlab import *
set_fortran_indexing(False)
view("traj.nc")
As we saw :ref:`earlier <latticeconstant>`, this will open an AtomEye viewer
window containing a visual representation of your crack system (as before
``fortran_indexing=False`` is used to number the atoms starting from zero). You
can use the `Insert` and `Delete` keys to move forwards or backwards through the
trajectory, or `Ctrl+Insert` and `Ctrl+Delete` to jump to the first or last
frame -- note that the focus must be on the AtomEye viewer window when you use
any keyboard shortcuts. The current frame number is shown in the title bar of
the window.
The function :func:`~qlab.gcat` (short for "get current atoms") returns a
reference to the :class:`~.Atoms` object currently being visualised (i.e. to the
current frame from the trajectory file). Similarly, the :func:`~qlab.gcv`
function returns a reference to the entire trajectory currently being viewed as
an :class:`~qlab.AtomsReaderViewer` object.
You can change the frame increment rate by setting
the :attr:`~atomeye.AtomEyeViewer.delta` attribute of the viewer, e.g. to
advance by ten frames at a time::
set_delta(10)
Or, to jump directly to frame 100::
set_frame(100)
You can repeat the ``view("traj.nc")``
command as your simulation progresses to reload the file (you can use `Ctrl+R`
in the `ipython` console to search backwards in the session history to save
typing).
.. _stress_analysis:
Stress field analysis
^^^^^^^^^^^^^^^^^^^^^
To compute and display the instantaneous principal per-atom stress
:math:`\sigma_{yy}` as computed by the SW potential for a
configuration near the beginning of your dynamical simulation::
mm_pot = Potential('IP SW', param_filename='params.xml')
at = gcat()
at.set_calculator(mm_pot)
mm_sigma = at.get_stresses()
sigma_yy = mm_sigma[:,1,1]
aux_property_coloring(sigma_yy)
The `mm_sigma` array has shape `(len(atoms), 3, 3)`, i.e. it is
made up of a :math:`3 \times 3` stress tensor :math:`\sigma_{ij}` for
each atom. The `sigma_yy` array is the ``[1, 1]`` component of each of
these arrays, i.e. :math:`\sigma_{yy}`. To read off the value of the
stress on a particular atom, just `right click` on it. As before, this
prints various information in the `ipython` console. The `_show`
property corresponds to the values currently being used to colour the
atoms. You will see that :math:`\sigma_{yy}` is very strongly peaked
near the crack tip. If you prefer to see the values in GPa, you could
do ::
aux_property_coloring(sigma_yy/units.GPa)
.. image:: sigma_yy.png
:align: center
:width: 600
The concept of per-atom stresses is a little arbitrary. The values we
are plotting here were obtained from partitioning the total virial
stress tensor, which is given by
.. math::
\tau_{ij} = \frac{1}{\Omega} \sum_{k \in \Omega} (-m^{(k)} (u_i^{(k)}-
\bar{u}_i) (u_j^{(k)}- \bar{u}_j) %\\
+ \frac{1}{2} \sum_{\ell \in \Omega} ( x_i^{(\ell)} - x_i^{(k)}) f_j^{(k\ell)}
)
where :math:`k` and :math:`l` are atom indices, :math:`ijk` are Cartesian
indicies, :math:`\Omega` is the cell volume, :math:`m^{(k)}`,
:math:`u^{(k)}`, :math:`x^{(k)}` and :math:`f^{(k)}` are respectively the
mass, velocity, position of atom :math:`k` and :math:`f^{kl}_j` is
the :math:`j`\ th component of the force between atoms :math:`k` and
:math:`l`. The first term is a kinetic contribution which vanishes at
near zero temperature, and it is common to use the second term to
define a per-atom stress tensor.
Note, however, that this requires a definition of the atomic volume. By default
the :meth:`~quippy.potential.Potential.get_stresses` function simply divides the
total cell volume by the number of atoms to get the volume per atom. This is
not a very good approximation for our cell, which contains a lot of empty
vacuum, so the volume per atom comes out much too large, and the stress
components much too small, e.g. the peak stress, which you can print in units of
GPa with::
print mm_sigma.max()/units.GPa
is around 4 GPa. Values of stress in better agreement with linear
elastic theory can be obtained by assuming all atoms occupy the same
volume as they would in the equilibrium bulk structure::
mm_pot.set(vol_per_atom=si_bulk.get_volume()/len(si_bulk))
mm_sigma = at.get_stresses()
print mm_sigma.max()/units.GPa
gives a value of around 25 GPa. As this is only a simple rescaling,
the unscaled virial stress values are perfectly adequate for locating
the crack tip.
Use values from the `sigma_yy` array to plot the :math:`\sigma_{yy}` virial
stress along the line :math:`y=0` ahead of the crack tip, and verify the stress
obeys the expected :math:`1/\sqrt{r}` divergence near the crack tip, and tends
to a constant value ahead of the crack, due to the thin strip loading. *Hint:*
use a mask to select the relevant atoms, as we did when fixing the edge atoms
above. You can use the matplotlib :func:`~matplotlib.pyplot.plot` function to
produce a plot.
.. _time_avg_stress:
Time-averaged stress field
^^^^^^^^^^^^^^^^^^^^^^^^^^
By now, you should have a few picoseconds of dynamics in your trajectory file.
Reload with ``view("traj.nc")`` to see what is happening. You can jump to the
end with `Ctrl+Delete`, or by typing `last()` into the `ipython` console. Here
is what the instantaneous :math:`\sigma_{yy}` looks like after 5 ps of dynamics:
.. image:: classical-crack-sigma-yy.png
:align: center
:width: 600
As you can see, the stress field is rather noisy because of
contributions made by the random thermal motion of atoms. The
:func:`~quippy.crack.find_crack_tip_stress_field` uses an exponential
moving average of the stress field when finding the tip. This average
is stored in the ``avg_sigma`` :attr:`array entry
<~quippy.atoms.Atoms.arrays>` inside the Atoms object, which is saved
with each frame in the trajectory. For techical reasons this is stored
as a reshaped array of shape ``(len(atoms), 9)`` rather than
``(len(atoms), 3, 3)`` array, so you can find the :math:`sigma_{yy}`
components in the 5th column (counting from zero as usual in Python),
i.e. ::
aux_property_coloring(gcat().arrays['avg_sigma'][:, 4])
You should find that the crack tip is more well defined in the average stress:
.. image:: classical-crack-sigma-yy-average.png
:align: center
:width: 600
.. _coordination:
Geometry and coordination analysis
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Press `k` to colour the atoms by coordination. This is based on the
:attr:`~.Atoms.nneightol` attribute of the Atoms object, which we set
to a value of 1.3 in the ``make_crack.py`` script. This factor acts as
a multipler for the covalent radii of the atomic species, taken from
the :attr:`quippy.periodictable.ElementCovRad` array. You can check
the maximum Si--Si bond-length this corresponds to with::
print 1.3*2*ElementCovRad[14]
Note that ``14`` is the atomic number of silicon. After the simulation has run
for a little while, you should be able to see both under-coordinated (green) and
over-coordinated (red) atoms near the crack tip.
Here is a typical snapshot at the end of 10 ps of dynamics. Note the
large number of defects, indicating that the fracture surface is not
atomically smooth as we find it to be in experiments. In your
simulation you may be able to spot signs of energy dissipation
mechanisms, such as dislocation emission from the crack tip.
.. image:: classical-crack-coordination.png
:align: center
:width: 600
.. _render_movie:
Rendering a movie of the simulation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
If you would like to make a movie of your simulation, you can use
the :func:`~qlab.render_movie` function. Arrange the AtomEye window so that the
crack is on the left hand side of the window at the beginning of the simulation
and near the right hand side at the end, then run the command::
render_movie('movie.mp4')
This function renders each frame to a ``.jpg`` file, before combining the
snapshots with the `ffmpeg <http://www.ffmpeg.org/>`_ tool to make a movie like
this one:
.. video:: classical-111 720 360
The example movie above makes the ductile nature of the fracture propagation
much clearer. We see local amorphisation, the formation of
strange *sp*\ :superscript:`2` tendrils, and temporary crack arrest. Comparing
again with the :ref:`experimental TEM images <si_tem_images>` makes it clear
that, as a description of fracture in real silicon, the SW potential falls some
way short.
.. _position_crack_tip:
Position of the crack tip
^^^^^^^^^^^^^^^^^^^^^^^^^
The :func:`~quippy.crack.find_crack_tip_stress_field` function works by
fitting per-atom stresses calculated with the SW potential (the
concept of per-atom stresses will be discussed in more detail below)
in the region near the crack tip to the Irwin solution for a singular
crack tip under Mode I loading, which is of the form
.. math::
\sigma_{ij}(r, \theta) = \frac{K_I}{2\pi r} f_{ij}(\theta)
where :math:`K_I` is the Mode I stress intensity factor, and the
angular dependence is given by the set of universal functions
:math:`f_{ij}(\theta)`.
You can verify this by comparing the position detected by
:func:`~quippy.crack.find_crack_tip_stress_field`, stored in the
`crack_pos` attribute, with the positions of atoms that visually look
to be near the tip --- `right click` on atoms in the AtomEye
viewer window to print information about them, including their
positions.
Compare the automatically detected crack position (printed as the
`crack_pos_x` parameter when you change frames in the AtomEye viewer,
or available via ``gcat().info['crack_pos_x']``) with what a visual
inspection of the crack system would tell you. Do you think it's
accurate enough to use as the basis for selecting a region around the
crack tip to be treated at the QM level?
.. _plot_G_and_crack_pos_x:
Evolution of energy release rate and crack position
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
For :ref:`netcdf` trajectories,
the :attr:`AtomsReaderViewer.reader.netcdf_file` attribute of the current
viewer object :func:`~qlab.gcv` provides direct access to the underlying NetCDF
file using the Python `netCDF4 module
<http://code.google.com/p/netcdf4-python/>`_::
traj = gcv()
dataset = traj.reader.netcdf_file
You can list the variables stored in `dataset` with::
print dataset.variables.keys()
To plot the energy release rate `G` as a function of simulation time,
you could do::
plot(dataset.variables['time'], dataset.variables['G'])
You should see that the energy release rate increases at a roughly
constant rate before stopping at constant value when the crack starts
to move (the increase is not linear since is is actually the `strain`
that we increment at a constant rate).
The following plot shows the evolution of `G` (blue) and of the
position of the crack (red; stored as `crack_pos_x`). Note that a
second vertical axis can be produced with the
:func:`~matplotlib.pyplot.twinx` function.
.. image:: energy-release-rate-crack-position.png
:align: center
:width: 600
In this case the crack actually arrests for a while at around :math:`t
= 6` ps. This is another characteristic feature of non-brittle
fracture, indicating that our simulation is failing to match well
with experiment. According to Griffith's criterion, fracture should
initiate at :math:`2\gamma \sim 2.7` J/m\ :superscript:`2`, whereas we
don't see any motion of the crack tip until :math:`G \sim 11` J/m\
:superscript:`2`. How much of this difference do you think is due to
the high strain rate and small system used here, and how much to the
choice of interatomic potential? How would you check this?
.. _plot_temperature:
Temperature and velocity analysis
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Using the method above, plot the evolution of the temperature during
your simulation. Here is another example plot, with the temperature
shown in blue and the crack position in red.
.. image:: temperature-crack-position.png
:align: center
:width: 600
You will see that lots of heat is produced once the crack starts to
move, indicating that the system is far from equilibrium. This is
another sign that our system is rather small and our strain rate is
rather high. How could this be addressed? Do you think an NVT
simulation would be more realistic? What problems could adding a
thermostat introduce?
If you have time, you could compare how well the atomic velocities
match the expected Maxwell-Boltzmann distribution of atomic
velocities, given by
.. math::
f(v)\,\mathrm{d}v = 4 \pi \left( \frac{m}{2 \pi k_B T} \right)^{3/2} v^2 \exp \left[ -\frac{mv^2}{2 k_B T} \right] \mathrm{d}v
Here's a Python function which implements this formula::
def max_bolt(m,T,v):
"Maxwell-Boltmann distribution of speeds at temperature T for particles of mass m"
return 4*pi*(m/(2*pi*units.kB*T))**(3.0/2.0)*(v**2)*exp(-m*v**2/(2*units.kB*T))
We can average the atomic speeds in the last 50 frames in our
trajectory and use the speeds data to produce a histogram::
m = traj[-1].get_masses()[0] # Mass of a Si atom
T = traj[-1].info['temperature'] # Temperature at end of simulation
v = traj.reader.netcdf_file.variables['momenta'][-50:,:,:]/m # Get velocities
s = sqrt((v**2).sum(axis=2)) # Speeds are magnitude of velocities
hist(s.reshape(-1), 20, normed=True, alpha=0.5) # Draw a histogram
ss = linspace(0., s.max(), 100) # Compare with Maxwell-Boltzmann distrib
plot(ss, max_bolt(m,T,ss), lw=2)
.. image:: crack-max-bolt-distrib.png
:align: center
:width: 600
.. _arsf:
Atom-resolved strain tensor
^^^^^^^^^^^^^^^^^^^^^^^^^^^
The virial stress expression above is only valid when averaged over
time and space, so this method of calculating per-atom stresses can
lead to unphysical oscillations [Zimmerman2004]_. One alternative is the
atom-resolved strain tensor, which allows the strain, and hence stress,
fields to be evaluated at the atomistic scale facilitating direct
comparisons with elasticity theory results [Moras2010]_.
A definition of the atom-resolved strain tensor can be obtained for
all the four-fold coordinated atoms in the tetrahedral structure (all
other atoms are assigned zero strain) by comparing the atomic
positions with the unstrained crystal. The neighbours of each atom are
used to define a local set of cubic axes, and the deformations along
each of these axes are combined into a matrix :math:`E` describing the
local deformation:
.. math::
E = \left(\begin{array}{ccc}
| & | & | \\
\mathbf{e}_{1} & \mathbf{e}_{2} & \mathbf{e}_{3} \\
| & | & |
\end{array}\right)
where, for example :math:`\mathbf{e}_{1}` is the relative deformation
along the first cubic axis. To compute the local strain of the atom,
we need to separate this deformation into a contribution due to
rotation and one due to strain. This can be done by finding the polar
decomposition of :math:`E`, by writing :math:`E` in the form :math:`E
= SR` with :math:`R` a pure rotation and :math:`S` a symmetric matrix.
Diagonalising the product :math:`EE^T` allows :math:`R` and :math:`S`
to be calculated. The strain components :math:`\epsilon_{xx}`,
:math:`\epsilon_{yy}`, :math:`\epsilon_{zz}`, :math:`\epsilon_{xy}`,
:math:`\epsilon_{xz}` and :math:`\epsilon_{yz}` can then be calculated
by rotating :math:`S` to align the local cubic axes with the Cartesian
axes:
.. math::
R^T S R = I + \epsilon = \left(\begin{array}{ccc}
1 + \epsilon_{xx} & \frac{1}{2}\epsilon_{xy} & \frac{1}{2}\epsilon_{xz} \\
\frac{1}{2}\epsilon_{xy} & 1 + \epsilon_{yy} & \frac{1}{2}\epsilon_{yz} \\
\frac{1}{2}\epsilon_{xz} & \frac{1}{2}\epsilon_{yz} & 1 + \epsilon_{zz}
\end{array}\right).
Finally if we assume linear elasticity applies, the atomistic stress
can be computed simply as :math:`\bm\sigma = C \bm\epsilon` where
:math:`C` is the :math:`6\times6` matrix of elastic constants.
The :class:`~quippy.elasticity.AtomResolvedStressField` class
implements this approach. To use it to calculate the stress in your
`crack_slab` Atoms object, you can use the following code::
arsf = AtomResolvedStressField(bulk=si_bulk)
crack_slab.set_calculator(arsf)
ar_stress = crack_slab.get_stresses()
Colour your atoms by the :math:`\sigma_{yy}` component of the
atom-resolved stress field, and compare with the local virial stress
results. Add the atom resolved :math:`\sigma_{yy}` values along
:math:`y = 0` to your plot. Do you notice any significant differences?
Repeat the minimisation of the crack slab with a lower value of
`relax_fmax` (e.g. :math:`1 \times 10^{-3}` eV/A). Do the stress
components computed using the two methods change much?
.. You can also use the :func:`~quippy.crack.plot_stress_fields` function
.. to plot the atom-resolved and Irwin continuum near-tip stress fields,
.. and the residual error between them after fitting.
When you are ready, proceed to :ref:`step3`.
3rdparty/quip/doc/Tutorials/adaptive-qmmm-step3.rst deleted 100644 → 0
View file @ 23bcf97f
.. _step3:
Step 3: LOTF hybrid MD simulation of fracture in Si
===================================================
In the final part of this tutorial, we will be extending our previous script for
classical molecular dynamics to carry out an adaptive QM/MM simulation of
fracture using the 'Learn on the Fly' (LOTF) scheme.
You will need the ``run_crack_classical.py`` script from :ref:`step2`. If you
don't have it, you can :download:`download it here <run_crack_classical.py>`,
and the ``crack.xyz`` input file from :ref:`step1`, which you
can :download:`also download here <crack.xyz>`.
3.1 Initialisation of the QM/MM atomic system (20 minutes)
----------------------------------------------------------
Import the relevant modules
^^^^^^^^^^^^^^^^^^^^^^^^^^^
Make a copy of your ``run_crack_classical.py`` script and name it
``run_crack_lotf.py``. Add the following extra import statements after those
that are already there::
from quippy.potential import ForceMixingPotential
from quippy.lotf import LOTFDynamics, update_hysteretic_qm_region
Definition of the simulation parameters
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. _parameters3:
Next, we need to add some additional parameters specifically for the
QM/MM simulation. Again, insert them in ``run_crack_lotf.py``, below the
existing parameters ::
qm_init_args = 'TB DFTB' # Initialisation arguments for QM potential
qm_inner_radius = 8.0*units.Ang # Inner hysteretic radius for QM region
qm_outer_radius = 10.0*units.Ang # Inner hysteretic radius for QM region
extrapolate_steps = 10 # Number of steps for predictor-corrector
# interpolation and extrapolation
The setup of the atomic structure and of the constraints is exactly the same as
before, so leave these sections of your script unchanged.
Setup of the QM and QM/MM potentials
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For the QM/MM simulation, we first need to initialise the classical SW potential
(`mm_pot`) and the quantum-mechanical one (`qm_pot`). The two Hamiltonians then need
to be combined into a hybrid QM/MM potential (`qmmm_pot`), which mixes the QM
and MM forces.
Leave the initialisiton of the SW classical potential as it is. After this, we
want to add some lines of code to setup the QM potential. Using the same
:class:`~quippy.potential.Potential` class, we initialise now the Density
functional tight binding (DFTB) potential. This is done by passing the new QM
`qm_init_args` as the `init_args` parameter and the same XML file as before for
the `param_filename` to the Potential constructor (note that the single file
``params.xml`` contains parameters for both the SW and DFTB potentials)::
qm_pot = ... # Initialise DFTB potential
The QM/MM potential is constructed using quippy's
:class:`quippy.potential.ForceMixingPotential` class. Here, `pot1` is
the low precision, i.e. MM potential, and `pot2` is the high
precision, i.e. QM potential. `fit_hops` is the number of hops used to
define the fitting region, `lotf_spring_hops` defines the length of
the springs in the LOTF *adjustable potential*, while the hysteretic
buffer options control the size of the :ref:`hysteretic <hysteretic>` :ref:`buffer
region <buffer>` region used for the embedded QM force calculation. ::
qmmm_pot = ForceMixingPotential(pot1=mm_pot,
pot2=qm_pot,
atoms=atoms,
qm_args_str='single_cluster cluster_periodic_z carve_cluster '+
'terminate cluster_hopping=F randomise_buffer=F',
fit_hops=4,
lotf_spring_hops=3,
hysteretic_buffer=True,
hysteretic_buffer_inner_radius=7.0,
hysteretic_buffer_outer_radius=9.0,
cluster_hopping_nneighb_only=False,
min_images_only=True)
The `qm_args_str` argument defines some parameters which control how
the QM calculation is carried out: we use a single cluster, periodic
in the `z` direction and terminated with hydrogen atoms. The positions
of the outer layer of buffer atoms are not randomised.
Change the line which sets the Atoms calculator to use the new
`qmmm_pot` Potential::
atoms. ... # Set the calculator
Set up the initial QM region
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Now, we can set up the list of atoms in the initial QM region using
the :func:`~quippy.lotf.update_hysteretic_qm_region` function, defined
in quippy. Here we need to provide the :class:`~.Atoms` system, the
centre of the QM region (i.e. the position of the crack tip), and the
the inner and outer radius of the :ref:`hysteretic <hysteretic>` QM
region. Note that the `old_qm_list` attribute must be an empty list
(``[]``) in this initial case::
qm_list = ... # Define the list of atoms in the QM region
The list needs to be attached to the `qmmm_pot` using the
:meth:`~quippy.potential.ForceMixingPotential.set_qm_atoms` method::
qmmm_pot. ... # Attach QM list to calculator
Milestone 3.1
^^^^^^^^^^^^^
Your ``run_crack_lotf.py`` script should look something
like :download:`run_crack_lotf_1.py`.
At this point you should run your script and check the initial QM region. For
testing, you should add a couple of temporary lines to force the script to
finish after setting the QM region and before repeating the classical MD::
import sys
sys.exit(0)
To visualise the initial QM region, you can type the following directly into
your `ipython` session (remember to do a ``from qlab import *`` first if you
haven't already)::
view(atoms)
aux_property_coloring(qmmm_pot.get_qm_atoms())
.. image:: crack-initial-qm-region.png
:align: center
:width: 600
In the image above, the red atoms are QM and the blue atom classical.
Internally, this list is actually saved as a :attr:`property
<~quippy.atoms.Atoms.properties>` inside the Atoms object named ``"hybrid"``,
which can also be displayed with ``aux_property_coloring("hybrid")``
3.2 Setup and run the adaptive QM/MM MD (20 minutes)
----------------------------------------------------
Initialising the Dynamics
^^^^^^^^^^^^^^^^^^^^^^^^^
The definition of the initial temperature of the system should be left as
in :ref:`Step 2 <step2>`. Don't forget to remove the temporary lines added above which
quit the script after setting up the initial QM region!
Instead of a traditional dynamics in the NVE ensemble, let's change the code to
use :ref:`LOTF predictor-corrector dynamics <lotf>`, using
the :class:`quippy.lotf.LOTFDynamics` class instead of
the :class:`~md.verlet.VelocityVerlet` class. We need to pass the following
arguments: `atoms`, `timestep`, `extrapolate_steps` (see :ref:`Parameters
section <parameters3>`)::
dynamics = ... # Initialise the dynamical system
The logger and crack tip movement detection functions can be left almost exactly
as before for now: we just need to make a small change to
the :func:`printstatus` function so to distinguish between extrapolation and
interpolation:
Change the line::
atoms.info['label'] = 'D' # Label for the status line
to::
atoms.info['label'] = dynamics.state_label # Label for the status line
This uses the :attr:`~quippy.lotf.LOTFDynamica.state_label` attribute to print
an ``"E"`` at the beginning of the logger lines for extrapolation and an ``"I"``
for interpolation.
Updating the QM region
^^^^^^^^^^^^^^^^^^^^^^
We need to define a function that updates the QM region at the
beginning of each extrapolation cycle. As before, we need to find the
position of the crack tip and then update the :ref:`hysteretic
<hysteretic>` QM region. Note that now a previous QM region exists and
its atoms should be passed to the
:func:`~quippy.lotf.update_hysteretic_qm_region` function. The current
QM atom list can be obtained with the
:meth:`quippy.potential.ForceMixingPotential.get_qm_atoms` method. To
find the crack position, use
:func:`~quippy.crack.find_crack_tip_stress_field` as before, but pass
the MM potential as the calculator used to calculated the stresses
(force mixing potentials can only calculate forces, not per-atom
stresses; we will check later that the classical stress is
sufficiently accurate for locating the crack tip)::
def update_qm_region(atoms):
crack_pos = ... # Find crack tip position
qm_list = ... # Get current QM atoms
qm_list = ... # Update hysteretic QM region
qmmm_pot. ... # Set QM atoms
dynamics.set_qm_update_func(update_qm_region)
Writing the trajectory
^^^^^^^^^^^^^^^^^^^^^^
Finally, we want to save frames to the trajectory every `traj_interval` time
steps but, this time, only during the interpolation phase of the
predictor-corrector cycle. To do this, we first initialise the trajectory file
(see :func:`~quippy.io.AtomsWriter`), and then define a function that only
writes to the trajectory file if the state of the dynamical systems is
`Interpolation`::
trajectory = ... # Initialise trajectory using traj_file
def traj_writer(dynamics):
if dynamics.state == LOTFDynamics.Interpolation:
trajectory.write(dynamics.atoms)
As before, we attach this function to the dynamical system, passing
`traj_interval` and and extra argument of `dynamics` which gets passed along to the
`traj_writer` function (see the :meth:`~quippy.lotf.LOTFDynamics.attach`
method)::
dynamics. ... # Attach traj_writer to dynamics
Now, we can simply run the dynamics for `nsteps` steps::
dynamics. ... # Run dynamics for nsteps
If you are interested in seeing how the LOTF predictor-corrector cycle is
implemented, look at the documentation and `source code
<_modules/quippy/lotf.html#LOTFDynamics.step>`_ for the
:meth:`quippy.lotf.LOTFDynamics.step` routine.
Milestone 3.2
^^^^^^^^^^^^^
The finished version of the ``run_crack_lotf.py`` script should look something
like :ref:`run_crack_lotf`. To clearly show the differences with respect to the
classical MD script, here is a :download:`patch
<run_crack_classical_lotf.patch>` which could be used to convert the classical
script into the LOTF one.
3.3 Visualisation and Analysis (as time permits)
------------------------------------------------
Predictor/corrector dynamics output file
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Let's first take a moment to look at the output of the script for the first
predictor/corrector cycle. Here we go through some example output, yours should
be similar. First there are a few lines about the initialisation of the system,
and then we get the results of the initial LOTF adjustable potential
optimisation::
Loading atoms from file crack.xyz
Fixed 240 atoms
25 atoms selected for quantum treatment
update_qm_region: QM region with 25 atoms centred on [-30.60517303 0.08401087 0. ]
Adding default springs
Got 1484 springs
Number of force components: 297
Number of parameters: 1484
Optimising 1484 adjustable parameters
RMS force component error before optimisation : .05630875465645784
Max force component error before optimisation : .34841292159055509
Using SVD for least squares fit, eigenvalue threshold = .00000000010000000
RMS force component error after optimisation : 0.27E-02
Max force component error after optimisation : 0.61E-02
Max abs spring constant after optimisation : 0.45E-01
You can see that before adjusting the parameters, the QM and classical potentials
differed by a maximum of 0.35 eV/A, with an RMS difference of 0.06 eV/A - in
this case the SW potential is actually doing a rather respectable job. After the
fit, which is this case involved 1484 spring parameters to fit 297 force
component, the force differences are of course much smaller.
Next we start the first predictor/corrector cycle. First we update the QM
region, and remap the adjustable potential to take account of any changes
since last time::
25 atoms selected for quantum treatment
update_qm_region: QM region with 25 atoms centred on [-30.6048418 0.08377744 0. ]
Adding default springs
Got 1484 springs
Number of force components: 297
Number of parameters: 1484
As this is the first step, there were no changes, so no re-optimisation is
required. Next we carry out 10 steps of extrapolation, with constant LOTF
adjustable parameters. During this time the strain is incremented as normal::
State Time/fs Temp/K Strain G/(J/m^2) CrackPos/A D(CrackPos)/A
---------------------------------------------------------------------------------
E 1.0 553.716406 0.08427 5.0012 -30.61 (-0.00)
E 2.0 547.749233 0.08428 5.0024 -30.61 (-0.01)
E 3.0 535.952151 0.08429 5.0036 -30.62 (-0.01)
E 4.0 518.731103 0.08430 5.0047 -30.63 (-0.02)
E 5.0 496.675925 0.08431 5.0059 -30.63 (-0.03)
E 6.0 470.538607 0.08432 5.0071 -30.64 (-0.04)
E 7.0 441.205418 0.08433 5.0083 -30.65 (-0.05)
E 8.0 409.663780 0.08434 5.0095 -30.66 (-0.06)
E 9.0 376.965040 0.08435 5.0107 -30.67 (-0.07)
E 10.0 344.184506 0.08436 5.0119 -30.69 (-0.08)
At the end of the extrapolation, it's time for a QM force evaluation
and another fit. Now the force errors before fitting are a little
larger, but the fit is still very good::
Optimising 1484 adjustable parameters
RMS force component error before optimisation : .10494977522791650
Max force component error before optimisation : .48515966905523733
Using SVD for least squares fit, eigenvalue threshold = .00000000010000000
RMS force component error after optimisation : 0.37E-02
Max force component error after optimisation : 0.96E-02
Max abs spring constant after optimisation : 0.83E-01
We next return to the initial dynamical state and re-run the dynamics,
interpolating between the optimised parameters at the two ends of the cycle.
Note that the strain is also returned to the initial value at :math:`t = 0`, and
that the temperature after one step exactly matches the interpolation phase
(since the forces and velocities at :math:`t = 0` are identical for
extrapolation and interpolation)::
State Time/fs Temp/K Strain G/(J/m^2) CrackPos/A D(CrackPos)/A
---------------------------------------------------------------------------------
I 1.0 553.716406 0.08427 5.0012 -30.65 (-0.04)
I 2.0 547.759567 0.08428 5.0024 -30.65 (-0.05)
I 3.0 535.982832 0.08429 5.0036 -30.66 (-0.05)
I 4.0 518.791314 0.08430 5.0047 -30.66 (-0.06)
I 5.0 496.773542 0.08431 5.0059 -30.67 (-0.07)
I 6.0 470.679783 0.08432 5.0071 -30.68 (-0.08)
I 7.0 441.394231 0.08433 5.0083 -30.69 (-0.09)
I 8.0 409.901969 0.08434 5.0095 -30.70 (-0.10)
I 9.0 377.251837 0.08435 5.0107 -30.71 (-0.11)
I 10.0 344.516566 0.08436 5.0119 -30.73 (-0.12)
To continue from here, we simply go back to the extrapolation phase and then
repeat the entire cycle.
QM active and buffer regions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Trajectory analysis
^^^^^^^^^^^^^^^^^^^
Open your new trajectory :ref:`as before <visualisation2>`, using the
:func:`~qlab.view` function from within a new `ipython` session, and
visualise the QM region by colouring the atoms using the
``hybrid_mark`` :attr:`property <quippy.atoms.Atoms.properties>` ::
aux_property_coloring("hybrid_mark")
This property is used internally to identify which atoms are used for the QM
active and buffer regions:
.. image:: crack-hybrid-mark.png
:align: center
:width: 600
The central green atoms have ``hybrid_mark == HYBRID_ACTIVE_MARK``, and they are
the atoms for which QM forces are used to propagate the dynamics. Classical
forces are used for all other atoms, including the red buffer region, where
``hybrid_mark == HYBRID_BUFFER_MARK``. As explained :ref:`above <buffer>`, the
purpose of the buffer region is to give accurate QM forces on the active atoms.
.. _cluster:
If you want to see the actual cluster used for carrying out the embedded DFTB
calculation, you could use the :func:`~quippy.clusters.create_cluster_simple`
function together with the same `args_str` cluster options defined above::
cluster = create_cluster_simple(gcat(),
args_str=("single_cluster cluster_periodic_z carve_cluster "
"terminate cluster_hopping=F randomise_buffer=F"))
view(cluster)
Colouring the cluster by coordination (press `k`) can be useful to check that
all cut bonds have been correctly passivated by hydrogen atoms:
.. image:: lotf-crack-cluster.png
:align: center
:width: 600
Comparison between classical and LOTF dynamics
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Step through your trajectory with the `Insert` and `Delete` keys to see what
happens in the LOTF dynamics. As before, you can jump to the end with
`Ctrl+Delete`. You should find that the dynamics is very different to the
classical case.
Check if the QM region is following the moving crack properly by looking at the
``hybrid_mark`` property. If you repeat the analysis of the :ref:`stress field
<stress_analysis>` carried out in :ref:`Step 2 <step2>`, you should find that
the :ref:`time averaged stress field <time_avg_stress>` is strongly concentrated
on the sharp crack tip. It is this stress field which is used
by :func:`~quippy.crack.find_crack_tip_stress_field` to follow the crack tip,
and hence to update the set of atoms in the QM region.
Here is a movie of a typical LOTF simulation on the :math:`(111)` cleavage
plane. To colour the QM atoms dark blue, we passed
the :func:`~qlab.highlight_qm_region` function as the `hook` argument
to :func:`~qlab.render_movie`:
.. video:: lotf-111 640 360
During the LOTF dynamics, the time-averaged stress field smoothly tracks the
crack tip, as can be seen in this movie, where atoms are coloured by
their :math:`\sigma_{yy}` component:
.. video:: elastic 640 360
And here is a head-to-head comparison of SW and LOTF dynamics:
.. video:: classical-vs-lotf 640 720
Fracture initiates much earlier in the LOTF case, i.e. at a much reduced energy
release rate, and is much more brittle, with none of the artificial plasticity
seen with the classical potential alone.
Note that if you continue the LOTF dynamics, however, we may see some defects in
the frature surface after the crack has propagated for a few nm. These are
associated with the relatively small system and high strain rate we are using
here, which leads to fracture at high energies and possibly to high speed
fracture instabilities [Fineberg1991]_. If you have time you can investigate
this in the :ref:`extension task on size and strain rate effects
<system_size_and_strain_rate>`.
.. video:: clas-vs-lotf 640 720
Although it is beyond the scope of this tutorial, you might be interested to
know that using an overall larger system, bigger QM region and lower strain rate,
as well as changing the Hamiltonian from DFTB to DFT-GGA, removes all of these
defects, recovering perfectly brittle fracture propagation. The DFT model also
gives an improved description of the fracture surfaces, which reconstruct to
form a Pandey :math:`\pi`\ -bonded chain, with it's characteristic alternating
pentagons and heptagons:
.. video:: silicon-111-dft-1400 640 360
.. _plot_G_and_crack_pos_x_lotf:
Evolution of energy release rate and crack position
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
If you follow the :ref:`previous approach <plot_G_and_crack_pos_x>` to plot the
energy release rate `G` and crack position `crack_pos_x` variables during your
LOTF simulation, you should find that the crack now advances monotonically, with
a constant crack velocity of around 2500 m/s, and at about half the energy
release rate of the classical case (6 J/m\
:superscript:`2` vs 12 J/m\ :superscript:`2`).
.. image:: lotf-energy-release-rate-crack-position.png
:align: center
:width: 600
For comparison, here is the classical plot again:
.. image:: energy-release-rate-crack-position.png
:align: center
:width: 600
You should find that the :ref:`temperature <plot_temperature>` still goes up,
but more gently than in the classical case, since the flow of energy to the
crack tip is closer to the energy consumed by creating the new surfaces. Some
heat is generated at the QM/MM border; usually this would be controlled with a
gentle Langevin thermostat, which we have omitted here in the interests of
simplicity.
.. _low_speed_instability:
Low speed instability on the (111) cleavage plane
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
If you are lucky, you may see the formation of a crack tip reconstruction
consisting of a 5 and a 7 membered ring on the lower fracture surface,
related to the Pandey surface reconstruction.
.. image:: lotf-crack-step-1.png
:width: 400
:align: center
This reconstruction can cause cracks to take a step down by one atomic layer,
which over time can build up via positive feedback mechanism into an
experimentally observable phenomena [Kermode2008]_.
.. image:: lotf-crack-step-2.png
:width: 400
:align: center
.. _extension_tasks:
.. _pred_corr_error:
3.4 Checking the predictor/corrector force errors (optional)
------------------------------------------------------------
Add `check_force_error=True` to the :class:`~quippy.lotf.LOTFDynamics`
constructor. This causes the LOTF routines to do a reference QM force evaluation
at every timestep (note that these extra QM forces are not used in the fitting,
so the dynamical trajectory followed is the same as before).
When checking the predictor/corrector errors, you need to disable the updating of
the QM region by commenting out the line::
dynamics.set_qm_update_func(update_qm_region)
Let's create a logfile to save the force errors at each step during
the interpolation and extrapolation. Add the following code before the
:meth:`dynamics.run()` call::
def log_pred_corr_errors(dynamics, logfile):
logfile.write('%s err %10.1f%12.6f%12.6f\n' % (dynamics.state_label,
dynamics.get_time()/units.fs,
dynamics.rms_force_error,
dynamics.max_force_error))
logfile = open('pred-corr-error.txt', 'w')
dynamics.attach(log_pred_corr_errors, 1, dynamics, logfile)
Finally, change the total number of steps (via the `nsteps` parameter) to a much
smaller number (e.g. 200 steps), close the logfile after the ``dynamics.run()``
line::
logfile.close()
Once the dynamics have run for a few LOTF cycles, you can plot the results with
a shell script called ``plot_pred_corr_errors.py``::
plot_pred_corr_errors.py -e 10 pred-corr-error.txt
The ``-e 10`` argument is used to specify the number of extrapolate steps. This
produces a set of four plots giving the RMS and maximum force errors during
extrapolation and interpolation:
.. image:: lotf_check_force_error.png
:align: center
:width: 600
Note that the scale is different on the extrapolation and interpolation plots!
Try varying the `extrapolate_steps` parameter and seeing what the effect on
force errors is. What is the largest acceptable value? You could also try
changing the `lotf_spring_hops` and `fit_hops` parameters, which control the
maximum length of the corrective springs added to the potential and the size of
the fit region, respectively.
Milestone 3.4
^^^^^^^^^^^^^
Here is a final version of the ``run_crack_lotf.py`` script including
checking of the force errors: :download:`run_crack_lotf.py`.
Further extension tasks
-----------------------
.. _qm_region_size:
QM region size
^^^^^^^^^^^^^^
Investigate the effect of increasing the QM region size, controlled by the
`qm_inner_radius` and `qm_outer_radius` parameters. When does the behaviour
converge qualitatively? What does this say about the size of the 'process zone'
in silicon?
.. _buffer_region_size:
Buffer region size
^^^^^^^^^^^^^^^^^^
We have used a hysteretic buffer region from 7 A to 9 A. How would you check if
this is sufficient? What criteria need to be satisfied for our results to be
considered to be converged with respect to buffer region size?
.. _freund:
Crack energy-speed relationship
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Try varying the flow of energy to the crack tip by changing the `initial_G`
parameter used when making the crack system in :ref:`step1`. How does this
affect the speed of the crack?
.. _other_orientations:
Other crack orientations
^^^^^^^^^^^^^^^^^^^^^^^^
Return to the beginning of :ref:`step1` and try classical and/or LOTF dynamics
(which will actually probably be faster!) on the :math:`(110)` surface. Do you
see any major differences? Can you find any dynamic fracture instabilities?
.. _system_size_and_strain_rate:
System size and strain rate effects
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
What is the effect of changing the system size on the critical energy
release rate for fracture? How would you converge with respect to this
parameter? Do you think experimental length scales can be reached? If
not, does it matter? Think about how the choice of loading geometry
helps here.
As well as finite size effects, and perhaps more severely, we are limited in the
time scales that can be accessed by our fracture simulations, especially when
using a QM method to describe the crack tip processes. Are there any scaling
relations that can help us out here? How would you estimate the effect of the
artificially high strain rate we have been forced to impose here.
3rdparty/quip/doc/Tutorials/adaptive-qmmm-theory.rst deleted 100644 → 0
View file @ 23bcf97f
Theoretical background
======================
.. _theory1:
Introduction to Atomic-scale Modelling of Fracture
--------------------------------------------------
Fracture is the main lifetime limiting failure mode of many materials, from
metals to ceramics and glasses. There are two broad classes of
fracture: *brittle* and *ductile*.
.. figure:: brittle-vs-ductile.png
:align: center
:width: 600
Brittle and ductile failure of steel, at low and high temperature
respectively. Image credit: `Internet Electron Microscopy website, UMIST
<http://pwatlas.mt.umist.ac.uk/internetmicroscope/micrographs/failure/ductile-steel.html>`_
In brittle fracture, there is no plastic deformation, and failure occurs along
particular low-energy cleavage planes. In contrast, ductile materials such as
steel tend to fail with a large amount of plastic deformation (or 'necking'),
pulling apart to leave rough fracture surfaces. The same patterns are evident
upon closer examination: brittle cracks are found to be atomically sharp,
proceeding bond by bond, while ductile cracks remain rough at the microscale,
driven by the formation and coalescence of microscopic voids.
.. figure:: brittle-vs-ductile-microscale.png
:align: center
:width: 600
Brittle and ductile failure mechanisms, for silica and aluminium respectively.
Image credits: C. Marlière and A. Weck.
Silicon is known to be an ideally brittle material, as shown by the image below.
Below the brittle to ductile transition, at around :math:`500^\circ~\mathrm{C}`, silicon
cleaves atomically smooth fracture surfaces (left and centre panel; sample in
central panel is tilted to show the crack front). At higher temperatures,
fracture is ductile, with the emission of multiple dislocations from the tip
(right panel).
.. _si_tem_images:
.. figure:: lawn-fracture-silicon.png
:align: center
:width: 600
Transmission electron microscopy (TEM) images of a crack in silicon.
Image reproduced from [Lawn1993]_.
.. _theory_griffith:
Fracture Mechanics
^^^^^^^^^^^^^^^^^^
The study of fracture mechanics dates back around 100 years, to
Griffith, who first proposed a thermodynamic energy balance criteria
to understand when cracks will propagate [Griffith1921]_. The key idea
is that stress concentrates at pre-existing flaws: this was motivated
by the observation that materials break at much lower loads than the
theoretical stress needed to break their chemical bonds For example,
in glass, the theoretical strength is given by
.. math::
\sigma_{theoretical} = \sqrt{\frac{E\gamma}{a}} \sim 10,000\;\mathrm{MPa}
where :math:`E`, :math:`\gamma` and :math:`a` are the Young's modulus, surface
energy and bond length, respectively.
For a slowly-moving crack of length `c` in an infinite plane, the
well-known Griffith criterion for fracture propagation is based on thermodynamic
energy balance between the release of elastic energy in an area
proportional to `c`\ :superscript:`2` and the cost of creating new surfaces,
which is proportional to `c`, as illustrated below.
.. image:: griffith-criterion.png
:align: center
:width: 600
This leads to a Griffith fracture strength for glass of
.. math::
\sigma_{fracture} = \sqrt{\frac{E\gamma\rho}{4ac}} \sim 100 \;\mathrm{MPa}
which is much lower than the theoretical strength. Here the additional
parameters are :math:`\rho`, the radius of curvature and the crack
length `c`. The effect of stress concentration increases for sharper and
longer cracks.
The Griffith criterion leads to a critical length :math:`c_0` below which it is
not energetically favourable for cracks to grow, since the elastic energy
released does not exceed the surface energy cost. Below :math:`c_0`, cracks
prefer to close up, meaning that not all flaws are unstable. This explains why
it makes sense to measure the length of cracks e.g. on aeroplanes, so that small
flaws can be identified and treated before they become critical.
In fracture mechanics it is common to use the energy release rate :math:`G` to
describe the flow of energy to a crack tip. :math:`G` is the driving force for
crack propagation. It is defined by
.. math::
G = - \frac{\partial U_E}{\partial c}
where :math:`U_E` is the total strain energy and `c` is the crack length.
The Griffith criterion can be reformulated in terms of :math:`G` to show that
crack propagation becomes favourable when
.. math::
G > 2\gamma
where :math:`\gamma` is the surface energy density, i.e. when the energy flow to
the crack tip exceeds the cost of creating two new surfaces.
.. _theory_atomic_fracture:
Atomic scale modelling of fracture
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Much work has been done to understand fracture at the continuum level
(see, e.g. [Freund1998]_ and [Lawn1993]_), but here we would like to
simulate fracture at the atomic level, to examine the combined effects
of stress and chemistry ('chemomechanics'). A first approach would be
to use classical interatomic potentials to carry out molecular
dynamics (MD). However, as we will see in :ref:`step2`, most classical
potentials fail to accurately reproduce brittle fracture. This is due
to stress concentration which has been shown to diverge as
:math:`\sigma \sim 1/\sqrt{r}` near a crack tip [Irwin1948]_, leading
to anharmonic stretching and rupture of bonds, which is typically not
well captured by simple interatomic potentials.
.. _irwin_sig_yy:
.. figure:: irwin-sig-yy.png
:width: 300
:align: center
Irwin near-field solution for :math:`\sigma_{yy}` for a singular crack.
Black is zero stress and yellow very high stress; note the divergence at the
crack tip.
Most potentials overestimate what is called the *lattice trapping
barrier*, the energy barrier to bond breaking at a crack tip that
arises from the periodicity of the crystalline lattice (in contrast to
continuum methods where the crack tip advances continuously). This
means that when fracture eventually does occur, there is too much
energy available, which is then dissipated by a variety of plasticity
mechanisms such as dislocation emission. This leads to results in
contrast with the expected brittle behaviour.
Interestingly, however, continuum theories and simple potentials agree
with one another until surprisingly close to the crack tip (~1nm),
where non-linear effects become important, as illustrated in the
figure below.
.. figure:: atomistic-vs-continuum.png
:width: 400
:align: center
Atomic and continuum calculations for the stress along the line ahead of a
crack tip in silicon. Agreement is excellent beyond ~ 2 nm from the tip.
Reproduced from G. Singh, J.R. Kermode, A. De Vita, R.W. Zimmerman,
*in prep.* (2013).
The region where atomistic and continuum theories disagree is the
non-linear *process zone*, where chemically interesting things are happening.
Here, we would like to treat this region at a quantum mechanical (QM) level.
.. _theory_multiscale:
Coupled multiscale approach
^^^^^^^^^^^^^^^^^^^^^^^^^^^
QM approaches such as density functional theory (DFT) do correctly
describe bond-breaking in silicon. However, the strong bidirectional
coupling between bond-breaking at the crack tip and the long-range
stress field driving fracture necessitates a multiscale approach. The
reason a fully DFT approach is not practical is that the boundaries of
the model system must be placed far enough away from the crack tip not
to affect the results, which means that large cells containing tens to
hundreds of thousands of atoms are needed. This exceeds the current
capabilities of most QM approaches. Fracture is perhaps the
archetypical coupled multiscale problem, with thousands of atoms
contributing to the elastic relaxation of the near-tip region. We will
describe how classical and QM descriptions can be coupled to study
problems in fracture using the 'Learn on the Fly' (LOTF) approach in
:ref:`more detail <theory3>` later in this tutorial.
.. figure:: multiscale-coupling.png
:width: 500
:align: center
Hierarchy of materials modelling techniques, showing simultaneous coupling
of QM methods and empirical interatomic potentials. Image source: G. Csányi.
.. _thin_strip:
Thin strip loading geometry and elasticity theory
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
We will use the thin strip fracture loading geometry illustrated below for our
fracture simulations, where the load is applied by displacing the top and bottom
surfaces.
.. image:: thin-strip.png
:align: center
:width: 500
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 :math:`\delta`, so the applied strain is
:math:`\epsilon_0 = \delta / h`. Far ahead of the crack, the strip is
in uniaxial tension: :math:`\epsilon_{yy} \to \epsilon_0` as :math:`x
\to \infty`.
The stress far ahead of the crack is given by :math:`\sigma_{0} = E'
\epsilon_{0}`, and therefore the elastic energy per unit length and
per unit thickness far ahead of the crack tip is
.. math::
W_\infty = \frac{1}{2}E'\epsilon_0^2\cdot 2h = \frac{\delta^2 E'}{h}
where :math:`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
:math:`W_\infty` is effectively replaced by a strip with energy
density zero.
The energy supplied to the crack tip is therefore equal to :math:`W_\infty`,
so the energy release rate is simply
.. math::
G = W_\infty = \frac{\delta^2 E'}{h}
In our simulations we will use periodic boundary conditions in the :math:`z`
direction, so we have plane strain loading (:math:`\epsilon_{zz} = 0`), which
means that the effective Young's modulus :math:`E'` is given by
:math:`E/(1-\nu^2)`, where :math:`E` is the Young's modulus in the :math:`y` relevant
direction and :math:`\nu` is the Poisson ratio, so finally we have
.. _thin_strip_equation:
.. math::
G = \frac{E \delta^2}{(1- \nu^2)h} = \frac{E \epsilon_0^2 h}{1 - \nu^2}
We can see that, in order to relate the strain we apply to the system
with the energy release rate, we will need to calculate the Young's
modulus and Poisson ratio for our model material. We will see how to
do this from the elastic constant matrix :math:`C_{ij}` :ref:`below
<youngs_modulus_and_poisson_ratio>`.
Classical interatomic potentials for silicon
--------------------------------------------
The Stillinger-Weber [Stillinger1985]_ interatomic potential provides a fairly good
description of many properties of crystalline and amorphous silicon. Its
functional form is a two- and three-body expansion suitable for the
representation of highly directional covalent bonds between Si atoms:
.. math::
V = \sum_{ij} V_2(r_{ij}) + \sum_{ijk} V_3(r_{ij}, r_{ik},
\theta_{jik}) \\
V_2(r_{ij}) = A\left(\frac{B}{r_{ij}^4} - 1\right) \exp
\frac{1}{r-a}, \; r < a \\
V_3(r_{ij},r_{jk},\theta_{jik}) = \lambda \exp \left[
\frac{\gamma}{r_{ij} - a} + \frac{\gamma}{r_{jk} - a} \right] \left(\cos
\theta_{jik} + \frac{1}{3}\right)
In particular, the three-body term stabilises the ideal tetrahedral
structure with respect to all the other possible structures. The
parameters of the SW potential were determined by fitting experimental
data with the constraint that the diamond-like structure must be the
most stable periodic arrangement of particles al low pressures.
Although this potential has not been fitted to the Si elastic
constants, it gives a reasonable description of all of them. As we
will see during this step, however, the SW potential fails to describe
the brittle fracture of silicon. A number of interatomic potentials
have been developed to go beyond the basic description of Si provided
by the SW potential (e.g. [Swadener2002]_, [Vink2001]_ [Buehler2006]_,
[Pizzagalli2013]_), however, they are usually not sufficiently
transferable to provide a general description of the inherently
quantum-mechanical processes occurring at the tip of a crack. As a
consequence, treating the tip at the QM level using hybrid QM/MM
methods is necessary to perform accurate MD simulations of brittle
fracture in Si. Here, we will use the SW potential because its simple
functional form and its speed make it a suitable choice for a
multiscale QM/MM approach, where only an accurate description of the
Si crystal far from the crack tip is required.
In this step, we will use the SW potential to perform a classical MD
simulation of the crack propagation in the NVE ensemble. The velocity
Verlet scheme [FrenkelSmit2001]_ will be used to integrate Newton's
equations of motion.
.. _theory3:
QM/MM Force Mixing and the 'Learn on the Fly' scheme
----------------------------------------------------
In this last part of the tutorial, we will perform an accurate MD simulation of
Si fracture using the "Learn on the fly" (LOTF) hybrid QM/MM technique
[DeVita1998]_, [Csanyi2004]_. In the present case, all atoms that are not suitably described by
the [Stillinger1985]_ potential (our MM scheme) will be treated
quantum-mechanically with the DFTB method [Elstner1998]_. These atoms are
those in the vicinity of the crack tip, where highly strained Si-Si bonds are
present, and where formation and rupture of chemical bonds occurs during
crack propagation.
Standard QM/MM techniques, usually developed for biological systems, adopt
energy-based approaches. The total energy of the system is written as a combination
of the QM energy, the MM energy and a QM/MM term, often specifically devised
for a particular system, that takes care of the interaction between the two regions.
While this approach allows the definition of a total energy which is
conserved during the dynamics, the forces used to propagate the MD are not accurate
enough because of the spurious effects due to the presence of the boundary
between the QM and the MM regions. Moreover, the necessity to suitably "terminate"
the QM region, does not allow the QM region to move during the simulation, which
is, however, required if we want to follow the motion of a crack tip.
The LOTF scheme adopts a force-based scheme instead, which allows the
QM region to move during the MD simulation and accurate forces to be
calculated even at the boundaries of the two regions. While the
details of the scheme have been thoroughly presented in a number of
articles [DeVita1998]_, [Csanyi2004]_, [Csanyi2005]_, [Moras2010]_,
[Kermode2008a]_, [Bernstein2009]_, we will here briefly explain the
basic concepts that will allow us to perform the crack simulation.
.. _buffer:
Calculation of the forces: buffered force mixing
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The first problem we have to face, when performing a QM/MM MD simulation is to
calculate accurate forces. While the calculation of the MM forces, using a
short-ranged potential (e.g. the SW potential) is trivially solved by
inexpensively computing the MM forces on the whole system, the problem of
calculating accurate QM forces is more complicated. The ultimate goal of any QM/MM
approach is to perform a simulation where all the system behaves instantaneously
as if it were fully QM. In a standard QM/MM approach, however, the QM forces
would be calculated by carving the QM region out of the system and chemically
terminating the resulting Si dangling bonds, for instance with H atoms. Even by
using more complex termination strategies, it is generally not possible to
mimic the presence of the surrounding MM system, and the forces on the atoms
close to the artificially created "QM surface" are not accurate.
To solve this problem, at the expense of an increased computational
cost, we increase the size of the QM region by adding a "buffer
region", as illustrated below. The calculation of the QM forces is
carried out on a cluster made from the QM + buffer regions, after
suitable termination with H atoms.
.. figure:: buffer-region.png
:width: 600
:align: center
The finite buffer termination strategy. Force on the atoms in the buffer
region (dark grey) are discarded to give accurate QM or MM forces on all atoms
(right hand panel). Reproduced from [Bernstein2009]_.
The buffer region must be large enough to minimize the errors on the QM forces
due to the presence of the outer artificial surface. The size of the buffer
region can be determined through some tests (see :ref:`this extension task
<buffer_region_size>`) , and is typically around 8 A, or 4 neighbour hops, for
Si. This near-sightedness of quantum mechanics is ultimately due to the locality
of the density matrix in covalent materials.
Once accurate QM forces have been obtained, only the QM forces on the atoms
belonging to the original QM region are used in the MD. The QM forces on the
atoms in the buffer region, which are strongly affected by the presence of the
outer QM surface, are discarded and replaced by the MM forces (as illustrated
above). In this way, we can obtain good forces on all atoms in the
system. These forces can be used in the MD simulation, provided that the
conservation of the total momentum is restored. This can be enforced by
subtracting the (typically small) mean force, so that the final
QM/MM forces sum to zero.
It is important to have good elastic matching between the QM and MM models, so
that there is no discontinuity at the boundary. For simple materials, this can
usually be achieved by scaling the classical model in space and energy to match
the lattice constant and bulk modulus of the QM model (for simplicity we omit
this step in this tutorial, but the mismatch here is not too big).
.. _hysteretic:
Hysteretic selection of the QM active and buffer regions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A key advantage of the buffered force mixing approach is that the QM region can
be redefined during a simulation. This works because new atoms first enter the
QM calculation at the outer boundary of the buffer region, where they do not
have a large effect on the forces at the centre, and therefore changing the QM
region does not introduce large inconsistencies.
In this tutorial, the QM region will be updated during the MD simulation in
order to follow the motion of the crack tip. When selecting the atoms that
belong to the buffered QM region, it is important to minimise fluctuations in
the set of QM atoms. This should be done in order to minimise the spurious heat
generation due to atoms whose description changes from MM to QM, or vice-versa
[Bernstein2009]_. This change of description is inevitable when the QM region
moves with the crack tip. However, solutions have to be found to minimise
fluctuations due to oscillation.
A robust way to minimise fluctuations is to employ a "hysteretic" algorithm for
the QM selection process [Bernstein2009]_. In the context of a fracture
simulation, atoms have to come within a certain *inner
radius* :math:`r_\mathrm{-}` from the crack tip to become selected as QM active
atoms (see the picture below). In our case, atoms have to come within 8 A from
the crack tip to become part of the QM region. However, using the hysteretic
algorithm, these atoms will remain QM until they move further than the
*outer radius* :math:`r_\mathrm{+}` (where :math:`r_\mathrm{+} > r_\mathrm{-}`)
away from the crack tip. In our simulation, this outer radius will be 10 A. We
refer to the current set of QM atoms in as the *QM active region*.
.. figure:: hysteresis.png
:height: 200
:align: center
Hysteretic QM selection algorithm. For fracture simulations the black *active region*
can be reduced to a single point at the crack tip. Reproduced from [Bernstein2009]_.
As well as using hysteresis to select the QM active atoms for which QM
forces will be used, we can also use the hysteretic selection
algorithm to minimise fluctuations in the buffer region. These radii
apply to the distance from any of the QM active atoms, so the buffer
takes the form of a shell of constant width around the QM atoms. Here
we will use inner and outer buffer radii of 7 A and 9 A, respectively.
This leads to overall QM active + buffer clusters with a radius of
around 15 A, containing around 150 atoms, including terminating
hydrogens (see the :ref:`example cluster <cluster>` below).
Further tricks which can be used to stabilise the QM and buffer region include
growing the regions using bond hopping rather than distance criteria, and using
time-averaged positions [Bernstein2009]_. For simplicity, in this tutorial we
will use only the hysteretic selection technique.
.. _lotf:
LOTF adjustable potential and predictor-corrector scheme
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The QM/MM forces calculated as just explained, using a buffered QM region, can
be directly used to perform an MD simulation, provided that the total sum of the
forces is constrained to be zero, as explained above. However, in many cases,
and in particular in this Si crack case, we can use yet another "trick" to speed
up our MD simulation, by at least a factor 10.
We first introduce a new Hamiltonian, which is the sum of the MM potential (the
SW potential here) and of a so called *adjustable potential*. This potential has
a general analytical expression and a set of parameters that can be adjusted. In
this case, we will use a simple two-body linear 'spring' potential with the form
.. math::
V_\mathrm{adj} = \sum_{\mathrm{pairs}\; ij} \alpha_{ij} r_{ij}
where :math:`i` and :math:`j` indicate a pair of atoms which are closer than an
arbitrary spring length (typically around 3 neighbour hops),
and :math:`\{\alpha_{ij}\}` are the set of adjustable parameters. We then define
a *fitting region*, typically slightly larger than the buffered QM region.
Our adjustable potential is then used to fit the difference between the QM and
the MM forces for each atom in this fitting region. (The functional form above is
particularly suitable as it can be optimised using linear algebra). Of course,
this difference will be zero for atoms in the fitting region that do not belong
to the QM region. In this way, after a successful force matching, the gradient
of our sum potential :math:`-\nabla (V_{SW}(\mathbf{r})+
V_{adj}(\mathbf{r},\{\alpha\}))` gives us the MM forces on the MM atoms and the QM
forces on the QM atoms. Of course a very small error on these forces
(typically :math:`\sim 10^{-3}` eV/A) is caused by requiring our forces to be
the gradient of a single Hamiltonian. This is however of the same order as the
error introduced by forcing the total sum of the forces to zero, and less than
the error associated with computing QM forces using clusters with a finite
buffer radius (which can be as much as :math:`\sim 0.1` eV/A compared to a
reference QM calculation in the full system with PBC). The figure below
illustrates the force fitting process.
.. image:: lotf-overview.png
:width: 500
:align: center
The definition of this new adjustable potential is very beneficial to
our simulation because it allows us to employ a predictor-corrector
algorithm that, as mentioned before, significantly speeds up our
simulation. This scheme is based on the possibility of varying the
values of the adjustable potential parameters :math:`\{\alpha_{ij}\}`
in both time and in space (i.e. each two-body term of the potential
has an independent parameter :math:`\alpha_{ij}`). Once these
parameters are adjusted to match the QM/MM forces in the fitting
region, we do not necessarily need to perform this fitting procedure
at the next time step. In many cases, in fact, the parameters would
change very slowly with respect to the changing atomic positions. This
means that the same parameters can be used for a small number of steps
(e.g. 10 in our case), after which the expensive QM forces are
recomputed, and then used to retune the parameters. This is the
philosophy behind the 'Learn on the Fly' approach.
The predictor-corrector scheme works as follows, with reference to the illustration below:
1. **Initialisation.** The QM region is selected, the buffered QM/MM forces are
calculated and the parameters of the adjustable potential are adjusted to
reproduce the difference between QM and MM forces in the fitting region.
2. **MD predictor: Extrapolation.** The
classical :math:`V_\mathrm{SW}+V_\mathrm{adj}` is used
with the fixed set of adjusted parameters :math:`\alpha_0` to evolve the
dynamics of the whole system starting from the positions :math:`R_0` for a small
number :math:`N` of steps of size :math:`\Delta t` for a total time
of :math:`\Delta T = N\Delta t`.
3. **QM/MM forces.** The buffered QM region method is used to compute the QM/MM
forces on the new configuration :math:`R_1'`.
4. **Force matching.** The parameters of the adjustable potential are tuned to
reproduce the difference between the new QM forces and the new MM forces to give a
new set of adjustable potential parameters :math:`\alpha_1`.
5. **MD corrector: Interpolation.** The formerly predicted MD steps are now
corrected by returning the system to :math:`R_0`, and
re-running the dynamics with the adjustable potential
parameters linear interpolated between :math:`\alpha_0` and :math:`\alpha_1`.
We then return to step 2. and iterate.
.. image:: lotf-pred-corr.png
:width: 500
:align: center
The number of extrapolation steps that can be made before the potential
parameters change too much can be checked by keeping an eye on the force errors
made by the LOTF scheme in comparison with doing a QM force evaluation at each
time step: there is an :ref:`extension task <pred_corr_error>` at the end of
this tutorial to do exactly that.
3rdparty/quip/doc/Tutorials/adaptive-qmmm.rst deleted 100644 → 0
View file @ 23bcf97f
.. currentmodule:: quippy
.. _qmmm_tutorial:
****************************************
Adaptive QM/MM MD of Fracture in Silicon
****************************************
This tutorial has been prepared for use at a hands-on session at the
`ADGLASS Winter School on Advanced Molecular Dynamics Simulations
<http://www.adglass.eu/adglass_news.html>`_, ICTP, Trieste, February
2013.
:Authors: James Kermode and Gianpietro Moras
:Date: February 2013
.. toctree::
:maxdepth: 2
adaptive-qmmm-step0.rst
adaptive-qmmm-theory.rst
adaptive-qmmm-step1.rst
adaptive-qmmm-step2.rst
adaptive-qmmm-step3.rst
adaptive-qmmm-solutions.rst
adaptive-qmmm-references.rst
3rdparty/quip/doc/Tutorials/adglass.rst deleted 100644 → 0
View file @ 23bcf97f
.. currentmodule:: quippy
Molecular Dynamics Simulation of Fracture in Quartz
***************************************************
This tutorial was prepared for use at a hands-on session at the
`Advanced Oxide Interfaces Workshop
<http://cdsagenda5.ictp.it/full_display.php?ida=a10139>`_, ICTP,
Trieste, May 2011.
Setting up your working environment
-----------------------------------
In this tutorial we are going to be using the `QUIP
<http://www.libatoms.org>`_ (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 <http://www.jrkermode.co.uk/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
:ref:`installation` section if you need some help on doing this.
.. To get started, you should set up your working environment with the
.. following shell command::
.. $ source /afs/ictp/public/shared/adglass/jkermode/start.sh
.. This sets up a few environment variables and gives you access to the
.. `QUIP` and `quippy` codes. The `$ADGLASS` environment variable should
.. now point to a directory containing the input files you'll need.
.. You will need to repeat this command in every new terminal you open
.. (or add it to the end of your `~/.cshrc` startup file).
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 :math:`\alpha`-quartz opening
on the basal plane :math:`(0001)`, with the crack front parallel to
the :math:`[\bar{1}10]` direction. The input structure is an
:download:`XYZ file <quartz_crack.xyz>` contains a :math:`150~\AA
\times 50 \AA` crack system with around 4000 atoms. The `y` axis is
aligned with the :math:`(0001)` plane, and the `z` axis with the
:math:`[\bar{1}10]` direction. The system contains a seed crack under
uniaxial tension in the `y` direction with a strain of 23%.
.. image:: quartz_seed_crack.png
:width: 600
:align: center
Download the following input files:
1. :download:`quartz_crack.xyz` - the input structure, pictured above
2. :download:`quartz_crack.xml` - parameters for potential and molecular dynamics
3. :download:`quartz_crack_bulk.xyz` - 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
$ cd $QUIP_ROOT
$ make Programs/crack
Then press `enter` to accept default settings.
.. Note::
The value you should assign to the variable :envvar:`QUIP_ARCH` depends on your
platform (see :ref:`installation` for more details). :envvar:`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 :envvar:`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 :envvar:`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 `quartz_crack_movie_1.xyz`.
You can monitor the progress of the simulation by opening the
trajectory file with my `modified version
<http://www.jrkermode.co.uk/AtomEye>`_ of `AtomEye`::
$ A quartz_crack_movie_1.xyz
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
<http://mt.seas.upenn.edu/Archive/Graphics/A>`_ or my notes on this
`modified version <http://www.jrkermode.co.uk/AtomEye>`_. In
particular, you might find this :download:`startup file <atomeyerc>`
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
<http://en.wikipedia.org/wiki/Moving_average#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).
.. image:: quartz_running_crack.png
:width: 600
:align: center
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 :math:`\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)
aq.write("quartz.xyz")
Here `a` and `c` are the lattice parameters (in :math:`\AA`), `u` is
the internal displacement of the silicon atoms and `x`, `y`, and `z`
the internal coordinates of the oxygens (see image below).
.. image:: quartz_unit_cell.png
:width: 300
:align: center
The variable `aq` is referred to as :class:`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 :attr:`~Atoms.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 :func:`alpha_quartz()` does. If
you are interested in learning more about using `QUIP` and `quippy`,
you could take a look at the online :ref:`tutorial-intro`.
Keep your `ipython` session open as we'll need it later. Open the
`quartz.xyz` file you created with `AtomEye`. You should see something
like this:
.. image:: quartz.png
:width: 300
:align: center
You can also download :download:`quartz.xyz` for comparison.
.. note::
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 :math:`E`, the Poisson
ratio :math:`\nu`, surface energy density :math:`\gamma` and the
Rayleigh Wave speed :math:`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
<http://jcp.aip.org/resource/1/jcpsa6/v133/i9/p094102_s1>`_ in
*J. Chem. Phys.* for more about this potential).
Apart from :math:`\gamma`, all of these quantities can be obtained
from the :math:`6 \times 6` matrix of elastic constants
:math:`C_{ij}`, so we will start by calculating this.
:math:`C` is defined by :math:`\bm\sigma = C\bm\epsilon` where :math:`\bm\sigma` and
:math:`\bm\epsilon` are six component stress and strain vectors following
the Voigt convention:
.. math::
\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)
The simplest way to calculate :math:`C` with `QUIP` is to use the
command line QUIP `eval` program. You will need the file `quartz.xyz`
you made earlier, as well as an :download:`XML file <TS_params.xml>`
containing the parameters of the classical potential::
$ eval init_args="IP TS" at_file=quartz.xyz 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 :math:`\alpha`-quartz, there
are only actually 6 independent values in the :math:`C` matrix:
.. math::
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)
where :math:`C_{66}` is given by :math:`\frac{1}{2}(C_{11} -
C_{12})`. We can exploit this symmetry to get all these values using
only two strain patterns: :math:`\epsilon_{xx}` and
:math:`\epsilon_{yy}+\epsilon_{zz}`.
If you like, you could try using the following `quippy` code to
evaluate :math:`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 :math:`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 :math:`(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 :math:`C` using the rotation matrix, and calculate
the compliance matrix :math:`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 :math:`\sigma_{ij} = c_{ijkl} \epsilon_{kl}`,
meaning that the transformed tensor is
.. math::
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)
aqr.write("quartz_rotated.xyz")
Open `quartz_rotated.xyz` in AtomEye and confirm that it is oriented
so that the :math:`(0001)` surface is aligned with the vertical (`y`)
axis like the fracture system or `quartz_0001.xyz`, 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 :download:`quartz_0001.xyz` contains a 54 atom unit cell for
the (0001) surface of :math:`\alpha`-quartz (shown below). Use either
the QUIP `eval` program or `quippy` to calculate the surface energy
density :math:`\gamma` predicted by our classical potential for this
surface.
.. image:: quartz_0001.png
:align: center
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 :class:`Potential` object using the
XML parameters from the file `TS_params.xml` as shown above. You can
then calculate the energy using the :meth:`~Potential.calc()`
function - see the tutorial section on :ref:`moleculardynamics` for
more details.
.. note::
The energy unit of QUIP are electron volt (eV), and the distance
units are Angstrom (:math:`\AA`).
**Hint**: you can use the expression
.. math::
\gamma = \frac{E_{surf} - E_{bulk}}{2A}
where :math:`E_{surf}` and :math:`E_{bulk}` are the energies of
surface and bulk configurations containing the same number of atoms,
and :math:`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 :math:`\gamma` of around 3.5 J/m\ :sup:`2`.
**Extension**: what effect does relaxing the atomic positions have on
the surface energy? (use the `relax` argument to the `eval`
program, or the :meth:`~Potential.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 \
atoms_in_file=quartz_0001.xyz 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
.. math::
G = - \frac{\partial U_E}{\partial c}
where :math:`U_E` is the total strain energy and :math:`c` is the crack length.
The well-known Griffith criteria uses an energy-balance argument to
equate the critical value of :math:`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 :math:`G > 2\gamma`, where :math:`\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.
.. image:: thin-strip.png
:align: center
:width: 600
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 :math:`\delta`, so the applied strain is
:math:`\epsilon_0 = \delta / h`. Far ahead of the crack, the strip is
in uniaxial tension: :math:`\epsilon_{yy} \to \epsilon_0` as :math:`x
\to \infty`.
The stress far ahead of the crack is given by :math:`\sigma_{0} = E'
\epsilon_{0}`, and therefore the elastic energy per unit length and
per unit thickness far ahead of the crack tip is
.. math::
W_\infty = \frac{1}{2}E'\epsilon_0^2\cdot 2h = \frac{\delta^2 E'}{h}
where :math:`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
:math:`W_\infty` is effectively replaced by a strip with energy
density zero.
The energy supplied to the crack tip is therefore equal to :math:`W_\infty`,
so the energy release rate is simply
.. math::
G = W_\infty = \frac{\delta^2 E'}{h}
In our simulations we will use periodic boundary conditions in the :math:`z`
direction, so we have plane strain loading (:math:`\epsilon_{zz} = 0`),
which means that the effective Young's modulus :math:`E'` is given by
:math:`E/(1-\nu^2)`, where :math:`E` is the Young's modulus in the :math:`y` relevant
direction and :math:`\nu` is the Poisson ratio, so finally we have
.. math::
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 :math:`G` (in units of J/m\ :sup:`2`)
and strain :math:`\epsilon` at which a sample of :math:`\alpha`-quartz
with a height of 50 :math:`\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 `quartz_crack_check.xyz` 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('quartz_crack_check.xyz')
params = CrackParams('quartz_crack.xml')
b = crack_rescale_homogeneous_xy(a, params, new_strain)
b.write('quartz_crack_rescaled.xyz')
Replace `new_strain` with the target strain which should be between
0.15 and 0.30. If you inspect the new file `quartz_crack_rescaled.xyz`
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 quartz_crack_bulk.xyz quartz_crack_rescaled_bulk.xyz
$ 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 :download:`crack-velo.sh` 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 :math:`G` and a universal velocity-dependent
function which he showed can be approximated as a linear
function of the crack speed :math:`V`
.. math::
2\gamma \sim G \left( 1 - \frac{V}{c_R} \right)
The Rayleigh surface wave speed :math:`c_R` sets the ultimate limit for the
crack propagation speed. The expected velocity as a function of the
loading :math:`G/2\gamma` is then
.. math::
\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 :math:`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 :math:`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 :math:`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 :math:`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
-------------------------
* :download:`quartz_crack.xyz`, :download:`quartz_crack.xml`, :download:`quartz_crack_bulk.xyz` - input files for crack simulation
* :download:`quartz.xyz` - quartz primitive cell
* :download:`quartz_0001.xyz` - quartz (0001) surface cell
* :download:`TS_params.xml` - potential parameters
* :download:`elastic.sh` - script to compute elastic constants
* :download:`cij.dat` - elastic constants obtained with `eval` program
* :download:`surface-energy.sh` - script to compute surface energy
* :download:`surface-energy-relaxed.sh` - script to compute relaxed surface energy
* :download:`surface-energy-annealed.sh` - script to compute annealed surface energy
* :download:`crack-velo.sh` - script to crack velocity
3rdparty/quip/doc/Tutorials/alphaquartz.png deleted 100644 → 0
View file @ 23bcf97f
3rdparty/quip/doc/Tutorials/alphaquartz.png

74.6 KiB

3rdparty/quip/doc/Tutorials/alphaquartz2.png deleted 100644 → 0
View file @ 23bcf97f
3rdparty/quip/doc/Tutorials/alphaquartz2.png

45.8 KiB