simulation methods for the calculation of surface

Document technical information

Format pdf
Size 8.1 MB
First found Nov 13, 2015

Document content analisys

Language
English
Type
not defined
Concepts
no text concepts found

Persons

Dennis K. Burke
Dennis K. Burke

wikipedia, lookup

James C. F. Huang
James C. F. Huang

wikipedia, lookup

William H. Fox
William H. Fox

wikipedia, lookup

Ernst A. Lehmann
Ernst A. Lehmann

wikipedia, lookup

Margaret B. Laird
Margaret B. Laird

wikipedia, lookup

Sol H. Weiss
Sol H. Weiss

wikipedia, lookup

Norman D'Amours
Norman D'Amours

wikipedia, lookup

Wang E
Wang E

wikipedia, lookup

Lee C. Teng
Lee C. Teng

wikipedia, lookup

Organizations

Places

Transcript

1
N e w st a t is t ic a l m e c h a n ic a l
SIMULATION METHODS FOR THE
CALCULATION OF SURFACE
PROPERTIES
Thesis for PhD. - Hannah Fox
First supervisor: Mike Gillan
Second Supervisor: Andrew Horsfield
D epartm ent of Physics and Astronom y
UCL
[UONDOff
UMI Number: U591209
All rights reserved
INFORMATION TO ALL USERS
The quality of this reproduction is dependent upon the quality of the copy submitted.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, th ese will be noted. Also, if material had to be removed,
a note will indicate the deletion.
Dissertation Publishing
UMI U591209
Published by ProQuest LLC 2013. Copyright in the Dissertation held by the Author.
Microform Edition © ProQuest LLC.
All rights reserved. This work is protected against
unauthorized copying under Title 17, United States Code.
ProQuest LLC
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, Ml 48106-1346
2
I, Hannah Fox, confirm th a t the work presented in this thesis is my own.
W here information has been derived from other sources, I confirm th at this has
been indicated in the thesis.
Abstract
I present two new methods for the calculation of surface properties. Firstly,
a m ethod of thermodynamic integration to calculate surface free energies. A
strain is applied to a unit cell of the bulk m aterial, th at opens up a vacuum
gap and creates two surfaces. A param eter s describes this process, from 5 = 0
(the bulk material) to s = si (large vacuum gap). The difference in free energy
between these two systems is then calculated by the integration of the stress on
the unit cell over s. I use this general theory to find the surface free energy of
the titanium dioxide ( 110 ) surface using density functional theory.
The second p art of the thesis gives a general transition state theory method
for the calculation of the desorption rate of a molecule from a surface, at any
coverage and tem perature. This approach depends on the density of molecules
as a function of the distance from the surface, and I show th a t this can be found
from the potential of mean force. This is especially useful a t low tem peratures,
where experiments are conducted but brute force simulation is computationally
unfeasible. I use this theory to calculate the desorption rate of water from
the (0 0 1 ) surface of magnesium oxide at 100 — 1200K and 0 —2/3 coverage,
with classical potentials. An im portant outcome of these calculations is that
the frequency prefactor (from the Polanyi-Wigner equation) is dependent on
tem perature.
3
Contents
1
In trod uction
16
2
M odeling o f m aterials
21
2.1
Density Functional Theory
...............................................................
22
2.1.1
The variational principle.........................................................
24
2 .1.2
The energy fu n c tio n a l............................................................
25
2.1.3
The Euler eq u atio n ..................................................................
27
2.1.4
The Kohn-Sham e q u a tio n ......................................................
28
2.1.5
Self-consistency.........................................................................
29
2.1.6
Accounting for our electron-electron in te r a c tio n s .............
30
2.1.6.1
The local density approximation (L D A ).............
31
2.1.6.2
Generalised gradient approximations (GGAs) . .
31
Practical application of D F T ................................................
32
2.1.7
2.1.7.1
Periodic boundary conditions and the plane wave
basis set
..................................................................
32
Pseudopotentials.....................................................
34
Interionic p o te n tia ls ...........................................................................
38
2 .2.1
Coulomb c o n trib u tio n .............................................................
39
2.2.2
Overlap repulsion......................................................................
39
2.1.7.2
2.2
4
CO NTENTS
5
2.2.3 Dispersion f o rc e s ......................................................................
40
2.2.4 Form for the rigid ion p o te n tia l.............................................
41
2.2.5 Fitting the p aram eters.............................................................
41
2.2.5.1
Empirical parametrisation
...................................
42
2.2.5.2
Ab initio param etrisatio n .......................................
42
2.2.5.3
Successes and failures of rigid ion potentials
. .
43
2.2.6 The shell m odel.........................................................................
43
3 F in ite tem perature m odeling
3.1
3.2
Molecular dynamics
...........................................................................
47
3.1.1 Simulation p a ra m ete rs.............................................................
47
3.1.2 In itia lis a tio n ............................................................................
49
3.1.2.1
Initial p o s it i o n s ......................................................
49
3.1.2.2
Initial v e lo c itie s ......................................................
50
3.1.3 Calculation of the f o r c e ..........................................................
51
3.1.4 Integrating the equations of motion
....................................
52
3.1.4.1
Verlet a lg o rith m ......................................................
53
3.1.4.2
Velocity Verlet a lg o rith m ......................................
53
3.1.4.3
Predictor-corrector alg o rith m ................................
54
3.1.5 E q u ilib ra tio n ............................................................................
55
Statistical m echanics...........................................................................
56
3.2.1 Microcanonical ensem ble..........................................................
56
3.2.2 Canonical ensemble
................................................................
57
.............................................................................
58
3.2.3.1
Nose th e r m o s ta t......................................................
58
3.2.3.2
Berendsen th e rm o s ta t.............................................
60
3.2.3 Thermostats
3.2.4 Ergodicity
3.3
45
................................................................................
61
Thermodynamic in te g ratio n ..............................................................
62
CO NTENTS
4
6
S urface free e n e rg y o f T i 0 2 (110)
66
4.1
Introduction..........................................................................................
68
4.1.1
Free energy ca lc u la tio n s.........................................................
68
4.1.2 Titanium dioxide and its (110)su rfa c e ..................................
70
4.1.2.1
4.2
4.3
4.4
B u lk .............................................................................
71
4.1.2.2 The (110) s u r f a c e ....................................................
72
4.1.2.3 Adsorbed species on TiC>2 ( 110 ) ..........................
76
Method of in v e stig a tio n ....................................................................
80
4.2.1 P seudopotentials......................................................................
80
4.2.2 Exchange-correlation functional.............................................
80
4.2.3 k-point m e s h ............................................................................
81
4.2.4 The simulation c e ll ...................................................................
82
4.2.5 Technical d e t a i l s ......................................................................
83
Ti (>2 and the structure of the ( 110 ) s u r f a c e ..................................
84
4.3.1
Bulk calculations......................................................................
84
4.3.2 (110) surface s tru c tu re .............................................................
85
4.3.3 Convergence of surface properties
.......................................
87
4.3.3.1
...with respect to vacuum gap L ...........................
87
4.3.3.2
...with respect to slab thickness n
89
4.3.3.3
Comparison with isostructural tin dioxide . . . .
.......................
91
4.3.4 Lattice dynamics of bulk TiC> 2 .............................................
93
4.3.5
97
C o n c lu s io n ...............................................................................
General theory for the calculation of surface free energies
4.4.1
....
99
Thermodynamic Integration...................................................
99
4.4.1.1
Zero T e m p e ra tu re ....................................................
99
4.4.1.2
Control interaction...................................................... 102
4.4.1.3
Finite T e m p e ra tu re ................................................... 103
CO NTENTS
7
4.4.2 Temperature Integration............................................................ 103
4.5 Application to the TiC>2 ( 110 ) s u r f a c e ................................................106
4.6 R esu lts....................................................................................................... 109
4.6.1 Finding an appropriate form for A(s)
................................... 109
4.6.2 Average stre s s ............................................................................... I l l
4.6.3 Calculation of FSWf ...................................................................113
4.6.4 Temperature integration m e t h o d .............................................114
4.7 Summary and conclusions......................................................................117
5 D esorp tion o f water from MgO(OOl)
119
5.1 Introduction.............................................................................................. 121
5.1.1 Adsorption of water on M gO (O O l)..........................................121
5.1.2 Desorption of water and other molecules from surfaces . . 122
5.1.2.1
Transition state t h e o r y .............................................123
5.1.2.2
Temperature programmed d e s o rp tio n ....................124
5.1.2.3
Desorption of water from M gO(OO l).......................126
5.2 Method of in v e stig a tio n .........................................................................127
5.2.1 Technical d e t a i l s ......................................................................... 128
5.2.2 The potential m o d e l...................................................................129
5.2.3 The unit cell of the s y s t e m ...................................................... 131
5.3 Theory of d eso rp tio n ...............................................................................133
5.3.1 Definition of the s y s te m ............................................................ 133
5.3.2 Formula for the desorption rate 7 ............................................. 134
5.3.3 Potential of mean f o r c e ............................................................ 136
5.4 Zero temperature adsorption e n e r g y ...................................................139
5.5 Equilibration of molecules on the surface.............................................141
5.5.1 Motion of molecules on the s u rfa c e ..........................................141
5.5.2 Surface diffusion coefficien t...................................................... 142
CO NTENTS
8
5.5.3 Orientation correlation
5.5.4 Randomness
5.6
............................................................ 145
...............................................................................145
The isolated m o le c u le ........................................................................... 147
5.6.1 Direct calculation of the desorption r a t e ................................147
5.6.2 Sticking coefficient......................................................................149
5.6.3 Potential of mean force m e th o d ................................................150
5.7
5.6.3.1
Average force (Fz)z ..................................................150
5.6.3.2
Potential of mean force<f>(z).................................... 151
5.6.3.3
y{z) and the desorption rate 7 ............................... 154
5.6.3.4
The frequency prefactor / ......................................... 154
5.6.3.5
Origin of the non-linearity of <j>mi n { T ) ................... 156
Higher C overages.................................................................................... 158
5.7.1 Sticking coefficient......................................................................158
5.7.2 Direct calculation of the desorption r a t e ................................159
5.7.3 Potential of mean force m e th o d ................................................162
5.7.3.1
Average force (Fz)z ...................................................163
5 .7.3.2
Potential of mean force <f>(z)...................................165
5.7.3.3
y(z) and the desorption rate 7 ................................166
5.7.3.4
Critical tem p e ra tu re...................................................167
5.8 Summary and conclusions..................................................................... 173
6
D iscussion
6.1
176
Surface free energy calcu latio n s............................................................ 176
6.2 Desorption calculations
.........................................................................178
A E w ald su m m a tio n
182
B E stim atin g errors using block averages
185
List of Tables
4.1 Lattice parameters and energy of the TiCh primitive unit cell
using different pseudopotentials and functionals. Experimental
values from Diebold’s review [18]........................................................
4.2 The surface energy and structure for every system studied.
...
85
88
4.3 Surface atom displacements after relaxation, compared with pre­
vious DFT studies and experiment.....................................................
89
4.4 Surface energy computed by different studies....................................
91
4.5 Frequencies found for the modes A \g, B 2g, and A 2u using different
methods...................................................................................................
4.6 Surface free energy at the three temperatures investigated
96
114
5.1 Parameters for the rigid ion potential describing interactions be­
tween: (a) Mg and O ions in the bulk or slab; (b) water molecules
and water/slab interactions.....................................................................130
5.2 Mean time between site hops Thop at 200 —500K and four coverages 143
5.3 Decay time rrot at 200 —500K and four coverages.............................. 146
9
L IST OF TABLES
10
5.4 The desorption rate 7 for the isolated molecule on the surface;
second column, as found from the number of crossings observed
during the simulation; third column, found from y{z)\ fourth col­
umn, the error between the two approaches........................................ 148
5.5 The desorption rate 7 for a molecule on the surface at various
tem peratures..............................................................................................150
5.6 Sticking coefficient S for three tem peratures at different cover­
ages. 5 corresponds to the probability th at a molecule approach­
ing a surface will not bounce off.............................................................159
5.7 Actual average number of molecules adsorbed on each surface at
each temperature, for a certain number of molecules in the system. 162
5.8 (a) The desorption rate (in ps-1 ) and (b) the frequency prefactor
/ (in s-1 ) at all tem peratures and coverages studied.........................166
List of Figures
3.1 Flowchart showing the basic operation of a program to perform
molecular dynamics................................................................................
48
4.1 The primitive unit cell of bulk titanium dioxide...............................
73
4.2 The (110) surface of titanium dioxide.................................................
73
4.3 Simulation cell used in molecular dynamics simulations..................
82
4.4 Surface energy U8urf(n = 4, L) with increasing vacuum gap L.
.
87
b) the PBE functional............................................................................
90
4.5 Relaxations of surface atoms with varying n, for a) the LDA, and
4.6 Surface energy with increasing n for different functionals....................92
4.7 Comparison of surface displacements for titanium dioxide and tin
dioxide, using PAW and PB E...............................................................
92
4.8 The A \g, B 2g and A 2u vibrational modes of bulk titanium dioxide. 94
4.9 Relative total energy vs. do for the A 2u mode, using three dif­
ferent functionals....................................................................................
96
4.10 Diagram showing the strain placed on the system as a function
of parameter s ............................................................................................ 100
4.11 Two slabs of titanium dioxide separated by some vacuum gap. . . 107
4.12 Form for control interaction A(s)............................................................ 110
11
12
L IST OF FIGURES
4.13 The fluctuating stress during a simulation and its running average.Ill
4.14 {atotal) vs. s at each tem perature......................................
112
4.15 Projection of the surface free energy up to 2000K using temper­
ature integration, and assuming the harmonic approximation.
. 116
5.1 a) Crystal structure of MgO, and b) the (0 0 1 ) surface................. 121
5.2 How the system looks, using periodic boundary conditions.
. . . 127
5.3
An example of the form of y (z).............................................................136
5.4
Configurations in which the water molecule adsorbs to the sur­
face: a) the energetically favourable flat configuration (FC); b)
the perpendicular configuration (PC); c) the bridging configura­
tion ( B C ) ................................................................................................ 139
5.5
The trajectory of a single molecule over the surface for lOOpsat
(a) 500K and (b) 200K............................................................................ 141
5.6 The trajectories for (a) sixth ML and (b) third ML for 15ps. . . 142
5.7 For the half coverage system at 400K: (a) mean square displace­
ment with time; (b) projection P (t) with tim e
5.8
144
The distribution y (z) at T — 800, 1000 and 1200K, calculated as
a histogram. The surfaces are a t z = 0 and 25.8A............................. 147
5.9
The z -coordinate of a molecule inbet ween the two surfaces, at
1200K. Bottom of slab is at z = 0..........................................................149
5.10 (a) the average force in the z-direction on the water oxygen as
a function of temperature, from 1200K down to 100K. (b) the
potentials of mean force <p(z) in the same temperature range. . . 152
5.11 Minimum of (f>{z)..................................................................................... 153
5.12 Detail of y(z) at all coverages at a tem perature of 1200K............... 160
LIST OF FIGURES
13
5.13 Desorption rate 7 found for 1—18 molecule systems: (a) assuming
a sticking coefficient 5 = 1 , and (b) using the values of 5 reported
in table 5.6.................................................................................................161
5.14 (a) the average force on the constrained molecule at 800K, for
zero, sixth, third and half coverage, and (b) the equivalent po­
tentials of mean force
164
5.15 Minimum of <f>(z) asa function of coverage a t 300, 400 and 800K. 165
5.16 Desorption rate at 800K, found by direct calculation and the
PM F method.............................................................................................168
5.17 Pressure vs. coverage a t 300, 400 and 800K.......................................168
5.18 Representation of the curve of po vs. a at 300K............................... 170
5.19 Density vs. coverage in the tem perature range 350-385K, in in­
tervals of 5K
171
List of Publications
1. H Fox, AP Horsfield and MJ Gillan: Density functional calculations of
surface free energies, J. Chem. P h y s 124:134709 (2006)
2 . H Fox, MJ Gillan and AP Horsfield: Methods for calculating the
desorption rate of an isolated molecule from a surface: Water on
MgO(OOl), Surf. Sci., 601:5061 (2007)
3. H Fox, MJ Gillan and AP Horsfield: Calculations of the desorption rate
of Water from MgO(OOl) at non-zero coverages, to be published 2008
14
Acknowledgments
First thanks must go to my two supervisors Mike Gillan and Andrew Horsfield.
Both have been understanding and patient during my (hopefully successful)
transformation from theory graduate to scientific researcher. Also, despite a
couple of physical relocations and a collegial reassignation, they have taken me
through the PhD. without any of the m ajor non-scientific dramas of some of my
fellow students, for which I am grateful.
Thanks also to other members of the Condensed M atter and Material Physics
group, for their help understanding and reasoning with the various computer
codes and systems I have used for my work. I have immensely enjoyed my three
plus years at UCL and more generally in London. Thanks to all the people in
the university th at have helped make th at happen.
15
Chapter 1
Introduction
JW Gibbs said th at “The whole is much simpler than the sum of its parts”, a
quote of wisdom in good company with “the simplest explanation is probably
the correct one”, and “to find the murderer, look for the one who stands to
benefit the most”. Many scientists to some extent dream of a grand unifying
theory, which will be simple, yet capable of explaining all the phenomena in the
universe. However, despite the great strides science has taken in the last century,
the distance to this goal cannot be ascertained, and indeed it is probably beyond
our reach. The problem is th at answers beget more questions, just like a child’s
innocent inquiry starts a neverending series of ‘why’s.
The proof of the substructure of the atom, the development of the theory of
quantum mechanics, these things were prompted by the results of experiments.
They answered the questions posed by experiments using new laboratory appa­
ratus, and in turn predicted more properties th at could not yet be experimen­
tally accessed. The two different strands of science compliment each other and
make progress together, answering the new question of ‘why’. They may be out
of step at any given moment, but time sees them resolve to the same conclusion.
16
CHAPTER 1. INTRODUCTION
17
The revolutionary new technology of the twentieth century, the computer,
introduced a third sibling to the scientific family. Its physicality suggests its use
as a piece of laboratory equipment, yet its logical capabilities suggest its suit­
ability to theory. The computer has forged a third way, which uses theoretical
ideas about how things work to perform a virtual experiment. As a result, we
can perform ‘experiments’ in regions which are not experimentally accessible,
or feasible.
In this thesis, I use computer simulation to investigate the properties of
a couple of systems, using two different approaches to the theory of molecu­
lar interaction. Comparison of these properties with experiment was not always
possible, but where it was, there is usually some quantitative discrepancy. How­
ever, parts of my work point towards ideas that experiment has not yet touched
on, just as I have used the accumulated experimental experience of the systems
of interest to further my work.
The free energy of a system is a very im portant quantity. Tied up in this ther­
modynamic potential is one of the keys to understanding the dynamic behaviour
of physical m atter at finite temperature. This is because of its dependence on
internal energy and entropy, with the balance between these two quantities de­
termining for example whether a chemical compound is solid, liquid or gas at
a given temperature. In my work, I have investigated two properties of sys­
tems involving surfaces, and the free energy comes into both the new methods
I propose for their calculation.
Firstly, I give a method for calculating the surface free energy Fsurf , of any
crystal surface, by a method of thermodynamic integration. The probability
of finding a system in one state or another is determined by the difference in
free energy between these states, so we can determine which surface of a given
C H A P TE R 1. INTRODUCTION
18
crystal is most likely to occur by calculate the free energy of each surface. The
equilibrium shape of a crystal is thus found by minimising the surface free energy
with respect to these possible surfaces. As well as this predictive power, they
can tell us how important surface entropy effects axe at a given temperature.
Despite this, the temperature dependence of Fsurf is generally not known for
most systems. This is because it cannot be directly calculated (unlike the surface
energy), or calculated by some average of the known properties of the system
(like the stress). Instead, some trick is needed, such as the method I use.
Secondly, I give a method for calculating the desorption rate of molecules
from surfaces. Experiments of desorption (of water molecules) are conducted at
low temperatures o f80—300K, but this region is difficult to simulate directly, be­
cause the rates are very low. An experiment is typically conducted over minutes
or seconds, but computers are currently only able to simulate nanoseconds of
molecular dynamics in reasonable time frames. This type of ‘rare event’ problem
- where interesting processes occur a t slow rates - is common in the simulation
of materials, and a number of tricks have been developed to get around it. In
my method, I use the potential of mean force method to find the desorption
rate. Here, I pull a molecule off the surface slowly, and the potential of mean
force is calculated a t a series of distances from the surface. The PM F is a free
energy function, describing how tightly adsorbed a molecule is to the surface
and the other adsorbed molecules, as a function of tem perature and distance.
Below I give an outline of the thesis, including the theory th at is used in the
simulations.
C h a p te r 2 introduces two different methods for the modeling of materials.
Firstly density functional theory, a method for finding the quantum mechanical
CHAPTER 1. INTRODUCTION
19
ground state of a system. Secondly interionic potentials, which find a simple
analytical form for the interaction between two ions in a system, treating the
ions classically.
C h a p te r 3 discusses the modeling of materials at finite temperature. Molec­
ular dynamics, the method used in my work, simulates the dynamics of a system
by moving all constituent atoms according Newton’s equations of motion. The
application of molecular dynamics requires the consideration of many aspects
of statistical mechanics, which I also discuss. The calculation of free energy
differences by the method of thermodynamic integration is also considered, due
to its importance in my work.
C h a p te r 4 contains my work on the calculation of the surface free energy
of titanium dioxide (110) using density functional theory. I firstly investigate
the particular system, examining the effect of various simulational details upon
the accuracy of its modeling. Then I use this knowledge to find the surface
free energy by a method of thermodynamic integration. Although the surface
free energy is a fundamental quantity, few studies have been performed to find
it, none of which I am aware using first principles modeling. I show that it is
feasible to calculate accurately with density functional theory.
C h a p te r 5 holds my work on the desorption of water from magnesium ox­
ide (001 ) using interionic potentials. The desorption rate was calculated at high
temperatures (through ’brute force’) and at all temperatures using the poten­
tial of mean force method. The methods can be checked against each other at
high temperature, showing the PM F to be accurate. As well as calculating the
desorption rate of a single molecule from the surface, I have investigated higher
CH APTER 1. INTRO D U CTIO N
20
coverages of |M L to f ML, which I don’t believe has been done before for any
system. This is not feasible using density functional theory, but I show th at the
statistical errors are manageable using classical potentials.
C h a p te r 6 contains a summary and discussion of the important results in
the thesis.
Chapter 2
Modeling of materials
The power of modeling, th at is the use of mathematical equations to simulate
and predict real events and processes, has been utilised by all from theoretical
and experimental scientists to financiers. There are very few systems whose
physics and equilibrium properties are known exactly; for example the ideal gas
and harmonic crystal. For all other materials, we make some kind of approxi­
mation, or assumption, about the true physics underlying the properties of the
material, and use this model to predict these properties, which we can then
compare with experiment.
Materials modeling in the computer age is dominated by two different meth­
ods. The first, density functional theory (DFT) which I discuss below in section
2 . 1 , aims to calculate the exact energy of a system by finding its electronic
groundstate. The second is classical modeling which treats the constituent par­
ticles of the material classically. For ionic systems, interionic potentials (section
2 .2 ) aim to model the short and long range forces that are responsible for the
bonding and interaction between ions, accounting for the consequences of quan­
tum mechanics without any explicit quantum calculation.
21
CHAPTER 2. MODELING OF M ATERIALS
2.1
22
D ensity Functional Theory
The large increase in computing power over the past twenty years has led to great
advances in the simulation of solids and liquids. Some of the chief beneficiaries
have been the various ab initio electronic structure methods, which require much
more computing time than their non-quantum counterparts. Of these, density
functional theory (DFT) is the grounding for many. The classic papers of DFT
are those published in the mid sixties by Hohenberg and Kohn [51], and Kohn
and Sham [59], although the earliest form of density functional theory was the
Thomas-Fermi model proposed in the 1920s [98, 24].
The aim of the following summary of DFT is to show how the ground state
energy of a system of interacting electrons in an external potential can be found
without making any approximations, in theory. The external potential is usually
defined as the potential felt by the electrons due to the nuclei of the material.
Of course we wish to be able to simulate a full material, both ions and electrons,
such that in principle we would want to solve the many-body time-independent
Schrodinger equation treating all particles in the system quantum-mechanically.
However this would be completely impossible and unnecessary. We make a
number of decisions about the system to simplify the problem:
• the electrons are non-relativistic.
• there are equal numbers of electrons with positive and negative spin.
• the nuclei are treated as classical particles, and we neglect their motion.
The last decision is justified by the Born-Oppenheimer Approximation. The
mass of an electron is many orders of magnitude smaller than that of an atomic
nucleus, and the velocity of a nucleus is much smaller than that of the electrons.
This means that we can consider the movement of the nucleus independently
CH APTER 2. MODELING OF M ATERIALS
23
of the movement of the electrons: that for a given movement of the nuclei the
electrons will follow them adiabatically.
Given these simplifications, we can now write the many-body time inde­
pendent Schrodinger equation (TISE) for a system of N e electrons at positions
n (i = 1 ...N e) in the potential from N j ions at positions R / (I = l...N j):
(2 .1)
where m is the mass of the electron, and Z i is the charge of ion I. The prime on
the second sum in the second term indicates th a t all sums where i ~ j should
be neglected. $ is the many-body wavefunction describing all the electrons in
the system. This is an eigenvalue equation for the energy E , which is a function
of solely the positions of the nuclei. We are interested only in the ground state
of the electrons, so we ignore all but the lowest energy solution to equation 2 . 1 .
Of course, it is impossible to solve equation 2.1 exactly, mostly due to the
interactions between the electrons described by the second term. If they didn’t
interact, the TISE could be split up into N e equivalent differential equations and
thus the resulting many-body wavefunction
would be an anti-symmetrised
product of single-particle wavefunctions. However the motions of the electrons
are correlated, and the particles exchange information about their motion, so
they cannot be treated independently. If we choose to neglect correlations, we
make the Hartree-Fock approximation, and then the TISE can be solved to find
the Hartree-Fock orbitals. The Hartree-Fock approximation, despite success
with some systems, is not accurate in systems where electron correlation is
significant. In DFT, we can treat electron correlation.
CHAPTER 2. MODELING OF M ATERIALS
2.1.1
24
The variational principle
A crucial building block of DFT is the variational principle. This describes a
method for the calculation of an upper bound on the ground state energy E 0 of
a system whose hamiltonian H is known.
Suppose we have, or can guess at, a wavefunction ip which is normalised
and has a form suitable for the true wavefunction of the system. Then we
can calculate the expectation value of the hamiltonian (ip\H\ip^ — E '. Now,
the lowest possible value of E ' is the ground state energy f?o, which would be
returned if ipwas the true wavefunction of the system. It is not likely we would
be th a t lucky, but E ' is still useful because it gives an upper bound on the
true ground state energy. Furthermore, we can repeat the calculation of E ' for
different trial wavefunctions ip, and the true ground state wavefunction of the
system will be the one th at minimises E '.
This can be proven by expanding the wavefunction ip in terms of the eigen­
states of the hamiltonian:
OO
w
= £ c ,in > ,
(2 -2)
n=0
where the Cn = (n\ip) are a set of complex coefficients subject to the constraint
£ „ \Cn\2 = 1. Then, using the normalisation condition th at (ip\ip) — 1, we find
that:
WHW = £
771
£
C^(m \E„C„\n) = £
n
\C„\2E n,
(2.3)
71
where E n is the energy eigenvalue of state n. The values of {Cn} th at give
the minimum possible value of (ip\H\ip) are C0 = e x p ia and Cn^o = 0, where
exp ia is an arbitrary phase factor. This means that:
(iP\H\iP) > E 0.
(2.4)
CHAPTER 2. MODELING OF M ATERIALS
2.1.2
25
The energy functional
Our hamiltonian for the system of N e interacting electrons in an external po­
tential ti(r) (from the atomic nuclei) can be expressed as:
H = T + U + V.
(2.5)
T is the kinetic energy of the electrons, and U describes all the electron-electron
interactions. V is the contribution to the energy from the external potential
such th at V ~
v b"*)- The ground state energy is given by:
Eo = < *oiH |*o>,
(2-6)
where 'to is the ground state wavefunction, assumed to be non-degenerate. The
electron density n(r) is given by:
n(r) = ( ¥ 0 |n(r)|¥o>;
JVe
n{r) = ^ <5(r - r 4).
i—l
(2.7)
DFT rests on two theorems that relate the external potential u(r), the ground
state energy Eo and the electron density n(r). Central is the idea that Eo can
be expressed as a functional of n(r); th at is for a given electron density one
value of the ground state energy can be calculated. We denote this functional
by Eo [n]. The validity of Eo [n] rests on the first theorem of DFT:
Two different external potentials cannot give the same ground state electron
density distribution.
We can prove this by contradiction by considering two potentials v(r) and
v'(r) which differ by more than a constant, and which give the same electron
density n(r). We then have two hamiltonians H and H ', and two different
CHAPTER 2. MODELING OF M ATERIALS
26
ground state wavefunctions >&o and ^q. The variational principle says that any
wavefunction can be used as a trial wavefunction and will give an energy larger
than the true ground state energy. So if we use
as a trial wavefunction for
H ', we find that:
( $ o | H #j * 0 ) = < * o | £ | * o > + < * o | H ’ - H \* o )
= E0 +
> K
J
d rn (r) [v '(r) - v(r)]
( 2 .8 )
,
where E 0 and E'0 are the corresponding ground state energies . Similarly;
+ ( * i|H - H #|*i>
= K -
J
drn(r) [v'(r) - t>(r)]
(2.9)
> Eo.
If we add equation 2.8 to 2.9, we find that:
Eq
+
Eq
>
Eq
+
E o,
(2 . 10)
which is a contradiction.
We have shown th at the external potential acting on the system of electrons
is uniquely specified by the ground state electron density, so th at it is possible to
express the E q as a functional of the density, namely Eo [n]. The next theorem
states that we can use a form of the variational principle to find the ground
state energy:
The ground state energy E q associated with a given external potential t>(r) is
CHAPTER 2. MODELING OF M ATERIALS
27
found by minimising the energy functional Eo [n] with respect to n(r) with v(r)
held fixed, and the n(r) that yields this minimum is the ground state density
distribution of electrons. n(r) is also subject to the constraint f dr n(r) — N e,
the number of electrons in the system.
The proof of this theorem rests on the fact th at any wavefunction $ uniquely
determines the electron density via equation 2.7. To find the ground state energy
for a given potential, we would need to find the
for which E is at a minimum,
by the variational principle. This is then equivalent to minimising the energy
with respect to the density.
2.1.3
The Euler equation
We now have a possible method for calculating the ground state energy of a sys­
tem of interacting electrons in an external field u(r). However, we need a clearer
definition of what is included in our energy functional Eo [n], and what it means
for it to be at a minimum. From our definition of the hamiltonian in equation
2.5, we can identify a number of different contributions to the functional using
equation 2 .6 :
Eq [n] = F [n] +
J
dr v(r) n(r),
(2 .11 )
where the functional F [n] is the expectation value of T + U:
F[n] = (¥0\T + U\V0)-
(2-12)
We can further split up F[n]:
F [n] = TV/ [n] + G [n]
(2-13)
CHAPTER 2. MODELING OF M ATERIALS
28
T n i [n] is defined to be the kinetic energy of a system of non-interacting elec­
trons in its ground state, whose density is n(r). This is different from just
the expectation value of the kinetic energy operator, but the reason for this
distinction will become clear in the next section when talking about the KohnSham equations. G [n] is defined precisely by equation 2.13, and includes all the
electron-electron interaction contributions to the energy.
We now want to find a condition for describing when E q [n] is at a minimum.
Taking the functional derivative of equation 2.11, we find th at the variation of
E is given by:
(2.14)
We are subject to the constraint th at the number of electrons in the system is
fixed, so J drSn(r) = 0. Using the method of Lagrange undetermined multipliers, we find the condition th at E be a minimum (that is, for SE to vanish) to
be:
(2.15)
where /x is the undetermined multiplier, defined such that we get the correct
number of electrons. Equation 2.15 is the Euler equation for the problem.
2.1.4
The Kohn-Sham equation
Equation 2.15 gives us a condition for finding the minimum of the functional
Eo [n], and thus the ground state energy of the system. However, we have no
practiced way of using it, as we do not know the form of G [n].
If we define an effective potential Vef f { r):
(2.16)
CHAPTER 2. MODELING OF M ATERIALS
29
we can rewrite equation 2.15 as:
(2.17)
This equation has exactly the same form as for a system of non-interacting
electrons in an external field Ve/ / . This is just a reclassification, transferring all
our problems into Vef f ; but assuming we can find some way of approximating
Vef f as a function of position and density, we can reduce the problem to solving
the Schrodinger equation for a system of non-interacting electrons in an external
potential Vef f :
(2.18)
This is called the Kohn-Sham equation, and its first N e/ 2 solutions of lowest
energy
{-0*} are the Kohn-Sham orbitals. They can be used to find the ground
state density distribution:
NJ2
"(r ) = 2
i * i 2’
(2.19)
however this is their only purpose; they are not the real wavefunctions of the
electrons in the interacting system.
2.1.5
Self-consistency
We do not yet have a form for our functional G [n], and therefore our effective
potential Ve//( r ) . However assuming that we can find some approximation to
it, we can proceed with calculations using the energy functional. We also do not
know the density of electrons n(r) in the ground state. Therefore, we must use
some iterative procedure to find the true density from some initial guess. This
procedure for calculating the ground state energy of our system of interacting
CHAPTER 2. MODELING OF M ATERIALS
30
electrons in an external potential runs as follows:
• we make an initial guess of the electron density distribution riin(r);
• we construct G[n] and hence K //( r ) ;
• we solve the Kohn-Sham equation;
• we calculate nout (r) from the Kohn-Sham orbitals.
At this point, we compare our initial density n*n(r) with our calculated density
n ou*(r ), and they will in general not agree. We now create a new density riin(r),
which we hope will be more self-consistent. This can be done for example by
linear mixing, using some proportion a of the initial density and the proportion
(1 —a) from the out density.
This process is repeated until self consistency between riin(r) and nOUf(r)
is achieved, and at this point we have found the n(r) for which the energy
functional Eo [n] is stationary.
2.1.6
Accounting for our electron-electron interactions
We have seen that the functional G [n] is responsible for the part of the total
energy functional treating the interactions between electrons in the system, as
defined in equations 2.12 and 2.13. We can split G\n\ into two contributions.
The first is the electrostatic energy of the electrons, called the Hartree energy:
( 2 .20 )
The second part, which we call the exchange-correlation energy E xc, is defined
only by E xc = G[n] — E h • The Hartree energy includes the unphysical self­
energy for each electron, and in exact DFT we would subtract this contribution
CH APTER 2. MODELING OF M ATERIALS
31
in our form for E xc. The total energy functional is now:
Eq [n] =
2 .1.6.1
J
dr V (r) n(r) + T n i [n] + E H [n] + E xc [n].
(2.21)
T h e local d e n sity a p p ro x im a tio n (LD A )
The local density approximation gives us a basic method for finding E xc. There
is one system for which the exchange and correlation energy can be accurately
calculated - that of a uniform electron gas in a uniform positive background
(known as jellium). The LDA posits that the exchange-correlation energy per
electron of the real system of interacting electrons at a position r where the
density is n(r), is equal to the exchange-correlation energy per electron exc of
jellium of the same density. Mathematically we can write:
(2 .22)
An exact formula for the exchange part of exc was reported by Dirac in 1930
[19]. No formula can be found for the correlation part of exc, so best estimates
are used from quantum Monte-Carlo calculations, where the energy is tabulated
as a function of the density of the electron gas.
The LDA can be derived, but only in materials where the electron density
varies slowly. Of course in practical calculations the density distribution always
changes rapidly with r. Another problem with the LDA is th at it fails to subtract
the unphysical self-energy term in the Hartree energy E j j • Nevertheless, it is
often satisfactory for a large number of basic materials, and is still widely used.
2.1.6.2
G e n e ralise d g ra d ie n t a p p ro x im a tio n s (G G A s)
The various different generalised gradient approximations attem pt to redress the
problems of the LDA by incorporating some consideration of the rate of change
CHAPTER 2. MODELING OF M ATERIALS
32
of electron density. A way of doing this is to replace exc by some function of
both the electron density and its gradient:
(2.23)
where:
(2.24)
There are a wide number of GGAs available, each of which have been shown
to model certain systems well, but do not necessarily work for other systems.
In many cases a particular functional has been seen to model the system less
accurately than the LDA. Thus it is important to test a system at the beginning
with both the LDA and various GGAs, and pick the one which best reproduces
experimental data, or perhaps previous theoretical results.
2.1.7
Practical application of D F T
2.1.7.1
P e rio d ic b o u n d a ry c o n d itio n s a n d th e p la n e wave b asis se t
We are faced with solving a Schrodinger equation for an electron in an effective
potential Vef f (r), to find the N e/2 wavefunctions with the lowest energy. As­
suming we want to simulate a large material containing the order of 1023atoms,
this is still an impossible task. This problem is addressed by using periodic
boundary conditions (PBC), where the material is built from an infinite num­
ber of small, repeating unit cells. This is not only used with DFT, but is rather
a general technique for the simulation of large systems. These unit cells are par­
allelepipeds, described by the vectors a ls a 2 and a 3 - The point r in the cell is
then equivalent to all the points r + R , where R is any vector of the Bravais lat­
tice R = li3ii +123.2 + / 3a 3 , where h are integers. This means that we need only
consider the number of electrons found in each cell, and th at the density n(r)
CHAPTER 2. MODELING OF M ATERIALS
33
and effective potential Vef / (r) also have the periodicity of the Bravais lattice.
Bloch’s theorem tells us that the eigenfunctions ipi(r) of any hamiltonian in
which the potential is periodic can always be expressed as:
^i(r) = Ui(r) exp ( ik .r ) ,
(2.25)
where k is a vector in the reciprocal lattice space. u*(r) is a function with the
periodicity of the Bravais lattice, and can be expanded in terms of plane waves:
=
E Ci5G exp ( iG .r ) ,
(2.26)
G
where the sum is over all reciprocal lattice vectors G defined by G.a* = 27rm
for all integers m and any lattice vector a». Putting this form for u*(r) into
equation 2.25, we have that:
V'tW = ^ 2 c*-g exp
+
•
(2-27)
G
So to find our set of V’t(r) we need to calculate the coefficients c^ g at every
point in k-space. However, a point k + G is equivalent to the point k (just as
in real space, r is equivalent to r + R ), so we need only find the coefficients for
values of k in the first Brillouin Zone (the k-space equivalent of the unit cell).
There are an infinite number of k points in the first Brillouin zone, so we
need a method of picking a finite number of points which will do a good job of
representing the whole cell. In my work I have used the scheme suggested by
Monkhorst and Pack [77]. It is worth noting that the larger the unit cell used,
the smaller the first Brillouin zone is, and therefore the less k points will be
needed. Whilst most primitive unit cells of crystal lattices are quite small and
will need a large number of k points, there are many situations in which it is
CHAPTER 2. MODELING OF M ATERIALS
34
necessary to use a supercell, which can be described by as few as one k point. A
supercell is a number of primitive unit cells, grouped together and used as the
repeating block of the Bravais lattice. They are useful for modeling properties of
systems which require symmetry breaking, and essential for examining surface
or interface behaviour.
Plane waves are not the only basis set we can use to express the wavefunctions ipi. Another class of basis sets that can be used in DFT are local orbitals,
which are centered on atomic nuclei, have an angular dependence and fall off
to zero well within half a nanometer from the nuclei. The most commonly used
set of local orbitals are the local atomic-like orbitals. These kinds of basis sets
have some advantages over plane waves: namely, th a t they already look similar
to the wavefunctions of an electron around a nucleus, and therefore less of them
may be needed to describe the true wavefunction. This will reduce the com­
putational time compared with plane waves. Also, they lend themselves easily
to physical interpretation. However, they bring their own disadvantages too.
Where plane waves are unbiased and treat every area of the unit cell equally,
the local orbitals have a clear bias that may affect their description of the elec­
tronic structure. There is also no systematic way of reducing the error due to
the incompleteness of the basis set. With plane waves, we can simply increase
the number of plane waves included in the expansion in equation 2.27.
2.1.7.2
P s e u d o p o te n tia ls
We are going to represent the wavefunctions i/>j(r) by a plane wave expansion.
This is difficult however, as tpi (r) behaves very differently near the atomic nuclei
from how it behaves inbet ween the nuclei. The electron density varies very
quickly near the nuclei, so lots of plane waves of high energy are needed in the
expansion. We can redress this by using the pseudopotential method.
Principle to this method is to adopt the frozen core approximation, assuming
35
CH APTER 2. MODELING OF M ATERIALS
the core electrons of an atom to be in exactly the same state as for a free atom.
Core electrons generally have much lower energy than valence electrons, and
it is only the valence electrons which are involved in determining the bonding
properties and density of states. Thus we do not need to treat them explic­
itly. Of course the core and valence electrons do interact via the Hartree and
exchange-correlation parts of the potential, so we need to modify our ionic po­
tential V (r) to include this. Another complication is th at the wavefunctions of
the valence electrons must remain orthogonal to those of the core orbitals during
the minimisation process. The energy eigenvalues of the valence electrons and
their nearby energy states must be reproduced by the pseudopotential. This
requirement is the equivalent of saying th at the pseudopotential has the same
phase shift upon valence electron scattering as the true potential.
We call our modified potential - pseudopotential - Vp8(r), and define a cutoff
radius r c. We require th at V^a(r) is equal to Vcore at distances greater than r c,
where Vcore is the potential from the ions and core electrons; and usually that its
first and second derivatives are equal at r c. We also require that Vps be smooth
within r < r c. These conditions mean th at the modified wavefunctions of the
valence electrons (the pseudo-wavefunctions) are equal to the real wavefunctions
at distances greater than r c, but are smooth and slowly varying within the cutoff
radius. This smoothness means that the pseudo-wavefunctions can be much
more easily expressed as plane wave expansions. In general, the cutoff radius
and pseudopotential will be different for valence electrons with different angular
momentum.
N o rm -c o n se rv in g p se u d o p o te n tia ls [41]
One method of ensuring that,
outside the cutoff radius, the true wavefunction if) and pseudo-wavefunction if)
CHAPTER 2. MODELING OF M ATERIALS
36
coincide is to apply the norm-conservation (NC) condition:
(2.28)
i.e. that the charge inside the region r < r c is conserved. This ensures th at the
pseudopotential is transferable to a number of different systems and environ­
ments.
U ltra s o ft p s e u d o p o te n tia ls
Norm-conserving pseudopotentials are success­
ful for many systems but are not suitable for systems containing highly localised
valence orbitals, for example transition metal atoms. This is because the NC
condition makes it impossible to construct a pseudowavefunction th at is much
smoother than the real wavefunction.
Vanderbilt [101] suggested a scheme that removes the NC condition, so that
a pseudowavefunction can be built th at optimises smoothness, and introduces
two or more reference energies which must be satisfied in order to ensure trans­
ferability between systems.
P r o je c to r a u g m e n te d waves
Blochl [8] suggested th at the real wavefunc­
tion and a well-behaved pseudowavefunction could be linked by a linear trans­
formation t , such that:
if)
=
Tlj).
(2.29)
Physical quantities can then be computed from the pseudowavefunction, instead
of being approximated. The tp and tp differ only within an augmentation region
around each atomic core, so the transformation can be broken into a sum:
The wavefunctions can then be represented as a linear sum over partial
waves & and 0*, which can be generated from numerical solutions of the radial
CHAPTER 2. MODELING OF M ATERIALS
37
Schrodinger equation. Associated with each partial wave <f>i is an orthogonal
projection function (p<|, so that the transformation is given by:
r = 1+
(l& > -
(P il-
( 2 *3 0 )
i
We are then free to choose our pi such as to maximise the smoothing of the
pseudowavefunctions and the transferability.
Kresse and Joubert [62] showed th at the ultrasoft pseudopotentials can be
formally derived from the PAW energy functional. They tested both potentials
with a range of solids and molecules, and found th at the PAW method is ex­
ceptionally precise, and th at the US pseudopotentials offer the same level of
precision except in magnetic systems.
CHAPTER 2. MODELING OF M ATERIALS
2.2
38
Interionic potentials
We have seen th a t density functional theory is potentially an extremely pow­
erful tool in the realistic simulation of materials. However, quite apart from
questions over the accuracy of the LDA and GGAs, DFT may not sometimes
be tractable due to its large demand of computer power. Classical potentials on
the other hand are cheap to use, and in situations where the electronic struc­
ture of the material is not directly important, they are capable of reproducing
various physical and chemical properties accurately.
A functional form for a classical potential, considering ionic materials, is
commonly found by consideration of the different short-range forces acting be­
tween ions in the material. The ions are assumed to be rigid ions, i.e. they have
no internal structure, such that the interionic interactions depend only on their
positions, i.e. V = V({ri}) where i = 1...7V, with N ions in the system. We
call V in this case a rigid-ion potential. In general, we can write V as a sum of
terms describing pair interactions, three-body interactions etc.;
V{{ri} ) = ' t v ( n ) + \ i ' t v ( r i,rj ) + 1-j : j : £
i
i
j^ i
i
V f o .r ,,* ) + ...
j^ i
(2.31)
There are N such terms, but the series is truncated quickly for practical rea­
sons. If we restrict V to pair-potentials, the computer time required will be
proportional to N 2, whereas including three-body terms would increase th at to
N 3. For this reason, three-body terms and higher are normally not included.
Although many-body interactions are often important, as the total energy of
the system must be influenced by certain groupings of ions larger than pairings,
this is often a suitable approximation. Instead of a normal pair-potential, we
can use an effective pair potential, which includes the many-body effects in some
averaged way.
CHAPTER 2. MODELING OF M ATERIALS
2.2.1
39
Coulomb contribution
The coulomb force acting between two ions is responsible for the vast majority
of the cohesive energy of an ionic crystal or molecule. We can write the potential
due to the electrostatic interaction between atoms i and j as:
V g ^ ( r) = ^
(2.32)
where Zi, zj are the associated charges on the ions and r = |r^ — |. So the
total contribution to the potential is:
(2.33)
These sums are over all ions in the system. Assuming we are using periodic
boundary conditions, as discussed in section 2.1.7.1, this means a sum over
an infinite lattice. This sum is only conditionally convergent, so we use Ewald
summation to calculate V coul. This breaks V coul into short range and long range
contributions: V coul = V ^ ul + V£ou*. The short range part can be calculated
quickly in real space, and the long range part converges quickly in Fourier space.
2.2.2
Overlap repulsion
At short range, there must be some force th at prevents a positive ion and a
negative ion collapsing into each other, and the origin of this repulsive force is
quantum mechanical in nature. As two atoms approach each other, there will
come a point when their electron clouds start to overlap. This pushes some of
the electrons into higher energy states to keep their wavefunctions orthogonal.
Also, this overlap leaves the nuclei of the ions incompletely screened, so they
repel each other coulombically.
The Born ionic model proposed the following functional form for this repul­
CHAPTER 2. MODELING OF M ATERIALS
40
sion:
Vi,jP(r ) = A iJ exP (~ r /Pi,j) i
(2‘34)
justifying it by the observation th at the wavefunction of a free ion falls off
approximately exponentially at large distances. A i j and pij are parameters
defined for each pair of species in the system.
2.2.3
Dispersion forces
Another contribution to the potential must come from dispersion forces, or Lon­
don forces, first described by London [73]. The origin of these forces is in the
fluctuation of an ion’s electronic structure. If we imagine an ion stationary in
a vacuum, we can make the statement th a t the ion is at rest, but this is only
accurate classically. Quantum mechanically, we know th at the electronic struc­
ture of the ion must be changing, even though on average, it will be spherically
symmetric. Thus, even an ion which is non-polar must have an instantaneous
dipole moment. Two ions which are within a short distance of each other can
induce moments in the other ion, so there will be a correlation between the
charge distributions of each ion. The energy of this correlation does not average
to zero, and indeed London showed th at the contribution to the energy from
these type of interactions was large.
To find an analytical expression for this contribution, one can follow the
method of Drude [20], by considering a model where an ion consists of a positive
charge Q, that remains stationary, and negative charge —Q, which oscillates
about the positive charge in the z-direction. If there are two such ions a and b,
separated by a distance r in the z-direction, then one can write the Schrodinger
equation for the two ion system as:
h2 d 2V
2m d zl
h2 d2$ /
+ 2m dz2+ (
l l2
2 0
1, 2
2 6
2zazbQ2\ _
47T€0r 3 J “
’
(
^
41
CHAPTER 2. MODELING OF M A TE R IA LS
where za and Zb are the instantaneous displacements of the negative charges.
The part of the potential proportional to r -3 represents the interaction between
the instantaneous dipoles. This equation can be transformed to a Schrodinger
equation for two independent simple harmonic oscillators by using reduced co­
ordinates, and then the energy of the interaction can be found to be:
i
rdisp _
ij
~
B i,j
r6
’
(2.36)
where B i j is related to the polarisability of the ions in question. In practice,
the Bi j are parameters to be found.
2.2.4
Form for the rigid ion potential
We can put all this together to give a simple functional form for a rigid ion
effective pair potential:
(2.37)
This is called the Born-Mayer-Huggins form, and the parameters A i,j, Bi j and
Pij need to be found for each pair of species in the system.
2.2.5
Fitting the parameters
Once we have decided on a functional form for our potential, such as equation
2.37, we need to find the parameters th at reproduce the properties of the ma­
terial to be studied as closely as possible. W hat exactly this involves depends
on what we intend to fit the potential to. The obvious answer is to reproduce
experimental observables; however with the widespread use of first-principles
modeling, we can also fit a potential to ab initio calculations. In some cases,
the best potential is found by fitting with both empirical and ab initio data.
CHAPTER 2. MODELING OF M ATERIALS
2.2.5.1
42
E m p iric a l p a ra m e tris a tio n
For a system containing n species of atoms, we need to fit (3 x Cg) parameters
(where C is the binomial coefficient). We therefore need this many independent
material properties, such as the cohesive energy, defect energies, lattice and
bonding distances, elastic and dielectric constants, phonon spectra. To reduce
this number, we can perhaps ignore the cation-anion dispersion term, or cationcation interaction, because cations are less polarisable and thus the dispersion
energies are usually smaller.
The principal benefit of fitting to experiment is th at it gives good results.
You can also engineer what the potential is good at modeling by fitting it with
the properties which are directly relevant. However this is also its main disad­
vantage, as although the potential may work well in the equilibrium area where
it is fit, there is no guarantee it will work away from this area. In other words,
the potential is not transferable away from its fitting locality, a large factor of
this being th at we have forced the complexity of the system to be represented
by a simple functional form. Experiment observes only macroscopic properties,
which are averaged over all the different microscopic possibilities. Modeling is
a microscopic process, so something may be lost in the fitting procedure.
2 .2.5.2
A b in itio p a ra m e tris a tio n
With ab initio methods, we can directly calculate microscopic properties, such
as the energy of interaction between two atoms as a function of the distance
between them, which is essentially the potential itself. We also do not need
to know in advance the functional form of our potential. So we avoid the po­
tential pitfalls of empirical parametrisation, finding a model that accurately
reproduces the microscopic behaviour of the atoms in the system in any ther­
modynamic state. The downside to ab initio fitting is of course the ab initio
CHAPTER 2. MODELING OF M ATERIALS
43
calculations themselves, which often do not agree with equivalent experimen­
tal data. Neither do the different ab initio methods always agree with each
other, different exchange-correlation functionals able to come up with signifi­
cantly different binding energies for example. So although this fitting method
is potentially more qualitatively accurate, it is quantitatively more suspicious
and therefore not necessarily as useful as empirical parametrisation.
2.2.5.3
Successes and failures o f rigid ion p oten tials
Through the use of rigid ion potentials, a lot of interesting systems have been
simulated without the time restrictions imposed by ab initio. These potentials
are widely successful for describing a crystal because most ions in a solid look like
free ions, and it is therefore reasonable to assume they interact like free ions. For
other systems like molecules, the electronic structure of the ions becomes more
important, but a potential will still be successful if the effective pair-interactions
are still a good approximation.
Despite this success, there are some im portant physical features of various
systems that they cannot reproduce, most notably the dielectric behaviour of an
ionic crystal. Our ions as described by the rigid ion model are not polarisable.
Although we include in an averaged way the interactions of the instantaneous
dipole moments of each pair of ions, we haven’t explicitly given them the ability
to have a dipole moment, so they obviously will not interact with an exter­
nal electric field correctly. The shell model presents one way of incorporating
polarisability into our model.
2.2.6
The shell model
To include polarisability, we could just assign each ion a dipole moment based
on its position and the electric field at th at position. However this ignores the
CHAPTER 2. MODELING OF M ATERIALS
44
change in the shape of the electron cloud of each ion, a change which will have
an affect on our description of the short range interaction. The shell model
includes polarisability by considering each ion to be composed of a massive core
and a massless shell. The core and shell interact only by a harmonic spring,
and the total charge of the ion is split by some appropriate fraction between
the two. The displacement of the shell from the core thus gives the ion a dipole
moment. The short range potential acts only between the shells of the ions, and
all coulombic interactions (except between core and shell of the same ion) are
included. The polarisability
of ion i is then:
yshelt
Oi =
,
(2.38)
where Z f hel1 is the charge of the shell, and k{ is the spring constant, which
can be found by fitting to dielectric data. Giving the shell zero mass means
that the shell follows the core adiabatically; this is called the static shell model.
It is possible to give the shell a small fraction of the mass; the dynamic shell
model. Given th at the fraction is small enough, the artificial shell oscillation
frequency will be much greater than the core oscillation frequencies and thus in
our classical movement of the ions the shell motion is still adiabatic.
Chapter 3
Finite temperature modeling
Chapter 2 gave us two methods for finding the energy of a system: the first
quantum mechanical, the second classical. If we are interested in the zero tem­
perature behaviour of the system, this can give us a full picture of its equilibrium
state. However, if we are interested in the finite tem perature behaviour, such
as phase transitions, the energy will not be enough to give us the necessary
information. In addition, experiments are commonly conducted at tempera­
tures of 100K upwards, so to properly compare with experiment we may need
to simulate at similar temperatures.
There are two commonly used methods for the simulation of materials at
finite temperature. Both aim to calculate the average of some variable (for
example a pressure or force, a distance), but use two very different approaches to
do so. Molecular dynamics, as discussed in section 3.1, simulates the dynamics of
the system in a life-like way, moving the particles in the system around according
to Newton’s equations of motion.
Averages are calculated by accumulating
information over time and then taking their time average. The second method,
called Monte Carlo, on the other hand, does not attem pt to reproduce the
45
CHAPTER 3. FINITE TEM PERATU RE MODELING
46
dynamics of the system. Instead, it picks a series of possible configurations
of the system randomly, ascribes each configuration a probability based upon
its energy, and then the average is found by the statistical mechanical average.
Statistical mechanics provides the link between the microscopic and macroscopic
properties of a system, and is important in molecular dynamics as well as Monte
Carlo. Section 3.2 gives a brief introduction to statistical mechanics.
I use molecular dynamics in all my work, but not Monte Carlo, For a
discussion of the Monte Carlo method, the reader is directed to Frenkel and
Smit [32].
CHAPTER 3. FINITE TEM PERATURE MODELING
3.1
47
M olecular dynam ics
Molecular dynamics (MD) is a technique for simulating the mechanics of classical
many-body systems, that can be used within DFT or using interionic potentials.
It uses Newton’s equations of motion to evolve a microsystem forward in time,
and from this can be extracted many equilibrium and transport properties of the
system. MD is a classical approximation to atomic mechanics, so when solving
Newton’s equations of motion we only consider the classical motion of the atoms.
This is a reasonable approximation at most temperatures, as quantum effects
become significant only somewhat below the Debye temperature, which is a few
hundred kelvin in most materials.
The operation of a basic MD program is shown in figure 3.1. The rest of
this section explains each step.
3.1.1
Simulation parameters
For an MD program to simulate the behaviour of a system of particles, it needs
to know many things about the system. Foremost, all the different species of
particle in the system, and how many there are of each, plus their mass and
charge. On the most elemental level, a particle is simply an ion.
When MD is performed using an interionic potential (such as that in equa­
tion 2.37) the related parameters are required to calculate the forces between
atoms in step 3 of figure 3.1. It is normal to specify a cutoff length lcut on
the calculation of all the short range components of the force (which typically
means all but the Coulomb term, which is calculated by Ewald summation).
This means that for two ions which are more than a distance of lcut apart,
the force is not calculated and is assumed to be zero. This can greatly reduce
simulation time without affecting the dynamics of the ions, as it is nearest and
next-nearest neighbours that have the greatest influence on what an atom does
CHAPTER 3. FINITE TEM PERATURE MODELING
48
START
1. read in parameters
for the simulation (e.g.
T, timestep, potential)
2. initialise the system
(give particles initial
positions and velocities)
3. compute the forces
on all particles
repeat
N mv
'
times
4. integrate Newton's
equations of motion
5. print averages of
measured quantities
STOP
Figure 3.1: Flowchart showing the basic operation of a program to perform
molecular dynamics.
CHAPTER 3. FINITE TEM PERATURE MODELING
49
next. If a cutoff length is used, it needs to be read in with the declaration of
potential.
The program also needs to know how long you wish to simulate for, and the
timestep of the motion. The timestep is the interval at which we move the ions
and calculate the forces. At the end of one interval and beginning of another,
the force is calculated and the equations of motion integrated, and each ion is
moved to a new position. The length of the simulation t8im is then equal to:
tsim — N mv At
(3-1)
where A t is the timestep, and N mv is the number of movements of the ions to
be completed. The timestep should be chosen according to two considerations.
If we choose a timestep th at is too small, we will waste a large amount of
computer time. But if we choose a timestep th at is too large, the system may
become unstable. We want to model the motion of the atoms so th at a plot
of some coordinate against time will be smooth and continuous. Therefore, an
approximate rule for calculating an appropriate timestep is to take the period
of the highest frequency vibrational mode in the materials of the system, and
divide by 20. A good choice of timestep for many systems is lfs, but not all
systems: in particular, those which include hydrogen, whose small mass lends
itself to very high frequency vibration.
3.1.2
Initialisation
3.1.2.1
In itia l p o sitio n s
The decision of where to position all the ions in the system at the beginning of
a simulation should not be influential on the subsequent equilibrium behaviour
of the system, but this does not mean we can simply assign random positions.
CHAPTER 3. FINITE TEM PERATU RE MODELING
50
For example, two ions within a molecule, placed at a distance th at is small
enough for their atomic cores to overlap, will repel each other strongly, and the
initial force between them may be enough to blast the molecule apart. So at
the least, all ions want to be positioned a minimum distance away from each
other, something of the order of 1.5A.
The best decisions made depend on the particular type of material under
investigation. For a crystal, the most sensible initial configuration is to have
all ions on their zero temperature lattice sites. For a liquid, it is also sensible
to place each constituent molecule on, for example the lattice site of a cubic
lattice, with a lattice parameter equal to the nearest neighbour distance in the
real liquid at th at temperature. Assuming you do not wish to simulate near the
freezing curve, the structure should melt quickly as it is not stable.
If you have already successfully simulated the system at a nearby state point,
it is a good idea to use the final positions of the ions in this simulation as the
initial positions, provided you can adjust the state variables to those you are
interested in. For example, if you have a liquid equilibrated at 300K, the final
configuration is a good choice of the initial configuration for the same liquid at
400K.
3.1.2.2
Initial velocities
The other side of initialisation is the assignment of initial velocities to the ions.
We want to give each ion a random velocity, but one appropriate to the tem ­
perature at which the simulation will proceed. One way to do this is to assign
each velocity component a random number drawn from the uniform distribution
[—0.5,0.5]. Then we can shift the velocities such th at the total momentum is
zero. To give the system the correct desired tem perature T, we can find the
current temperature T' using the equipartition theorem, which states that each
CHAPTER 3. FINITE TEM PERATURE MODELING
51
degree of freedom has associated with it an energy k T /2. So T ' is given by:
N
3 mv2
(M )
t = l 01=1
where N is the number of ions in the system, m* is the mass of ion i, and v,
is the velocity of ion t. Then we can scale all the velocities by (T j T ' ) 1/2 to get
the system at the appropriate temperature.
However this method gives us a velocity distribution th at is on average
uniform, which is non-physical. At a given tem perature, the distribution of
velocities (in one direction) will be Gaussian, and can be derived from the
Maxwell-Boltzmann distribution. This gives the probability distribution for
one component of the velocity as:
(
777 \ 1/2
J
exp ( ~ m v 2J 2 k T ) .
(3.3)
We can approximate this distribution by generating a large number N of random
numbers from the uniform distribution [—1,1], and adding them together to find
a number d in the range [—iV, iV], If a large number of d are calculated, the
resulting distribution is approximately gaussian, and exactly guassian in the
limit of infinite N and d. The values of d can then simply be scaled, and this
number attributed to a component of the velocity of an ion. In practice, N = 12
is generally considered to be large enough.
3.1.3
Calculation of the force
The program now has all it needs to begin the simulation. Its first step is to
calculate the overall force on each ion. For systems where periodic boundary
conditions (PBC) are to be used (discussed in section 2.1.7.1), the short range
components of the force F ij are only calculated between ion i and the nearest
CHAPTER 3. FINITE TEM PERATURE MODELING
52
periodic image of j : all other ions j in the infinite lattice are assumed not to
contribute. To calculate the force, the first thing is to calculate all the interionic
distances r i j that apply. Then the force is calculated from the derivative of the
potential.
There are techniques to speed up the calculation of the r, j , such as the Verlet
list [102] which can be used when a cutoff length is specified. In this method, a
second cutoff length ly is defined, with ly > lcut , and at the beginning of the
calculation a list of all ions within a distance ly of ion i is made, for each ion.
Then for subsequent calculations of the force, only those ions on the list need
be considered in the force calculation. Of course this list needs to be updated
throughout the simulation, either at regular intervals or when any ion moves a
distance greater than (ly —lCut)~
Another technique for increasing efficiency is to use cell lists [50]. This
approach divides all the space of the system into cells, and a force between two
ions is only calculated if they are in the same cell or neighbouring cells (of which
there are 26 in 3D). The success of the technique thus depends on the size of
the cells. It is possible to combine the Verlet and cell list approaches, by using
a cell list to construct a Verlet list.
3.1.4
Integrating the equations o f m otion
Integration of the equations of motion - that is, calculating where the ions will
move to - is done by one of a number of algorithms developed for the purpose.
The simplest, which is also accurate and stable, is the Verlet algorithm.
CHAPTER 3. FINITE TEM PERATURE MODELING
3.1.4.1
53
V erlet a lg o rith m
This is derived from the following Taylor expansion of the coordinate r(t) of a
particular particle around time t:
r (t + At) = r(t) + v (t)A t + ^ a (t)A t2 + ib ( t) A t3 + 0 (A t4),
1
o
(3.4)
where v(f) is the velocity, a(f) the acceleration, and b(t) the third derivative
of r(t), at time f, and At is the timestep. Writing equation 3.4 again with
At —►(—A t), and adding to the above, we get:
r(t —At) = 2r(t) —r(t —At) + a(t)A f2 -+- 0 (A t4).
(3.5)
a(t) is found by dividing the force on the ion in question, calculated as in section
3.1.3, by its mass. Truncating the series at the fourth term means an error in
the positions of order A t4. We will also want to calculate the velocities of the
ions. These may be needed to calculate the kinetic energy of the system, or as a
measure of instantaneous temperature, both properties that are crucial to check
whether a simulation is progressing reasonably. The velocities can be estimated
by subtracting equation 3.4 from its (—At) equivalent, giving;
v{t) = *
+ A \-A f -
^
+0 (A n
(3.6)
These will only be accurate to order A t2.
3.1.4.2
V elocity V erlet a lg o rith m
The velocity Verlet algorithm provides a more accurate way of finding the ve­
locities, without changing the trajectories of the ions as calculated by the Verlet
CHAPTER 3. FINITE TEM PERATURE MODELING
54
algorithm. It can be summarised as follows:
r(t + At) = r(t) + v (t)A t + ^ a (i)A t2;
v(t + At) = v(t) + ^ (a(t) + a (t + At)) At.
(3.7)
Firstly, the positions r(t 4- At) are found as according to the Verlet algorithm,
equation 3.5. Then the new acceleration a(i + A t) and velocity v (t + At) are
calculated.
The velocities calculated by equation 3.7 are accurate to order A t4. The
velocity-Verlet algorithm is then in theory better than the Verlet algorithm.
However in modern computing simulation, where the position and its derivatives
can be stored to 13 significant figures, they will give the same results.
3.1.4.3
P re d ic to r-c o rre c to r a lg o rith m
Another algorithm that can calculate r(t + At) and in principle all of its deriva­
tives is the predictor-corrector algorithm. This uses the same taylor expansion
as in equation 3.4, usually truncated at the A t3 term. The evaluation of this
equation gives a projected position, rp, and the equivalents vp, ap and bP be
found in the same way. Then, the forces are calculated using r p, to give a new
acceleration ac which is in general different from ap . The difference between
the two is called the error signal Aa, the result of truncating the series. A term
proportional to A a is added to all the projections, for example for the position:
r(t + At) = I-** + cAa.
(3.8)
The coefficients {c} are constant during the run, and are usually chosen such as
to maximise the stability of the system.
The predictor-corrector algorithm does not in general conserve the energy
CHAPTER 3. FINITE TEM PERATU RE MODELING
55
of the system. This can lead to the total energy of the system drifting over a
simulation run. The Verlet algorithm is preferred, because with this the energy
of the system will not drift provided the timestep is not too large.
3.1.5
Equilibration
Whenever we are starting a new simulation on a system, there will be a period
at the beginning when it is not in thermal equilibrium. We may prepare the
positions and velocities of each ion in the system well, but nevertheless, at the
beginning of the simulation the properties of the system may change rapidly.
For example, in the case where we start a simulation of a liquid by placing
each molecule on the site of a cubic lattice, some time will pass before this
ordered arrangement has disappeared and the radial distribution function (i.e.
the probability of finding a molecule at a certain distance from another molecule)
takes the correct form. Once the system has reached therm al equilibrium - that
is, the properties of the system no longer change with time (on average) - the
production run can begin and real measurements can be taken.
CHAPTER 3. FINITE TEM PERATU RE MODELING
3.2
56
Statistical mechanics
Molecular dynamics can give us the trajectory of every ion in a system of ions at
a finite temperature. This allows us to observe explicitly how different materials
behave on very small timescales. However, this kind of detail is useless when
we come to compare our insights with experiment, which measures macroscopic
properties averaged over many of these microscopic configurations. Statistical
mechanics provides the bridge between these micro and macro pictures.
The state of our system of N classical ions is specified by the positions r*
and momenta pi of all ions. Together these constitute a 6iV-dimensional phase
space, and the particular point in phase space R = [{r»} , {pt}] at which we find
the system at a given time determines the total energy of the system E = E (R)
at that time. We have to assume th at there are a number of points in this phase
space th at give the same energy. In lieu of any other information, we can assume
th at each point R is equally likely.
There are a set of points {R} in phase space th at it is possible for the
system to be in. We can imagine th at we have large number of mental copies
of the system, one copy for each possible point in phase space. Then this
set of copies is called a statistical ensemble.
There are a number of types
of ensemble which I will consider when performing molecular dynamics: for
example, the microcanonical ensemble. This corresponds to a set of system
copies th at all have the same energy, volume and number of particles contained
in them. Simulation of the system will then find it identical to one of these
copies at any point in time. Below I detail this and other possible ensembles.
3.2.1
M icrocanonical ensemble
In this ensemble, we keep the number of ions in the system, the volume of
the system, and its energy all constant. When we simulate the system in the
CHAPTER 3. FINITE TEM PERATU RE MODELING
57
microcanonical ensemble, we are sampling a constant energy surface in phase
space. The probability density of finding the system at a point R ' on the surface
is equal to th at of finding it at any other point R also on the surface.
It is this ensemble that the Verlet algorithm of section 3.1.4 samples in an
MD program. The use of any other ensemble requires the use of some additional
trick: see for example, the section on thermostats below.
3.2.2
Canonical ensemble
This is like the microcanonical ensemble, except we allow the energy to vary but
keep the temperature T fixed. This is done by considering the system to be in
thermal equilibrium with a large reservoir of tem perature T. Particles inside the
system can therefore exchange energy with the bath, although the total energy
of system and bath combined is fixed. We can express the probability density
that at any moment in time the system has a position in phase space R ' and
energy E by the Boltzmann distribution in the classical limit:
P(E) = ±exp(-l3E),
where
(3.9)
= 1 /k T (k being the boltzmann constant), and Z is the partition
function, defined as:
J dpdr exp(—8 E)
(3.10)
We can use this to calculate the thermal average (or ensemble average) of some
observable A , which we may be able to compare with an experimental measure­
ment:
(A) = j d pdtA P (E ) — ^ J d p d r A e x p ( —0E )
CHAPTER 3. FINITE TEM PERATURE MODELING
58
Taking as our observable the energy, we can see that:
IPS -
}~
=
I dPd r E e x P ( ~ E / k T )
J dpdrexp(—E / k T )
1
J
d In [/d p d re x p (—E / k T ) \
d { l/k T )
(3.13)
d\n Z
d { lfk T )
(3.14)
The Helmholtz free energy F is related to the energy of a system by the relation
F — E — T S . Dividing by T and taking the derivative with respect to (1 / T )
gives us the relation:
Therefore, we can see that F is related to the partition function Z:
F = -k T \n Z .
3.2.3
(3.16)
Therm ostats
As mentioned above, our algorithm for integrating the equations of motion in
molecular dynamics is only capable of sampling states in the microcanonical
ensemble. However, a number of methods have been devised to work around
this obstacle. Here I look at the extended Lagrangian formulation used by
the therm ostat of Nose [80], which can also be applied to other ensembles, for
example with fixed pressure. I also briefly look at the alternative Berendsen
therm ostat.
3.2.3.1
N ose th e rm o s ta t
We want to pretend that our system is in thermal equilibrium with a large reser­
voir of tem perature T, and do MD within the canonical ensemble. Therefore,
59
CHAPTER 3. FINITE TEM PERATU RE MODELING
we introduce an additional degree of freedom into the system, called s, which
acts on behalf of this reservoir. The total energy of the system + reservoir will
be constant, but the energy of each fluctuates with time.
If we take {rj} to be the real coordinates of the particles in the system,
with conjugate momenta {p^}, we can introduce virtual variables r* = r ' with
conjugate momenta pi = spj. The additional degree of freedom also acts on the
timestep of the motion such th at the virtual tim e step is At = sA t. We then
postulate an extended lagrangian:
(3.17)
where s is our additional degree of freedom, Q is an effective mass associated
with s, and g is a parameter to be fixed. The momentum conjugate to s is:
(3.18)
This gives us a Hamiltonian of the extended system of:
N
(3.19)
The extended system generates a microcanonical ensemble of (6N + 2) degrees
of freedom, and the partition function of this ensemble is:
H-Nose)
(3.20)
If we define the function in the square brackets to be /(s ), we can evaluate
CHAPTER 3. FINITE TEM PERATURE MODELING
60
the delta function by using the property 6 [/(a)] = <S(s —s o ) /|/,(s)| where so is
the single root of f(s). If we define the Hamiltonian of the real system to be
J /( r ',p ') = E p ' 2/2 mi + V (r,Ar), then:
.3 J V + 1
(3.21)
where:
(3.22)
Equation 3.21 can then be simplified by integrating the delta function over s,
to give:
(3.23)
where we have chosen g = 3N -+- 1. This partition function looks identical to
that for the canonical ensemble.
The Nose thermostat is usually now implemented in the formulation of
Hoover [52, 53], which simplifies the equations of motion.
3.2.3.2
B e re n d se n th e rm o s ta t
The thermostat of Berendsen [7] pushes the instantaneous tem perature Tjnst
towards the desired temperature T by scaling the velocities at each timestep,
by a factor %'
(3.24)
CHAPTER 3. FINITE TEM PERATURE MODELING
where
tt
61
is a relaxation constant. This is a straightforward way of simulating
at a given temperature, but it does not generate genuine trajectories within a
canonical ensemble. It has an error associated with averaged quantities of order
1/N.
3.2.4
Ergodicity
To calculate the average value of an observable in an ensemble where each
possible value A has an associated probability density p, we have to find:
dpdrAp.
( 3 .2 5 )
(A) is called the ensemble average of the system.
In molecular dynamics simulation, by its nature we cannot perform ensemble
averages, only averages in time. If we take for example the average density of
atoms at a distance r from atom i over a time t, this is given by:
(3 .2 6 )
We must assume th at provided t is long enough, pi(r) does not depend on the
initial conditions of the system, that is the system’s initial configuration in phase
space. If this is true, we can average over a number of initial conditions and get
the same result:
£ /< ? ( lim t->oo H o
Pi { r ) =
P t ( r , r ;v ( 0 ) , p N ( 0 ) , t ,) d t ')
number o f IC s
(3.27)
where (rAr(0),piV(0)) denotes the position in phase space at time t' = 0. If we
consider the case where we want to average only over ICs compatible with the
chosen ensemble, for example the microcanonical ensemble with fixed N , V , E ,
CHAPTER 3. FINITE TEM PERATURE MODELING
62
then:
Pilr)
= fl(JV,V,B) / E drA' dP'V
*(r,rw(0),pw(0),O<*').
(3.28)
where the integral is over all ICs with energy E, and ft(N, V, E) is the associated
partition function. If we swap around the two integrals, we can recognise the
ensemble average (p*), so that:
1 f*'
pi^ =t¥SolJ0
W r>fJV(0)’PJV(°)»*,)d*')tfV£-
(3-29)
The ensemble average does not depend on the length of the simulation t , so this
reduces to the ergodic hypothesis:
Pi(r) - {pi(r)) -
(3.30)
So we can calculate the average of an observable either by computing its time
average, or its ensemble average. This holds true for a vast number of systems,
but it should be noted th at many systems are not ergodic, for example nearly
harmonic solids. In a harmonic system, different degrees of freedom do not
exchange energy. The normal modes of the material obey the equation of motion
r + u ;2r = 0, and have constant total energy. A system which is nearly harmonic
will have transfer of energy between degrees of freedom, but the sampling of
phase space will be slow.
3.3
Therm odynam ic integration
In chapter four, I present a method for the calculation of surface free energies.
It uses thermodynamic integration [32], which is a technique for calculating
differences in free energy.
CHAPTER 3. FINITE TEM PERATURE MODELING
63
The free energy of a system, or some state of a system, is a very important
quantity. The second law of thermodynamics states that, for a closed system
with fixed energy, particle number and volume, the entropy 5 is at a maximum
when the system is in equilibrium. The free energy is defined as F = E — T S ,
so F (at constant volume and temperature) is at a minimum in equilibrium.
If we wanted to know which of two phases of a material was stable at a given
state point (i.e. temperature, volume), we would compare the free energy of
each phase.
So F is an important quantity, but its calculation is not as straightforward
as for other variables, such as pressure or tem perature. These are found by the
ensemble average of some function f = /(R (£)) where R (t)) is the trajectory
of the system in phase space. F cannot be expressed in this form. If we look
at equation 3.16, we can see th at the free energy is determined by the partition
function of the system, which is related to the volume in phase space th at is
accessible to the system. The entropy is another variable like the free energy,
that cannot be calculated directly in a simulation.
In order to calculate a free energy then, we need to relate it to what we can
calculate, for example the pressure p and energy E . From the definition of the
free energy, we can see that:
dF = dE — T d S - S d T
= - SdT -pdV ,
using the first law of thermodynamics dE = T d S — pdV. The derivative of the
free energy with respect to volume at constant tem perature is:
(3.31)
64
CHAPTER 3. FINITE TEM PERATURE MODELING
We also have equation 3.15 relating the free energy to the energy:
(3.32)
So to find the difference in free energy A F between two phases a and /?, we
can find a reversible and continuous path in the V — T plane from a to /?,
and integrate equations 3.31 and 3.32. In practice, this means calculating the
pressure and energy at a series of points on the path.
It is possible to calculate the absolute free energy of some state if we can link
it to another state whose free energy is known. For example, the free energy of
the ideal gas is exactly known, so the free energy of a liquid can be found by
calculating A F between the liquid and the ideal gas. Another state whose free
energy is known is the low tem perature harmonic crystal.
The thermodynamic integration over volume and/or temperature above is
essential when the integration path must be over an experimental variable. In
a simulation, we can use this technique too, and in fact we have an advantage,
because the path we choose to integrate over does not have to be physical. If
we can express the potential energy as a function of some variable A, we can
calculate the free energy difference between a state with A = 0, and another with
A = 1. For example, A can switch on some additional interaction, or external
field.
So we can write the energy as a function of A: E = E ( A), and also the
partition function: Z = Z ( A). Rewriting equation 3.16 for the free energy:
F = - k T \ n Z = -Jfcrin
J
dpdr exp (—/3E(A))
(3.33)
65
CHAPTER 3. FINITE TEM PERATURE MODELING
where
= 1/kT. We can then differentiate F with respect to A:
VT
-kT
f dpdrexp (—0E(X))
J
dpdr(-/3)
exp (—
A)) (3.34)
(3.35)
So the derivative of the free energy with respect to our path parameter A is
equal to the ensemble average of the derivative of the energy with respect to
A. The free energy difference between a system corresponding to A = 0, and
another with A = 1 is then equal to:
(3.36)
This means we can calculate a free energy difference by the calculation of an
ensemble average.
I have stated that the path over which we integrate needs to be continuous
and reversible. This is an important condition because in general, the simplest
path between two states of a system does not fulfill it. When a solid for example
is heated up, it will melt some time after the point at which the solid-liquid phase
transition should have occurred. And when you cool the liquid back down, it
will crystallise at a lower temperature than it melted at. These types of phase
transitions and many others exhibit this hysteresis, because there is a large free
energy barrier separating the two phases at or near coexistence. In such cases,
some trick is needed to ensure that the integration path is reversible, and one
such trick is used in my calculation of surface free energies in chapter 4.
Chapter 4
Surface free energy of
TiO2(110)
In this chapter, I present a general method of thermodynamic integration to
calculate the free energies of surfaces. Surface free energies are very important
in determining the equilibrium shape of crystals, as they govern the relative
stability of one possible surface over another (as recognised as early as 1901 by
Wulff [105, 49]). Away from zero temperature, entropic effects can become too
large to ignore and therefore the surface energy is no longer useful in this respect.
Despite this fact, there has been little work reported on the calculation of surface
free energies in comparison to surface energies. In section 4.4 I outline a general
method of calculating them, which can be applied using many simulational
methods.
This involves a variant of thermodynamic integration, as discussed in section
3.3. In this case, the parameter A over which we integrate is related to the strain
on a repeating unit cell caused by the presence of a vacuum gap. Suppose we
have a unit cell of bulk crystal, and then we stretch the cell to open up a
66
CHAPTER 4. SURFACE FREE EN E R G Y OF TIO2(110)
67
vacuum gap in a continuous way. The size of the gap is linearly related to a
strain parameter we call s, and we can find the stress on the system at any
value of s. By picking an appropriate range of values of s we can integrate the
stress to give the reversible work required to form the surface. Temperature
integration can also be used to check the results and give the dependence of the
surface free energy on temperature.
In this chapter I apply the method to the titanium dioxide (110) surface
using density functional theory. I aim to show th at the method is feasible even
for implementation with computationally demanding ab initio codes, and indeed
is accurate at zero temperature, where the free energy is equal to the surface
energy.
TiC>2 is a prototypical metal oxide and industrially important substance,
much studied experimentally and theoretically. Its (110) surface is the most
energetically favourable, and despite the large amounts of lab and computer
time given over to it there has not been universal agreement as to its structure
in the past. It has properties not shared by other isostructural dioxides, and
its ab initio modeling is particularly sensitive to the way electron exchange
and correlation are handled, as I shall show in section 4.3. I have spent some
time addressing the effect of functionals and other simulational details on the
modeling of the material, in preparation of the surface free energy calculations.
Section 4.4 gives the general theory of the thermodynamic integration method;
section 4.5 applies the general theory to the particular case of the TiC>2(110) sur­
face, and section 4.6 gives the results of this application using density functional
theory. First, I shall give an overview of previous studies that have calculated
surface free energies, and then the material I shall study in more detail. Section
4.2 will then talk about how all the simulations were performed.
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2 (110)
4.1
Introduction
4.1.1
Free energy calculations
68
As already noted, the equilibrium shape of crystals is governed by the relative
free energies of the exposed surfaces. The crystal will have a number of different
surfaces of varying size, and the equilibrium shape will be that which minimises
the free energy. At high temperatures, the role of entropy becomes more im­
portant, find if surface entropy effects are significant, the difference between the
free energy and internal energy of a system may be great, so th at the latter will
be a poor guide when comparing the forms of crystals. An example of this is in
the growth of natural crystals and the sintering of ceramics, both of which are
controlled by surface free energies near the melting tem perature of the material.
There has been a large amount of work on the calculation of surface energies for
a wide variety of materials, but comparatively little work has been reported on
the calculation of thermodynamic surface free energies. At present, very little
is known about the temperature dependence of surface free energies.
There are a number of known methods for calculating free energy differences
[40]. For certain systems, it is suitable to calculate the free energy difference
between two states from direct counting of the number of configurations in each
state. This is because the free energy is related to the partition function Z by the
Helmholtz relation F = —k T l n Z . The partition function counts the number of
possible microstates of the system in a given state, so a free energy difference
between states a and b can be found by AF = —k T In (ZajZ^). This method
is only useful when the two states are similar and can be observed within one
simulation, and so is obviously not applicable to the calculation of a surface free
energy.
A perturbation approach can also be used. Here, we express the change in
the hamiltonian between the two systems of interest as a series of small pertur­
CHAPTER 4. SURFACE FREE E N E R G Y OF T I 0 2(110)
69
bations, and then adding together the effects on the free energy of each of these
perturbations. Given the Helmholtz relation, a change A A in the hamiltonian
will give a derivative of the free energy:
dF(A) _ F (A + AA) - F(X)
d\
AA
kT
AA
’J d p d re x p { —f3H(\ + AA)}"
f dpdr exp {—jOH(X)}
kT
- — In[{exp{ - /? (H (A + AA) - #(A ))})a]
(4.1)
where within the logarithm is the thermal average of a ratio of Boltzmann factors
at a given A. For some path of A, a free energy difference can be given by the sum
of these small perturbations AA. If the two systems of interest axe sufficiently
similar, it is possible to take AA to describe the full difference between the two
systems, and so equation 4.1 needs to be calculated only once.
Another method is to calculate the potential of mean force, where the free
energy is expressed as the logarithm of the probability of finding the system at
some reaction coordinate J2, which is some coordinate of the system, typically
a spatial coordinate. The probability P ( R ) can then be found by letting the
system explore the possible configurations corresponding to possible values of
the reaction coordinate. To improve the sampling of the configurational space,
an umbrella potential can be added to the Hamiltonian to ensure th at certain
areas of the space are sampled, for example high-energy configurations.
Another method used is the quasiharmonic (QH) approach. Here, the full
interatomic potential is replaced by its quadratic expansion about the atomic
equilibrium positions. The vibrational modes of the material are calculated, and
then the free energy is a sum of the total energy of the equilibrium crystal, plus
a summation of terms contributed by each vibrational frequency [45, 29, 65].
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
70
A widely used general method for calculating free energy differences is the
lambda-integration or thermodynamic integration method. Here a system is
switched reversibly to another system by adjusting a continuous parameter A.
The derivative of the internal energy with respect to A can then be integrated to
find the difference in free energy between the two systems. Grochola et al. [38]
and Foiles [29] use this method to find the surface free energy of iron b c c (lll)
and copper (100) respectively. Hansen et al. [45] use A-integration to find the
surface free energy at a reference temperature, and then use the method of
tem perature integration to find its temperature dependence. Davidchack and
Laird cleave a surface in a hard-sphere crystal by introducing two walls which
force two surfaces to form, and then integrate the pressure on the walls over the
distance between them [15]. The vast majority of studies use either LennardJones [15], empirical potentials [45] or the embedded-atom method [38, 29, 76].
Only one study that I am aware of has solely used first principles modeling to
calculate surface free energies (using the quasi-harmonic method) [65].
4.1.2
Titanium dioxide and its (110) surface
Titanium dioxide is popular among experimentalists and theoreticians alike,
thought of as a prototypical metal oxide. It has a number of phases, the most
common of which is the rutile phase, which I investigate here. It has many
industrial uses, principally as a white pigment and in heterogeneous catalysis,
which means it is cheap and widely available. It also has a simple structure in
comparison with many more complex oxides, whilst sharing their physical and
chemical properties, so it can be modeled by a large range of empirical and first
principles theories to give information on itself and other materials.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
4.1.2.1
71
B u lk
TiC>2 has a tetragonal primitive unit cell, containing two titanium and four
oxygen atoms, and this structure is known by the name rutile (shown in figure
4.1). The structure and electronic properties of bulk titanium dioxide have
been examined many times, see for example Glassford and Chelikowsky [35],
and I will not attem pt to review th at here. One aspect of bulk rutile is however
relevant to my work on the (110) surface; the vibrational modes of the material
may have an effect on the calculation of surface atom displacements. Thus, I
shall review here the material’s lattice dynamics.
The lattice dynamics of bulk rutile were investigated by Traylor et al. [99]
using coherent inelastic neutron scattering in 1971. Prom this early stage, it
was clear th at the lattice dynamics were very important in determining the
technological properties of the material. For example, rutile has an exception­
ally high static dielectric constant along the c-direction, which increases as the
temperature is lowered. This is explained in terms of the transverse optic (TO)
A 2u vibrational mode, which at room tem perature is soft (173cm-1 [99]), and
which gets softer at lower temperatures. It never becomes unstable and so the
material does not undergo a ferroelectric phase transition, and is therefore clas­
sified an incipient ferroelectric. The TO A 2u mode has been often noted since
this was identified, its softness being peculiar to this material and not shared by
other rutile oxides like tin dioxide and germanium dioxide (with frequencies of
465cm-1 and 455cm-1 respectively). A number of ab initio studies have been
performed, calculating the vibrational modes at the gammapoint [66, 79, 92],
usually using LDA. This is one of many areas where the use of GGAs must be
checked against LDA results, and in this instance they do not compare favorably.
Montanari and Harrison [79] did calculations using LDA, and the functionals
PW91 and PBE. They found that the frequencies predicted by LDA were in
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
72
excellent agreement with experiment and previous calculations, but th at the
GGAs consistently fared worse. The TO A 2u mode frequency was predicted to
be imaginary using PBE, which would mean the mode was unstable and a phase
transition would occur. The mode found with PW91 was very soft. They con­
cluded th at these results were due to the GGA’s overestimation of the lattice
parameters, and demonstrates the discrepancy between a GGA and an LDA
description of a system.
4.1.2.2
The (110) surface
The (110) surface of titanium dioxide is shown in figure 4.2. There are six atoms
in the surface layer per unit cell; one six-fold titanium , one five-fold titanium,
two oxygens in the surface plane, one oxygen sitting below it and another oxygen
above it. This last atom is called the bridging oxygen. The sideways view shows
the primitive unit cell one would use to model the surface, for a slab of three
layers.
A comprehensive review of the surface science of titanium dioxide was writ­
ten by U. Diebold in 2003 [18], which pulls together all the experimental and
theoretical studies performed on the (110) surface. The first ab-initio studies of
the Ti (>2 (HO) surface were performed by Ramamoorthy et al. in 1994 [88, 89].
They found that the relaxations of the surface atoms are substantial, and that
they are responsible for a large reduction in the calculated surface energies.
Since then, the surface’s energetics and structure have been investigated within
density functional theory using various plane wave [46, 5, 69] and atomic or­
bital [96, 9] codes. Bates et al [5] found that the magnitudes of the surface
relaxations oscillate substantially with the number of layers that make up the
crystal slab. Each layer consists of six atoms in the simulation cell, and they
are stacked in an ABAB manner, so that each layer is displaced by half the
lattice parameter in the [-11 0]-direction relative to the layers above and below
73
CH APTER 4. SURFACE FREE E N E R G Y OF T I 0 2(UO)
a
Figure 4.1: The primitive unit cell of bulk titanium dioxide. Rutile structure,
with two titanium and 4 oxygen atoms per unit cell. Grey atoms are titaniums,
red oxygen. Dashed lines indicate the plane in which the oxygens lie. The
distance of an oxygen in this plane from the nearest titanium is given by u a^/2,
where u is a dimensionless fraction.
top view
[110)
[001]
top layer
/
bridging oxygen
/
0(B)
five-fold titanium
/
TifS)
middle layer
in-plane oxygens
om
/
\
\
•
-
bottom layer
s
six-fold titanium Ti(6)
Figure 4.2: Looking down on the (110) surface, and sideways at a thrcc-laycr
slab.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
74
it. The properties of the surface can then take different values, depending on
whether you use an even or odd number of layers. When an even-layered slab
is used, there is no symmetry plane in the middle of the slab and so the atoms
are able to displace a large amount. For odd-layered slabs, there is a symmetry
plane which restricts cooperative displacement through the slab, so such large
displacements aren’t possible.
Harrison et al
[46] compare results from previous plane wave LDA and
GGA studies, with their own results using full-potential linear augmented plane
waves (FP-LAPW) and linear combination of atomic orbitals (LCAO) methods,
with the PBE functional. They found poor agreement between the theoretical
results, even between the two methods they used themselves, in which they
took care to control the effects of all numerical tolerances. The comparison
of the position of the bridging oxygen in particular was unfavorable. They
postulated that the energy surface with respect to vertical displacement of the
bridging oxygens is flat, and found a very soft, anisotropic and anharmonic
surface rigid-unit mode which involves surface ion displacements of 0.15A for
thermal vibrations corresponding to room temperature. Bredow et al [9] also
found a large oscillation of the interlayer distances, surface energy and electronic
structure with increasing number of layers. They explained the effect by surfaceinduced hybridization of Ti 3d and O 2p orbitals among the layers, leading to
strong bonding between first and second layers, and weaker bonding between
second and third layers. Thus, in even-layered systems, they describe the slab
as a stacking of pairs of layers. When they removed the 3d orbitals from the
Ti basis set they found that Ti02 behaved in the same way as the isostructural
Sn0 2 (110) surface.
Hameeuw et al [42] recently modeled three different systems in an attem pt
to eliminate the large oscillation in displacements between even and odd num­
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U0)
75
bered layers. They relaxed all the atoms in an n-atom slab, and then repeated
the calculation with the two middle layers fixed for even layered slabs (to ensure
the slab was bulk-like in the centre), and then again for all slabs holding two
external layers fixed, so there is only one surface of interest. They found that
this last system was the optimal one for converging all the surface structural
and electronic properties in the quickest time, because it eliminates the symme­
try differences between odd and even layered slabs. They compare their results
for the surface atom displacements with the experimental values of Lindsay et
al. [72] (see next paragraph) very favorably. The bridging oxygen, of particular
interest, found good agreement. This shows that the large oscillation of surface
structure with slab thickness is most likely an artefact caused by the different
symmetry conditions of even and odd layers.
The benchmark experimental results of the (110) surface have, until recently,
been those of Charlton et al. [13]. Here, surface x-ray diffraction (SXRD)
was performed to determine the structural relaxations in the (lx l) surface unit
cell. They found that the six-fold titanium atoms and the in-plane oxygen
atoms relaxed outwards, away from the surface, and that the five-fold titanium
atoms and bridging oxygens relaxed inwards. These relaxations are noted as
creating a rumpling of the surface layer, which is still present in the second layer
but at half the magnitude. In the past year, a new study has been published
by Lindsay et al. [72], which studies the surface with quantitative low-energy
electron diffraction (LEED). In contrary to the previous results, they have the
bridging oxygen relaxing outwards. The magnitudes of the relaxations of the
other surface atoms are also non-negligible. The difference between the two
studies is marked, and an explanation for this was offered, by way of Harrison
et a/.’s work finding a soft anharmonic surface mode. The SXRD study was
performed at room temperature, where this mode could displace surface atoms
CHAPTER 4. SURFACE FREE E N E R G Y OF TlCh(HO)
76
by 0.15A . The effect of vibrations on the LEED interpretation however were
found to be small.
4.1.2.3
A d so rb e d species on T i 0 2 (110)
Many experimental studies have looked at adsorbed species on the surface, and
have largely concluded that water adsorbs molecularly, dissociating only at de­
fect sites and therefore at low coverages, using STM [11, 91], HREELS [47],
TPD [54, 47, 48], These defect sites are typically explained as missing bridging
oxygens. However, there is also a large amount of literature which concludes
that water adsorbs dissociatively on the surface; using Hartree-Fock methods
[23], semi-empirical cluster modeling [10], and DFT [36, 68]. This discrepancy
has repeatedly been attempted to be resolved.
Lindan et al [71] performed a series of calculations with various surface
reconstructions and coverages, paying particular attention to the introduction
of possible computational errors, in a hope to put an end to the debate. They
found that water bound more strongly to the surface in dissociated form over
a range of coverages and cell sizes, and th at the energetics of adsorption varied
considerable with coverage. When two water molecules were present on the
surface, a mixed state (i.e. one molecule and one dissociated molecule) was the
most favorable state at all but 1ML (full coverage), because of the stabilisation
offered by hydrogen bonding between them. They also found a large barrier
to dissociation, that increases as coverage decreases, due to deformation of the
substrate. On the contradiction of their results with experiment, they thought
that, at a temperature large enough to overcome the dissociation barrier, the
water molecules would be mobile enough to find a defect site. Thus, experiment
would expect to never see dissociation anywhere but here. Lindan also looked
at this earlier [70], investigating the adsorption of water to the surface using
first principles molecular dynamics. They firstly simulated the bare surface,
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U 0)
77
and then after sufficient equilibration, they placed a water molecule 3A above
one of the surface five-fold titanium atoms. The molecule was drawn in towards
the surface, and then it dissociated. The bond in the bridging hydroxyl pointed
towards the oxygen in the terminal hydroxyl, to maximise the hydrogen bonding.
They also investigated the vibrational spectra of the hydroxylated surface at
120K. A sharp peak around the water stretch mode frequency was associated
with the terminal hydroxyl, and a broad range of frequencies in this region
were assigned to the bridging hydroxyl. This was explained by the hydrogen
bond between the terminal oxygen and the bridging hydrogen: there is coupling
between the hydrogen’s motion and the broad range of low frequency modes of
the oxygen.
There have been a number of other MD studies. Predota et al. [86, 87] have
investigated the water/crystal interface using MD, but with classical force fields
for the interactions of the surface with water molecules and using the SPC /E
model for bulk water. SPC /E does not allow for dissociation of water molecules,
so they investigated the non-hydroxylated and hydroxylated surface. Water
adsorbed to the surface molecularly, for the former case, and they found above
the first layer of oxygen atoms, a second layer of adsorbed water molecules which
occupied distinct positions relating to the underlying crystal surface structure
in both cases. After this point, the structure quickly decayed. Increasing the
temperature from 298K to 448K had minimal effect on this interface structure.
The choice of which exchange-correlation functional is naturally as impor­
tant in this area as in surface structure. Casarin et al. [12] cluster-modeled
both water and carbon dioxide on the surface, and found that the use of GGA
as opposed to LDA strongly reduces the adsorption energies of the molecules
to the surface, bringing them closer to available experimental values. Most
of the studies previously mentioned favour GGA over LDA, because has been
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
78
found that GGA provides a more accurate description of molecular dissocia­
tion energies, adsorption energies and of hydrogen bonding in particular. The
Perdew-Burke-Ernzerhof functional and its revised version are commonly used.
Zhang and Lindan [110, 111] have investigated water on the surface also at
more than 1ML coverage, using MD and static relaxation calculations (and a
(2x1) cell reconstruction). Firstly, they identify four states possible at 1ML; two
mixed states, which were found to be the most stable, and two states with disso­
ciative adsorption only. To these surfaces they added another water molecule in
a variety of possible positions, and find th at there are many adsorption sites for
this third water molecule, due to a large number of hydrogen bonding possibil­
ities. Proton transfer occurred in some systems, leading to chainlike structures.
Secondly, they continue to add molecules up to a coverage of 3ML. Here they
find that the water near the surface tends to become more molecular with in­
creasing coverage, so the more layers of water that are on the surface, the less
favorable the partially dissociated structures at the interface become. They also
note that there is a struggle in the near-surface adsorbed molecules (in the 2nd
layer) between the determinate influence of the surface and first layer, and the
hydrogen-bonded network in the outer layer. It is possible that this struggle
could be tipped in the favour of the latter when more layers are added.
Kornherr [60] attempts to bridge the gap between a few molecules and bulk
water, using force fields which simulated separately the oxide slab, the adsorbed
water molecules, and the amorphous bulk water. They found that molecules
in the third and fourth layer interacted with the slab at about 10% of the
interaction of the first and second layers with the slab, and are thus very loosely
attached to it. They investigated also a surface reconstruction with defects
present, and found that adsorption at the defect was much more exothermic
than adsorption above a five-fold coordinated titanium.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U 0)
79
The study by Stefanovich et al. [93] bucks the trend, concluding from ab
initio embedded cluster calculations th at an isolated water molecule adsorbs on
the surface in the molecular form. They offer as reason for the discrepancy with
other theoretical results, the fact th at no periodic boundary conditions were
used in their calculations, so they modeled a true isolated molecule. Bandura et
al. [4] back up these results, finding that dissociation is favoured with periodic
DFT and molecular adsorption is much more favourable with embedded cluster
HF. Langel [63] study the system using Car-Parrinello MD, and find that a
water molecule placed on the surface adsorbs associatively, whereas two hydroxyl
groups will quickly recombine (at 320K).
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(UO)
4.2
80
M ethod of investigation
All calculations were done with the ab initio MD program VASP [61].
4.2.1
Pseudopotentials
We saw in chapter two th at it is necessary to replace the real potential of
the system by some sort of pseudopotential when using a plane wave basis
set to describe the wavefunction. This is because the true wavefunction will
oscillate rapidly near the atomic cores, and therefore a large number of plane
waves would be needed to represent it, which is very computationally expensive.
Three different approaches to constructing a pseudopotential were offered. In
my work, I have used ultrasoft pseudopotentials (USPP, Vanderbilt [101]) and
the projector augmented-wave approach (PAW, Blochl [8]). I expect them to
give very similar results, as they have been shown to give similar results for this
type of system [62].
4.2.2
Exchange-correlation functional
Also in chapter two, I noted how the choice of how to deal with exchange
and correlation could make a non-negligible difference to the properties of the
material simulated. We have loaded all our theoretical approximations and
uncertainties into this part of the energy functional, so the quality and suitability
of the exchange-correlation functional plays a crucial role in the accuracy of our
simulation.
There exist now a number of functionals, and although some are reported
to model certain systems better than others, the only way to assess which one
is best is to try them all (usually, for some test case). For my calculations,
I have used four different functionals: the local density approximation (LDA)
of Ceperley-Alder [84], and the generalised gradient approximations (GGAs) of
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(H0)
81
Perdew and Wang (PW91) [83]; Perdew, Burke and Ernzerhof (PBE) [82]; and
the revised PBE (RPBE) [43].
4.2.3
k-point mesh
The wavefunction describing the system will be calculated at a number of points
in k-space, the reciprocal of real space. I use the scheme of Monkhorst and Pack
[78], which uses a mesh of fc-points. The user inputs the size of mesh they want
to use: for example, a 2 x 2 x 2 mesh gives a total of 8 fc-points, with two equally
spaced in each direction. A mesh of 1 x 1 x 1 gives just one A:-point, at the center
of the first Brillouin zone (the gamma-point).
I have used a variety of meshes in my calculations. For modeling bulk tita­
nium dioxide, with six atoms per primitive unit cell, I used a mesh of 3 x 3 x 5.
The ratio between the number of fc-points in each direction is a result of the
shape of the primitive unit cell. The crystal structure of rutile has lattice param­
eters a = b > c, so more Appoints are needed in the 2 -direction as the reciprocal
cell is larger in this direction. Using a mesh of 4 x 4 x 6 leads to a difference in
energy of order < 4meV, so I conclude that 3 x 3 x 5 is dense enough.
For the simulation of the surface at zero temperature, I used a mesh of
4 x 4 x 2 . There is no advantage to using more than one fc-point in the direction
perpendicular to the surface. This is because we want our slab to be sufficiently
far away from its periodic image that any interaction between the two will be
negligible, and so a large fc-point density will not be required. Also, the lattice
parameter in this direction will be very large due to the large vacuum gap, so
the length of the reciprocal unit cell in this direction will be small. I have used
a mesh with two points in the direction perpendicular to the slab; however,
this did not waste computer time, because VASP reduces the problem to the
calculation for one fc-point in this direction, due to symmetry.
CHArTER 4. SURFACE FREE ENERGY O r TIG, (110)
82
' "1001]
Figure 4.3: Simulation cell used in molecular dynamics simulations. Primitive
unit cell has been doubled in the [001] direction.
Using a larger mesh of 5 x 5 x 2 leads to a difference in the surface energy of
order < 3meV. Using a mesh of 2 x 2 x 2 increased the error a negligible amount,
with the pressure on the cell changing by less than a kbar. This implies that
such a dense mesh as 1 x 4 x 2 is not necessary, an observation that is most
useful at finite tem perature when MD is performed and therefore efficiency is
more of an issue. Therefore I used a mesh of 2 x 2 x 2 at finite temperatures.
4.2.4
The simulation cell
At zero temperature, the smallest possible simulation cell is used. In the MD
rims, however, I could not do this because periodic boundary conditions are
imposed, therefore sill atoms of one type on the surfaces would have identical
motion, which is unphysical. So I doubled the cell in the [OOlj-direction, as
shown in figure 4.3, the simulation rell having double the number of atoms
that were in the cell at zero temperature. Static relaxation calculations at OK
show' th at this arrangement reproduces the surface structure significantly more
accurately than doubling the cell in the [110]-diicction.
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(U0)
4.2.5
83
Technical details
A plane wave cutoff of 500eV was used in all calculations. This means that only
waves of kinetic energy 500eV and less are used in the plane wave expansion of
the wavefunction. Static relaxations were performed until the forces on all the
atoms were less than 10-4 eV/ A .
A timestep of lfs was used in molecular dynamics simulations. The maxi­
mum frequency of vibrations at zero wave vector in titanium dioxide is 850cm
[66], or 26THz, so that our chosen timestep is approximately 40 times smaller
them the minimum vibrational period. The Nose-Hoover thermostat was used
to regulate the temperature of the system within the canonical ensemble.
CHAPTER 4. SURFACE FREE ENERGY OF T I02(110)
4.3
84
T i0 2 and the structure of the (110) surface
In this section I detail the work I did on th e m aterial titanium dioxide. As
reviewed in section 4.1.2, the structure of the (110) surface has been the cause
of some debate, with a lack of agreem ent between ab initio theoretical results,
and w ith experim ent. On the theory side, one explanation for this has been
the different exchange-correlation functionals used to approxim ate the electronelectron interactions. In preparation for using the m aterial to test my m ethod of
therm odynam ic integration, I have therefore thoroughly investigated the surface
structure with a num ber of functionals, and pseudopotentials to check whether
this may also make a difference.
Firstly, I have modeled the bulk unit cell. This will provide th e first test of
the functionals, as GGAs are widely known to overestim ate lattice param eters.
4.3.1
Bulk calculations
To find the equilibrium bulk lattice param eters o, 6 , c of an orthorhom bic crystal,
and in the case of TiC>2 also the internal param eter u (where ua-^/2 is the nearest
O-Ti distance), we minimise the energy function U — U (a, b,c, u) with respect
to each of the four param eters. (This becomes even more involved if the unit cell
is not orthorhom bic.) VASP performs this m inim isation process upon request,
with one caveat. I use a plane wave basis set, containing plane waves of energy
500eV and less. W hen the program changes the lattice param eters, this doesn’t
change the num ber of waves in the set, but it does change the length of the
wavevectors. This can lead to an error in the calculation of the stress tensor,
called the Pulay stress. However, if a large enough energy cutoff is used, this
error will be negligible. I used an energy cutoff 1.3 times the default value set
by the potentials.
Table 4.1 shows the lattice param eters and energy using different pseudopo-
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(110)
85
potential - functional used
a (A )
c (A )
u
energy (eV)
USPP - LDA
USPP - PW91
PAW - LDA
PAW - PW91
PAW - PB E
PAW - R PB E
4.569
4.657
4.571
4.655
4.659
4.691
2.933
2.978
2.928
2.969
2.970
2.977
0.304
0.305
0.304
0.305
0.305
0.305
-58.531
-53.565
-58.071
-53.124
-52.812
-52.263
experim ent
4.584
2.953
0.305
-
Table 4.1: Lattice param eters and energy of the TiC>2 primitive unit cell using
different pseudopotentials and functionals. Experim ental values from Diebold’s
review [18].
tentials and functionals. It shows th a t using a GGA makes a large difference
to th e calculated energy of the unit cell and to the lattice param eters. There
is very little difference between the results using ultrasoft pseudopotentials and
PAW. The largest difference in a between the systems tested is 0 . 1 2 A , which
is large, b u t not far from the largest difference w ith experim ent ( O . l l A ) . The
LDA results here seem to give the best agreem ent with experiment. All GGAs
overestim ate the lattice param eters.
4.3.2
(110) surface structure
I have investigated the (110) surface, also with different pseudopotentials and
functionals. To model the surface correctly, a slab geom etry is used, where the
slab should be thick enough for the inside to behave like the bulk m aterial.
Above and below the slab is vacuum.
Due to periodic boundary conditions (see section 2.1.7.1), what I am actually
modeling is a stack of slabs, separated by a vacuum gap L. The slab is considered
to be made up of a num ber of layers n; in this system each layer consists of
2
titanium atom s and 4 oxygen atom s. In theory, the more layers in your slab, the
more accurately the surface will be modeled. L is defined such th at the length
of the simulation cell in the direction perpendicular to the surface is equal to
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(110)
86
(ay/2 x n /2 ) + L. For an even layered slab, when L = 0 the system returns to
being bulk titanium dioxide. L should be large enough for the atom s in each
slab not to feel the atoms of the slabs above and below it. Thus, n and L are
two quantities which must be sufficiently large to converge the surface energy
and structure.
W hen a slab is modeled, the atom s are allowed to relax, but the shape and
size of the simulation cell (i.e. slab + vacuum) are kept fixed. Because the (110)
surface is found from cleaving the prim itive unit cell along a diagonal, the slab
does not have the same shape as the prim itive unit cell. This means th a t the kpoint mesh used to relax the primitive unit cell is no longer suitable to model the
simulational cell, for which I use a 4 x 4 x 2 mesh. The k-point density (number of
k-points per reciprocal cell volume) for each type of calculation is different. This
means th a t, if I construct a slab using the a and c found above, there will be a
small pressure on it which could affect the surface calculations. The way around
this is to relax the new simulation cell w ith vacuum w idth L = 0, and then fix
the size of this cell, add a vacuum gap and relax the atoms. However, doing
this would lead to each calculation of fixed n being performed with different
lattice constants, and only be possible when n is even (due to the alternating
position of subsequent layers). To test the effect of this discrepancy in fc-point
density, I performed a slab calculation in two ways: using the param eters found
from bulk calculations, and relaxing this cell before adding a vacuum gap, then
relaxing the atoms. Comparing the two m ethods, I found a difference in the
energy of ImeV and a maximum difference to the surface atom displacements
of 5mA. Thus I conclude th at using the first m ethod gives perfectly acceptable
accuracy with respect to this issue.
87
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(110)
0.0514
>CD
gj 0.0513
5>
c
a>
©
oO
•tC
30
<
0.0512
0.0511
vacuum gap L (A)
Figure 4.4: Surface energy Usurf(n = 4,L ) with increasing vacuum gap L.
4.3.3
Convergence of surface properties
4.3.3.1
...w ith respect to vacuum gap L
I performed a series of calculations to determ ine what value of L would be
sufficient to converge th e surface form ation energy t/sur/ , defined by:
tt
ti _
U s u r f v 1-! L>) —
Ut ot ( ni L ) — U t 0t ( n ,
2^4
0)
’
(A
(4.2)
where 17f0t(n, L) is the total energy of a slab of n layers, separated by its periodic
image by L , and A is the area of each created surface. Figure 4.4 shows Ulot{n =
4, L) vs. L for a slab four layers thick. The energy climbs and then levels off,
with the difference between the surface energies w ith L = 6A and L = 9A only
0.1%, so I conclude th a t
calculations.
6A
is sufficient. L =
6A
is used in all subsequent
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(110)
U S P P -L D A
en e rg y
su rfa c e en e rg y
11(6)
Ti (5)
OflP)
0(B )
uspp-pmi
e n e rg y
su rfa c e en e rg y
T»(6)
Ti<5)
OflP)
0(B )
en e rg y
su rfa c e e n e rg y
H(6»
H<5>
D (IP)
0(B)
_________________________number of (ay ra
11
7
4
6
9
10
5
6
-232.175 -29&36T -349.174 -407576 -466.234 -524719 •563-320 -641.639
0-9449
0.9226 0.9515 09309
04593
1.0105 0.9040 0.9677
0-274
0.264
0203
0276
0.296
0255
0.235
0251
-0164
-0.167
-0.166
-0.163
■0.168
-0.199
-0.156
-0.176
0.177
0.169
0.169
0.156
0.165
0.170
0.166
0.175
0.006
0.074
0.094
0.055
0.063
0.036
0.076
0.051
7
6
4
5
6
-213.126 -266322 -320.196 -373.602 -427319
06695
05066
05460 0.6203
05629
0241
0-325
0205
0326
0313
-0.177
■0.146
-0.151
-0.176
-0.165
0204
0.194
0.170
0.213
0.165
0.037
0002
0.106
0.115
0.119
4
7
5
6
4
-230397 -288.151 -346.474 -404.427 -462.614
0.9716 0.8737 0.9321
0.8915
05299
0262
0271
0221
0264
0.166
-0.177
-0.135
-0.151
-0.140
-0.140
0.162
0.174
0.165
0.178
0.168
0.089
0.042
0.099
0.009
0563
P A tV - P H & f
e n e rg y
su rfa c e en e rg y
n<6»
11(5)
OflP)
0(B )
4
7
5
4
4
-211562 -264318 -317.774 -370234 -424.030
0-4642
0.4333 05985
0.4660 05472
0254
0325
0203
0335
0336
-0.112
-0.132
-0.149
-0.139
-0595
0211
0.167
0241
0206
0255
0517
0.144
0.130
0.065
0.148
ftW-Pflg
e n e rg y
su rfa c e en e rg y
11(6)
H{5>
OflP)
0(B )
4
7
5
6
6
9
-210366 -262.805 -315.957 -366500 -421591 -474325
0.4120
05796
0.4456 05272
0.4615 05045
0257
0331
0205
0343
0349
0290
-0567
-0.145
-0.139
-5.106
-0.130
-0-116
0217
0210
0.190
0265
0229
0.249
0.139
0.152
0519
0.069
0.159
0.101
P A W -P P B t
en e rg y
su rfa c e en e rg y
H(6)
R(5)
OflP)
0(B )
4
7
5
6
6
9
-206-790 -260-600 -313329 -365339 -417.694 •470500
02057
0-1331
0.1388 02513
03301
0.1338
0.417
0217
0.473
0295
0516
0353
•0579
-0-002
-0.111
-0.138
0.046
-0.076
0.367
0295
0200
0240
0.415
0262
0.027
0-102
0219
0.276
0320
0.159
energies in eV, surface energies in Joiies per metres squared, and <Ssplacements in angstroms.
Ti(6) = sa-fold coordinated titanium atoms
TH5) = tve-Wd coordinated ttenium atoms
0(1p) = tvpiane oxygen atoms
0(B) * bridging oxygens
Table 4.2: The surface energy and structure for every system studied.
CHAPTER 4. SURFACE FREE ENERGY OF TlCh(UO)
Surface atom:
T i( 6 )
Ti(5)
USPP - LDA
PAW - LDA
LCAO - LDA [96]
LEED-IV [72]
SXRD [13]
0.27
0.26
-0 .1 7
-0 .1 4
-0 .1 7
-0 .1 9 ± 0 .0 3
-0 .1 6 ± 0 .0 5
0 .2 2
0.25 ± 0 .0 3
0.12 ± 0 .0 5
O (IP)
0.18
0.19
0.13
0.27 ± 0 .0 8
0.05 ± 0.05
89
O(B)
0.08
0.08
0 .0 1
0.10 ± 0.05
-0 .2 7 ± 0 .0 8
Table 4.3: Surface atom displacements after relaxation, compared with previous
D FT studies and experiment.
4.3.3.2
...w ith respect to slab thickness n
Finding a value of n which converges the surface properties is a more com­
plicated problem. Many studies before have found th a t th e surface energy and
atom displacements oscillate with the num ber of layers in the slab, taking a long
time to converge to a definite value. This has been widely explained, as previ­
ously mentioned, by the sym m etry breaking in even-layer slabs. An odd-layer
slab has a n atural sym m etry plane at th e centre, the atom s of which therefore
will not move from their initial positions. I also found strong oscillations of
surface properties, but the severity of these oscillations was strongly dependent
on the pseudopotential and particularly exchange-correlation functional used.
Table 4.2 has my full results, varying these two details and detailing the surface
properties for various n-layered slabs, but I make a suitable comparison in figure
4.5, which shows the oscillation of the surface relaxations with the number of
layers, for a system modeled with a) the LDA (on the left), and b) the PBE
functional, both with the PAW potentials.
The displacements in a) take some time to converge, as expected, but one
would be confident about relying on a slab of thickness
8
layers or more. W ith
system b) however, no such statem ent can be made at a thickness of
8
layers
or indeed 9. The displacements are converging very slowly, if indeed we can be
sure they are converging at all, as the displacements of even and odd-layer slab
atoms appear to be increasing.
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(U0)
90
0.4
<
0.3
0.3
0.2
y 0.2
c<D
E
s
Ti(6)
• Tl(5)
O(IP)
x—* 0(B)
jg O.H
Q.
eg
T5
ni
- 0.1
-
-
0.1
0.2
number of layers
number of layers
Figure 4.5: Relaxations of surface atom s with varying n, for a) the LDA, and
b) the PB E functional.
Comparing the oscillations of the six-fold and five-fold titanium atom s, we
see th a t they are out of phase for a) and in phase for b). This is mirrored in the
other systems I studied; LDA finds them out of phase and all three GGAs find
them in phase. These results illustrate the im portance in choosing the correct
functional and potential when modeling this system. The PW 91 displacements
behave in a similar m anner but w ith less severe oscillations, but for the R PB E
functional the oscillations are even more severe, having a m agnitude of 0.15A
between even the eighth and ninth layers.
Table 4.3 shows my results for an eight layer slab (LDA, both USPP and
PAW) compared with previous D FT study using linear combination of atomic
orbitals and the LDA, and with the two main experimental data. My results
are in good agreement with the previous study, and the LEED-IV study by
Lindsay et al. [72]. Comparison of both my results and the other D FT results
with the previous experim ental SXRD study is very poor, lending support to
the superiority of the more recent study.
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(U0)
91
surface energy (Jm 2)
US - LDA n = 4 (this study)
US - LDA n = l l (this study)
PAW - LDA n = 8 (this study)
US - LDA n = 4 [8 8 , 89]
LCAO - LDA [96]
PW - LDA [42]
0 .8 6
0.94
0.89
0.76
0.90
0.90
Table 4.4: Surface energy computed by different studies. LCAO - LDA results
from Swamy et al., with convergence w ith respect to slab thickness to within
0.01 Jm ~ 2, P W (plane wave) - LDA results from Hameeuw et al., using a four
layer slab with the two bottom layers held fixed.
The surface energy oscillates with slab thickness in the same way as the
surface displacements. Figure 4.6 shows the surface energies from table 4.2.
T he convergence is good in most cases, with U8urf calculated with the LDA
converging quickly. The GGA of R PB E does not give good convergence however,
as for the surface displacements. Table 4.4 compares some of my values of Usurf
with those from previous studies. Ram am oorthy et al. [8 8 , 89] found a surface
energy of 0.76Jm - 2 , which is a O .lJm - 2 difference with my result; however,
my calculations are probably more accurate, as I used a larger energy cutoff,
double the number of k-points, and relaxed the atoms much more (down to a
force of 1 0 - 4 e V / A as opposed to 1 0 - 1 e V r/ A ) . More recent studies put Usurf a t
0.90Jm ~2, which agrees reasonably with my results.
4.3.3.3
Com parison w ith isostructural tin dioxide
Tin dioxide has the same rutile structure as titanium dioxide, but does not
share some of its more peculiar properties, for example the very high dielectric
constant in the c-direction. To offer a counterpoint to the complicated situation
with TiC>2 , I performed some slab calculations on Sn (>2 using PAW and the
PBE functional. The results are shown in figure 4.7. Firstly, we can see th at
the displacements vary with less m agnitude for tin dioxide, and th a t for the six-
CHAPTER 4. SURFACE FREE ENERGY OF TJO2(U0)
32
0.8
o —o US-LDA
□—□ US-PW91
0 —0 PAW-LDA
0 — 0 PAW-PW91
0 —0 PAW-PBE
0—0 PAW-RPBE
^ 0.6
2 0.4
4
5
6
7
8
9
10
11
number of layers n
Figure 4.6: Surface energy with increasing n for different functionals.
<
I Ti(6) titanium dioxide
I T.(5)
►CKIP)
t 0(8)
I Ti(6) bn dioxide
I Ti(5)
►0(1P)
I 0(B)
I
- 0.1
-0.24i-
_L
5
J
6
7
number of layers n
Figure 4.7: Comparison of surface displacements for titanium dioxide and tin
dioxide; using PAW and PBE.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
93
fold titaniums and in -plane oxygens there is a large difference in the amount
of relaxation. For the latter, as well as for the Ti(5), the oscillations are out of
phase between the two systems. Overall, the displacements in the tin dioxide
system appear to converge more rapidly, with even a five layer slab appearing
to reproduce the surface structure well.
4.3.4
Lattice dynamics of bulk T i02
A previous study [46] has implied that the slow convergence of (110) surface
properties could be related to the vibrational modes of the surface, and perhaps
bulk titanium dioxide. They showed that the energy surface was very flat with
vertical displacement of the bridging oxygen, and th at this was related to a pos­
sible soft vibrational mode. As noted in the literature review, TiC>2 has a very
soft transverse optical mode A 2u, which is peculiar to this particular material
and responsible for some of its properties, in particular its high dielectric con­
stant. Although this particular mode cannot be explicitly responsible for the
behaviour of the BOs (because it involves motion only in the [001] direction),
similar such modes and combination of modes at high temperatures could affect
the outcomes of experiments, and thus explain discrepancy between experiment
and theory.
For example, from table 4.3, the displacement of the BO is seen to be 0.10A
in a LEED study at 140K, and —0.27A in a SXRD study at room temperature.
The latter result disagrees with almost all OK theoretical investigations, which
find that the BO moves away from the surface slightly.
I have performed a number of static minimum energy calculations to find
the frequency of a number of modes at the gamma-point (i.e. center of the first
Brillouin zone). Figure 4.8 shows the three modes I have looked at. In the A ig
mode, each oxygen moves in the x — y plane, away from and then toward their
CH APTER 4. SURFACE TREE E N E R G Y O P TIO2(U 0 )
B mode
94
A, mode
mode
Figure 4.8: The A\g, B2g and A2u vibrational modes of bulk titanium dioxide.
For the first two modes, only the oxygens move, and this movement is in the
plane containing them, perpendicular to the z-direction. In the A2xi mode,
all movement is in the z-direction, with the oxygen-titanium movement out of
phase.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U0)
95
opposite. In the B 29 mode, the motion is similar to the A ig mode, except that
when oxygens in the x —y plane at 2 —0 move towards each other, the oxygens
in the plane at z = | c move away from each other. The Am mode involves the
motion of both the titanium and oxygen atoms, the two species moving up and
down the z —axis respectively. Both species are involved, so the magnitude of
their motion is different because the two species have different mass, and there
is a 2:1 ratio of oxygens to titaniums. We can find the relation between the
magnitude of the oxygen displacement do and of the titanium displacement dm
by recognising that the center of mass of the unit cell must be conserved. This
gives us the formula dm = do x 2m o /m m Lattice dynamics are defined by the dynamical matrix, which is the second
derivative of the total energy with respect to the displacements of atoms in the
cell, and in general is a (n x n) matrix, where n is the number of degrees of
freedom, equal to 3xno. of atoms in the primitive unit cell. To calculate the
frequency of a particular mode, we do not need to calculate the full dynamical
matrix; if we know how this mode moves the atoms, we can displace them by
a small amount d and calculate the energy of the unit cell. By repeating this
for different displacements, we can plot a graph of energy vs. d. Provided the
displacements are kept small, the graph should be quadratic, as the motion will
be harmonic. The second derivative of the energy at d = 0 is the spring constant
k of the motion, from which the frequency of the mode can be found:
The mass m and the displacement need to be considered carefully to get the
correct frequency. For the A ig and B 2g modes, four oxygen atoms per unit cell
are involved in the movement, so the mass is equal to 4mo- For the A 2u mode,
both atom species move, but in opposite directions. If we shift our point of
9C
CH ARTER 4. SURFACE FREE E N E R G Y OF TIO2(U 0 )
Method
LDA
PW91
PBE
Neutron scattering
IR & Raman
Ai9mode
613
600
580
610
612
#2 j mode
811
787
783
825
827
A2Umode
156
173
167
Tabic 4.5: Frequencies found for the modes Ais , B 2g, and A2u using different
methods, given in cm-1. Neutron scattering results from ref. [99], IR & Raman
from [21],[85].
view to that of the titanium, this is mathematically equivalent to the titaniums
remaining stationary and the oxygens moving a displacement [do + d n ) —
do { 1 + 2mo / n i T i ) - So if we obtain a graph of energy vs. do , and then scale
the displacement by
(1
I 2m o/m n ), the second derivative at d — 0 will be the
spring constant k' for which the mass ra' = 4mo- Then the frequency w of the
mode will be equal to (fc'/m')1/2.
0.07
0.06
| 0.04
USPP-LDA
PAW - LDA
PAW - PW91
USPP - PW91
PAW-PBE
S 0.03
0.02
0.01
-
0.1
-0.05
0
0.05
0.1
displacement of oxygens (A)
Figure 4.9: Relative total energy vs. do for the A2u mode, using three different
functionals.
Table 4.5 shows the frequencies found for the three modes. In every case,
C H APTER 4. SURFACE FREE E N E R G Y OF TIO-2(110)
107
Figure 4.11: Two slabs of titanium dioxide separated by some vacuum gap. The
two exposed surfaces are the (110) surface. Atj is the distance between the
bridging oxygens in the z-direction, which stand above the plane of titanium
atoms.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
97
LDA outperforms the GGAs, showing at worst a 10% difference with exper­
iment. PW91 and PBE underestimate the frequencies for the A ig and B 2g
modes, and in the A 2u case a it is unreasonable to try and extract a number.
Figure 4.9 shows the plot of energy vs. displacement of the oxygen for the A 2u
mode, for each of the functionals. The LDA curve looks like a good quadratic.
For this mode however, the energy is not seen to be quadratic near the origin
for either of the GGAs. Both curves approaching from large displacement are
seen to fall, then rise, then fall again before reaching zero displacement. In
addition, the PBE curve dips below zero. More than 5 calculations were done
in the region d = [0,0.05], to check the behaviour. The mode is thus seen to
be unstable here, as found previously by Montanari and Harrison [79], and the
GGAs seem to fail to reproduce an interesting and particular property of the
material.
4.3.5
Conclusion
On the basis of all these results, the best functional to use is definitely the
local density approximation. It reproduces the lattice parameters to less than
0.015A, the surface properties converge quicker and more convincingly, and it
reproduces the vibrational modes of the material well. All the GGAs produce
less accurate results across all these areas. It should be noted however that if
I were to introduce some molecules into the system, such as water, then new
tests would need to be done to ascertain the best functional. Some studies
have concluded that the GGAs, in particular PBE and its revised edition, are
better at modeling small molecules and their interactions with surfaces than
LDA (see for example [12]). As regards the use of ultrasoft pseudopotentials or
the projector augmented-wave approach, there is little difference between the
two. The property which displays the biggest difference between the two is the
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
98
surface energy, which based upon comparison with other ab initio calculations,
the PAWs reproduce the best.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
4.4
99
General theory for th e calculation o f surface
free energies
Firstly in this section, I outline a general method for the calculation of surface
free energies by thermodynamic integration. This method can be applied at zero
temperature, where the surface free energy reduces to the surface energy, and
also at higher temperatures. Here, the results can be checked by temperature
integration, which I discuss in section 4.4.2.
4.4.1
Thermodynamic Integration
In section 3.3, I outlined the basic principle of thermodynamic integration for
finding the difference in free energy between two states of a system. I apply this
principle by considering the two states as being two simulation cells illustrated
in figure 4.10. The first cell contains the bulk crystal phase of the material. In
the second cell we have created two surfaces that are separated by some distance
L, the vacuum gap. In doing this, we have increased the length of the cell in the
direction perpendicular to the surfaces by L. We can then link these two cells
by a dimensionless parameter s, which is proportional to L. This parameter is
what we shall integrate over to find the difference in free energy between the
two cells.
4.4.1.1
Zero T e m p e ra tu re
Let us call the three vectors describing a simulation cell a a (a = 1,2,3), with
each vector having components aa\ (A = 1,2,3) . We can cause the cell to
change size and shape by applying a strain to it. This is like some force pushing
or pulling the cell. Mathematically, it will be a 3 x 3 tensor. When we apply
a strain to a cell, we are putting it under pressure. As this pressure is not in
C H A P T E R 4.
S U R F A C E F R E E E N E R G Y O F T I 0 2 (110)
continuous.
reversible path
A
bulk unit cell
100
sA
■
s=0
s=s.
Figure 4.10: Diagram showing the strain placed on the system as a function of
param eter s. A is some small distance relating s to the size of the vacuum gap
L, such th at L = sA.
general isotropic, it must be defined also as a 3 x 3 tensor, called the stress, and
is denoted cr\
To find the stress on the cell, we can exert some small strains e\^ to explore
the space around a a :
O'aX =
^
An infinitesimal change in the strain
^ i t dgf f
(4 -4 )
causes the following change to a'aX:
(4.5)
do-ax —
The components of the stress tensor are then defined to be :
I ( d u \
« J „ o ’
(4.6)
where V is the volume and U is the total energy of the repeating cell. U depends
on the vectors a Q and thus the strain
and the positions of the atoms in
the cell. When we come to calculate the stress tensor, we will relax all the
atoms in the cell to their equilibrium positions. Therefore, this dependence of
U disappears, so th a t the energy is a unique function of E\^ and the derivative
in equation 4.6 becomes an absolute derivative.
Equation 4.6 can then be rearranged to leave us with a formula for the
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
101
infinitesimal change in energy U:
dU — V
ov\d £ \v.
(4.7)
Av
We want to express the energy change in terms of our parameter s, which defines
our path for the thermodynamic integration. The vectors describing the cell can
be expressed as:
“ I- ® A / i ( ^ ) ] ® a r / i (0)•>
® aA (^ ) — ^
(^ * ® )
with eA^(O) = 0. The change in aa\ due to an infinitesimal change in s is:
daQa = ^
(0)ds.
(4.9)
t*
Expressing a aA(0) in terms of aQ^(s) from equation 4.8 and substituting into
equation 4.9 we get:
ddaA = ^ 2 ~foTds
fl
V
+ e (s ))-1 ] ^ ( s)-
(4-10)
Comparing equations 4.5 and 4.10 we infer that the infinitesimal strain is:
< f e ^ = * 2 ^ [ ( I + e(*))-V -
(4-H)
We can now put this into equation 4.7 to give us the variation of U with s:
dU = V (s)ds'52a„ x ^ - [ ( l + e ( s ) ) - V .
Xfiv
(4.12)
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
102
Finally, the change of energy as s goes from 0 to si is:
U(S l) - U ( 0 ) = r V ( S) y '< T „ A ^ [ ( I + e ( s ) ) - V d s Jo
&
ds
4.4.1.2
(4.13)
Control interaction
We need to introduce an additional interaction to ensure that a vacuum gap
opens in a continuous and reversible manner as the cell is strained. This is
because, if one starts with the unstrained simulation cell and strains it, a vacuum
gap will not open: the cell will simply be deformed and the atoms will spread
out to compensate. As the strain is increased there will come a point at which
it will be energetically favorable to form surfaces, and the material will rapidly
break into two, so there will be a discontinuity in the stress. In addition, the
process will not be reversible: the point at which the surfaces spontaneously
form will not be the same as that when the surfaces snap together if we were to
run the process backward.
I introduce then an additional control interaction, called Ucon(s). Its role
is to push the surfaces apart in a continuous and reversible way. It must be
continuous, and satisfy the conditions Ucon(0) = I/con(si) = 0, so th at it does
not contribute to our value of surface free energy.
The total energy is now Utotai = U + Ucon• We need to account for the
interaction’s effect on the stress by adding a term <j^n:
U M - 1/(0) = f " dsV (S) V K x + < x" ) ^ [ ( I + e ( s ) ) - V .
J0
\»V
^
(4.14)
The surface energy is then:
(4, 5,
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(110)
103
where A is the area of each created surface.
4.4.1.3
F inite Tem perature
To extend the above approach to cases of finite temperature, we define the stress
tensor as:
^
4
G
£ L
’
(4l6)
where F is the Helmholtz free energy, and the derivative is taken at constant
temperature T. The change in free energy resulting from a finite strain is then
given by the analogue of equation 4.14:
F(Sl) - F(0) = P dsV(s) £
(<r„x + <?*) ^ [ ( 1 +
(4.17)
Here, {.) denotes that the thermal average of the total stress must be taken. To
find the true surface free energy, the integration in equation 4.17 should be taken
out to infinity. At a large enough value of
however the error is negligible.
The surface free energy is then:
(4, 8 )
4.4.2
Temperature Integration
If we have Faurf at one temperature, we can calculate it at another temperature
by the method of temperature integration.
In the canonical ensemble, free energy is given by F ( N ,V ,T ) = —k T \n Z ,
where k is Boltzmann’s constant, Z is the partition function, V is the volume,
and N is the number of particles. Dividing by T and taking the partial derivative
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2 (110)
104
with respect to T gives:
8_
8T
= g t [ - * ln Z ]
= -Tfji f-fcln J dpdrexp(-/3f/)|
= - | / d p £ir e x p ( - ^ ) A ( - ^ )
Taking this final derivative gives the Helmholtz relation:
d_
dT
(4.19)
where {U) is the thermal average of the total energy. Then, integrating over T :
(4.20)
We have two systems: the simulation cell at zero pressure, s = 0, and the
strained cell at s = s i . Writing the equation 4.20 for both systems, subtracting
and dividing by 2A:
Fm r/(T2)
T2
i2 d
dTr
Feurf(T \) _ J _ [f T‘
T yT'1
j « £ ! . = „ > ( 4 . 2 1 )
T,
2 A Jj Tt
This means that we need only calculate Fsurf at one temperature, and with
an appropriate number of further average energy calculations we can find it at
any temperature (as is done in [45]). In practice, it is unwise to use just one
temperature to find F (T) at a very different temperature, where the behaviour
of the material may differ.
We can simplify equation 4.21 if we assume that the motion of the atoms is
harmonic. This is often a good assumption for some materials away from high
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(UO)
105
temperatures close to the melting point of the material (and its surfaces), and
away from very low temperatures where quantum effects become important.
Since the number of degrees of freedom is the same for the two systems, by the
equipartition theorem, the difference in energy ((Us=8l) —(Us- 0)) is constant
and equal to 2AUsurf, where Usurf is the zero temperature equilibrium surface
energy. Then, equation 4.21 simplifies to:
F,u rf(Tl) = F ,ur/(T i) + ( l - ^ ) (Usurf ~ Faurf(Tl)) .
(4.22)
106
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2 (110)
4.5
Application to the T i02 (110) surface
The previous section outlined a general theory for calculating surface free ener­
gies, which can be applied to any material. In order just to test its validity, I
have chosen titania, and here adapt the theory specifically for this system.
When considering the (110) surface of titanium dioxide, I model a simulation
cell which is orthorhombic, described by the vectors ai = ai, a 2 = 6j, a 3 = c(s)k.
We can express c(s) as:
(4.23)
c(s) = c(0) + s A,
where s is our strain parameter and A is a small increment. Thus the strain
tensor e(s) can be seen to be (from equation 4.8):
e A <t(3 ) — S x z S f t z ^
8.
(4.24)
The volume is given by F (s) = a6[c(0) + sA]. Putting these into equation 4.17,
we get a form for the change in free energy:
ds ab A (ozz + acon)
(4.25)
It just remains to find an appropriate form for Ucon(s), the interaction that
pushes the opposing surfaces apart. It could have many forms, including a single
particle potential that depends on the positions of selected atoms in the surface
region, or a two-particle potential depending on the relative positions of selected
atoms. For this particular surface, see figure 4.11, the simplest way to ensure
that a surface is formed is to include a repulsive force between the bridging
oxygens (BOs) on opposing surfaces.
These BOs stand above the crystal planes on the surface of a slab; ensuring
they do not move towards each other will prevent the movement of the atoms
C H A P T E R 4.
S U R F A C E F R E E E N E R G Y O F T I 0 2(1 10)
107
Figure 4.11: Two slabs of titanium dioxide separated by some vacuum gap. The
two exposed surfaces are the (110) surface. A 77 is the distance between the
bridging oxygens in the z-direction, which stand above the plane of titanium
atoms.
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(110)
108
beneath them and thus maintain the presence of surfaces. The simplest form
for
U con
then is:
Ucon —N b o A ( s )At/,
(4.26)
where N bo is the number of bridging oxygens on each surface in the repeating
unit cell. Ar\ represents the distance in the k direction between BOs on op­
posite sides. A(s) is a force curve that will need to be found with appropriate
calculations.
We can now find acon from:
acon =
_1
V
N bo
V(s)
fd U c o n \
V
fezz )
The derivative of Arj with respect to ezz is At/. This is because the strain is
homogeneous; any vector R in the cell is transformed by R ' = (1 + e)R. And
for the second term:
ds _ ds dc(s) _ c(s)
deZz
dc(s) dezz
A
(4.28)
Using these results in equation 4.27, it can be shown that:
a'.con
_
N b q A ( s )A tj
V (s)
N b o * 9A
abA ^ ds
(4.29)
CHAPTER 4. SURFACE FREE ENERGY OF TIO2(U0)
4.6
109
Results
Here I present my results for the calculation of the surface free energy of titanium
dioxide (110) using density functional theory. To begin, I needed to make some
technical decisions based on my experience of modeling the material, as set
out in section 4.3. Here, I investigated the accuracy of the modeling of the
titania (110) surface, taking into account the pseudopotentials and exchangecorrelation functional used, and considerations such as the slab thickness, and
vacuum gap. I concluded that the choice of pseudopotential had little effect,
but th at the use of generalised gradient approximation over the LDA reproduced
surface properties less accurately. For the forthcoming simulations I therefore
decided to use the LDA, with ultrasoft pseudopotentials.
As regard the slab thickness, I found th at a slab of seven or eight layers
reproduced the surface properties well. These are however 42 and 48 atom
systems, if the smallest possible cell is used, and at finite temperature this
would be doubled (see section 4.2.4). Therefore, I have used only a four layer
slab, to reduce the amount of computer time needed. This is sufficient as I
only aim to demonstrate the method works, so the absolute value of surface free
energy obtained for this particular system is not the main point of interest.
4.6.1
Finding an appropriate form for A(s)
A(s) defines the force that needs to be exerted upon the surface bridging oxygens
to ensure that a vacuum gap is opened continuously as the cell is strained. We
have the requirement that A(0) = A (si) = 0, so th at the interaction is switched
on and then off. To construct a suitable A(s) I did a series of calculations where
all the atoms in the cell were allowed to relax except the BOs, which were
’pinned’, such that the distance between them through the slab was held fixed
at its s = 0 value. Once the other atoms are fully relaxed, there were forces left
CH APTER 4. SURFACE FREE E N E R G Y OF T I 0 2 (110)
110
1.50
—
—
1.25
Force on BOs
Adjusted curve for A(s) _
1.00
0.50
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
vacuum gap (A)
Figure 4.12: Form for A(s). A value of s = 40 corresponds to a vacuum gap of
2A.
on the bridging oxygens, which wanted to bring them closer and close the gap.
Figure 4.12 is a graph of the force on each BO for various vacuum gaps (i.e.
various values of s). This curve does not decay to zero, as even at very large
vacuum gaps the BOs will relax from their bulk term inated positions like the
other surface atoms. Section 4.3.3 showed th a t the BOs like to move away from
the surface by some 0.08A, increasing the distance in the k direction between
themselves and the titanium atoms below. Thus, the curve decays to a force in
the region of 0.6 —0.8eV/A. Figure 4.12 also shows an adjusted curve which
falls to zero at s = 40. My small increment A is equal to 0.05A, so this is
equal to a vacuum gap of 2A. This was used as A(s); it satisfies the endpoint
and continuity criteria, and ensures th a t a vacuum gap opens up in a reversible
manner. To find a real surface energy, the cell must be strained until each
surface no longer interacts with the other. A vacuum gap of 5A converges the
surface energy to within 0.5%, so I have taken a corresponding value of si = 100.
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2 (110)
111
-100
-150,
1.5
0.5
time(ps)
Figure 4.13: The fluctuating stress during a simulation and the running average.
Simulation performed for s = 15.
A(s) drops to zero at s ~ 40, and is zero in the range s = [40,100].
4.6.2
Average stress
I performed calculations at three temperatures: OK, 295K and 1000 K. The static
relaxations at OK were relatively quick, so I was able to simulate at 46 different
s values. However the molecular dynamics simulations took much longer, and so
1 have only 6 points at the higher temperatures, and these simulations lasted for
2 —3ps. The length of the simulation was determined by observing the running
average of the stress, an example of which is shown in figure 4.13, so that the
errors in the average stress were less than 5 kbar.
Figure 4.14 shows a graph of (crtotal) at each temperature. The black line
shows the fit of the average total stress at zero temperature. The curve does
not start at the origin, which we would expect for a system in equilibrium. This
is because there is a large contribution to the stress at s = 0 from the term in
<jcon that depends on the derivative of ^4(s) (equation 4.29). It can be seen from
C H APTER 4. SURFACE TREE E N E R G Y OF TIO2(U 0 )
112
200
— OK
• 295K
• 1000K
100
jo -100
-200
05
1.0
1.5
2 .0
2 .5
3 .0
3 .5
40
4.5
5.0
vacuum gap (A)
Figure 4.14: {(Ttotai) vs- s at each temperature. The lines are fitted curves,
and the symbols are calculated points. 295K results are shown by the circular
symbols, 1000K results by the square symbols, and OK results are shown by the
plain line. The sharp discontinuity at s = 40 is because dA/ds ^ 0 for s < 40,
but dA/ds = 0 for s > 40.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(110)
113
figure 4.12 that this derivative is non-zero both here and at $ = 40, which is the
reason for the discontinuity at this point.
Given our freedom to choose a form for Ucon, I could have chosen a curve for
A(s) which had zero derivative at both these points. In this case, the stress curve
would start from the origin and there would be no discontinuity at s = 40. There
is however, despite the strange nature of the curve, no problem with integrating
the stress if the curve is fit correctly.
To perform this integration I split the curve into appropriate regions and
treated each separately. The first rapidly climbing region (s = [0,13]) was fit
with the fermi function, and the second flat region (s = [13,40]) was fit with a
parabola. The third region (s — [40,100]) was fit with an exponential curve. In
the last region, the stress was assumed to fall off exponentially in the same way
for finite T as for zero T.
The error in a time average can be found using the sample variance and
estimating a correlation length r. r has the meaning that, a system at a time
t+ r wall have forgotten its state at the previous time t. In this case, the error in a
sample is equal to a2r / t 8im, where o \ev is the variance of the sample, and t8im is
the number of timesteps. I calculated r to be approximately 250 timesteps using
blocking averages (see the appendix for an explanation of blocking averages).
The errors found by this method were (for all bar one point) smaller than Skbar,
or 8 eV/A3.
4.6.3
Calculation of Fsurj
To find the surface free energy, I integrated the stress over s from s — 0 to
s = 100, and used equations 4.25 and 4.18. At zero temperature, this is simply
the surface energy, equal to the energy required to create each surface per unit
area. I did not need to go to all this trouble to calculate the surface energy, as
CHAPTER 4. SURFACE FREE E N E R G Y OF T I 0 2 (110)
temperature (K)
Faurf (Jm-2 )
0
0.86
295
0.80
114
1000
0.68
Table 4.6: Surface free energy at the three temperatures investigated,
this can be found simply by:
as was done in section 4.3.3. But it was good to go to this trouble as I can check
that the theory works.
The comparison of surface energy is then as follows:
FBUrf(0K) = U8urf = 0.858Jm~2 calculated by integration over s
U8UTf = 0.859Jm-2 calculated from equation 4.30
The values compare very favourably and show th at the method is accurate. I
found values at finite temperature of Faur/(295K) = 0.80Jm-2 and F*ur/(1000K)
0 .68 Jm~2. Table 4.6 summarises the surface free energies at every tem perature
found by thermodynamic integration.
4.6.4
Temperature integration m ethod
Equation 4.22 gives an equation for the surface free energy at a tem perature T2
as a function of the surface free energy at a temperature 7 \ . So if we have the
surface free energy at one temperature, we can in theory calculate it any other
temperature.
However its foremost use within my work is to check that the thermody­
namic integration method is still accurate away from zero temperature, where
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U0)
115
its accuracy has been proven. So, the following gives the surface free energy
calculated taking firstly, T\ = 295K, and then T\ = 1000K:
Fsurf(Ti = 295K) = 0.80Jm-2 =* Fsurf(T2 = 1000K) = 0.66 Jm “ 2
FsurfiTi = 1000K) = 0 .68 Jm - 2 =* F8urf(T2 - 295K) = 0.81JnT2
Comparing the result for T2 = 1000K, the two integration methods agree well,
with a difference of only 0 .02 Jm - 2 . Comparing the result for T2 = 295K,
the difference is smaller, equal to 0.01 Jm - 2 . These differences are not the
same because of the weighting of the surface free energy by the temperature
ratio T2 /T 1 in equation 4.22. The difference in the energy calculated by the
integration over temperature between 295K and 1000K is constant, but the
actual difference between Fgurf at each tem perature is different, depending on
whether you calculate the 295K result from the 1000K result, or vice versa.
We might have expected the two methods to agree exactly. However there
are a number of sources of error which can explain why they do not. Firstly,
there is the possible error in the calculation of Fsurf at the finite temperatures
with the thermodynamic integration. At zero temperature the method has been
shown to be accurate to within 0.005Jm - 2, but at finite T less points on the
stress curve were calculated, and there is a larger possible error attached to each
one. This is unavoidable, because of the large computational effort required to
perform ab initio finite temperature MD. So, given th at there is an error in
Fsurf (Ti), and this error will be different for different values of Ti, I would
expect there to be some discrepancy between the results.
Secondly, there is an error associated with my adoption of the harmonic
approximation in the temperature integration. There is perhaps a large degree
of anharmonicity in the material, as previous studies have found anharmonic
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(110)
116
0.9
0.8
0.5 -
0.4
500
1000
1500
2000
temperature (K)
Figure 4.15: Projection of the surface free energy up to 2000K using temperature
integration, and assuming the harmonic approximation.
modes of vibration on the surface layers a t room tem perature [46].
I can use my results to estimate the surface free energy at any temperature.
Figure 4.15 shows the projection of Faurf up to a temperature of 2000K, which
is near the melting point of the material. Here, it is estimated to be 0.46Jm - 2 ,
which is nearly half the surface energy. The reliability of the projection at this
kind of temperature is uncertain, as there may be large anharmonicity here.
Also, we are dealing with a surface not a bulk material, so the melting point
will be lower and there may be a large number of anharmonic surface modes
explored at high temperatures.
CHAPTER 4. SURFACE FREE E N E R G Y OF TIO2(U0)
4.7
117
Summary and conclusions
The calculation of surface free energies is non-trivial but important for the
comparison of equilibrium crystal surfaces. I have used a method of thermody­
namic integration to calculate the surface free energy of titanium dioxide ( 110 )
at temperatures of OK, 295K and 1000K using density functional theory. At zero
temperature, where Fsurj — Ugurf is calculated both by thermodynamic inte­
gration and by static relaxation of the surface, the correct value was found to be
0.86Jm-2 , with agreement between the two methods being within O.OOlJm-2 .
At finite temperature, I found values of F8urf{T = 295K) = 0.80Jm-2 and
Fsurf (T = 1000K) = 0.68Jm-2 . This shows th at the surface free energy has a
strong temperature dependence, dropping by 0.18Jm-2 from the surface energy
at 1000K. This implies that there is a large increase in entropy as the tempera­
ture of the system is raised. Previous studies [79, 46] and my own work (section
4.3.4) have shown that the bulk material and its surfaces have a number of soft
vibrational modes. The energy surface is very flat with respect to movement
of the bridging oxygens, for example. If the atoms of the surface layers are
vibrating large amounts at relatively low temperatures, this would increase the
entropy and explain the rapid drop of the surface free energy.
The method of temperature integration was used to check the accuracy of
Fsurf at the finite temperatures studied. They agreed to within 0.02Jm - 2 ,
showing that the thermodynamic integration method is accurate at high tem­
peratures as well as zero temperature. The small error is down to two concerns:
firstly, the error associated with the finite length of the MD simulations per­
formed, which are computationally expensive. This expense is not prohibitive
however, as at zero temperature the static relaxations are quick and many points
on the integration path can be found. This provides a good template for the
curve of stress against strain at finite temperature.
CHAPTER 4. SURFACE FREE E N E R G Y OF T I0 2(UO)
118
The second source of error is anharmonicity in the motion of the atoms.
I have assumed the harmonic approximation in the temperature integration,
which simplifies the relationship of surface free energy to temperature to a linear
dependence. This approximation does not seem unreasonable, as an error of
0.02Jm - 2 is only 3% Faurf at 1000K. Projecting the linear relation between
F8Urf and T up to temperatures of 2000K implies that it will drop to nearly half
its OK value by the melting point of the bulk (2143K). However, the projection
is likely to be unreliable in this temperature region, because the anharmonicity
in the system will increase with increasing tem perature, and the surface will
start to exhibit melting behaviour at some earlier point.
In summary, I have demonstrated that the calculation of surface free energies
with density functional theory is possible, and accurate.
Chapter 5
Desorption of water from
MgO(OOl)
In this chapter, I present and discuss my investigation into the adsorption of
water on the Magnesium Oxide (001) surface, and its desorption. Understanding
the way simple molecules such as water and other larger molecules interact
with crystal surfaces is a basic and important piece of knowledge, and this is
reflected by the large number of studies due to the wide number of applications
to catalysis, crystal growth, corrosion, nanotechnology, biophysics, geology etc.
This particular system has been much studied, both experimentally (by tem­
perature programmed desorption) and theoretically (both ab initio and classi­
cally). Experiments provide useful information about the desorption, such as
rates as a function of temperature. The process of translating the results into a
microscopic understanding of the mechanisms is made much more robust by the
use of computer simulation, and some progress has been made in this regard.
In this chapter I detail my work on the system, using interionic potentials to
model water on the surface from zero to two-thirds coverage. My aim was to use
119
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
120
a number of methods to calculate the desorption rate over a wide tem perature
range; this would take in high temperatures where molecules desorb easily, with
lifetimes on the surface of lOps and less, down to low temperatures of 100—300K,
where a molecule can stay adsorbed for hours. At these lower temperatures
brute force calculation of the desorption rate becomes impossible, due to the
long simulation times necessary. However, it is also at these low temperatures
that TPD experiments are conducted, and therefore theoretical access to this
thermal region is highly desirable. In order to gain this access, I have used the
potential of mean force method at both low and high temperatures; in the case
of the latter, to compare with the brute force method and check its validity.
To begin, I introduce the different methods of investigating desorption. Sec­
tion 5.2 gives all the details of the system, and how the molecular dynamics
was performed. Section 5.3 outlines a general theory of the desorption of any
molecule species from a surface, including the use of the potential of mean force
to this end. Section 5.4 presents results for the adsorption energy of a water
molecule at zero temperature, and the following section 5.5 investigates the mo­
tion and equilibration of molecules on the surface over a range of coverages.
Finally, sections 5.6 and 5.7 contain the theory of desorption as applied to the
H2 0 /Mg 0 (001 ) system for first, the isolated molecule, and second for higher
coverages.
C H A P T E R 5. D E S O R T T IO X OF W A TE R F R O M MGO(OOl)
a) unit cell
b) (001) surface
Figure 5.1: a) Crystal structure of MgO, and b) the (001) surface. Atoms arc
magnesium in red and oxygen in blue.
5.1
In tr o d u c tio n
In this section I firstly introduce tire MgO(OOl) surface, and discuss the previous
work that has been conducted on the adsorption of water molecules. Then I
give a general background of the most common theoretical and experimental
methods for the investigation of desorption of the molecules from the surface,
and other studies which have applied them to this particular system.
5 .1 .1
A d s o r p tio n o f w ater o n M gO (O O l)
Magnesium oxide holds a status as a model oxide for surface modeling, so a
great number of studies have been conducted on this material. It has the rocksalt structure, as shown in figure 5.1, with a nearest cation-anion distance of
a /2 = 2.105A (Wyckoff [106]). The most energetically favourable surface is the
(001) surface.
The adsorption of water on MgO was once considered to be solely molecular,
with dissociation only occurring at steps or defect sites: for example [90, 64]
for theory and [27, 95, 107] for experiment. One of the first theoretical stud­
ies to suggest otherwise was that of Giordano et al. [33, 34], who found using
DFT that on a perfect surface with one monolayer coverage, a mixed state
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
122
with one third of the molecules dissociated was the most stable system. This
was attributed to the interactions between the adsorbed molecules, with some
molecules bonding to each other like a dimer. They concluded th at this dimerisation of adsorbed molecules aided dissociation, considerably increasing the ad­
sorption energy. Odelius [81] came to the same conclusion using ab initio MD,
adding th at at a temperature of 200K some of the dissociated molecules re­
combine. This suggests that temperature has to be taken into account when
comparing with experiment. Other studies have since come to the same conclu­
sions [14, 17, 74].
One of the first experimental studies to suggest dissociation on the surface
was Kim et al. [56, 57], using metastable impact electron spectroscopy (MIES)
and ultraviolet photoelectron spectroscopy (UPS). They performed ab initio
calculations of the density of states (DOS) at one monolayer coverage for three
systems: no molecules dissociated, one third dissociated and one half dissoci­
ated. Comparison of the DOS with th at from the experiment showed th at a
state with one third of molecules dissociated was the most favourable. Further
adsorbed water layers were seen to be solely molecularly adsorbed, and they
were not able to discern whether dissociation still occurred in the first layer
when a number of monolayers were present. Other studies have agreed that
dissociation in the first monolayer is observed [31, 55, 108].
5.1.2
Desorption of water and other m olecules from sur­
faces
The adsorption and desorption of water and other molecules from surfaces is a
much studied process, due to the many applications to catalysis, crystal growth,
corrosion, nanotechnology, biophysics etc. Experiments provide useful informa­
tion about the desorption, such as rates as a function of temperature. The
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
123
process of translating the results into a microscopic understanding of the mech­
anisms is made much more robust by the use of computer simulation. However,
straightforward simulation of desorption from a surface is not practicable at low
temperatures, as desorption events are rare and thus require unaffordably long
simulation times to reproduce. Therefore, more sophisticated techniques are
needed for calculating desorption rates at a variety of temperatures, and within
a reasonable amount of computer time.
Firstly I shall talk about transition state theory, which provides a theoretical
framework for investigating desorption, and some of the techniques used to
reduce the computational time required for its implementation. Then I shall
look at the most common experimental method used for observing desorption.
Then I will summarise the previous work done on the particular system of
interest.
5.1.2.1
Transition state theory
Desorption from a surface is commonly treated theoretically by transition state
theory (TST) [44, 100, 3, 103, 104], which uses statistical mechanics to find the
rate of crossing of a surface separating the initial state from the final state. In
the case of desorption, the initial state is a molecule adsorbed on the surface,
and the final state is the start of its escape into the gas phase. The rate of
crossing of this surface at a given temperature T in the canonical ensemble, can
be expressed as an average of the crossing rate for systems of energy E = [0, oo]
in the microcanonical ensemble with temperature T:
k (T) =
1
f°°
/
dE exp (-J3E) M E ) k (E) = (k(E ))T .
(5.1)
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
124
Zi and fa are the partition function and density of states for the initial state.
Evaluating k(E) leads to the standard TST relation:
k{T) = j K^ - e x p ( - / J A £ ) ,
(5.2)
where Z t s is the partition function for the transition state and AE is the
energy barrier separating the initial and transition states. So the crossing rate
is proportional to the ratio of transition state and initial partition functions,
and a Boltzmann factor.
When TST is applied to the process of desorption, its use alone still requires
the running of very long simulations at low temperatures so that a desorption
event is observed, so some method is usually employed to reduce the compu­
tational time required. Becker and Fichthom [6 , 28] studied the desorption of
single alkane molecules from graphite, using two methods to speed up the run.
Firstly, they identified some transformation which linked the high and low tem­
perature systems, ran the simulations at a high temperature and then scaled
down to their target low temperature. Secondly, they added a compensating po­
tential to interactions between molecules and surface, to reduce the adsorption
energy and allow the molecule to desorb more easily. A similar compensating
potential was used by Grimmelmann et al. [37] to study the desorption of xenon
from platinum.
5.1.2.2
Tem perature program m ed desorption
A very important way of measuring the desorption rate directly is temperature
programmed desorption (TPD). In this method, all the molecules are initially
stably adsorbed on the surface of interest, and then the temperature is slowly
increased. The number of desorption events per unit time can be measured as
the experiment progresses, with the time elapsed related to the tem perature of
CHAPTER 5. DESORPTION OF W ATER FROM MGO(001)
125
the system. In simple cases, there will be one peak where most of the molecules
desorb, and the temperature at which this occurs will tell us the adsorption
energy of the molecule. The time distribution is usually analysed by fitting
to the Polanyi-Wigner equation [58], in which the rate of change of adsorbate
coverage 9 (the number of molecules per adsorption site) is given by:
(5.3)
where / is a frequency prefactor, n is the reaction order and A E the activation
energy for desorption. This equation is obviously related to equation 5.2.
The frequency prefactor is often thought of physically as an “attem pt fre­
quency”, given by the number of times the molecule tries to escape from the
surface per unit time. Therefore, when an estimate is required to analyse TPD
experiments, a common value selected is / = 1013s- 1 , which is in the range
of typical vibrational frequencies in solids. Another interpretation is th at /
is related to the ratio of configurational space available for the adsorbed and
free states of the molecule. This is reflected in the form of equation 5.2, where
the analogy of / is the ratio of the partition functions of desorbed to adsorbed
molecules [100, 109]. In theoretical calculations, / emerges naturally from cal­
culation of the desorption rate, and these values can potentially be used in TPD
analysis.
Ahmed et at. [1] showed using TPD and x-ray photoelectron spectroscopy
that the desorption spectra observed from such experiments is highly sensitive
to the way the surface is prepared. Stirniman et at. and Xu et at. [95, 107]
investigated the desorption of D 2O from MgO(OOl), using equation 5.3. The
former used a prefactor of / = 1013s-1 , finding an adsorption energy of 0.65eV.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
5.1.2.3
126
D esorption o f water from MgO(OOl)
Although there is a great deal of experimental data for the H 2O/M gO( 0 0 1 ) sys­
tem, the methods needed to make comparisons between experiment and theory
have not been fully developed, particularly in the area of desorption. Previ­
ous theoretical work investigating H2 0 /M g 0 (001 ) was done by McCarthy et
al. [75], who modeled the system using potentials fit to post-Hartree-Fock cal­
culations. They used TST coupled to Monte Carlo simulations to find surface
hopping and desorption rates of the isolated molecule. More recently, Alfe and
Gillan [2 ] conducted a first principles investigation into the desorption of an
isolated water molecule from MgO(OOl) using TST arguments and a potential
of mean force (PMF) method, enabling them to calculate the desorption rate
7 down to 100K. They found that they seriously overestimated 7 compared to
experiment [107, 95], and suggested th at this might be due to inaccuracies of
the exchange-correlation functional used.
CH APTER o. D ESO R TTIO N OF WATER FR O M MGO(OOl)
MgO slab
I
repeating unit
in z-direction
•mu
Uni
I
A
Figure 5.2: How the system looks, using periodic boundary conditions. There
are two water molecules in the system (in grey), and the ions of the slab are
magnesium (in red) and oxygen (in blue). The distance between the two surfaces
in the unit cell is not to scale: typically in simulations, a distance of over 15.4.
is used.
5.2
M eth o d o f in v estig a tio n
In a TPD experiment, water is deposited onto one surface of a thin magnesium
oxide film. This film is usually mounted on a substrate of some other material.
There will typically be a few monolayers of water covering the surface. The
system is heated up slowly at a constant rate, and the number of desorptions is
counted with the time elapsed. A typical experiment lasts a number of minutes.
We cannot simulate such a long time in a molecular dynamics simulation (as
the simulation of just a nanosecond takes hours of computer time), so we cannot
conduct a theoretical investigation like we would a TPD experiment. Instead,
we can simulate the interactions of the water molecules with the surface at a
constant temperature. It would be possible to simulate a thin film of magnesium
oxide mounted on a substrate, or kept still by other means, and a number of
molecules interacting with the surface. We could place a virtual wall some
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
128
distance away from and parallel to the surface, to make the system closed.
However, it is more convenient to use periodic boundary conditions (see section
2 .1 .7 . 1 ), so that any molecule which desorbs from the top surface will readsorb
on the bottom surface. The system is illustrated in figure 5.2.
I perform two types of simulations in my work. In the first, all molecules are
free to move, the system is simulated at high temperatures and the molecules
move back and forth between the two surfaces within the repeating unit, and how
often they do this depends on the desorption rate. In the second, the potential
of mean force calculations, all the molecules in the system are adsorbed to just
one of the surfaces.
In the remainder of this section, I outline all the particular details of how
the simulations were performed.
5.2.1
Technical details
MD simulations were performed with the code DL_POLY [30] which can be
used to simulate a wide variety of systems; simple atom mixtures, point ions,
rigid molecules, polymers, macromolecules, metals and covalent systems. In
addition to
p e r fo r m in g
MD, DL_POLY has a zero-kelvin energy minimisation
algorithm, used to relax a system before an MD run. This algorithm searches
for an energy minimum by giving ions a starting velocity, letting them move
until the force exerted on them is opposite to their velocity. At this point all
forces are set to zero, and the process begins again.
The Nose-Hoover therm ostat (see section 3.2.3.1) was used to regulate tem­
perature within the canonical ensemble. When attempting to model the isolated
molecule at low temperatures (< 100 /if), it was found that this therm ostat reg­
ulated the temperature in an unphysical way: T was zero for 90% of the simula­
tion, but would regularly spike to large temperatures for very short periods, so
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
129
th at the average temperature was correct. To get around this problem, I used
the Berendsen thermostat to equilibrate the system, which gave a tem perature
oscillating normally around the desired temperature. Then I would switch to
Nose-Hoover for the production run.
A time-step of 0.5 fs was used, to properly simulate the movement of the
hydrogen atoms. The fastest vibrational mode of the water molecule has a
frequency of 3756cm-1 , corresponding to a time period of ~ 9fe. A timestep of
0.5 fs is then a suitable timestep, as it is close to l/2 0 th of this period.
5 .2 .2
T he p oten tial m odel
In my investigation I have used the rigid ion pair potential whose form was given
in equation 2.37 to describe ion-ion interactions:
vi
eXP (~ r / p i J ) -
-
(5-4)
Values for the parameters describing interactions between Mg and O ions were
taken from Lewis and Catlow [67], who fitted to experiment for a large number
of oxides, successfully applying them to the calculation of perfect lattice, defect
lattice and surface properties. They also included the shell model; however I
did not use the shell model. There will be no external electric field in all my
simulations, and the dispersion term in the potential includes the effect from in­
stantaneous dipole moments. I calculated the adsorption energy of a molecule on
a slab that was described by the Lewis and Catlow model (i) including the shell
model and (ii) not including it. The difference in energy was less than 25meV,
so I am confident that the model does not misrepresent the interaction between
molecule and slab. The short range interactions between the magnesium atoms
were set to zero (to aid fitting by reducing the number of parameters), as was
the Mg-0 dispersion parameter. The other parameters are given in table 5.2.2
CHAPTER 5. DESORPTION OF W ATER FROM MGO(001)
130
(a)
-
0.0
0.0
1
Pi,j (-M
B iyj (eVA°)
0
Mg - Mg
o
parameter
(eV)
22764.4
0.1490
20.37
Mg - O
821.6
0.3242
0 .0
(b)
parameter
(eV)
pi (A)
B i (eVA6)
charge qi
Mgatirf
982.0
0.2358
53.1
1.966
Gsurf
4154.8
0.2294
10.9
-1.966
OW
14378.5
0.2299
6.7
-0.820
HW
30.4
0.2967
0.8
0.410
Table 5.1: Parameters for the rigid ion potential describing interactions between:
(a) Mg and O ions in the bulk or slab; (b) water molecules and water/slab
interactions. For (b), the parameters are subject to the combination rules A i j =
(AiAj)'/*, PiJ = (pi + Pj)/2 and B tJ = (B .B ,)1/ 2.
(a).
For the interactions between water molecules, and of water molecules with
the MgO(OOl) surface, I used parameters from McCarthy et al. [75]. These
were fitted to ab initio Hartree-Fock calculations with correlation corrections
to the energetics. They specifically investigated the structure and energetics of
this system, and found th at little charge redistribution took place upon the ad­
sorption of a molecule, so a classical description of the interaction is reasonable.
The parameters are given in table 5.2.2 (b). Note that the model is such that
the molecules interact only with the surface ions of the slab, Mgsur/ and 0 8urf.
The intramolecular interactions were described separately by harmonic po­
tentials, with equilibrium geometry corresponding to the SPC representation of
water [94]. The spring constants used were k\ = 31.2eV/A 2 for oscillation of
the O-H bond and ka = 5.0eV/rad2 for oscillation of the H-O-H angle. This
model does not allow for dissociation. As noted in section 5.1.1, dissociation is
not the normal method of adsorption to the surface, but is observed in a mixed
state of both physisorbed and chemisorbed molecules at non-zero coverages in
C H APTE RS. DESORPTION OF W ATER FROM MGO(OOl)
131
certain proportions. Therefore we might expect dissociation to occur when we
investigate coverages of more than one molecule. For the present however, I
have assumed that adsorption is only physical, and for what I aim to achieve
this is a reasonable assumption.
5.2.3
The unit cell of the system
The slab will be composed of a number of layers niayere: each layer contains
a plane of ions, with an equal number of Mg and 0 ions. Two planes one on
top of the other will be identical, but displaced by half the lattice constant in
the x and y-directions, so the interior of the slab looks like bulk MgO. The two
surfaces of the slab will be separated by some distance L. W hat we choose for
L and niayer8 will affect the structure of the surface layers and their interaction
with the molecules; we want to have enough layers th at the interior of the slab
behaves like bulk MgO, and a large enough L th at molecules adsorbed on one
surface do not interact with the other surface across the gap. This needs to be
balanced with the associated increased computational cost.
In practice, I calculated the surface formation energy as a function of niayers
at OK using the energy minimisation algorithm (see section 5.2.1). The (001)
surface is very flat, with only very small displacements of surface ions from bulk
positions. I found that a three-layer slab gave a surface formation energy within
0.03eV compared with a five layer slab, and have therefore used a three layer
slab in my calculations. For L, I have used a distance of 25.8 A , at which point
the long range coulomb force has negligible effect. Decreasing L by 10A results
in a change of adsorption energy of less than an meV.
The number of ions in each layer is also a very im portant factor. If the
size of the cell in the x and y-directions is small, the distance between two
equivalent molecules adsorbed on a surface will be small. This will impact on
C H A P T E R S. DESORPTION OF W ATER FROM MGO(OOl)
132
the nature of their motion: molecule movement will be correlated and unphysical
because of the restrictions imposed by the periodic boundary conditions. Again,
any increase in the cell size would need to be weighed against an increase in
computational cost. For my investigation I have used 36 atoms per layer. The
adsorption energy was found to be converged to within lOmeV compared to 100
atoms per layer. This size of cell means th at the distance between a molecule
and its image is 12.657A, which compares very favourably to the second nearest
neighbour distance at 1ML coverage of ~ 4.219A.
The water coverage on the surface is defined in terms of monolayers, such
th at a 1ML system has one molecule adsorbed for every Mg in the surface
layer. With a slab th at has 36 atoms per layer, this corresponds to 18 ad­
sorbed molecules. Lower coverage is defined by fractions: half ML = 9 adsorbed
molecules, third ML = 6 adsorbed molecules etc.
C H A P TE R S. DESORPTION OF W ATER FROM MGO(OOl)
5.3
133
Theory of desorption
The following theory offers a number of methods for calculating the desorption
rate of a molecule from a surface, and is completely general regardless of cover­
age, molecule and surface type. I start by defining the system and its properties
of interest which will be required in the theory.
5.3.1
Definition of the system
As discussed in section 5.2, my system consists of a pair of parallel surfaces
perpendicular to the z-direction, and a number of molecules between them that
can be either adsorbed on a surface, or in the gas phase (see figure 5.2). We
define a molecule to be in the gas phase when its interaction with the two
surfaces is weak, this typically being when it is more than a few A away from
both surfaces.
The molecules adsorbed on a surface will move around it - adsorb at differ­
ent sites, in different patterns, try out different orientations - until they have
equilibrated. At this point, they have forgotten how they first hit the surface.
These translation, combination and orientation processes will occur at different
rates, and the theory to be presented requires that the rate of desorption is
slower than the slowest surface equilibration rate; th at is, th at the molecules
are equilibrated before they have a chance to desorb.
We have v identical molecules in the system, with /i atoms in each molecule.
The position of atom k (k = l.../i) in a molecule j (j ~
is given by Tj,k-
We assume that the position of the molecule can be described by the position
of one of its component atoms, say atom k = 1 . This atom would commonly
be the one with the largest mass, or the one in the centre of the molecule. In
the case of H2 O, the oxygen in the molecule is the natural choice. Then, the
CHAPTER 5. DESORPTION OF W ATER FROM M GO(001)
134
density of water molecules p(r) at a point r is:
p(r) = J 2 (S (r - ry.i)) = v (S (r - n ,i) > ,
(5.5)
i - i
where (...) signifies a thermal average.
We now define the average of p(r) = p(x, y, z) over x and y for a given 2 to
be p(z), and a quantity y(z), which is analogous but dimensionless:
where po is the average number density of molecules in the gas phase (any
molecule having z <
zq
is considered to be adsorbed on the surface) and is
constant. We define the number of molecules adsorbed on a surface per unit
area, d, by:
(5.7)
The following two subsections describe the methods we have used to calculate
the desorption rate for a given number of molecules v in the system.
5.3.2
Formula for the desorption rate 7
If we perform an MD simulation on the system in which no ions are constrained,
the molecules in the system will desorb and adsorb freely. Each of them will
approach a surface, stay on the surface for some amount of time, and then
desorb, to cross back to the other surface. On average, the number of molecules
adsorbed on each surface is constant, and equal. There is a detailed balance
between the number of molecules desorbing from and re-adsorbing to a surface,
so the rate at which molecules cross a plane z = z ‘ in the gas phase in one
direction is equal to the rate at which they cross it in the other direction. This
CHAPTER 5. DESORPTION OF W ATER FROM MGO (001)
crossing rate
k
135
we define as the number of molecules crossing z — z' in the
positive direction per unit area per time.
The number of molecules per unit area in the gas phase with ^-coordinate
between z' and z* + dz' is equal to podz1. Using the Maxwell-Boltzmann distri­
bution, the number of molecules per unit area with a velocity between vz and
vz + dvz that cross the plane with z = z* in a time dt is:
(
m
V /2
(5.8)
where m is the molecular mass. Integrating with respect to vz (when vz is
positive only) and dividing by dt gives the crossing rate
k:
(5.9)
It is possible that a molecule arriving from the gas phase will not stick to
the surface. It may just bounce off, and then cross the plane z — z ' , so that
not all of the rate
k
can be attributed to spontaneously desorbing molecules.
Unless one wishes to assume that all arriving molecules stick to the surface,
k
can be found by multiplying the right hand side of equation 5.9 by the sticking
coefficient 5, with 0 < S < 1 , which can be found from simulations.
The crossing rate
k
is the product of a (the number of molecules adsorbed on
the surface per unit area), and the desorption rate 7 , the number of molecules
leaving the surface per unit time per adsorbed molecule. Using equations 5.7 &
5.9, we have:
(5.10)
This simple equation shows that the key quantity for determining the spon­
taneous desorption rate is the distribution y{z). An example of y(z) is given
in figure 5.3. y(z) can be found at high temperatures by performing a long
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
136
peeks corresponding to time
molecules spend on the surface
3*
position
of the slab
Figure 5.3: An example of the form of y(z).
simulation, assembling a histogram from the ^-coordinates of all the molecules
at each time-step, and then normalising so th a t it is equal to 1 out in the gas
phase. We expect y(z) to have two peaks, corresponding to the time spent by
the molecules on each surface, and be flat between these peaks, corresponding
to the gas phase.
The accuracy of this method can be checked by counting the number of
times, N crose, a molecule crosses a plane z = z' in the gas phase. N croa8 can
then be converted to a desorption rate using the relation
k
= N croasf A tSim,
and equation 5.7, where A is the area of the plane.
5.3.3
Potential of mean force
y(z) is our key to finding the desorption rate. At high temperatures, the desorp­
tion rate is high and there are plenty of desorption events in a relatively short
space of time, and so it can be found by the histogram method. However at low
temperatures we might simulate a few seconds (of simulated time, not computer
time) and still not see a desorption event, so another approach is needed.
First we need to express our density p(r) in terms of the statistical mechanics
of thermal equilibrium. Desorption is not an equilibrium process. However, any
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
137
molecules adsorbed on a surface at a given time may be in thermal equilibrium
with the atoms of the surface, and with each other. In practical terms, this
means that the statistical mechanics we use will be valid, providing the rate of
desorption is slower than the slowest surface equilibration rate. Assuming this
is the case, and using equation 5.5:
p(r) = v f J J d R j
i
j,k
Z - 1t f ( r - r i , i ) e x p ( - ^ t 7 ) ,
(5.11)
where we denote the positions of the N atoms in the solid by R» (* = 1...N), and
let U ({ R i,rjtfc}) be the potential energy of the toted system, and /? = 1 /k s T .
The partition function Z = J n » I I j,k dTj,k exP
Therefore p(z), the
average of p(r) over x and y, is given by:
p(z) = A ~ l v
J
J I d t li d x d y
n
drjjkZ 1J(r - r M ) exp (-/? ?/).
(5.12)
Returning to y(z); we can express y(z) in terms of the potential of mean
force (f>(z), by:
y(z) = exp {-P4>{z)).
(5.13)
To determine a way to compute (f>(z), we rearrange equation 5.13 and differen­
tiate with respect to z:
Tz = - k B T T z lDy(z) =
(5.14)
The derivative d in p{z)jdz can be obtained from equation 5.12, giving:
d(f> _ f U i d R i d x d y H j k drjtkS(r - r i ti)(dU /dz) exp (- 0 U )
dz
J
dRt dx dy J],-,* drj<k6(r - n .i) exp (~j3U)
(5.15)
This is in the form of a thermal average where the system is constrained, so that
CHAPTER 5. DESORPTION OF W ATER FRO M MGO(OOl)
138
the ^-coordinate of atom 1 in molecule 1 has the fixed value 2 . The thermal
average is over d U fd z , so that:
(5.16)
where Fz is the force acting on atom 1 in molecule 1, in the z-direction. The
ergodic hypothesis (see section 3.2.4) states th a t an ensemble average is equiva­
lent to a time average for a sufficiently long simulation, so we can compute this
as a time average using molecular dynamics. The potential of mean force <f>(z)
is then found by integrating {Fz)z;
CH APTER J. D ESO R PTIO N OF WATER FR O M MGO(OOl)
133
Figure 5.4: Configurations in which the water molecule adsorbs to the surface:
a) the energetically favourable flat configuration (FC); b) the perpendicular con­
figuration (PC); c) the bridging configuration (BC), seen at particular heights
of the oxygen atom in PMF calculations. Side and plan views are shown for
each configuration. Red and blue circles show Mg and O atoms respectively.
Numbers 1-4 indicate equivalent desorption subsites.
5 .4
Zero tem p era tu re a d so rp tio n en erg y
To find Eads, I performed a series of calculations with the DL_POLY zero-kelvin
energy minimisation algorithm. The final orientation of the molecule proved to
be highly dependent on its starting position, so that it was possible to identify
two different configurations: the “flat configuration” (FC) where the molecule lies
roughly parallel to the surface and near a surface Mg, with each of its OH bonds
pointing towards a neighboring surface 0; and the “perpendicular configuration”
(PC) where the molecule lies in a plane perpendicular to the surface, with one
OH bond pointing up, and the other lying over a Mg-0 neighbour bond (figure
5.4). When the molecule is in the FC, it lies a distance of ~ 2.lA above the
surface magnesium. I calculated the corresponding adsorption energies to be
Eads = 0.874eV and E ££ = 0.785eV. This is a difference of nearly a tenth of
CHAPTER 5. DESORPTION OF W ATER FRO M MGO(OOl)
140
an eV, showing that the PC is a local minimum, and that there is an energy
barrier between the two states. McCarthy et al [75] find an adsorption energy
of 0.759eV in the FC, but they use a rigid array to represent the slab and a rigid
water molecule. Alfe and Gillan [2], using DFT, find th at the adsorption energy
can vary by a factor of 2 depending on the exchange-correlation functional
employed. They found E ads = 0.46eV using the PBE functional, or E ads =
0 .95 eV using the LDA, in both cases for molecules in the flat configuration.
Stimiman et al. [95] found an adsorption energy of 0,65eV using TPD, although I
would not expect experiment to agree with my result, because the water-surface
interaction potential that I use was fit to ab initio calculations, not experiment.
Each Mg is a possible site for FC or PC adsorption, but in addition there
are four possible subsites (labeled in figure 1 ), which due to the symmetry of
the surface are energetically equivalent.
There is also a third configuration which is not a stable equilibrium orienta­
tion but is important during PMF calculations. When constrained at z — 2.89A,
the molecule can also be found in the tcbridging configuration” (BC), where the
plane of the molecule cuts a surface unit cell diagonally, with OH bonds pointing
at surface oxygens. It is similar to the FC, but with the oxygen pulled upwards.
The energy of this arrangement (at this height) is virtually identical to th at of
the PC. We shall see however that at very low temperatures, the distinction
between these two configurations is very important.
141
C H A P T E R 5. D E S O R P T IO N O F W A T E R F R O M M G O (O O l)
5 .5
E q u ilib r a tio n o f m o le c u le s o n t h e su r fa c e
When laying out the theory in section 5.3, I stated th a t it was only valid pro­
viding th at the desorption rate was much slower than the slowest rate of equili­
bration. This is because the theory presupposes th at the molecules adsorbed on
the surface are in equilibrium with it. Therefore I have investigated the motion
of molecules on the surface in detail.
5 .5 .1
M o t io n o f m o le c u le s o n t h e s u r fa c e
We have seen th at a single water molecule on the (001) surface of magnesium
oxide prefers to lie flat, some 2.1.4 above a surface Mg ion, at zero tem perature.
At non-zero temperatures, the molecule will move around the surface, sampling
lots of different configurations, hopping from site to site, and subsite to subsite.
The same will be true at higher coverages, but with the freedom of a given
molecule now restricted by the positions of the other molecules.
OW positioi
Mg atoms
O atoms
4
2
0
•2
-4
•6
-8
-6
-4
-2
0
2
4
-
-8
-6
-4
■2
o
2
4
Figure 5.5- The trajectory of a single molecule over the surface for lOOps at (a)
500K and (b) 200K. Mg ions are represented by black circles, O ions by black
crosses, and the trajectory by the blue line.
Figure 5.5 shows the trajectory of a single molecule over the surface at two
142
CH APTER 5. D ESO R PTIO N OP W ATER PR O M MGO(OOl)
temperatures for lOOps. At 500K, the molecule makes a total of seven site hops
in this time, easily overcoming the energy barriers on the surface to diffuse
across it quickly. At 200K however, the molecule fails to hop from its initial site
during the time. It makes lots of hops between subsites on this particular Mg
ion, but stays away from a position directly over the Mg.
o
m
t
«
-e
-4
-2
0
2
4
Figure 5.6: The trajectories for (a) sixth ML and (b) third ML for 15ps.
Figure 5.6 shows the trajectory at 400K for 15ps of three molecules and six
molecules, corresponding to sixth ML and third ML. In both cases the molecules
arc still able to move freely across the surface. At all coverages, this movement
has stopped (for simulations of Ins) at a temperature of 200K, as it did for the
isolated molecule.
5 .5 .2
Surface d iffu sion co efficien t
The surface diffusion coefficient D gives a number that describes this type of
behaviour. It uses the fact that, for a particle on a random walk, the average
displacement of the particle after a given time is zero, but its square displace­
ment will always be positive. If we store the positions of the particle, we can
CHAPTER 5. DESORPTION OF W ATER FROM MGO(001)
coverage =
200K
300K
400K
500K
OML
371
127
9
iM L
368
65
13
iM L
364
70
18
143
iM L
467
58
13
Table 5.2: Mean time between site hops Thop at 200 —500K and four coverages.
Thop is given in ps. At 200K, the diffusion coefficient D, which is proportional
to l/T/j 0p, was zero at all coverages.
calculate how this square displacement varies with time. For a system with
many particles, we can average over each molecule, giving us the mean square
displacement, which has the form:
(Ar(*)2) ->4X?£ + C,
(5.18)
where (A r(i)2) = (\r(t — t0) - r(£0)|2^ is the mean square displacement of a
molecule after a time t , and C is a constant. The accuracy can be increased
by averaging over multiple time origins to within a given simulation. Figure 5.7
(a) shows the mean square displacement for the 0.5 ML system at 400K. At the
origin, the graph curves up quadratically, before flattening off to an approximate
straight line. Because we are averaging over multiple time origins, there is less
averaging done at high t, so the line falters and is completely lost by 2ns. D is
obtained by fitting a straight line to the curve at small t.
From D we obtain the mean time between hops Thop from D — cP/Arhop-,
where d = 2.98A is the distance between neighbouring Mg sites. Table 5.2 gives
Thop for a range of temperatures and four coverages, calculated from simulations
lasting 2.5ns. At 200K, the diffusion coefficient was found to be zero at all
coverages. This is because, as noted above, no inter-site hopping was observed.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
144
(a)
(b)
Cl
Figure 5.7: For the half coverage system at 400K: (a) mean square displacement
with time; (b) projection P(t) with time
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
5.5.3
145
Orientation correlation
Another way to characterise the motion of the molecules is to assess the corre­
lation of their orientations. If a molecule is moving rapidly across the surface,
its orientation will be constantly changing. We can assess this change using
the unit vector n which bisects the H-O-H angle of the molecule, if we store
the positions of all the ions in each water molecule. Then we can calculate the
average projection P (t), given by:
P(t) - {hi(t0).ni(t0 + 1))
(5.19)
As for the diffusion coefficient, P(t) is averaged over each molecule in the system
and over multiple time origins. Figure 5.7 (b) shows P(t) at one third ML
coverage at 400K. It must be equal to 1 at t = 0, but it drops to zero after
approximately 50ps, and stays around zero for higher t. The time at which
it can be said P(t) is zero (which I label rrot) is an indication of how fast a
molecule will lose memory of its orientation, and is therefore a measure of how
quickly the molecule equilibrates.
Table 5.3 shows the decay time rrot a t the four temperatures and coverages.
It is difficult to find an exact number at which P(t) drops to zero at non-zero
coverages, presumably because the runs were not long enough. Extracting a
number for Trot from figure 5.7 (b) is obviously not an exact science, but the
figures give an idea of the rate at which a molecule changes configuration.
5.5.4
Randomness
I also wanted to discern whether the direction in which the molecule hopped
from one Mg site to another was random. There are four possible directions the
CHAPTER 5. DESORPTION OF W ATER FRO M M GO(001)
coverage -200
300
400
500
0ML
50
7
2
2
£ML
50
25
5 -1 0
£ML
250 - 500
50
1 0 -1 5
146
±ML
500
50
1 0 -1 5
Table 5.3: Decay time Trot at 200 —500K and four coverages. rrot is given in
ps.
molecule can hop in. The fraction of hops th a t were in the same direction as
the previous hop over a 1.8ns simulation at 400K was found to be 0.25 as we
would expect for a completely random process, which confirms th at memory of
direction is lost.
147
CHAPTER 5. DESORPTION OP WATER PROM MCO(GOl)
60
— 800K
— 1000K
— 1200K
50
40
N
30
20
10
5
15
10
20
z(A)
i
25
Figure 5.8: The distribution y(z) at T = 800, 1000 and 1200K, calculated as a
histogram. The surfaces arc at z = 0 and 25.SA.
5.6
T h e iso la ted m olecule
The simplest system I can investigate is that of an isolated molecule on the sur­
face; the limit of zero coverage. In this case there are effectively no inter-molecule
interactions, so complex behaviour such as molecules clustering together will not
be observed. The molecule will be free to sample all the configurations available
and move unobstructed across the surface.
Below, I use the two methods outlined in section 5.3 to calculate the des­
orption rate of an isolated molecule horn the MgO(OOl) surface.
5.6.1
D ir e c t c a lc u la tio n o f th e d e so r p tio n rate
I performed molecular dynamics at four temperatures from 600K to 1200K, for
a period of 20ns. During this length of simulated time, the number of times the
molecule desorbed from a surface ranged from just
12
at the lowest temperature,
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
temperature (K)
600
800
1000
1200
148
7 (observed, ps 1)
7 (predicted, ps J)
error (%)
6.02 x 1 0 ~ 4
6.68 x 1 0 ~ 4
11.02
1.32 x 10“ *
8.28 x 10 “ *
2.01 x 1 0 - 1
1.67 x 10-*
7.77 x 10-*
1.91 x IQ" 1
-9.61
-6.18
-4.99
Table 5.4: The desorption rate 7 for the isolated molecule on the surface; second
column, as found from the number of crossings observed during the simulation;
third column, found from y(z); fourth column, the error between the two ap­
proaches.
to over 2000 at the highest temperature. At 600K and lower temperatures, the
computational time required to observe a statistically reasonable number of
desorption events is too large for this ‘brute force’ method to be feasible. This
is why we develop other methods, such as the potential of mean force method,
to take us into the regions th at brute force won’t allow.
The ^-coordinate of the water oxygen was recorded a t every time-step, and
from this a histogram was accumulated, and normalised by po to give the func­
tion y{z). Figure 5.8 shows y(z) at 800, 1000 and 1200K. All curves average
at 1 in the centre as required, and there are two peaks corresponding to the
cumulative amount of time spent on each surface by the molecule. The desorp­
tion rate increases with increasing T, so we see th at the peaks are smaller at
high temperatures, indicating that for a large proportion of the simulation the
molecule is traveling from one surface to the other.
Table 5.4 shows the desorption rate 7 calculated from equation 5.10, com­
pared with that found from the number of observed crossings, as explained in
section 5.3.2. The two methods should give the same result in the limit of in­
finite simulation time, but they even agree well at 1000 and 1200K after 20ns
of simulation. By studying sub-sections of the 20ns runs, I saw that the er­
ror between the two methods does decrease with increasing simulation times at
T = 1000 and 1200K. At 600K there are only 6 crossings observed, so that we
149
CHAPTER 5. DESORPTION OF W ATER FRO M MGO(001)
20
<
So
■cg
N
-10
-2 0,
20
time (ps)
60
80
100
Figure 5.9: The z -coordinate of a molecule inbetween the two surfaces, at
1200K. Bottom of slab is at z = 0.
cannot expect the statistics to be accurate.
5.6.2
Sticking coefficient
I also used these MD simulations to estimate the sticking coefficient S, by study­
ing the ^-coordinate of the oxygen atom over a series of simulations. Figure 5.9
shows the 2 -coordinate of the oxygen for a section of the simulation at 1200K.
The molecule is seen to bounce off the surface once in this section, at around
60ps. Over the range of temperatures 800 — 1200K, the molecule was seen to
stick 95 —97% of the time, corresponding to sticking coefficients between 0.95
and 0.97.
This does not affect our comparison in table 5.4, but can affect comparison
with the below potential of mean force results. At a sticking probability of 95%
however we can assume a sticking coefficient of 1, as I only aim to show the
method is feasible.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
temperature (K)
7 (ps *) : LR
100
200
400
600
800
6.68 x 10“ 4
1.67 x 10" 2
1000
1200
7.77 x 10_!i
1.91 x 10" 1
7 (ps-1 ) : PMF
1.13 x 10~a&
8.29 x 10"1B
1.33 x 10"*
8.87 x 10” 4
1.55 x 10" 2
7.62 x 10“ 2
1.92 x 10" 1
150
prefactor / (s 1)
9x
1x
2 x
5x
2 x
9x
101S
1017
1016
1015
1015
1014
Table 5.5: The desorption rate 7 for a molecule on the surface at various temper­
atures. The long run numbers are found by integrating y(z) from the histogram,
and then using equation 6 . The PMF numbers are found using y(z) calculated
from equation 8 .
5.6.3
Potential o f mean force m ethod
1 have attempted to use the PMF method at a large range of temperatures: from
200 to 1200K in intervals of 200 K, from 0 to 100K in intervals of 25K. For each
7\ I performed a set of constrained MD simulations, in which the ^-coordinate
of the water oxygen was constrained at 18 different values. The duration of the
simulation at each 2 -coordinate was at least 200 ps, depending on how long it
took to take an average. The molecule was free to translate in the x —y plane,
and to rotate and vibrate. A Mg atom in the middle layer of the slab was also
constrained in the 2 -direction, so th at the slab would not drift up towards the
molecule.
5.6.3 .1
A verage force {Fz)z
Figure 5.10(a) shows the average force {Fz)z at temperatures from 100 to 1200K,
where 2 is the distance from the middle of the slab.
The running average
of the force with increasing simulation length was observed, so th at at most
temperatures and 2 -values the error in (Fz)z was less than 5meV. However, for
T < 100K, I found th at the sampling of Fz became very slow in the region
2 ~ 5A. The reason is that in this region the molecule switches between the
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
151
BC and PC geometries described in section 5.4. The switching between the two
is rapid at high T, but for T < 100K the molecule tends to become trapped
in one or the other, presumably because of an energy barrier between them.
Although the energies of the two geometries are similar, the mean values of Fz
in the two differ by ~ 0 .2 eV/A, which is a large fraction of {Fz)z itself. The
statistical error on (Fz) thus becomes large in this region. We can see th at (Fz)z
depends strongly on T, and it is striking th at at low T it becomes rather flat
in the region 5 < z < 5.5 for T < 200K, even developing a weak subsidiary
minimum at z ~ 5.5A at 100 K. By studying movies of the simulations, I saw
that the typical orientation of the molecule and average positions of the H atoms
undergo large changes as z varies in this region, and believe th at this is why
(Fz)z depends on z as it does.
5.6.3.2
P o te n tia l o f m ean force <j>(z)
The PMF <f>(z) is now obtained from equation 5.17. The integral was performed
by fitting a cubic spline to (Fz)z. Figure 5.10 (b) shows the PM F at tem­
peratures of 100 to 1200K. Given the strong T- dependence of (Fz)z, it is no
surprise to see th at <p{z) also depends strongly on T. In particular, the value
4>min
— <f>(zmin) at the position zmin of its minimum, plotted as a function of T
in figure 5.11, varies by nearly a factor of three between T = 0 and T = 1200K.
The statistical sampling problems that affect (Fz )z at low T also cause a signifi­
cant uncertainty in 4>m in - This is illustrated in figure 5.11, where two alternative
values of 4>m in are plotted for T
<
100K. These values were obtained by assum­
ing that (Fz)z in the region z ~ 5A is dominated by configurations associated
with either the BC or PC geometries. Which of the alternative low-T curves for
<t>min we should choose becomes clear when we consider consistency with the ad­
sorption energy E adS obtained from our static relaxation calculations in section
5.4. I found that the most stable adsorbed geometry (FC) gave E ^ £ = 0.874eV.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
152
(a)
(b)
Figure 5.10: (a) the average force in the ^-direction on the water oxygen as a
function of temperature, from 1200 K down to 100K. (b) the potentials of mean
force (f>(z) in the same temperature range, z is the distance from the middle of
the slab. The arrows show the direction of increasing temperature.
153
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
-0.3
-0.4
^ -0.5
•—• assuming BC where possbte
»--• assuming PC where possbie
-0.7
-
0.8
-0.9
0
200
400
600
800
1000
1200
T (K)
Figure 5.11: Minimum of
At low temperatures, there are two possible
configurations, as explained in the text.
Now, as T —y 0, the dominant configurations for any given constrained value of
z must be those of lowest energy for that 2 -value. Integration of the force Fz
in those dominant configurations must give the adsorption energy of the most
stable adsorbed geometry 1. From this, we see that the lower of the two <f>min
curves at T < 100K is closer to being correct. Taking this curve, I find that
4>min obtained by integration of the mean force is 0.871eV, which is within 5meV
of Eads- We discuss in detail the physics underlying the strong T-dependence
of (ffmin below.
The PMF has a clear minimum, and at low temperatures this will dominate
the expression for the desorption rate. We can make a harmonic approximation
1A water molecule approaching the surface may travel down one of several unconnected
valleys, leading to different states of adsorption, in this case the molecule being adsorbed in
either the BC or PC. Each state will have a PMF minimum associated with it, <£,. Then the
true PMF minimum is given by (j>min & —/^ T ln
exp(—/?<&)] - At low temperatures the
sum will be dominated by the deepest well, and so we can neglect all the other contributions,
which is why we can ignore the PC curve in figure 5.11.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
154
about this minimum, namely:
<f>{z) ~
_
TTlLiJ
(z
ZTn{n')
(5.20)
where m is the mass of a water molecule, zmin is the position of the minimum
and ur is an effective vibrational frequency th at gives the correct curvature of
the PM F at the minimum. If we substitute this harmonic approximation into
equations 5.10 & 5.13 we get:
7 = ^ e x p / t y ,rm%n
(5.21)
which has the form of an attem pt frequency times a Boltzmann factor. This
shows how important 4>min is in determining the desorption rate.
5.6 .3 .3
y(z) and the desorption rate 7
I then used <j)(z) to calculate y(z) at all the temperatures studied, and from
equation 5.21 obtained the desorption rate 7 . These PMF results for 7 are
compared in table 5.5 with the results of direct calculations at T > 600K. The
two sets of values are in quite close agreement for 800, 1000 and 1200K, and
in somewhat less good agreement at 600K; however, we recall th a t the direct
results at 600K are not expected to be accurate, because of the small number
of observed events.
5.6.3.4
The frequency prefactor /
The results for 7 can also be used to deduce the frequency prefactor of the
Polanyi-Wigner equation (equation 5.3). Desorption of isolated molecules at
CHAPTER 5. DESORPTION OF WATER FROM MGO(001)
155
low coverage is a first-order process. In this n = 1 case;
7 = —Q~xdBfdt = / exp ( - A E j k BT ) .
(5 .22 )
If we set the activation energy A E equal to the zero-temperature adsorption
energy (AE — E ^ a = 0.874eV in the present case), we can then deduce the
value that / must have at each temperature in order to reproduce the calculated
7 . The values of / obtained in this way are also reported in table 5.5.
We note two important features of these results: First, / is much greater
than the value of 1013s - 1 often assumed in the analysis of experimental results;
second, / is not a constant, but varies by a factor of 104 over the range 200 1200K. It is not unexpected that / is greatly enhanced above the typical “at­
tem pt frequency” of ~ 1013s- 1 , as several previous experimental and theoretical
studies on a variety of systems have found prefactors as high as 1018s- 1 [75]. In
the first-principles calculations of the desorption rate of H2 O from MgO(OOl)
by Alfe and Gillan [2 ] the value / = 2.7 x 1015s - 1 was found, and it was pointed
out th at the enhancement of / is closely related to the temperature dependence
of (f>(z).
We saw in section 5.6.3.2 that the desorption rate can be written as an at­
tempt frequency w/27r multiplied by a Boltzmann factor exp
5.21). If we now write
deviation of
<f>m in
proximately / =
<j>m in { T ) ~
—E ads
+
A ^mfn(T), so th at A <f>m i n is the
from its zero-temperature value
(u r /2 ,K )expl3A < f> rnin ( T ) .
(equation
Since
(u j/
—E ad s,
then we have ap­
2 t t ) ~ 1013s- 1 , the factor
exp(/?A0min(T)) describes the enhancement of the frequency prefactor. Note
that if A 4>min (T) were linear in temperature, then the enhancement factor would
be a constant. The strong T-dependence of / therefore arises from the curvature
of A (j>min (T) as a function of T.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
5.6.3.5
156
O rig in o f th e n o n -lin e arity o f <pmin(T)
To see why a strong non-linear dependence of A <pmin(T) on T might be expected,
we return to our definition of p(z) in equation 5.12, and relate this to <p(z)
through equation 5.13. We can simplify the statistical mechanics by ignoring
the degrees of freedom of the MgO substrate and considering only those of
the molecule. We can assume also that the bond-stretching and bond-bending
vibrations are not significantly affected by adsorption of the molecule, so that
we need consider only the three translational degrees of freedom x , y and z
and the three angles 0 , (p and tp describing the orientation of the molecule.
Here, 0 and <p are the polar and azimuthal angles specifying the direction of the
molecular bisector, and ip is the angle describing rotation of the molecule about
the bisector. Then we write the energy of the adsorbed molecule relative to its
energy in the gas phase as:
U (x, y , z, 0, <p, ip) = - E ads + i m u 2 (z - z min)2 + v (x, y, z, 0, <p, ip).
(5.23)
Then we have:
i f
i
rJr
f2n
f 2jt
<P{z) — —fcsTln{— j dxdy-^-^
<20 sin 0 J
d<p
dip
j
exp - 8 ( - E a d , +
( z - Zm in) 2
J
+ V (x,y,Z,9,<j>, ^ ) ) }.
(5-24)
where the x —y integral covers a large surface area A. Now, using our approxi­
mation that (p{z) = —E ads + \m<jj2{z — Zmin)2 + A<pmin{T) we obtain:
&<pmin{T) = - k BT h i j i
J
d x d y
J
<20 sin
oj
exp [ - p v ( x , y :z,0,(p,ip)]}.
dtp
J
dip(5.25)
(5.26)
C H APTE RS. DESORPTION OF W ATER FROM MGO(OOl)
157
This means that A0mmCO could only be linear in T if the argument of the
logarithm were independent of T, which is clearly not the case. In fact, if we go
to low T, the molecule will become confined to small oscillations of x, y, d, <f>and
ip about the most stable geometry. In the classical statistical mechanics we are
using, as T —¥ 0, each of these five variables becomes confined to a region whose
width is proportional to T 1/ 2, so that A <f>(T) becomes proportional to —T ln T .
This means that (j>(T) approaches —E a^s with an infinite slope as T —> 0, which
is consistent with the curvature of <f>min(T') seen in figure 5.11.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
5.7
158
Higher Coverages
Considering a single water molecule on the surface led to some interesting re­
sults, showing that the frequency prefactor which is commonly assumed to be
a constant does in fact depend on temperature. The next natural step was to
increase the number of water molecules in the system. Although the isolated
molecule case is interesting, looking at higher coverages is more relevant to com­
parison with experiment, and more complex behaviour will be observed due to
intermolecular interactions.
The methods outlined in section 5.3 to calculate the desorption rate apply at
all coverages and shall be used again here. I have looked at systems containing 3,
6 , 9,12 and 18 molecules. When all molecules are adsorbed on one surface, this
corresponds to a sixth, a third, a half, two-thirds and one monolayer coverage. I
start with a discussion of the sticking coefficient, which is more important than
for the isolated molecule.
5.7.1
Sticking coefficient
When there was just one molecule in the system, I found that the probability
of the molecule bouncing off the surface upon approach was less than 0.05.
With more molecules in the system, I would expect this probability to increase,
due to the incoming molecule interacting with a molecule already adsorbed on
the surface. It is highly probable that at any moment in time, an adsorbed
molecule will have one O-H bond pointing in the air, and also th at an incoming
molecule will have an O-H bond pointing downwards. Therefore if an incoming
molecule approaches a region of the surface already occupied by a molecule, it
will be pushed back. Another factor is th at the incoming molecule may ’run
into 1 another molecule that has just desorbed, which will also push it back into
the gas phase.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(001)
coverage =
400K
800K
1200K
£ML
0.975
0.900
0.575
IML
0.725
0.625
0.425
±ML
0.650
0.325
0.250
159
|M L
0.0
0.0
Table 5.6: Sticking coefficient S for three temperatures at different coverages.
S corresponds to the probability th at a molecule approaching a surface will not
bounce off.
S was calculated by equilibrating a number of molecules on one surface at a
given temperature. Another molecule was then placed in the gas phase, given a
random velocity (appropriate to the tem perature) towards the surface in ques­
tion, and its behaviour observed. This procedure was repeated 40 times. In
the case of half coverage for example, 8 molecules were placed on the surface,
with another molecule in the gas phase; I have labelled this case half coverage.
Table 5.6 shows S at three temperatures and four coverages. At all tempera­
tures, S decreases with increasing coverage, as expected. S also decreases with
increasing temperature, due to the larger velocities of the molecules on average.
5.7.2
Direct calculation o f th e desorption rate
I performed molecular dynamics again at 600, 800,1000 and 1200K, for a period
of 4ns. y(z) was found in the same way as for the isolated molecule. Figure 5.12
shows y(z) for systems containing 1-18 molecules for 1200K. Only one of the
two peaks (that you see for example in figure 5.8) is shown. This clearly shows
that, as the number of molecules in the system is increased, the peak of y(z)
decreases in size, and hence the desorption rate increases. The figure also shows
that y(z) dips below 1 at a distance of around 4 —5A from the surface at high
coverages. Based on my observations of the behaviour of incoming molecules
when finding the sticking coefficient, I believe this to be because of molecules
which do not stick to the surface. My y(z) only takes account of the position of
the oxygen in each water molecule, and not of the position of the hydrogens. As
CH APTER 5. D ESO R PTIO N OF WATER FRO M MGO(OOl)
1G0
1mol
Figure 5.12: Detail of y(z) at all coverages at a temperature of 1200K. z is the
distance from the surface.
I have already mentioned, there is a high a probability that a molecule adsorbed
on the surface will have an O-H bond pointing at the gas phase, so this region of
4 - oA is where the hydrogens of adsorbed and incoming molecules will interact
unfavourably, resulting in the reverse of the latter. Therefore this region acts
as a sort of ‘no-go zone’, and makes a qualitative difference to y(z) due to the
low sticking coefficient at this temperature. At lower temperatures the effect is
smaller due to an increased 5.
Figure 5.13 (a) shows the desorption rate over all coverages, before we include
the reduction due to molecules bouncing off the surface instead of sticking to
it. 7 increases with coverage and temperature. As for the isolated molecule, I
compared these 7 with the desorption rate found from the number of crossings.
The number of crossings do not distinguish between a desorbing molecule and
one that has bounced off, so this comparison is viable. Not including the 600K
results, the two agreed within 13% for all coverages.
C H A P T E R S . D ESORTTIOH OF W ATER FRO M M G 0 (001)
1G1
(a)
•—• 1mol
•—• 3moIs
*—- 6mols
•—>9mols
•—* 12mols
•—• Ittmois
S
f?
i
g0)
O
0.1
(b)
5
I
0 DS
1200
Figure 5.13: Desorption rate 7 found for 1 —18 molecule systems: (a) assuming
a sticking coefficient 5 = 1 , and (b) using the values of 5 reported in table 5.6.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
temperature (K)
600
800
1000
1200
3 mols
1.49
1.39
1.17
0.94
6 mols
2.98
2.74
2.22
1.75
9 mols
4.46
4.00
3.12
2.49
12 mols
5.92
5.11
3.99
3.11
162
18 mols
8.64
6.88
5.31
4.30
Table 5.7: Actual average number of molecules adsorbed on each surface a t each
temperature, for a certain number of molecules in the system.
Figure 5.13 (b) shows 7 with the effect of S included. This graph is the same
as in figure 5.13 (a), but with each data point scaled by the appropriate sticking
coefficient. At temperatures where S was not calculated, I have interpolated
between the nearest higher and lower temperatures. This shows th at the true
desorption rate in fact decreases with increasing coverage at high temperatures
(T > 1000K). The inset figure however shows that, at 600K, this is reversed. At
this temperature, 7 increases with increasing coverage. Whether this is a trend
that continues to lower temperatures I cannot tell from these results.
The coverage of molecules on each surface is not easily defined by the number
of molecules in the system. During a long MD simulation, there are at any
moment a number of molecules on one surface, a number on the other, and
a number in between. We define the average number of molecules on each
surface a through equation 5.7, and table shows a as a function of temperature
and coverage. At 600K, where there were few desorption events, the average
number of molecules on each surface (equal to a x A) is very close to half the
number of molecules in the system. As the temperature is increased, there are
more desorption events and therefore a greater proportion of a molecule’s time
is spent inbet ween the surfaces, and the average coverage decreases.
5.7.3
P otential of mean force m ethod
I have used the PM F method for higher coverages in the same way as for the
isolated molecule. For example, in the half coverage case, there are 8 molecules
CHAPTER 5. DESORPTION OF W ATER FROM M GO(001)
163
adsorbed on the surface, and one more molecule which is constrained at various
heights above the surface. There is however a complication in taking the method
to higher coverages: at high temperatures, the unconstrained molecules will
desorb from the surface, and therefore the coverage will not be constant, not
for the length of the simulation or over all the simulations. To get around this
problem, I wrote some extra code into DL_POLY th at reversed the velocity
of any unconstrained molecule attem pting to desorb (i.e. with a height above
the surface above a certain value, equal to 3.8A). This kind of intervention was
necessary at high temperatures (800K), and high coverages (§ML at 400 and
800K).
I studied the system at 300, 400 and 800K, at coverages of a sixth to two
thirds of a monolayer. For a given tem perature and coverage, I performed
simulations at 16 different z-heights for 1 —4ns each. The length of simulation
needed to be longer than for one molecule, as the surface equilibration processes
happen at a slower rate, and therefore it takes longer for the system to accurately
sample the available phase space. This is also why it is not feasible to go to
temperatures lower than 300K.
5.7.3 .1
Average force (Fz)z
Figure 5.14(a) shows the average force at 800K for four coverages. The error in
(Fz)z was kept to within 5meV/ A in the same way as for the isolated molecule,
by monitoring the running averages. The force has the same basic shape as for
the isolated molecule, with a second subsidiary minimum again observed at 300
and 400K.
The figure does show something new however; a maximum around z = 6.75A,
where the force becomes positive and then dips below zero again. This means
that, in this region, the force acting on the molecule would be such as to push
it away, further into the gas phase. At low temperatures, this feature is still
CHAPTER 5. DESO RPTIO N OF WATER r R O M MG0(001)
1C4
j5 -0.1
- 0.2
(b)
— 1mol
— 3mols
01
Z(A)
Figure 5.14: (a) the average force on the constrained molecule at 800K, for zero,
sixth, third and half coverage, and (b) the equivalent potentials of mean force,
z is the distance from the middle of the slab.
1G5
CH APTER 5. D ESO R PTIO N OF WATER FRO M MGO(OOl)
i
'-------- 1
<
1
<
1
1-------- 1-------- 1-------- 1
-0.3
••
•
-0.4
1
■ ■ 300K
♦ ♦400K
• • 8<XK
r_
>
}«
■
- 0.6
0
2
4
6
8
10
12
number of molecules
Figure 5.15: Minimum of
as a function of coverage at 300, 400 and 800K.
observed to a lesser degree. The origin of this force may come from the kind of
H-H interactions considered when exam ining th e sticking coefficient' molecules
adsorbed on the surface are not arranged favorably to accept another molecule
onto the surface.
5.7.3.2
Potential of mean force <£(*)
Figure 5.14(b) shows the potential of mean force at 800K at four coverages. We
can see that the small maximum in {Fz)z has affected the form of 4>{z), with a
clear maximum at 6.4. We saw for the isolated molecule that the minimum of
4>(z) is what determines the desorption rate. Figure 5.15 shows this quantity at
all temperatures and coverages studied. It shows that our conclusion about how
the desorption rate varies with coverage in section 5.7.2 - that
7
increases with
increasing coverage - only applies at high temperatures. According to figure
5.15, at SOCK, the desorption rate will be at a minimum for a coverage of 5 or
6
molecules.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
166
(a)_____________________________________________________
_
coverage -OML
§ML
iM L
3ML
300K
1.41 x 10“ 9 2.19 x 10“ 10 1.42 x 10“ 10 5.10 x 10~ i0 3.49 x 10“ **
7.92 x 10- 7
2.24 x 10"°
400K
4.03 x 10- 5
1.33 x 10" 6 7.15 x 1(T 7
2.92
x
10
~
2
800K
1.55 x 10~ 2 1.92 x i c r 2
7.47 x 10_i
(b)
OML
£ML
coverage ±ML
iM L
|M L
300K
1.1 X 1017 6 .8 x 10 ie 2.5 x 10l7 1.7 x 10Itf
6.8 x 1017
400K
1.4 x 1017 7.4 x 10ie 8.2 x 10 ie 2.3 x 1017 4.2 x 1018
800K
5.0 x 1015 6.1 x 10 1& 9.4 x 1015 2.4 x 101S
Table 5.8: (a) The desorption rate (in ps x) and (b) the frequency prefactor /
(in s-1) at all temperatures and coverages studied.
5.7.3.3
y(z) a n d th e d e s o rp tio n r a te 7
From the PMF, y(z) was found from equation 5.13 and the desorption rate
from equation 5.10, as for the single molecule. As implied by figure 5.15, at
low temperatures 7 initially decreases when more molecules are added to the
system, but then starts to increase again at around a third coverage. The
frequency prefactor varies in the same way, first decreasing with coverage then
increasing, as shown in table 5.8.
For our theory to be valid, we need the desorption rate to be much slower
than all surface equilibration rates, as calculated in section 5.5. Comparison
between the two shows th at 7 is always at least 100 times slower than either
the hopping rate or reorientation rate. For example, at half coverage at 400K,
the hopping rate is 1,7 x 10-2 , compared with a desorption rate of 2.24 x 10-6 .
At 800K, results for 7 as a function of coverage can be compared with those
from direct calculation in section 5.7.2, and figure 5.16 shows this comparison.
The agreement is not as good as for the isolated molecule, with a factor of two
difference at large coverages. This error could be explained by the lack of a gas
phase in the PMF calculations. For very low coverages, the density of molecules
in the gas phase po is very low, at any temperature. As you put more molecules
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
167
in the system however, po will increase and become more important. At high
temperatures, we have already seen that there are a large number of molecules
in the gas phase at high coverages. For example, at 1200K for a nine molecule
system, there are on average 2.5 molecules on each surface, and 4 molecules in
the gas phase.
At 300 —400K, the density po is very small, and therefore the PM F calcula­
tions at these temperatures should not be affected. But at 800K, it is possible
that the lack of a gas phase will cause an error in 7 .
5.7.3.4
Critical tem perature
At coverages of larger than one molecule, intermolecular interactions will become
important in the system. For example, if molecules begin to cluster on the
surface, so that their positions are no longer random but correlated, there will
be some reduction in free energy associated with this behaviour, and so there
will be a larger barrier to desorption. We may expect to find some critical
temperature Tc at which clustering on the surface starts to occur.
The desorption rate is calculated in my method with the assumption that
there are a number of molecules in the gas phase, and a number adsorbed on
the surface, with the desorption rate equal to the adsorption rate. Therefore 7
is determined by the relationship between po and <r, given by equations 5.9 and
5.10:
(5.27)
I have found the phase diagram of po vs. a. I could do this using the above
equation 5.27, but it will be more convenient to rearrange equation 5.7 to get
the density:
(5.28)
Using the integrals of y{z) which were evaluated to find the desorption rate
CH APTER 5. D ESO R PTIO N O r WATER FROM M G 0(001)
008
£
ICS
direct calculation
PMF method
0 04
0.02
average number of molecules adsorbed
Figure 5.16: Desorption rate at 800K, found by direct calculation and the PMF
method.
£C
3
|
O
300K exact
300K approximate
400K exact
400K approximate
600K exact
800K approximate
i
Cl
(0
re
at
c
£
1
I
no. of adsorbed molecules
Figure 5.17: Pressure vs. coverage at 300, 400 and 800K. The exact curve finds
the pressure using no approximations. The approximate curves assume that
the potential of mean force is harmonic around its minimum (please refer to
the text). Units of pressure arc arbitrary, and results at both temperatures arc
placed on the same graph only for convenience.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
169
within PMF, I have calculated po- Figure 5.17 shows density as a function of
coverage at 300,400 and 800K. The units of pressure are arbitrary, and the graph
does not represent the true relation between the pressure at each temperature,
rather the results have been scaled simply to put them on the same graph.
The density is approximately 1000 times larger at 400K compared with 300K,
and approximately 10000 times larger at 800K compared with 400K. The figure
shows two curves for each temperature. In the exact case, the density is found
from equation 5.28 making no approximations. In the approximate case, I have
substituted y(z) with its potential of mean force definition (equation 5.13),
and substituted <f>(z) with the harmonic approximation used previously for the
isolated molecule (equation 5.20):
rZo
pa = a / I
dz exp [-fi<f>(z)\
<f>min{cr,T) + ~muj2 (z - z,'min.
/ ^®Xp [ ^ ^ mjn((T, T)] £
dz exp
<xexp \p4>min{v-, T)] x constant
■mm
(5.29)
So here we have the density as a function of the minimum of <f>(z). The figure
shows that this is a good approximation at all coverages.
The reason why I want to make this approximation is because I want to find
the critical temperature. Tc will be somewhere in the range [300K,800K], so
if I can find some way of modeling the behaviour between the temperatures I
have studied, I can find 7 and po at any temperature and coverage. I cannot
model the integral of y(z) in equation 5.28 with temperature, or the desorption
rate, because they scale exponentially with temperature. The minimum of the
potential of mean force 4>min, however, as shown in figure 5.15, will be relatively
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
170
o
Figure 5.18: Representation of the curve of po vs. a at 300K.
easy to model with temperature. This is where the approximation becomes
useful.
Figure 5.17 shows three different types of curve for po vs. <x. At 800K, the
curve is always rising, and the gradient is always increasing. At 400K, the curve
always rises but the gradient dips at coverages of 1 —4 molecules. At 300K,
the curve has a local m inim um and maximum. The gradient d of the curve
defines two important points in the behaviour of the system. The first is the
point at which 7 ceases to decrease when molecules are added to the system. At
800K, 7 always increases with a, but at 400K, at first 7 decreases then increases.
Through equation 5.27, we can see that the temperature at which this behaviour
changes is given by the condition d(d + Sd) > d(d).
The most important point is the critical point, where we have the condition
that, at one point on the curve only, d = 0. This point occurs somewhere
between 300 and 400K. We can decipher its meaning by considering the meaning
of the 300K curve, represented in figure 5.18. Suppose we have a surface in a
box, with vacuum above it. Then we start to add molecules into the vacuum,
letting the system equilibrate after each addition. As the number of molecules
in the gas phase increases, the coverage on the surface increases. At the point
where there is density p'0 and coverage <t', a change will occur in the system. A
CH APTER 5. D ESO K TTIO N OF WATER FR O M MGO(OOl)
171
number of molecules from the gas phase will adsorb to the surface, increasing
the coverage to a" and correspondingly decreasing the pressure from the gas
phase to pg. Interactions between molecules - i.e. clustering of molecules on
the surface - are responsible for this behaviour. Above the critical temperature
where the local maxima and minima disappear, no clustering will be observed
on the surface.
350K
385K
2e-10
num ber of adsorbed m olecules
Figure 5.19: Density vs. coverage in the temperature range 350-385K, in inter­
vals of 5K.
To find the critical temperature Tc, I went back to the minimum of the PMF.
To find 4>mi„ = 4>min (d,T) between the temperatures 300 and 400K, where Tc
is to be found, I fit a polynomial to the curves of <j>mi„ against coverage at both
temperatures, not including the coverage of 12molecules (two-thirds coverage).
This gave four parameters n. b. r and d at two temperatures. I then assumed
that each parameter had a linear dependence on temperature, find constructed
figure 5.19, which shows p0 vs. d in intervals of 5K from 350-385K. I find
Tc1
—
—
I/I UXv.
C H A P TE R S. DESORPTION OF W ATER FROM MGO(OOl)
172
Ferry et al [26] investigated the structure of adsorbed water molecules on
the surface using LEED and HAS (Helium atom scattering). They found that
upon the heating of a water layer, a phase transition from a 2 D solid to a
2D gas occurred at “ 210K. The solid phase corresponds to a large interaction
between adsorbed molecules, i.e. clustering; the gas phase to no clustering.
My value of the critical temperature is then over 150K larger than experiment.
This implies that the interaction between water molecules is too strong. This
kind of overbinding and overstructuring is common in ab initio representations
of water (see for example [39, 25]), and the potential I use has been fit to ab
initio calculations, so this is perhaps the cause of the lack of agreement with
experiment. The use of an interionic potential fit to experiment would perhaps
give a more accurate critical temperature.
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
5.8
173
Summary and conclusions
The adsorption and desorption of molecules on and from surfaces is a very
important area of science, and as a result a large number of experimental and
theoretical studies have been conducted upon it. However the link between the
two has not been fully realised, due to differences in the way an experiment
and a simulation are conducted. In order to contribute to the understanding of
desorption processes, I have calculated the desorption rate of a water molecule
from the MgO(OOl) surface at a number of coverages, using a general method
that can be applied to any gas/surface system. A classical interionic potential
was used, so that all statistical averages could be well converged.
Two methods for calculating the desorption rate 7 were used; one being based
on the direct counting of desorption events, and the other on the calculation
of the potential of mean force (PMF). The first method can only be used at
high temperatures (T > 600K), where the 7 is of the order of 10- 2ps- 1 , but the
second method can also be used in the experimentally interesting region a t lower
T where 7 is of order Is- 1 . The two methods give results for 7 at coverages of
0 —| ML th at are in close agreement in the high-T region where both methods
can be used, which gives confidence in the correctness of the PMF method.
One of the most important conclusions to come out of the calculations is
about the nature of the frequency prefactor / in the Polanyi-Wigner formula
commonly used to analyse experimental measurements of the desorption rate.
It is often assumed th at / is of order 1013s - 1 and that it is independent of T,
yet I have found that for the present system neither assumption is correct. I
find that / increases from 1014 to 1018s ” 1 as T decreases from 1200 to 200K
for the isolated molecule, and th at it has a similar dependence and magnitude
at higher coverages. The reason for this strong dependence on temperature
lies in a consideration of equation 5.30, which gives the frequency prefactor as
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
174
proportional to the ratio of partition functions in the transition state, to the
adsorbed state:
7 = / ex p ( 0E ada) =
1 Z rr q
exP ( - P E *ds) ,
(5.30)
Whilst Z t s will change little with temperature, Z ad will be strongly temperature
dependent. At low temperatures, the movement of the adsorbed molecules is
severely restricted, with the diffusion rate and reorientation rate dropping to
zero at 200K. This means that adsorbed molecules explore less configurational
space, and Zad is decreased from its higher temperature value, leading to an
increase in / . This explanation for the enhancement of the prefactor is not
new: for example, in a detailed experimental study of desorption of a series of
alkane molecules from MgO(OOl), it was argued th at confinement of rotational
degrees of freedom is mainly responsible for the observed variation of f from
~ 1013 for methane to ~ 1019s -1 for decane [97]. In MD simulations of alkane
desorption from A u (lll), Fichthom and Miron [28] found qualitatively similar
results. W hat is new is the explicit dependence of / on temperature, which I
don’t believe has received much attention thus far. Even though the interaction
model used is highly simplified, the results for / are instructive, because they
suggest that the strong T-dependence of / might also be found in real systems.
Alfe and Gillan [2] investigated the isolated molecule on the MgO(OOl) sur­
face with density functional theory, using the same methods as used here. They
also found enhancement of the prefactor - 1015s -1 - but did not calculate / as
a function of temperature. Their results suggest some dependence on temper­
ature, but much weaker than found in my results, perhaps meaning that the
molecular degrees of freedom are less strongly confined in the DFT calculations.
Where Alfe and Gillan were only able to simulate the isolated molecule in
contact with the surface, because of the high computational cost of DFT MD,
CHAPTER 5. DESORPTION OF W ATER FROM MGO(OOl)
175
I was able to simulate higher coverages. I found that the desorption rate was
dependent on coverage, with 7 at half coverage being half th at at zero coverage
at 800K, once the sticking coefficient is taken into account. The prefactor was
also dependent on coverage, and was dependent on temperature in the same
way as the isolated molecule.
The presence of a number of molecules on the surface means th at intermolecular interactions will be important. At high temperatures, these interactions
have little effect on the desorption, because of the large energy of the molecules.
But at low temperatures, clustering of molecules on the surface becomes en­
ergetically favourable. I found a critical tem perature of 375K from the phase
diagram of density in the gas phase (i.e. pressure) vs. surface coverage, above
which clustering will not be observed.
Chapter 6
Discussion
6.1
Surface free energy calculations
The free energy of surfaces is a fundamental quantity, but there is no one
straightforward method to calculate it. In chapter four I proposed a method
of thermodynamic integration, and demonstrated its use on the titanium diox­
ide (110) surface. Thermodynamic integration requires the performance of a
number of long molecular dynamics simulations to get just one value of the
surface free energy. In my calculations, I performed six MD simulations at each
temperature studied, each of length 2 —3ps.
I did this using density functional theory, so in total the twelve MD simula­
tions performed entailed a large computational expense. Each simulation took
between 6 and 12 hours on sixteen processors, using an SGI Altix, amounting to
100 —200 computer hours in total. Although this is a large effort, a number of
points can be made about the suitability for the method using DFT. Firstly, at
zero temperature, calculations along the integration curve are cheap, so the way
that the stress varies with the applied strain can be pinned down quite finely:
176
CHAPTER 6. DISCUSSION
177
see figure 4.14 for the form of this curve. This provides a good template for the
curve of stress against strain at higher temperatures too: the figure shows that
there is no change in the basic features of the curve even at a temperature of
1000K. Of course, I cannot say the same for any system other than titanium
dioxide (110), but I do not see any reason why this could not be a general obser­
vation. Therefore, with a knowledge of the stress-strain curve at OK, a number
of special values of strain s could be picked, which would give the shape of the
curve with minimal computational effort.
Secondly, although a large amount of computer time is needed, the resulting
value of surface free energy calculated will be accurate. The method finds the
surface energy at OK to within 0.001 Jm - 2 . At finite temperatures, temperature
integration shows th at values of the surface free energy are also accurate. The
error has increased to around 0 .02 Jm - 2 , but this may be explained by defi­
ciencies in the temperature integration method, and in any case is still quite
accurate considering all the time averages which are required in the method.
Thirdly, although perhaps not completely accurate, the method of temper­
ature integration will give us a projection for the surface free energy a t every
temperature with the calculation of Fsurf at only one temperature. This po­
tentially reduces the work required down to 4 —8 long MD simulations, and
this to find Fsurf as a function of T. The accuracy of such an approach would
obviously depend on the particular system - the presence of soft vibration modes
and anharmonicities for example - and perhaps what value of T you choose to
simulate. Projections of Fsurf over large temperature ranges are likely to be
less accurate than those over small ranges around your chosen temperature.
The projection of Faurf for TiO 2 (110 ) implies that we may expect the surface
free energy to drop to half its OK value before it melts. Even if this is unreliable,
the drop from 0 — 1000K is large at 0.18Jm- 2 . This implies that there is a
CHAPTER 6. DISCUSSION
178
large increase in entropy as the temperature is raised- This is compatible with
previous studies and my own work which have found soft vibrational modes and
flat energy surfaces.
Early in chapter four I modeled the (110) surface with static calculations,
for different exchange-correlation functionals, and converging surface properties
with respect to slab thickness etc. I firstly found that the three generalised
gradient approximations I used ([83, 82, 43]) reproduced all properties of the
material studied less well than did the local density approximation ([84]), with
PBE and RPBE being particularly bad. I also found th at surface properties
oscillated strongly with the number of layers composing the material (a slab
in periodic boundary conditions). I concluded th at a slab with seven or eight
layers converges the surface properties within reasonable limits, when using the
LDA. However, in my surface free energy calculations I used a slab of only four
layers. Using an eight layer slab would quadruple the length of the simulations
required, which would have required too many computer hours. My results
for F8urj therefore will not be benchmark calculations. For the application
of the method to other materials however, a large number of layers may not
necessarily need to be used. Magnesium oxide (001) for example, is commonly
modeled using only a three layer slab. This is a result of the simplicity of the
rocksalt crystal structure in comparison with rutile.
6.2
D esorption calculations
The second part of my thesis looked at the desorption of water molecules from
the magnesium oxide (001) surface. Temperature programmed desorption ex­
periments can find rates of desorption as a function of temperature, but com­
parison of these rates with those calculated by theory is often inappropriate
because of coverage considerations. TPD experiments have one or two mono­
CHAPTER 6. DISCUSSION
179
layers of molecules adsorbed on the surface, whereas the majority of theoretical
studies consider desorption in the limit of zero coverage. In my work I aimed
to simulate desorption at large coverages, up to two-thirds of a monolayer, us­
ing classical potentials to try and reduce statistical errors. Although the one
molecule system has been studied in first principles calculations ([2]), ab initio
is currently too expensive to use for large coverages.
Firstly however, I studied the single molecule on the surface. This itself
threw up some interesting conclusions. The frequency prefactor (as described by
the Polanyi-Wigner equation) is typically assumed to be of the order of 1013s- 1 .
I however found it to be enhanced above this by a factor of 100 at 1000K.
Previous studies have found similar enhancements, but I furthermore found
that the prefactor depends strongly on tem perature, which has not previously
been reported. This dependence increases it to 1018s - 1 at a temperature of
200K. I have shown th at the cause of this strong variation is the decrease in
configurational space explored as the temperature of the system is decreased.
An indication of this was given by consideration of diffusion of molecules. A
single molecule on the surface will hop to a nearby adsorption site on average
every 9ps at 500K, but by 200K one hop is not observed in Ins. The first
principles calculations of Alfe and Gillan [2] did perhaps find some temperature
dependence, but it was not as pronounced as I found. This may mean that
the method of classical potentials used may overstate the variation of / with
temperature, as after all the interaction between molecule and surface has been
highly simplified.
My work on the desorption from higher coverages also finds that the prefactor
varies with temperature, and in fact varies with coverage too. This mirrors
the way the desorption rate 7 varies with coverage. 7 is largely affected by
the sticking coefficient S at high coverages, which defines what percentage of
CHAPTER 6. DISCUSSION
180
molecules which approach the surface will stick to it, as opposed to bouncing
off. For the isolated molecule, S is over 0.95 a t all temperatures studied, but
for higher coverages S is found to decrease with both increasing coverage and
temperature, decreasing from 0.65 at 400K for half coverage to 0.25 at 1200K.
At low temperatures, clustering will occur on the surface, and I found the
critical temperature below which clustering will occur. Although comparison
with experiment shows that I overestimate Tc by 165K, my results show it is
feasible to obtain the phase diagram of pressure, temperature and coverage.
There is a huge amount of computational work required in investigating
desorption at higher coverages, compared to the limit of zero coverage. Firstly,
there are obviously more ions in the system, and therefore more calculations
of the force need to be performed. Secondly and most importantly, the time
needed to reduce the statistical errors to acceptable levels is massively increased.
For example, for the calculation of the average force a t one point on the PMF
curve, 200ps fire needed for the isolated molecule and 4ns are needed for twothirds coverage. The latter requires 16 hours of computer time on one processor
of a Pentium PC. Then consider th at over ten of these simulations need to be
performed to find the desorption rate at just one tem perature and coverage.
The implication of this for performing higher coverage simulations with density
functional theory is th at it is currently out of reach. A large increase in computer
power will be required before it is feasible to calculate the desorption rate in
this manner at coverages similar to those found in TPD experiments.
Therefore, the use of classical potentials to explore these type of processes
is valuable, even if one cannot hope to come up with definitive quantitative
answers to these questions. It has already posed a question - about the tem­
perature dependence of the frequency prefactor - and offered a quantity in the
critical temperature where comparison with experiment is straightforward and
CHAPTER 6. DISCUSSION
181
independent of many technical particulars and assumptions.
If I wanted to try and get better agreement with experiment, it would be
worthwhile to repeat the investigation using an interionic potential fit to experi­
mental data, as opposed to ab initio calculations. This would become important
at high coverages and low temperatures, which is the region of interest, where
intermolecular interactions become important. Clustering on the surface not
only determines how adsorbed molecules arrange themselves, but also how large
a free energy barrier there is to desorbing from the surface. A molecule will
desorb less easily if it has become ’friendly’ with its neighbouring molecules. A
potential fit to experiment would hopefully model the water-water interactions
accurately.
Appendix A
Ewald summation
Ewald summation is a technique used to compute the contribution to the to­
tal energy of a system from Coulombic interactions when periodic boundary
condition s are used. This contribution is given by:
(A.1)
t= l
where the potential 0 (r») is:
(A.2 )
The prime on the summation denotes that the sum is over all periodic images
n and particles j , except for j = i when n = 0. The 1 /r dependency makes the
sum is only conditionally convergent, meaning that the evaluation of the sum
depends on what order the summations axe performed.
The method proposed by Ewald [22, 16] deals with this problem by screening
each of the point charges with a gaussian charge distribution of opposite sign.
182
183
APPENDIX A. EWALD SUM M ATION
This involves defining the potential in the following way:
4>coul = <t>Ewald = <i>R + <i>F - 0S-
(A.3)
<f)R is the potential felt at a point r* due to a set of screened charges, and
this can be calculated by direct summation, as the electrostatic potential from
a screened charge rapidly drops off with distance. We have added a set of
gaussian charge distributions to the potential, so we must now add a term 4>F ,
that is the potential from the same distributions, with opposite sign. <f>F is a
smoothly varying periodic function, so it can be represented by a Fourier series,
which also converges quickly. The final term <j>s is the interaction between ion
i and its compensating charge cloud. It is necessary to subtract this, because
we had to include the contribution to 0 F (r») from ion i in order for <j>F to be
periodic. Each part of the energy is given below; for a derivation of these results,
the reader is directed to Frenkel and Smit [32].
4>n is a real space sum, which converges rapidly, and its contribution to the
energy is given by:
=
(A.4)
2U
r»
where e r fc is the complementary error function:
pO O
e r fc (a r) =
I
Jx
exp (—t2) dt.
(A.5)
4>F is a sum in Fourier space, and its contribution given by:
UF = j
^ l ^ l 2exp ( - fc2/4a),
(A.6 )
APPENDIX A. EWALD SUM M ATION
where the fourier transform of the density is:
1 N
= v X ) z*exp
i=l
•
The contribution from the self-interaction term (j>s
Appendix B
Estimating errors using block
averages
The ergodic hypothesis states th at an ensemble average (i.e. a thermal average)
is equivalent to an average over time t as t —» oo (see section 3.2.4). When we
perform molecular dynamics, we have to simulate for some finite t, so there will
be an error in our value of the average of an observable A. We can estimate this
error by the method of block averages.
At each timestep t», we record a value of the observable Ai. Our distribution
{Ai} is our statistical sample, and it will be close to a gaussian distribution. If
all the elements Ai were independent, then the error in the mean A would be:
error =
N
(B.l)
where N is the number of elements in the sample, and cr2(A) is the sample
variance, given by:
(B.2)
185
APPENDIX B. ESTIM ATING ERR O RS USING BLOCK AVERAGES
186
However the Ai are obviously not independent, as where we find a particle at
a given time will have a very strong dependence on where it was the timestep
before. Each particle will only be able to move a certain small distance during
each timestep. Despite this, there will be a value of x, say, for which
will
not be correlated to A{. If we have this correlation length Xi the error in the
mean is given by:
error =
cr2(A)
.
(B.3)
x *.
To find Xi we can break the length of the simulation into a number of equally
sized blocks. There will beNb blocks, each containing
elements, such that
Nt,nb = N . We can calculate the average of A within each block, denoted A b.
We now have a new sample population of 7V&block averages. The mean of this
new distribution must be equal to A , and the variance must be:
,
Nb
**(.#) = n ' E W - t f i=l
(B-4>
The correlation length can then be found by:
X = lim n °
n b—nx>
a {A)
So x
is
(B-5)
foundbybreaking the sample into blocks for different values of n&,
plotting the right hand side of equation B.5.
Bibliography
[1] S. Ahmed, S.S. Perry, and O. El-Bjeirami. J. Phys. Chem. B, 104:3343,
2000.
[2] D. Alfe and M.J. Gillan. J. Phys. Condens. Matter, 18:L451-L457, 2006.
[3] O.L. Anderson and K. Zou. J. Phys. Chem. Ref. Data, 19:69, 1990.
[4] A.V. Bandura, D.G. Sykes, V. Shapovalov, T.N. Troung, J.D. Kubicki,
and R.A. Evarestov. J. Phys. Chem. B, 108:7844-7853, 2004.
[5] S.P. Bates, G. Kresse, and M.J, Gillan. Surface Science, 385:386-394,
1997.
[6 ] K.E. Becker and K.A. Fichthorn. J. Chem. Phys., 125:184706, 2006.
[7] H.J.C. Berendsen, J.P.M. Postma, W.F. Gunsteren, A. DiNola, and J.R.
Haak. J. Chem. Phys., 81:3684, 1984.
[8 ] P.E. Blochl. Phys. Rev. B, 50:17953, 1994.
[9] T. Bredow, L. Giordano, F. Cinquini, and G. Pacchioni. Phys. Rev. B,
70:035419, 2004.
[10] T. Bredow and K. Jug. Surf. Sci., 327:398-408, 1995.
187
BIBLIO G RAPH Y
188
[11] I.M. Brookes, C.A. Muryn, and G. Thornton. Phys. Rev. Letts, 87:266103,
2001 .
[12] M. Casarin, C. Maccato, and A. Vittadini. Appl. Surf. Sci., 142:196-199,
1999.
[13] G. Charlton, P.B. Howes, C.L. Nicklin, P. Steadman, J.S.G. Taylor, C.A.
Muryn, S.P. Haarte, J. Mercer, R. M cGrath, D. Norman, T.S. Turner,
and G. Thornton. Phys. Rev. Letts., 78:495-498, 1997.
[14] J-H Cho, J.M. Park, and K.S. Kim. Phys. Rev. B , 62:9981, 2000.
[15] R.L. Davidchack and B.B. Laird. Phys. Rev. Letts., 85:4751-4754, 2000.
[16] S.W. de Leeuw, J.W . Perram, and E.R. Smith. Proc. R. Soc. London A ,
373:27,1980.
[17] L. Delle Site, A. Alavi, and R.M. Lynden-Bell. J. Chem. Phys., 113:3344,
2000.
[18] U. Diebold. Surf. Sci. Reports, 48:53-229, 2002.
[19] P.A.M. Dirac. Proc. Cambridge Philos. Soc., 26:376, 1930.
[20] P.K.L. Drude. The theory of optics. Longman, London, 1933.
[21] D.M. Eagles. Phys. Chem. Solids, 27:1243, 1964.
[22] P.P. Ewald. Ann. Phys., 64:253, 1921.
[23] A. Fahmi and C. Minot. Surf. Sci., 304:343-359, 1994.
[24] E. Fermi. Rend. Accad. Lincei, 6:602, 1927.
[25] M.V. Fernandez-Serra, G. Ferlat, and E. Artacho. Molecular Simulation,
31:361-366, 2005.
BIBLIO G RAPH Y
189
[26] D. Ferry, A. Glebov, V. Senz, J. Suzanne, J.P. Toennies, and H. Weiss. J.
Chem. Phys., 105:1697, 1996.
[27] D. Ferry, S. Picaud, P.N.M. Hoang, C. Giradet, L. Giordano, B. Demirdjian, and J. Suzanne. Surf. Set., 409:101-116, 1998.
[28] K. Fichthom and R.A. Miron. Phys. Rev. Letts., 89:196103, 2002.
[29] S.M. Foiles. Phys. Rev. B, 49:14930-14938,1994.
[30] T.R. Forester and W. Smith, http://www.cse.scitech.ac.uk/ccg/software.
[31] M. Foster, D. Passno, and J. Rudberg. J. Vac. Sci. Technol. A, 22:1640,
2004.
[32] D. Frenkel and B. Smit. Understanding molecular simulation: from algo­
rithms to applications, chapter 7. Academic Press, 1996.
[33] L. Giordano, J. Goniakowski, and J. Suzanne. Phys. Rev. Letts., 81:1271,
1998.
[34] L. Giordano, J. Goniakowski, and J. Suzanne. Phys. Rev. B, 62:15406,
2000.
[35] K.M. Glassford and J.R. Chelikowsky. Phys. Rev. B , 46:1284-1298,1991.
[36] J. Goniakowski and M.J. Gillan. Surf. Sci., 350:145-158, 1996.
[37] E.K. Grimmelmann, J.C. Tully, and E. Helfand. J. Chem. Phys., 74:53005310, 1981.
[38] G. Grochola, S.P. Russo, I.K. Snook, and I. Yarovsky. J. Chem. Phys.,
117:7676-7684, 2002.
[39] J.C. Grossman, E. Schwegler, E.W. Draeger, F. Gygi, and G. Galli. J.
Chem. Phys., 120:300-311, 2003.
BIBLIOG RAPH Y
[40] W.F.van. Gunsteren, X. Daura, and A.E. Mark.
190
Helv. Chim. Acta,
85:3113-3129, 2002.
[41] D.R. Hamann, M. Schluter, and C. Chiang. Phys. Rev. Letts., 43:1494,
1979.
[42] K.J. Hameeuw, G. Cantele, D. Ninno, F. Trani, and G. Iadonsi. J. Chem.
Phys., 124:024708, 2006.
[43] B. Hammer, L.B. Hansen, and J.K. Norskov. Phys. Rev. B, 59:7413-7421,
1998.
[44] P. Hanggi, P. Talkner, and M. Borkovec. Rev. Mod. Phys., 62:251, 1990.
[45] U. Hansen and P. Vogl. Phys. Rev. B, 60:5055-5064, 1999.
[46] N.M. Harrison, X.-G. Wang, J. Muscat, and M. Scheffler. Faraday Dis­
cuss., 114:305-312,1999.
[47] M.A. Henderson. Surf. Sci., 355:151-166, 1996.
[48] M.A. Henderson. Langmuir, 12:5093-5098,1996.
[49] C, Herring. Phys. Rev., 82:87, 1951.
[50] R.W. Hockney and J.W. Eastwood. Computer simulations using particles.
McGraw-Hill, 1981.
[51] P. Hohenberg and W. Kohn. Phys. Rev., 136.B864, 1964.
[52] W.G. Hoover. Phys. Rev. A, 31:1695-1697,1985.
[53] W.G. Hoover. Phys. Rev. A, 34:2499-2500, 1986.
[54] M.B. Hugenschmidt, L. Gamble, and C.T. Campbell. Surf. Sci., 302:329340, 1994.
BIBLIO G RAPH Y
191
[55] M.A. Johnson, E.V. Stefanovich, and T.N. Truong. J. Phys. Chem. B,
103:3391, 1999.
[56] Y.D. Kim, R.M. Lynden-Bell, A. Alavi, J. Stulz, and D.W. Goodman.
Chem. Phys. Letts., 352:318-322, 2002.
[57] Y.D. Kim, J. Stultz, and D.W. Goodman. J. Phys. Chem. B, 106:15151517, 2002.
[58] D.A. King. Surf. Sci., 47:384, 1975.
[59] W. Kohn and L.J. Sham. Phys. Rev., 140:A l 133, 1965.
[60] A. Komherr, D. Vogtenhuber, M. Ruckenbauer, R. Podloucky, and G. Zifferer. J. Chem. Phys., 121:3722-3726, 2004.
[61] G. Kresse and J. Furthmuller. http://cm s.m pi.univie.ac.at/vasp.
[62] G. Kresse and D. Joubert. Phys. Rev. B, 59:1758, 1999.
[63] W. Langel. Surf. Sci., 496:141-150, 2 0 0 2 .
[64] W. Langel and M. Parrinello. J. Chem. Phys., 103:3240, 1995.
[65] M. Lazzeri and S.de. Gironcoli. Phys. Rev. Letts., 81:2096-2099, 1998.
[66 ] C. Lee, P. Ghosez, and X. Gonze. Phys. Rev. B, pages 13379-13387,1994.
[67] G.V. Lewis and C.R.A. Catlow. J. Phys. C: Solid State Phys., 18:11491161, 1985.
[68 ] P.J.D. Lindan, N.M. Harrison, and M.J. Gillan. Phys. Rev. Letts., 80:762765, 1998.
[69] P.J.D. Lindan, N.M. Harrison, M.J. Gillan, and J.A. White. Phys. Rev.
B, 55:15919-15927, 1997.
BIBLIOG RAPH Y
192
[70] P.J.D. Lindan, N.M. Harrison, J.M. Holender, and M.J. Gillan. Chem.
Phys. Letts., 261:246-252, 1996.
[71] P.J.D. Lindan and C. Zhang. Phys. Rev. B, 72:075439, 2005.
[72] R. Lindsay, A. Wander, A. Ernst, B. Montanari, G. Thornton, and N.M.
Harrison. Phys. Rev. Letts., 94:246102, 2005.
[73] F. London. Z. Phys., 63:245, 1930.
[74] R.M. Lynden-Bell, L. Delle Site, and A. Alavi. Surf. Sci. Letts., 496:L1L6 , 2002.
[75] M.L. McCarthy, G.K. Schenter, C.A. Scamehom, and J.B. Nicholas. J.
Phys. Chem., 100:13989-16995,1996.
[76] J. Mei and J.W. Davenport. Phys. Rev. B, 46:21-25,1992.
[77] H.J. Monkhorst and J.D. Pack. Phys. Rev. B, 13:5188-5192,1976.
[78] H.J. Monkhorst and J.D. Pack. Phys. Rev. B, 13:5188, 1976.
[79] B. Montanari and N.M. Harrison. Chem. Phys. Letts., 364:528-534,2002.
[80] S. Nose. J. Chem. Phys., 81:511-519, 1984.
[81] M. Odelius. Phys. Rev. Letts., 82:3919, 1999.
[82] J.P. Perdew, K. Burke, and M. Ernzerhof.
ACS-Symposium Series,
629:453, 1996.
[83] J.P. Perdew and Y. Wang. Phys. Rev. B, 45:13244, 1992.
[84] J.P. Perdew and A. Zunger. Phys. Rev. B, 23:5048-5079, 1981.
[85] S.P.S. Porto, P.A. Fleury, and T.C. Damen. Phys. Rev., 154:522-526,
1967.
BIBLIO G RAPH Y
193
[86 ] M. Predota, A.V. Bandura, P.T. Cummings, J.D. Kubicki, D.J.
Wesolowski, A.A. Chialvo, and N.L. Machesky.
J. Phys. Chem. B,
108:12049-12060, 2004.
[87] M. Predota, Z. Zhang, P. Fenter, D.J. Wesolowski, and P.T. Cummings.
J. Phys. Chem. B, 108:12061-12072, 2004.
[88 ] M. Ramamoorthy, R.D. King-Smith, and D. Vanderbilt. Phys. Rev. B,
49:7709-7715, 1994.
[89] M. Ramamoorthy, D. Vanderbilt, and R.D. King-Smith. Phys. Rev. B,
49:16721-16727,1994.
[90] C.A. Scamehom, A.C. Hess, and M.I. McCarthy. J. Chem. Phys., 99:2786,
1993.
[91] Schaub.R., P. Thostrup, N. Lopez, E. Laegsgaard, I. Stensgaard, J.K.
Norskov, and F. Besenbacher. Phys. Rev. Letts, 87:266104, 2001.
[92] R. Sikora. J. Phys. Chem. Solids, 66:1069-1073, 2005.
[93] E.V. Stefanovich and T.N. Truong. Chem. Phys. Letts, 299:623-629,1999.
[94] F.H. Stillinger. Science, 209:451, 1980.
[95] M.J. Stirniman, C. Huang, R. Scott Smith, S.A. Joyce, and B.D. Kay. J.
Chem. Phys., 105:1295, 1996.
[96] V. Swamy, J. Muscat, J.D. Gale, and N.M. Harrison. Surface Science,
504:115-124, 2002.
[97] S.L. Tait, Z. Dohnalek, C.T. Campbell, and B.D. Kay. J. Chem. Phys.,
122:164708, 2005.
[98] L.H. Thomas. Proc. Cambridge Philos. Soc., 23:542, 1927.
BIBLIO G RAPH Y
194
[99] J.G. Traylor, H.G. Smith, R.M. Nicklow, and M.K. Wilkinson. Phys. Rev.
B, 3:3457-3472, 1971.
[100] D.G. Truhlar and B.C. Garrett. Acc. Chem. Res., 13:440, 1980.
[101] D. Vanderbilt. Phys. Rev. B, 41:7892-7895,1990.
[102] L. Verlet. Phys. Rev., 159:98, 1967.
[103] A.F. Voter. Phys. Rev. Letts., 78:3908, 1997.
[104] A.F. Voter. J. Chem. Phys., 106:4665-4677, 1997.
[105] G. Wulff. Z. Krystallog. Mineral, 34:449, 1901.
[106] R.W.G. Wyckoff. Crystal Structures. New York: Wiley, 1963.
[107] C. Xu and D.W. Goodman. Chem. Phys. Letts., 265:341-346,1997.
[108] Y. Yu, Q. Guo, S. Liu, E. Wang, and P.J. Moller. Phys. Rev. B, 68:115414,
2003.
[109] A. Zangwill. Physics at Surfaces. Cambridge University Press, 1988.
[110] C. Zhang and P.J.D. Lindan. J. Chem. Phys., 118:4620-4630, 2003.
[111] C. Zhang and P.J.D. Lindan. J. Chem. Phys., 119:9183-9190, 2003.

Similar documents

×

Report this document