text: 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 188.8.131.52 The local density approximation (L D A )............. 31 184.108.40.206 Generalised gradient approximations (GGAs) . . 31 Practical application of D F T ................................................ 32 2.1.7 22.214.171.124 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 18.104.22.168 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 184.108.40.206 Empirical parametrisation ................................... 42 18.104.22.168 Ab initio param etrisatio n ....................................... 42 184.108.40.206 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 188.8.131.52 Initial p o s it i o n s ...................................................... 49 220.127.116.11 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 126.96.36.199 Verlet a lg o rith m ...................................................... 53 220.127.116.11 Velocity Verlet a lg o rith m ...................................... 53 22.214.171.124 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 184.108.40.206 Nose th e r m o s ta t...................................................... 58 22.214.171.124 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 18.104.22.168 4.2 4.3 4.4 B u lk ............................................................................. 71 188.8.131.52 The (110) s u r f a c e .................................................... 72 126.96.36.199 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 22.214.171.124 ...with respect to vacuum gap L ........................... 87 22.214.171.124 ...with respect to slab thickness n 89 18.104.22.168 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 18.104.22.168 Zero T e m p e ra tu re .................................................... 99 188.8.131.52 Control interaction...................................................... 102 126.96.36.199 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 188.8.131.52 Transition state t h e o r y .............................................123 188.8.131.52 Temperature programmed d e s o rp tio n ....................124 188.8.131.52 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 22.214.171.124 Average force (Fz)z ..................................................150 18.104.22.168 Potential of mean force<f>(z).................................... 151 220.127.116.11 y{z) and the desorption rate 7 ............................... 154 22.214.171.124 The frequency prefactor / ......................................... 154 22.214.171.124 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 184.108.40.206 Average force (Fz)z ...................................................163 5 .7.3.2 Potential of mean force <f>(z)...................................165 184.108.40.206 y(z) and the desorption rate 7 ................................166 126.96.36.199 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. 220.127.116.11 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 22.214.171.124 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. 188.8.131.52 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 184.108.40.206, 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 188.8.131.52 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. 188.8.131.52 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 22.214.171.124 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. 22.214.171.124 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 18.104.22.168), 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 188.8.131.52 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. 188.8.131.52 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. 18.104.22.168 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. 126.96.36.199 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. 184.108.40.206 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) 184.108.40.206 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. 22.214.171.124 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. 18.104.22.168 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 188.8.131.52), 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 184.108.40.206 ...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. 220.127.116.11 ...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. 184.108.40.206 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. 18.104.22.168 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 126.96.36.199 (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. 18.104.22.168 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. 22.214.171.124 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. 22.214.171.124 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) 188.8.131.52 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 188.8.131.52) 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. 126.96.36.199 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. 188.8.131.52 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 126.96.36.199 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) 22.214.171.124 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. 18.104.22.168 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. 184.108.40.206 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 . 22.214.171.124 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.