UFDC Home  myUFDC Home  Help 



Full Text  
THEORETICAL MODELING AND DESIGN OF COMPLEX MATERIALS By MAOHUA DU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003 Copyright 2003 by MaoHua Du This dissertation is dedicated to my late grandfather, MeiLin Du. ACKNOWLEDGMENTS I am truly grateful to many people who contributed to my work at the University of Florida. First of all, I would like to acknowledge my research advisor and committee chair, Professor HaiPing Cheng, for her guidance and encouragement throughout my graduate career. Her openminded attitude toward research creates a free environment that allows me to investigate scientific subjects that truly excite me. Her constant enthusiasm and optimism gave me great encouragement and are greatly appreciated. I would also like to thank the members of my supervisory committee (Professors Jeffrey L. Krause, David H. Reitze, Adrian E. Roitberg, and Samuel B. Trickey) for the time and the insight they provided. My gratitude also goes to current and former members of our group and many of my colleagues. I would especially like to thank Mr. LinLin Wang, who has been my officemate for the last few years. I enjoyed many conversations with him about research and life. My appreciation also goes to Mr. Chun Zhang, Ms. Zhihong Chen, Mr. Rongliang Liu, Dr. Magnus Hedstrom, Dr. Andrew Kolchin, Professor Frank Harris, Professor Rodney Bartlett, Professor Andrew Rinzler, Professor James Fry, and Professor Amlan Biswas for many useful discussions and the help they have offered. I would also like to thank Professor Lin Ye, who was my undergraduate advisor at Fudan University. I was lucky to have such a knowledgeable mentor. Her integrity is also deeply respectable. Finally, I would like to express my great appreciation to my parents for their unconditional love and support; to my parentsinlaw, for their help and wonderful food I have enjoyed during the time of writing this dissertation; and especially to my wife, Li, for her love, support, understanding and patience, and for the joyful life I have been so lucky to have with her. TABLE OF CONTENTS Page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ............................................. .. .. .... .............. viii LIST OF FIGURES ......... ........................................... ............ ix ABSTRACT .............. ..................... .......... .............. xii CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. D ensity Functional T heory ........................................ ....................................... 1 C om plex Sy stem s ........................................ .................... 2 ChemoMechanical Phenomena in Amorphous Silica.......................................3 Nanopeapods: M materials by Design.......................... ........... ............... 3 2 M E T H O D S ............................................................................. 5 TotalEnergy Pseudopotential Calculations ...................................... ............... 5 BornOppenheimer Approximation...................... ... ......................... 6 Density Functional Theory ...................... ......... ............................. 6 LocalDensity Approximation..................... ...... ............................ 9 G eneralized G radient Correction...................................... ........ ............... 10 Pseudopotential Approxim ation ....................................................................... 10 The H ellm annFeym an Theorem ............................................. ............... 13 M molecular Dynam ics ................. ..................... .......... ....................14 Density Functional Theory for Periodic Systems.................................................14 Bloch's Theorem ................... ................ .................... ..... ........ 14 k p o in t S am p lin g ........................................................................................... 1 5 M om entum Space Form alism ........................................ ......... ............... 15 D ualSpace Form alism .................... ........................ .. ............... ... 20 Density Functional Theory for Finite Systems: A DualSpace Formulism...............21 Combining the Density Functional Theory with Classical Potential Functions: A M ultiscale A pproach..............................................................................27 In tro d u ctio n ................................................................................................... 2 7 Partitioning of the System ............................................................................31 T total E n erg y ............................................................................................3 2 The QM /CM interface ................................. ....................... ...............35 T he linkatom m ethod ....................................................... .....................35 T he pseudoatom m ethod ............................................................................36 3 MULTISCALE SIMULATIONS OF AMORPHOUS SILICA.............................37 In tro d u ctio n ........................................................................................................... 3 7 Test of the LinkAtom Method ........... .. ......... ................... 42 Test of the PseudoAtom M ethod.........................................................................52 Hydrolysis of a TwoMembered Silica Ring on the Amorphous Silica Surface........54 Generation of the Amorphous Silica Surface........................... ............... 54 Hybrid QM /CM Approach ............................................................................ 59 Sim ulation Setup .......................... ............ ............... .... ....... 61 R results .................... .... ..... ...... ....... ........ ... ..................... .......... ...... 63 Adsorption of One Water Molecule on the TwoMembered Ring....................63 Adsorption of Two Water Molecules on the TwoMembered Ring ...................72 Adsorption of W ater M olecules on the Silanols .......................... ..................73 D iscu ssion ......... .. .. ..................................................... ...................... .. .76 C o n clu sio n s..................................................... ................ 7 8 4 ENERGETIC AND ELECTRONIC STRUCTURE OF THE CARBON NAN OTUBE COM PLEXES........................................... ............................... 79 Structure of a SingleW all Carbon Nanotube............................................... 80 Electronic Structure of the Carbon Nanopeapods ..............................................81 In tro d u ctio n .................................................................... .. 8 1 M ethod and Sim ulation Setup ........................................ ......... ............... 83 Structures and E nergetics .................................................. ...............................84 Manipulation of FullereneInduced Impurity States in M@C60@SWNTs .........85 Electronic Structure of (KxC60)n@ SW NTs ........................................................ 92 Bromine Doping of SingleWall Carbon Nanotubes: Separating Metallic from Sem conducting N anotubes ............................................................................95 C o n c lu sio n s......................................................................................................... 1 0 0 5 CONCLUSIONS ......................................... .....................102 APPENDIX A TOTAL ENERGY OF A PERIODIC SYSTEM IN MOMEMTUMSPACE R E P R E SE N T A T IO N ....................................................................... ..................104 B POPULATION ANALYSIS BASED ON ELECTROSTATIC POTENTIALS .....112 L IST O F R E F E R E N C E S ......... .. ............... ................. .............................................. 115 BIOGRAPHICAL SKETCH .............. ........... ... ............. 123 LIST OF TABLES Table page 31. Several SiO bond lengths in combined QM/CM model and cluster model (cluster B ) ....................................... ......................................... 6 3 32. The adsorption energies of one and two H20 molecules in the cluster model with or without fixing the positions of H atoms that saturate the dangling bonds of the clu ster. .............................................................................. 69 33. The adsorption energies of water physisorbed and chemisorbed on amorphous silica surface. .............................................................................7 1 41. The cohesive energy per C60 molecule for (10, 10), (16, 0), (17, 0), (17, 0)b, and (19, 0) peapod..................................................... ......... ... ........... 85 LIST OF FIGURES Figure pge 21. The division of a system into QM and CM regions. ................ .............. 31 31. The structure of the Si207H6 cluster. ........................................................................42 32. The structures of the cluster shown in Fig. 31, obtained from (a) allquantum calculation, (b) combined QM/CM calculation using the linkatom method, and (c) allclassical calculation. ......................... .................. .................. ........... 43 33. The structures of the training system with isodensity surfaces, obtained from (a) all quantum calculation, and (b) combined QM/CM calculation........................ 45 34. The total energy change due to the change of (a) the Si(1)O(b) bond length or (b) the Si(2)O (b) bond length. .............................................................................. 47 35. The force on (a) Si(1) and (b) O(b), when the Si(1)O(b) bond length varies. .........48 36. The force on (a) Si(1) and (b) O(b), when the Si(1)O(b) bond length varies. .........49 37. The structure of the Si309H6 cluster. ................................ .. ......................... 50 38. The structures of the cluster shown in Fig. 37, obtained from (a) allquantum calculation, (b) combined QM/CM calculation using the linkatom method, and (c) allclassical calculation (bottom panel). ................ .......................... ............... 51 39. Cluster structures obtained from combined QM/CM calculations using the pseudo atom m ethod ...................................................... ................. 53 310. Partial pairdistribution functions for amorphous silica at 300 K ..........................56 311. One layer of the amorphous silica surface. ................................... ............... 57 312. Number density of atoms as a function of the distance from the center of the slab. 58 313. (a) The average SiO bond length, and (b) the average SiOSi bond angle, as functions of the distance from the center of the slab. ............................................58 314. A schematic picture of a silica surface (side view). .............................................. 62 315. Snapshots and potential energy of QM atoms along the dissociation path of one H20 m olecule on a silica surface. ............................... .................................. 64 316. C lu ster m models .........................................................................65 317: Valence electron number density contours on the plane crossing Sil, 01 and 02 atoms in the combined QM/CM surface model (a) and the cluster model (b).........66 318. Contour plot of(a) excess electron density and (b) deficiency electron density in the watersurface system relative to superposed electron densities of the surface and H20 on the plane crossing Si(1), 0(1) and 0(2) atoms in the combined QM/CM surface m odel. ........................................................................67 319. Contour plot of (a) excess electron density and (b) deficiency electron density in the watercluster system relative to superposed electron densities of the cluster and H20 on the plane crossing Si(1), 0(1) and 0(2) atoms in the cluster model...........68 320. Snapshots and potential energy of the QM atoms along the water dissociation path ("hydrogen relay process") on a silica surface. ............................................. 74 321. Two H20 molecules on top of two silanols. ...................................................75 41. A graphene sheet. The chiral vector and the translational vector are also shown.....80 42. The structures of (a) (10, 10), (b) (17, 0)a, and (c) (17, 0)b nanopeapods ...............85 43. Energy band structures of empty and C60 encapsulated (17, 0) nanotube. ...............86 44. Energy band structures of(a) C60@(16, 0)a, (b) C60@(17, 0)a, (c) C60@(19, 0)a. ....87 45. (a) The density of states projected on the nanotube wall in slice (I), (II), (III) and (IV). The sliced nanotube is shown in (b). Contour maps of the charge distribution of the lowest tiuderived energy band are shown for (c) a (17, 0)a semiconducting peapod and (d) a (10, 10) m etallic peapod. ........................................ ..................89 46. Energy band structures of(a) K@C60@(17, 0), (b) Ca@C60@(17, 0), (c) Y@C60@(17, O)a (spinup), (d) Y@C60@(17, 0)a (spindown), and (e) K @ C 6o@ (16, 0). ......................................................................90 47. The geometry ofK3C6o@(17, )b. ................ ............................................. 93 48. The density of states of (a) C6o@(10, 10), (b) KiC60@(10, 10), (c) K2C60@(10, 10), (d) K3C60@(10, 10), (e) K4C60@(10, 10), and (f) K5C6o@(10, 10). .......................94 49. The density of states of(a) (C6o)2@(17, 0)b, (b) (KiC6o)2@(17, 0)b (c) (K2C6o)2@(17, 0)b, (d) (K3C6o)2@(17, 0)b ........... ..................................... 95 410. The band structure of (a) the (10,0) nanotube, and (b) the (10,0) nanotube with a bromine atom/unit cell adsorbed on the surface. ....................................... 98 411. The band structure of (a) the (9,0) nanotube, and (b) the (9,0) nanotube with a bromine atom/unit cell adsorbed on the surface. ............. ........................ ......... 98 412. The band structure of (a) the (6,6) nanotube, and (b) the (6,6) nanotube with a bromine atom/unit cell adsorbed on the surface. ............. ........................ ......... 99 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THEORETICAL MODELING AND DESIGN OF COMPLEX MATERIALS By MaoHua Du December 2003 Chair: HaiPing Cheng Major Department: Physics Theory has been taking a particularly important role in the study of complex materials. This dissertation presents studies on two important complex systems: amorphous silica and carbon nanopeapods. A distinctive approach that combines density functional theory with classical potential functions is formulated, developed, and implemented to study silicon dioxide systems. We test this method on small silica clusters and applied it to the study of hydrolysis of the amorphous silica surface. The water physisorption, chemisorption, and reaction pathway on a surface defect site, known as a twomembered ring, are studied in detail. A new reaction pathway that involves the cooperative events is discovered. A carbon nanopeapod is a carbon nanotube encapsulating a chain of fullerenes. The encapsulated fullerenes insert impurity states in the energy band gap of the host semiconducting nanotube. The type of metal atoms encaged in the fullerenes and the space between the nanotube and fullerenes are found to be two important factors that control the energy levels and electron occupation number of these impurity states, and thus control the type of the majority carrier (p or n) and the carrier density of a semiconducting peapod. The possibility of superconductivity of potassiumdoped peapods is also discussed. CHAPTER 1 INTRODUCTION Over the last few decades, the fundamental techniques of theory, modeling, and simulation of materials have undergone dramatic advances. These new techniques, along with the parallel revolutionary advances in computer hardware and experimental techniques, enable us to better understand material properties and to design materials with specific desired characteristics. Based on the knowledge of the energy band structure for many important materials accumulated in the 1960s and 1970s, the computational materials research has moved on to many more complex issues, such as surfaces, defects, liquids, clusters, amorphous solids, etc. More recently, nanoscale materials have been added to the list. Density Functional Theory In this exciting journey of understanding the nature of materials, researchers around the world have developed many theories, among which the density functional theory (DFT) has emerged as one of the most powerful tools in studying material properties. This theory dates back to Thomas (1927), Fermi (1927, 1928a, 1928b) and Slater (1951), but the modern version dates to the seminal work by Kohn and his collaborators in the mid1960s (Hohenberg and Kohn, 1964; Kohn and Sham, 1965). In that work, they proved that the groundstate manyelectron problem can be replaced by an exactly equivalent set of selfconsistent oneelectron equations (KohnSham equations); and that the manybody part of the problem can be included by an effective exchangecorrelation potential. Therefore, the complex manybody problem is essentially reduced to the problem of a noninteracting inhomogeneous electron gas in an effective potential. If one knows all of the appropriate electronic interactions, for example, the electronelectron electrostatic interaction, the exchangecorrelation, the ionion interaction, and the electronion interaction, one can add up all these interactions and find the total energy. Since many of physical properties are related to total energies or to differences between total energies, we can study a wide range of problems within the framework of DFT. With some appropriate approximations and optimized algorithms, density functional theory has become a costeffective firstprinciples method and has reproduced a variety of groundstate properties within a few percentage of error compared with experimental data. Complex Systems In the past century, the approach in condensed matter physics was often to obtain elegant solutions to problems reduced to their ultimate simplicity. Physicists have identified fundamental interactions in matter that govern the electronic and transport properties of the materials. On the other hand, in the real world, many problems are of high complexity. Some materials or molecules have complex structures. Some phenomena are driven by multiple physical or chemical mechanisms. These complex problems are not always easily reduced to problems solvable by simple models. We must develop new tools to achieve new understanding of the complex problems in nanoscience, life science, and other emerging new scientific areas. Complexity takes many forms. In this work, we identify two complex systems. One system is amorphous silica, in which the complexity lies in the combination of the amorphous nature of the system and the interaction between the chemical process at nanoscale and the elastic field at macroscopic scale. The other system is carbon nanotube and fullerenebased hybrid materials, in which simple structures with different properties are selfassembled to form complicated structures with enhanced properties. ChemoMechanical Phenomena in Amorphous Silica The framework of solid state physics was built based on Bloch's theorem for electronic states in a periodic system. The study of the properties of an amorphous system has been a challenge for many years. In this work, we address a classic problem, the description and prediction of socalled chemomechanical phenomena, in which mechanical stress and chemical reactivity amplify one another. Such phenomena are extremely important in technologies as diverse as environmental remediation, advanced aircraft structures, and innovative energy storage devices. However, the theoretical description and prediction of chemomechanical phenomena are still in a primitive state because the chemical processes at nanoscale and the mechanical properties at meso or macroscopicscale interplay and have been traditionally studied with different methods in different fields. In this work, we combine density functional theory with classical interatomic potential functions via a linkatom method (Du and Cheng, 2003a). This multiscale simulation method has been implemented and tested in small silica clusters (Du and Cheng, 2003a), and has been used to study hydrolytic weakening of the amorphous silica surface (Du, Kolchin, and Cheng, 2003a, 2003b). Nanopeapods: Materials by Design Central to nanoscience is the assembly and manipulation of fundamental building blocks on the nanoscale to create functional structures, materials, and devices. In nanoscale, zerodimensional dots or nanocrystals, onedimensional wires, and two dimensional films often have unusual properties distinctly different from those of the same material in its bulk form. Furthermore, these materials can be used as building blocks to design more complex materials. Unlike traditional materials, nanostructured materials are often highly complex and heterogeneous in shape and substance. The discovery of carbon fullerene (Kroto et al., 1985) and nanotube (Iijima, 1991) is another manifestation of Richard Feynman's 1959 prediction that there "is room at the bottom". Inserting metal atom(s) inside a fullerene cage creates a metallofullerene with different chemical properties from the empty fullerene. The selfassembly of two basic building blocks (nanotube and fullerene) results in a new class of hybrid materials: so called nanopeapods. In our work, electronic structures of several semiconducting and metallic carbon nanopeapods have been studied using the density functional theory (Du and Cheng, 2003b). We find that the fullereneinduced states inside the band gap of the hostsemiconducting nanotube can be treated as impurity states. The manipulation of these impurity states controls the type of the majority carrier (p or n) and the carrier density of a semiconducting peapod. In addition, we discuss the possibility of superconductivity of potassiumdoped peapods. CHAPTER 2 METHODS This chapter describes the methods that have been used in this work in studying material properties. The calculation of the electronic and geometric structure of a complex system is a formidable task. Such calculations can only be performed with many simplifications and approximations. Even within the framework of density functional theory, there is a large volume of literature describing numerous techniques and algorithms that optimize the accuracy and speed of the totalenergy calculation, and this literature is still growing. A thorough review of this field is certainly beyond the scope of this work and the author's knowledge. Instead, a short introduction to some of the important ingredients in the totalenergy calculations is given. In this chapter, I first introduce the BornOppenheimer approximation, density functional theory, the pseudopotential method, the HellmannFeynman theorem, and molecular dynamics. An excellent review of the totalenergy pseudopotential method can be found in the article by Payne et al. (1992). Next, I describe the details of the implementation of the DFT calculation in periodic and finite systems (Ihm et al., 1979; Yin and Cohen,1982b; Bamett and Landman, 1993). Only the planewave based implementation is discussed. Finally, I describe the hybrid method that combines DFT with classical potential functions (Du and Cheng, 2003a). TotalEnergy Pseudopotential Calculations To perform accurate and efficient totalenergy calculations, a series of simplifications and approximations is needed. These include the BornOppenheimer approximation to separate the electronic and nuclear wavefunctions, density functional theory to model the electronelectron interaction, and the pseodopotential method in conjunction with a planewave basis set to model the electronion interaction. BornOppenheimer Approximation The quantum mechanical total energy and the subsequent minimization of that energy are dependent on the electronic and nuclear coordinates. To simplify the problem, the nuclei can be treated adiabatically, leading to the separation of the electronic and nuclear coordinates in the manybody wavefunctions. This socalled BornOppenheimer approximation is justified by the fact that the nuclei are much heavier than electrons; and thus electrons respond essentially instantaneously to the motion of the nuclei. This approximation reduces the manybody problem to the solution of the dynamics of the electrons in a frozen configuration of nuclei. Density Functional Theory Hohenberg and Kohn (1964) proved that the total energy of an electron gas in the presence of an external potential is a unique functional of the electron density. The minimization of the totalenergy functional of the electron density yields the groundstate energy of the system, and the corresponding electron density is the exact groundstate electron density. Kohn and Sham (1965) then showed that the manyelectron problem can be replaced by an exactly equivalent set of selfconsistent oneelectron equations (Kohn Sham equations). To simplify the notation, I will present the total energy functional and KohnSham equations for the spinunpolarized case. They can be easily extended to the spinpolarized case. For more details about the density functional theory, see Parr and Yang (1989), Jones and Gunnarsson (1989), and Dreizler and Gross (1990). The total energy of a system of ions and electrons can be written as Eota[ {R }, p(r)] = E, ({R,})+Eeec[p (r)] (21) Here, ELo,({R, })is the Coulomb energy associated with the interactions among the ions at positions {R}) given by En JZ (22) I,j IR Rj I where Z1 and Z, are the charges of the I and Jth ions, and Hartree atomic units are used. Eeec[p(r)] is the groundstate energy of the electrons. The groundstate electronic energy is given by Ez,, = TE + E + E (23) where Te is the kinetic energy, EeL is the electronion interaction energy, and E,, is the electronelectron interaction energy. These terms are given by T=If _ V2 ), r (24) ,ZfJy( 2 (24) E = JV,, (r)p (r)d3r (25) E = p(r)p(')d rd r' + E [(r)] (26) 2 11 rr' p (r) = f ,(r)2 (27) where Vo is the static total electronion potential, p (r) is the electron density, f is the electron occupation number, and E.x[p(r)] is the exchangecorrelation energy functional. The set of oneelectron states y, that minimize the energy functional Eeec is given by the solutions to the KohnSham equations (Kohn and Sham, 1965): V2 + o,"(r) + VH(r) + (r) (r) = s (r) (28) where VH (r) is the Hartree potential of the electrons given by VH P r (29) V r rr' The exchangecorrelation potential, V (r), is given by the functional derivative V (r) (r) (210) 6p(r) The KohnSham equations represent a mapping of an interacting inhomogeneous electron gas in the presence of an external potential V(r) onto a noninteracting inhomogeneous electron gas moving in an effective potential Vf (r) = V(r) + VH (r) + Vxc (r). If the exchangecorrelation functional were known exactly, the KohnSham equations would give the exact groundstate p and Eiec . The KohnSham equations are nonlinear and must be solved selfconsistently. Once the eigenvalues s, and the wave functions y, are obtained, the total electronic energy Ee ec can be determined from Eq. (23) to (27), or alternatively from the formula c = c P(r)p(r )d3rd3r'+Ex [p(r)] ixc(r)p(r) (211) r2 rr' where N is the number of occupied KohnSham orbitals. Here N N S(212) The sum of the singleparticle KohnSham eigenvalues does not give the total energy, because it overcounts the electronelectron interaction and the exchangecorrelation energy. LocalDensity Approximation The KohnSham equations [Eq. (28)] are exact provided that an exact exchange correlation functional E, is known. However, the exact E, is usually not known, hence an approximation must be used. The simplest approximation is the localdensity approximation (LDA) proposed by Kohn and Sham (Kohn and Sham, 1965). In LDA, the exchangecorrelation functional can be written as Ec = c(r) p (r)dr (213) Thus, C Ec [p (r)] = [p(r)sx(r)] 6p(r) p (r) The local density approximation assumes that the exchangecorrelation energy functional is purely local. The exchangecorrelation energy per electron s, is equal to the exchangecorrelation energy per electron in a homogeneous electrongas that has the same density as the electron gas at r. Thus, a'(r) f'A[p(r)] =s [p(r)] (215) A commonly used parameterization for afLA is the VoskoWilkNusair parameterization (Vosko et al., 1980, 1982) of the CeperlyAlder data (Ceperly and Alder, 1980). In this parameterization, the quantum Monte Carlo results by Ceperley and Alder were interpolated by Vosko et al. to provide an analytic form for the exchangecorrelation energy. Considering the inexact nature of the local density approximation, it is remarkable that calculations using LDA have been so successful. Generally, totalenergy differences between related structures can be believed to within a few percent and structural parameters to at least within a tenth of an A. Cohesive energies, however, can be in error by more than 10%, sometimes much more. Generalized Gradient Correction For small density variation around the uniform electron density, LDA is an excellent approximation to the exchangecorrelation energy. However, it tends to introduce a larger error for systems having strong spatial variation of electron density. Thus, the generalized gradient approximation (GGA) is often used forE, In it E, is a functional of both the electron density and the gradient of the electron density (Perdew et al., 1996). In many applications, the GGA showed improved results over LDA in total energies (Perdew et al., 1992); atomization energies (Perdew et al., 1992, 1996; Becke, 1992; Proynov et al., 1995); energy barriers and structural energy differences (Hammer et al., 1993; Hammer and Scheffler, 1995; Hamann, 1996; Philipsen et al., 1994). The GGA usually expands and softens bonds, an effect that sometimes corrects and sometimes overcorrects LDA results. Typically, GGA favors the density inhomogeneity more than LDA does. Pseudopotential Approximation The totalenergy density functional method is often implemented in a planewave basis set, which is discussed later. Here, I introduce the pseudopotential approximation that greatly reduces the computational cost of the planewave based DFT method. A large number of planewaves are needed to expand tightly bound core orbitals and rapidly oscillating valence wavefunctions in the core region. The pseudopotential approximation (Phillips, 1958; Heine and Cohen, 1970; Yin and Cohen, 1982a) allows one to replace the core electrons and the strong ionic potential by a weaker pseudopotential that acts on pseudowavefunctions rather than the true valence wavefunctions. For a central field atom, the radial KohnSham equation is given by S 2+ t+ [p(r)] rRn,(r)=+ s[rR ,(r) (216) 2 dr2 2r2 The pseudowavefunctions are obtained from an allelectron atomic calculation that self consistently solves Eq. (216). Once the pseudowavefunction is obtained, the screened pseudopotential is given by the inversion of Eq. (216), PPcr / 1(1+ 1) 1 d2 , [I +rRP (r)] (217) 2r2 2rRPP(r) dr2 where RPP(r) is the radial pseudowavefunction. The ionic pseudopotential can be obtained by subtracting the Hartree VPP(r) and exchangecorrelation VP (r) potentials calculated from the valence pseudowavefunctions from the screened pseudopotential PPsc (r). Vr (r) = V PPs ( (r) VV4 (r) (218) The pseodopotential needs to meet a set of requirements. The standard requirement is "normconservation", which means that the pseudowavefunctions are identical to the real wavefunctions outside a core region, and that they generate same electron density inside the core region. A widely used normconserving pseudopotential is the Troullier Martins pseudopotential (Troullier and Martins, 199 la). An accurate pseudopotential should reproduce the scattering of the ion core. Obviously, an accurate description of such scattering requires a nonlocal pseudopotential, i.e. one that uses a different potential for each angular momentum component of the wavefunction. The pseudopotential operator can be decomposed into local and nonlocal parts, V~p (r) = VPPc(r) + VfPPnc (r)P (219) where VPPl'(r) is the local potential, and the semilocal potential JPPn"ic(r) is given by SPP'nic(r) = VPP (r) VPPloc (r) (220) Here, P, is the operator that projects out the /th angular momentum component from the wavefunction. The nonlocal potential is shortranged, since at large r, the ionic potential is independent of angular momentum and approaches the value of, where Z is the core r charge. In principle, the local potential can be arbitrarily chosen. However, since the summation in Eq. (219) can be only taken for a few angular momentum terms, the local potential should be chosen such that it adequately reproduces the atomic scattering for all the higher angular momentum channels. Kleinman and Bylander (1982) suggested a nonlocal potential to replace the semi local term in Eq. (220), KC .f/ lc(r)(DP (r)) (PP (r)V "c(r) V, Kc(r) =I (221) (O ((r) ) OPP(/)c ( (r) where PP'"f"(r) is given by Eq. (220), and OPP(r) is the atomic pseudowavefunction for the angular momentum 1. This form of the nonlocal potential substantially reduces the number of integrals that need to be evaluated. The normconserving condition requires that the total pseudocharge inside the core region match that of the allelectron wavefunction. Thus for those highly localized valence orbitals (e.g., 2p for firstrow atoms and 3d for transitionmetal atoms), a large number of planewaves are still needed. Vanderbilt relaxed the normconserving rule and constructed ultrasoftt" pseudopotentials with a lower cutoff energy for planewave expansion (for detailed discussion, see Vanderbilt, 1990; Laasonen et al., 1991; Kresse and Hafner, 1994). The advantage of such pseudopotentials is that they reduce the necessary energy cutoff for transition metals and firstrow elements by a factor of 2 to 4 (Kresse and Furthmuller, 1996a). The HellmannFeyman Theorem Once the total energy of the system is obtained via Eq. (21), the force on ion I, f, is given by dE OE OE Oy OE 8Oy df, C E C_ E (222) dRI ORI OJ R, ORI O 17R, aE If {u, } are eigenstates of the Hamiltonian, and given that E is just Hy,, the last two terms in (222) can be rewritten as R,Hy, + Y R, (RI / ORI/ f, (224) aRI This is the HellmannFeynman theorem (Hellmann, 1937; Feynman, 1939), which holds for any derivative of the total energy. Molecular Dynamics The motion of ions on the BornOppenheimer potential energy surface is governed by the forces obtained by Eq. (218) via integration of the Newtonian equations of motion m, RI = fI (225) The integration can be performed using the Verlet algorithm (Allen and Tildesley, 1987). In contrasted to empirical potentials, the forces determined by abinitio methods are much more accurate, and should be reliable in describing the dynamical evolution of a system. Density Functional Theory for Periodic Systems For a periodic solid with infinitely large size and an infinite number of electrons, we must make use of the periodicity of the system and apply Bloch's theorem in order to perform a realistic totalenergy calculation. Bloch's Theorem A wavefunction in a periodic system can be written in the following form (Ashcroft and Mermin, 1976) VY,k (r) = exp(ik.r)u,k (r) (226) Here, u1,k (r) is a periodic function that satisfies the condition uk (r) = uk (r + nR) (227) where R is any lattice vector and n is an integer. Such a function can be expanded in a discrete planewave basis set with wave vectors as reciprocal lattice vectors of the crystal: k (r)= cG exp(iG.r) (228) G Therefore the electronic wavefunction can be written as ,,k (r) ,kG Iexp [i(k + G).r] (229) G, k+G2 Ecu The summation is over all reciprocal lattice vectors up to a cutoff energy Ek,. But in principle Eut = o. kpoint Sampling The number of k points in the first Brillouin zone is equal to the number of unit cells, which is an infinite number. However, there are only a finite number of occupied states for each kpoint. Therefore, the Bloch theorem transforms the problem with an infinite number of electronic wavefunctions to one with a finite number of electronic wavefunctions at an infinite number of k points. To further simplify the problem, the infinite number of k points can be sampled using a special set of k points in the first Brillouin zone (Chadi and Cohen, 1973; Monkhorst and Pack, 1976). The accuracy of the calculation can always be improved by using a denser set of k points. With Bloch's theorem, kpoint sampling, and a kinetic energy cutoff, we reduce the problem with an infinite number of variables to one with a finite number of variables, and thus are ready to implement the density functional theory in a periodic system. MomentumSpace Formalism In this section, I show the derivation of the momentumspace formalism of DFT for periodic systems. For relevant references, see Ihm, Zunger, and Cohen (1979), Yin and Cohen (1982b). The total energy of a crystal in the pseudopotential approximation can be expressed as follows: E,,a = T, + E +EH Ec + E,, (230) where Te =:ZfJd3' (r)f V2 v(r) (231) E = f d'r~ (r)V (rRts)p,(r) (232) i,R,s EH = d3r p(r)VH(r) (233) Ex = p (r)sc (r)d3r (234) E =NY' zSS (235) o 2 R,s,s' It, t, + R p(r) = f k,(r)12 (236) The symbols / and y, (r) are, respectively, the electron occupation number and the valence pseudowavefunction for the oneelectron state i. Nis the number of unit cells in the crystal, and R is any lattice vector. VH (r) is the Hartree potential, given by Eq. (29). The symbols Zs and ts are, respectively, the core charge and the basis vector for atom s in the unit cell. The prime in Eq. (234) means that the s=s' term is excluded in the summation if the s and s'th ions are in the same unit cell. VPP is the angularmomentum dependent pseudopotential and is given by PP (r)= ~PP(r = V(PPlc()c+ PPPnlc (r) (237) 1 / The momemtumspace representation of each energy term per unit cell in Eq. (2 30) is given by T  IC kk Ck+G 2(k+G)2 (238) 2 k k G E p (G) EH 47. (239) 2G)G Exc = Q2 p*(G)Fxc (G) (240) G Ee = cZp (G) e t rPPlc (G) GO s + W .* c Yr(G'G)otsPPn (241 + ""cZ k ,kek+G Lk+G'Zl L s,l,k+G, (241) k,,l G,G' +QZl d'r VsPPic(r) s s irts 1n Ierfc( tt, t,, + R) 2f E, 2 Z tt, +R (242) 4xn 1 G \ + 4 exp G cos[G*(ts ts ) q G#o G2 4r1J2,1 q Q where p G)= p(r)e G*rdr (243) P (r)= ZZC kf, k Z Ck+GCk+G e(G' G)*r (244) 4 k i G,G' VPPlc (G)= Jd3rVfPP (r)e,G.r (245) ec (G)=' ex(r)e G*rd3r (246) V",. kG, (21+l)P(cosy)ji j(, (k+ GIr) (r) j,(k+G'r)rdr (247) (k + G).(k + G') cosy = (248) k+G k+G'l Here, qc is the volume of the unit cell, cok is the normalized weight for each k point, Qe is the total charge of electrons in the unit cell, j, and P, are spherical Bessel functions and Legendre polynomials, respectively, and is a parameter controlling the convergence of the Ewald summation. See Appendix A for details of the derivation. If the KleinmanBylander type of nonlocal pseudopotential is used, Eq. (247) becomes P 47t (2/+ 1)P (cosy) [ (k +( r)pp RP i(r)r'2dr V"PI +( R(r)l, PPnlc ) RPs I l) Ri (249) x j^i (k+G'r)(PPnc RIs 2dr where RPs (r) is the radial pseudowavefunction. If there are n different G vectors, for each I and each k point, one reduces the required number of integrals from n(n+1)/2 to n by using the KleinmanBylander nonlocal pseudopotential. The integrals in Eq. (249) can be evaluated in either reciprocal or real space. In reciprocal space, the number of operations scales with the size of the basis set, which increases proportionally with the size of the system. Since the nonlocal pseudopotential is shortranged and thus is only nonzero within a sphere around each ion, the number of operations necessary to evaluate the integral in realspace does not increase with the system size. Therefore, it is more efficient to evaluate the integrals in Eq. (249) in real space for large systems. In order to do so, the wavefunction for each angular momentum must be represented in real space. However, the discrete realspace Fourier grid introduces errors when projecting angular momentum components of the wavefunctions because of the nonuniqueness of the spherical harmonic projection in a discrete real space grid (Payne et al., 1992). To reduce these errors, KingSmith, Payne and Lin (1991) proposed a procedure for improving the accuracy of the realspace calculation of the nonlocal contribution to the total energy. Tests of this method showed good agreement between real and reciprocalspace calculations (KingSmith, Payne and Lin, 1991). The KohnSham equations in momentumspace representation can be derived variationally from the expression for the total energy in Eq. (230), ZI k+ G GG +k,G,G ,k+G =E,k C,k+G (250) G' 2 where VkGG' = k VH(G' G)+V,(G' G)+ el(GG) rlc(G ' G)+ sk+G,k+G'(251) Here, VH (G) is given by Eq. (A. 10). In practice, the divergent terms VH (0) and V'pp1' (0) are set equal to zero when solving Eq. (250). (See Appendix A for details of the discussion.) These contribute to the diagonal terms in the Hamiltonian matrix. Since the trace of a matrix is independent of the representation, setting these terms to zero corresponds to a constant shift of the total energy. A compensating term [the last term in Eq. (241)] should be added to the total energy. Once Eq. (250) is solved, the total energy can be obtained, E,,to = ZZokfkS,k EH + EXC q P* (G)Vxc (G) + E k i G + IQ dr VP(r) ) (252) + Q r Ve~o(r) y where EH, Ex and E,,, are given by Eq. (239), (240) and (242), respectively. This expression is simpler to evaluate than Eq. (230), (238)(242), since the double summation over reciprocal lattice vectors in Eq. (241) is absorbed in the simple summation of the eigenvalues of the occupied states in Eq. (252). The HellmannFeynman force on each ion is given by F Fs,+Fs, (253) where F, = V'E5 = iQ p* (G) Ge'G'rVPP'c (G) (254) G#O 7i GG )t.., PP~nl c i ZZ kf.kCk+G ,k+G (G'G)Z[ e(G' G)'t P ,k+G k,,l G,G' s and Fs,, = V E = Zs' 2 exp sin [G(t t)] (255) ss iC o GO 42) Rt t, )R r erfc( t(r tR t + R) 21R et,t,,+R2} SIt(t ts + RRI J I t ts, + e R t,t,+R t,t,+R Fs,, is the force contribution from the valence electrons, and Fs, is the force contribution from the ions. DualSpace Formalism The procedure of evaluating the kinetic energy in momentum space and the potential energy in real space can be thought of as a dualspace formalism for pseudopotential calculations (Martin and Cohen, 1988; Troullier and Martin, 1991b; Wentzcovitch and Martin, 1991) in contrast with momentumspace formalism (Ihm, Zunger and Cohen, 1979) in which all the calculations are carried on in momentum space. The dualspace formalism is more computationally efficient than the momentum space formalism for large systems. The results from both methods can converge to the exact solution with increasing energy cutoff with similar rates of convergence. Next, I will introduce the dualspace formalism for the pseudopotential calculations of finite systems. Density Functional Theory for Finite Systems: A DualSpace Formulism For the study of finite systems, the method described in Section (22) can be used in conjunction with a supercell containing a large volume of vacuum that reduces the interaction energy between repeated replica of the physical system to a negligible value. However, if the system has a net charge or dipole moment, the interaction between image 1 1 point charges decreases as, and the corresponding value for dipole interaction is, L L where L3 is the volume of the supercell. Thus, the energy converges very slowly with respect to the size of the supercell. Although methods have been devised to accelerate the convergence (Makov and Payne, 1995), it is still difficult to treat a net charge distribution in the supercell. Here, I will introduce a planewave based method that does not use a supercell, and thus does not involve repeated images of the system are involved (Barnett and Landman, 1993). This approach allows one to study systems possessing large multiple moments efficiently, though they may become prohibitive for a supercell method. Following the notation of Barnett and Landman (1993), the eigenfunctions of the KohnSham equations are denoted xy,, where o is the spin index. The corresponding density of spin o is given by p, (r) = f (r)2 (256) where {fj } are the occupation numbers, with ~ fj, = N, and N, is the number of valence electrons. In a finite system, the total energy is expressed as usual by Eq. (230). Each energy term in Eq. (230) is given by =T f= J (Jy (V2 jy) (257) E ZZ = (Z EH = d3r p (r)VH(r) (259) E, p (r);s (r)d3r (260) Son = (261) s' s' ts ts, where VH is the Hartree potential, given by Eq. (29), VPP is the pseudopotential, which can be decomposed into local and nonlocal parts, JVPP = (r, r') = VPP" (r) + VJ," (r, r') (262) Thus, E1 = Ec + En (263) where ELc = d3r p(r)VPP"c (r t,) ) (264) r I V^PPn (Ir^I VPn' (Ir tlP jo) EjV c (r Y /R, (r)Y, ()) R,/ (r)Y,() V nc r t j) S s dr[rR, (r)]v I nc (r) (265) ,j C s,l,m 1 F1, = (266) fdr rRs, (r) V flC (r) Ks,Rm (r)= V nc (r)R], (r) Ym (i) (267) In Eq. (265), the KleinmanBylander nonlocal pseudopotential is used (Kleinman and Bylander, 1982). The KohnSham Hamiltonian is given by 1 H= V2 + +VH + VX (268) 2 where V E a c(p)] (269) 6p ap Vat (r,r') = Vic (r) + Vfc (r, r ') = VPPc (rt, )+ j ,tKs,, (rts)K,,m (r'ts) (270) s s l,m and VH is the Hartree energy, given by Eq. (29). Barnett and Landman (1993) suggested a dualspace method in which the kinetic energy is evaluated in the momentum space and the remainder of the energy is evaluated in the real space. The Hamiltonian matrix element (jyo H lka) is also evaluated in the dualspace, (jc H kc)= (j ( iV2+V ((r,r')+ (r)+V, (r)) k) =Jd 3k(j I V2k)l)(kl ko) + 3d k d3rld r'(jo k)(k r')(r' (V, (r,r')+VH (r)+Vc(r))Ir)(rI k) (271) = Jd k(j V2 k)(k ko)+ d k(jcy k)Jd r(k r) x [r (VH (r)+ V (r)+ VePP"c (r))r)(r ky)+ d3r'(r rPPn"c r')(r' ky)] In practice, the integrals over real and momentum spaces are replaced by summations over the discrete real and momentumspace grids. The size of the realspace cell is defined as follows r (x,y,z); 0 The wavefunction \ ,, (r) = (r jo) is expanded in a planewave basis within the cell, i.e., 1 y (r)= 4 (g)eg'r (273) g9gmax for r in the cell, and vo (r) = 0 (274) for r outside the cell. In Eq. (273), Q is the volume of the cell, Q = LLyL and gmax is the cutoff momentum in momentum space. The momentumspace grid (gspace) is defined by g = 2 k, k (275) Lx Ly L, are integers satisfying where {k }a= ,y, and Nx N N ax (277) Lx Ly Lz 7x For a finite system, the wavefunction can be chosen to be real in realspace, i.e., o v(r) = J (r) (278) and thus, J, (g) =j*T (g) (279) The density p, (r) can also be expanded in a Fourier series but requires a doubled momentumcutoff with respect to the wavefunctions, i.e., p,,(r). 1 :I I^ (0: (0 (g)g').r g g (280) 1 SGD (G)eGr The Gspace grid in Eq. (280) is defined as G=27r m my m (281) [L L,) where {m }a=,,y,z are integers satisfying M < m < and M = 2Na (282) 2 2 Now, a realspace Sgrid can be defined as a dual of the Gspace grid. The density and the KohnSham potential can be evaluated on this Sspace grid, i.e., S= nx L nyL, nZ L M 'My ',M MMM A A x y z/M where n x,y (. acLx,v, are integers satisfying The wavefunction can be evaluated on the Sgrid as The wavefunction can be evaluated on the Sgrid as 0 (S)=JM (S) (284) (285) I (g)e'g.s Jr g where M = MXM The density on the Sgrid is given by D, (S)= If, (S)2 (286) Now the different energy terms in Eq. (259), (260), (264) and (265) can be directly evaluated on the Sgrid as follows 1 EH (S)VH (S) (287) Exc = p (S)Ec (S) (288) S s Ec = p (S) I VPP,c (S t, ) (289) S 2 Ecd = f Fs,, K, (S ts)yk, (S) (290) j,0 s,i,m S The kinetic energy is evaluated in momentum space. The Hamiltonian matrix element in Eq. (271) is written as (283) 1 2 (jo H ko) = 2 2 (g) (g) + (g) V, (S) + V, (S) + IP (S t) j (S)(291) + ZZF ,K (S ts)Z Km '(S ) t (S e g.s s Im S' The force on each ion is given by F V = VVE J  (292) Combining the Density Functional Theory with Classical Potential Functions: A Multiscale Approach Introduction Many issues of material behavior, such as corrosion, surface polishing, etching, and crack growth in solids, are of great interest for both fundamental and practical reasons. Studies of these material behaviors are essential to today's electronics, optics and construction industries. Many important phenomena of materials have been observed experimentally, but the underlying mechanisms are still not clear and require theoretical studies. However, the majority of material behaviors including those mentioned above involve complex chemical and mechanical processes. Theoretical studies of those complex processes require large scale calculations. Such computational studies of materials are highly interdisciplinary, requiring broad knowledge of physics, chemistry and computational techniques. Different material properties require different levels of refinement. If the mechanical and thermodynamical properties are of interest and trajectories of the particles are dictated by classical mechanics, classical methods normally can work well to model the system. Continuum mechanics methods and approximations, such as finite element, have been used for decades to describe mechanical properties of specific materials on the meso and macroscopic length scales. However, for simulations of systems that are greatly perturbed from equilibrium, continuum mechanics methods that are based on linear elasticity theory are not appropriate. Instead, atomistic statistical mechanical methods, such as molecular dynamics (MD), normally are used to study such systems, in particular to describe thermal fluctuations and pressure waves, etc. Many subjects, including surface collisions, crack propagation, and firstorder phase transitions, have been studied by classical molecular dynamics (CMD), since it offers good statistics with respective to trajectories of a large number of atoms. In classical MD, the atomistic interactions are described by potential energy functions. Normally they are derived from fitting parameters to match experimental data or quantum mechanical calculations or their combination. With modern computer technology, millions of atoms can be simulated by classical MD. However, the absence of any explicit description of the electronic structure makes CMD an untrustworthy method in describing chemical reactions, in which charge transfer is important. In most cases, quantum mechanics (QM) must be used to describe the chemical processes of materials accurately, since they involve bondbreaking and bondforming. But the computationally demanding nature of QM methods prohibits their use in large atomic systems with thousands of atoms. Therefore, a challenge arises when both chemical and mechanical properties are important in material behaviors (for example, the chemical environment can change the mechanical properties of materials). In fact, this challenge is part of a very large problem that stood unresolved for many years. Researchers in basic science often work in different fields and develop theories and methods for different types of problems at different time and length scales. Much research work has been dedicated to two extremes: extremely small and extremely large systems. The key links between them are often missing. This problem is central to many problems in biology and materials science in which information at different lengthscales must be integrated to analyze complex processes. The work we have been doing is to solve a particular aspect of this problem, in which the driving forces of material behaviors in microscopic and macroscopic levels interact with each other. For example, the chemical environment can influence the mechanical properties of many materials. Experiments show that water, ammonia, hydrazine and methanol can increase the rate of slow crack growth in silica and alumina (Michalske and Bunker, 1984). The corrosion or crack growth rate is influenced by the number of chemically active sites along the crack front and the rate of their reaction with the environment. Unfortunately, we lack good tools to study this kind of material behavior involving both chemical and mechanical processes. The simulation of chemical reactions at the crack tip requires the use of quantum mechanical methods, but the large number of atoms which are needed to simulate crack growth excludes the possibility of allquantum mechanical simulations. In many previous studies on surface corrosion and slow crack growth, large systems were replaced by small model systems. However, this approach causes errors, such as missing structure constraints and neglected longrange interactions. Therefore, new methods that can integrate different scales of calculations are urgently needed. In the last decade, combined quantum mechanical and molecular mechanical methods (QM/MM) have become very useful tools in the studies of the reaction mechanisms of enzyme and protein, solvent effects, catalysis and many other topics. These methods integrate quantum mechanical and classical mechanical methods and aim to treat large biological systems (Field et al., 1990; Zhang et al., 1999; Eichinger et al., 1999; Antes and Thiel, 1999; Reuter et al., 2000; Kairys and Jensen, 2000; Titmuss et al., 2000; Hall et al., 2000; Nicoll et al., 2001; Chaban and Gerber, 2001; Cui and Karplus, 2002). Broughton et al. (1999) developed a multiscale method, which couples finite element, molecular dynamics, and semiempirical tightbinding methods, to study crack propagation in silicon. However, the tightbinding method is oversimplified and not able to describe bond breaking properly. Sauer and Sierka (2000) have developed a combined quantum mechanicsinteratomic potential function approach to treat extended systems. A very useful application of this method is the prediction of reaction rates for an elementary process on the surface of solid catalysts and how these rates differ between different catalysts with the same active site. Based on recent progress in multiscale simulation techniques, we have developed a combined quantummechanical and classicalmechanical (QM/CM) method (Du and Cheng, 2003). We also realize that one of the important applications of such combined QM/CM method is the study of amorphous materials, which are difficult to treat in an all quantummechanical method since a large number of atoms must be included in the simulation to describe the amorphous nature of the system. Thus, we have implemented this method to study amorphous silica. The test of the method in small silica clusters and its application to the hydrolysis of amorphous silica will be described in the Chapter 3. First, I will introduce the combined QM/CM method. Partitioning of the System In many chemical processes, for instance, chemical reactions on the surface of an amorphous material, charge transfer only takes place in small areas where the chemical reactions occur. In other areas, the atoms only form an environment for the particles involved in the chemical reactions. The first step in any QM/CM approach is the division of the system into two regions: QM region and CM region, as shown in Figure 21. The particles inside the QM region are modeled as electrons and ions, and treated by QM methods. These particles are involved in the chemical reactions. The particles inside the CM region are modeled as atoms, and treated by CM methods. These atoms serve as an environment to the QM region. Long range interactions QM Figure 21. The division of a system into QM and CM regions. Total Energy The total energy of the system is divided into three components: the energy of the quantum subsystem, the energy of the classical subsystem, and the interaction between the quantum and classical subsystems, i.e., E = EQM (tr }, p (r)) + Ec ({r, + E c (r },{r ) (293) Here, the QM nuclei, CM atoms are labeled by indices n and s, respectively. In Eq. (2 93), only EM is explicitly dependent on the electron density p(r), which is obtained by solving the KohnSham equations selfconsistently V2 +V V,(r d' p(r' + 6E ,(r) = s ,(r) (294) L2 I r r' 6p(r) p (r) => f ),(r)12 (295) where V, (r) is the ionic pseudopotential, Ex is the exchangecorrelation energy, and {f } are the electron occupation number. The different energy components can be computed as following: +f fJ ~* (r)V, (rV, (r)dr +1 jdrdr' p(r)p(r')+E~[p(r)] I2 nr r' 1 1*2 Ecm Z R. +U({R }) (297) EQmCM = V({Rn},{Rj}) (298) where Z, is the core charge of the nth ion. U({Rs}) and V({Rn},{Rs}) are described by sums of classical interatomic potential functions in pairwise or threebody form. The particles inside the QM region are treated as ions and electrons when computing EQM. However, for computing E,,,M, they are simply treated as atoms and interact with classical atoms via the same potential functions that are used to describe the classical interatomic interactions. In other words, a classical atom does not distinguish an atom inside the QM region from one inside the CM region. Such treatment of EQMeCM ensures the continuity of the mechanical properties across the QM/CM boundary. If simple pairwise classical interatomic potential functions (, are used, U({Rs}) and V({Rn}, {Rs}) are expressed as U({R,})= ~(SP,,) (299) s s'>s V({R },{R }) = I I ,D, (r, ) (2100) n s One can also use threebody potential functions having the form D = C f2 (i, j)+ jY f (i, j, k), (2101) 1 J>1 i j>i k>j where s is an energy parameter,f2 is the pairwise interaction between the ith andjth atoms, andf3 is the threebody term as a function of the i,j, kth atoms. In this case, U({Rs}) and V({Rn}, {Rs}) are expressed as U({RJ})= II F f:Ls23 +3(, ,,, ,,,) (2102) Ss'>s s s'>s s">s' V({R },{R }J)= Y f2 (s) (2103) +ZZ j f,(r,, ^ ,,^ ) + ZZjZ 3(,r r,,r,s,) n S S >S n ns>N S n s s'>s n n'>n s The forces acting on the QM and CM atoms are respectively given by F = V,(EQM +EM/CM) (2104) and F, = V,(EM + EQM M) (2105) In many other hybrid QM/CM methods, the QMCM electrostatic interactions are described by the sum of the interactions between QM electrons and CM point charges and the interactions between QM ions with CM point charges, i.e., E electrotahc P(r)Cqd Zn q, EQMc r d r + Rq (2106) s ir R. n,s Rn Rsl where q, is the partial charge ofsth classical atom, and Z, is the nth quantum ion. In these methods, EQM cM includes the electrostatic interactions between the QM electrons and CM point charges and is explicitly dependent on the electron density p (r). However, for those CM atoms close to the QM/CM boundary, the point charges assigned on each of these CM atoms result in an unrealistically large electrostatic potential being applied to QM electrons. Also, Newton's third law breaks down, i.e., the force acting on the QM particles due to the classical particles is not equal to the force acting on the classical particles due to the QM particles. These problems are amplified if the CM atoms are included in the structure optimization procedure. Based on these considerations, we treat QMCM and CMCM interactions at an equal footing. Thus, there are no direct electrostatic interactions between QM electrons with CM atoms. Such treatment is a good approximation for the systems that do not contain strongly polar or charged groups. It can be applied to many covalently bonded materials, such as silicon and silica, but may be problematic for those ionic metals. The QM/CM interface If the QM/CM boundary does not cross any covalent bond in the system, the interface would not require special treatments. Thus, the hybrid QM/CM method can be applied easily to problems such as solvent effects. However, in covalent materials, the QM/CM boundary must cross covalent bonds. Without any special treatment of the boundary, the dangling bonds will introduce errors to energy states such as adding energy levels in the energy gap. There are several methods to treat the QM/CM boundary. I will introduce the linkatom and the pseudoatom methods. The linkatom method One way to treat the dangling bonds at the QM/CM boundary is to introduce link atoms. See an example in the article by Antes and Thiel (1998). The linkatoms generally are hydrogen atoms. They are attached to the broken covalent bonds on the QM/CM boundary to saturate the valence. Since the electronegativity ofH is between that of Si and O, H accepts electrons from Si (as O does), and donates electrons to O (as Si does). This explains the success of Htermination for clusters of silica and zeolites. In our method, the linkatoms are included in the QM calculation, but do not interact with the CM atoms, because they are introduced only to satisfy valence considerations and are not present in the real system. The total energy can be rewritten as follows, E= EE, ({r},{rL}, p (r))+E ({r})+C EQcm ({r},{r,}) (2107) where L stands for the linkatom. If the atoms A, e {QM} and B, e {CM}, we define the constraint such that the linkatom CL is at a constant distance from An along the A B bond. The An CL bond length is optimized to minimize the force on C, in some relevant small cluster calculations. The pseudoatom method The alternative to using hydrogen atoms to terminate the dangling bonds at the boundary is to use pseudoatoms to "absorb" these bonds. Taking the silicon dioxide system as an example, a pseudooxygenatom O* can be defined to have seven valence electrons, a nuclear charge of seven, and an effective core potential. Thus, O* has just one free valence to make a bond with the neighboring Si atom inside the QM region. One can either replace the pseudoatom by a real atom, an F atom in the case of the foregoing example, or use a pseudopotential as was done by Zhang, Lee, and Yang (1999). In the pseudoatom method, the dangling bonds are "absorbed" into the pseudoatoms. However, the orientation of the dangling bond is not taken into account if simply using a real atom as the pseudoatom. If one chooses to reparameterize the pseudopotential, the bondangle dependence is still difficult to include in the pseudopotential. Thus, it is one advantage of the linkatom method since such bondangle dependence is naturally included if the linkatoms are constrained to lie on the line of the bonds that connect the QM and CM parts. CHAPTER 3 MULTISCALE SIMULATIONS OF AMORPHOUS SILICA Introduction Silica is the most abundant chemical species in the earth's crust and has been used in microelectronics, optical devices, nuclear waste containment vessels, space shuttle insulation tiles, and many other technologies. With large abundance on earth, both water and silica play fundamental roles in a vast variety of human activities, and their interaction is one of the most important problems in materials and environmental sciences. The phenomenon of hydrolytic weakening was first addressed in geoscience in 1960s (Griggs, 1967). Since then, hydrolysis of silica surfaces has been a subject attracting much attention in mineral and semiconductorceramic research. Amorphous silica is abundant in the mineral and rock. The study of waterrock interaction has taken on increased importance now that toxic materials are to be stored in geologic aquifers. The suitability of such storage is determined, in part, by predicting the chemical reactions that will take place between hypothetically released waste and the rocks in the repository. The waterrock interaction is important because it creates a surface with chemical properties that are unique from the underlying mineral. Therefore, the study of the surface reactivity of the silica is fundamentally important to the materials and environmental sciences (Hochella and White, 1990). Although many experiments have been carried out to study watersilica interaction (Zhuravlev, 1993; Zhuravlev, 2000; Bunker, et al., 1989a, Bunker, et al. 1989b, Chuang and Maciel, 1996; Chuang and Maciel, 1997; Bolis, et al., 1991; Chiang, Zegarski, and Dubois, 1993; D'Souza and Pantano, 1999), the mechanisms of some phenomena are not very well understood. For example, water can increases crack propagation speed in stressed silica by many orders of magnitude (Michalske and Bunker, 1984). This type of behavior is also an example of a chemomechanical phenomenon, in which mechanical stress and chemical reactivity amplify one another. The basic question in these phenomena is how the strained chemical bonds break under the presence of water and other chemical species. In the absence of strain, SiO bonds are known to be inert to chemicals (H20, NH3, etc) that can cause stress corrosion in silica glass. However, strained SiO bonds greatly increase reactivity by creating acidic and basic adsorption sites on silicon and oxygen. These reactive sites are believed to play crucial roles in stress corrosion and crack growth. For silica exposed to the atmosphere, the surface is covered by silanols (SiOH) with coverage of about 4.9 silanols/nm2 (Zhuravlev, 1993; Zhuravlev, 2000). (The number of Si atoms on the surface is about 8 Si/nm2.) H20 can be adsorbed upon the surface silanols via hydrogen bonds with an energy of adsorption higher than 44 kJ/mol (latent heat of liquefaction of water) (Bolis, et al., 1991). Thus the silica surface in atmosphere is hydrophilic. On the other hand, hydrophobocity can be developed through heat treatment, which results in the condensation of silanols into siloxane bridges. The dehydroxylated silica surface (heat pretreatment at temperature higher than 4000C) is hydrophobic. Further rehydration of the dehydroxylated silica surface cannot restore the silanol concentration on the surface to its value before heat treatment (Zhuravlev, 1993; Zhuravlev, 2000; D' Souza and Pantano, 1999). On the dehydroxylated silica surface, some highly strained defects consisting of a twomembered ring can be identified by IR peaks at 888 and 908 cm1. These characteristic peaks disappear upon exposure to water, ammonia, methanol and organosilanes (Bunker, et al., 1989a; Bunker, et al., 1989b). The high strain in the twomembered ring structure creates a local reactive site that is vulnerable to chemisorption of gasphase species. The reaction between water and strained SiO bonds is an important aspect of surface hydroxylation. Reliable computational models for watersilica interaction are critical for a detailed microscopic description of the silica surface hydroxylation, and the prediction of reaction pathways that lead to stress corrosion. Garofalini (1990) reported classical molecular dynamics (MD) simulations of water adsorption on the amorphous silica surface. Reaction pathways for the bond rupture and formation on the surface were investigated in that work. However, the classical description of dissociative chemisorption is not reliable and no energetic information was obtained. This is because the interatomic potential functions normally do not take into account the details of electronic behavior. Clearly, high level quantum mechanical calculations are required. While ab initio methods have been successfully used in the calculations for large clusters and crystalline surfaces, amorphous surfaces remain a challenge. The disordered nature of such surfaces requires relatively large simulation cells, which makes the direct use of ab initio methods difficult. In order to overcome this difficulty, several models have been proposed. A frequently used approach is the cluster model, in which relevant clusters are chosen from the amorphous surface and treated by ab initio methods. Walsh and coworkers (2000) studied the activation barriers of water dissociative reactions with a number of clusters that represent defect sites on the amorphous surface. For the SiO2 system, the charge transfer is local and the electronic response is shortranged. The cluster model may be safe in this respect. However, the elastic field underneath the surface also has a significant effect on the surface reconstruction due to water adsorption. Penev and coworkers (1999) suggest that smaller clusters generally overestimate the activation barrier and reaction energy. Therefore, great caution must be taken when explaining surface phenomena based on the data from cluster calculations. An alternative approach is to use a fcristobalite surface to model a hydroxylated amorphous surface (Ceresoli, et al., 2000; larori, 2001). This approach was adopted because the density and refractive index of fcristobalite silica are closest of all the crystalline phases of silica to those of amorphous silica. Also, two types of the surface silanols (geminal and single silanols) exist on the (100) and (111) surfaces of / cristobalite (Chuang and Maciel, 1997). Iarlori and coworkers (2001) have used the f cristobalite model surface to study dehydroxylation and adhesion reactions, and compared results with experimental data on the amorphous surface. The calculated dehydroxylation energy is much lower and the activation barrier is higher than the respective experimental values. The amorphous silica surface has a variety of surface sites (two, three, four, five, sixmembered rings) but the fcristobalite model surface has only one type of surface ring structure, which limits the applicability of this approach. Moreover, crystalline silica is considered to have higher surface reactivity than amorphous silica. Cristobalite is the most pathogenic polymorph of silica due to its high surface reactivity (Fubini, 1998). The chemical response to amorphous silica is much less pronounced than to crystalline polymorphs. Crystalline silica is also more hydrophilic and more difficult to dehydroxylate than amorphous silica (Fubini, 1998). The longrange order and hence longrange surface response in the flcristobalite model do not exist in amorphous silica. In short, there is not an adequate rationalization for assuming that the surface reactivity and adhesion properties of the amorphous silica surface resemble those on the fcristobalite surface. In the last decade, hybrid methods that integrate quantum mechanics and classical interatomic force field have been investigated in great detail and used extensively in the studies of reaction mechanism of enzyme and protein, solvent effects, catalysis, surface chemistry, crack propagation, etc (Field, Bash, and Kurplus, 1990; Zhang, Lee, and Yang, 1999; Eichinger et al., 1999; Sauer and Sierka, 2000; Broughton et al., 1999; Ogata et al., 2001). The problems studied with multiscale schemes usually involve local chemical processes in systems that do not have longrange order. Apparently, multiscale methods are well suited to the studies of amorphous materials as well. Therefore, we have developed a combined quantum mechanical (QM) and classical mechanical (CM) method to study phenomena related to amorphous silica (see Chapter 2, or Du and Cheng, 2003a). The advantage of such a multiscale method is that it can more realistically simulate chemical reactions involving surface atoms (treated by QM method), which are in an elastic field provided by the surrounding bulk material (treated by simpler CM method). Its accuracy can be systematically improved by increasing the size of the QM region. We are interested in the chemomechanical phenomena in amorphous silica, in particular, stressinduced corrosion. As the first step, we chose a particular site (the two membered ring) on the amorphous silica surface, and studied both physisorption and chemisorption of water on this site. In the following part of this chapter, I will first show the results of test of the QM/CM interface (the linkatom and the pseudoatom methods) on small silica clusters (Du and Cheng, 2003a), and then apply the combined QM/CM method on the studies of hydrolysis of the twomembered ring on the amorphous silica surface (Du, Kolchin, and Cheng, 2003; Du, M.,H, Kolchin, A., and Cheng, H,P, 2004, J. Chem. Phys., in press). Test of the LinkAtom Method A small silica cluster (Si207H6, see Fig. 31) is chosen as a training system to test the combined QM/CM method. This cluster consists of two tetrahedrons, which are building blocks of amorphous silica. Amorphous silica can be thought of as a random network of such tetrahedral units. H 0 0 O o i o s i 0 SiH Figure 31. The structure of the Si207H6 cluster. The gray, dark gray, and light gray spheres represent the Si, 0, and H atoms, respectively. The system is then divided into two regions: the QM region on the left, including the bridging oxygen atom, and the CM region on the right (the middle panel of Fig. 32). The linkatom (H atom) is placed between the bridging O(b) atom and the Si(2) on the right, and is only allowed to move on the line of the O(b)Si(2) bond. The resulting O(b) H(link) bond length after relaxation is 1.82 bohr. In the calculations of other silica systems, the linkatoms will be constrained on the line of O( {QM})Si(e {CM} ) bonds with a fixed O(e {QM})H(link) bond length of 1.82 bohr. The BeestKramerSanten (BKS) potential is chosen to describe the classical interatomic interactions (van Beest, Kramer, van Santen, 1990). For atoms i andj, the BKS potential has the form of q,q exp( D A exp(b)ii r r (31) where r, is the distance between the vi andjth atoms. QM H 0 3.103 3.028 0 H 1 3.028 H o Si(1) O(b) Si(2) H 07.05 H 0 O H (a) IQM/CM I H O0 H O H O QM 3.103 3.046 3.154  Si(1) O(b) 0 .38 H S109.38 CM 3.059 Si(2) 1ink) 109.700 (b) (c) 3.060 X060 3.144 0 Si(1) O(b) Si(2) 09.680 Figure 32. The structures of the cluster shown in Fig. 31, obtained from (a) allquantum calculation, (b) combined QM/CM calculation using the linkatom method, and (c) allclassical calculation. The optimized structures obtained from the allquantum calculation, combined QM/CM calculation, and allclassical calculation, are shown in Fig. 32. The Si(1)O(b) bond length is 3.03, 3.05, and 3.14 bohr for the QM, combined QM/CM and CM methods, respectively. The corresponding Si(2)O(b) bond length is 3.03, 3.15, and 3.14 bohr for three different methods. These bond lengths, along with other structure parameters shown in Fig. 31, indicate that the QM part of the hybrid system has a structure very similar to that in the allquantum calculation, and the corresponding CM part has a structure nearly identical to that in the allclassical calculation. The difference between the structures of the QM and CM parts of the hybrid system is mainly due to the difference between the QM and CM methods. In other words, the interface is of a high quality that is largely transparent. The structure of the training system with isocharge density surfaces is shown in Fig. 33. If including the electrostatic interaction between electrons in the QM region and atoms in the CM region [see Eq. (2106)], the partial charges [q, and q, in Eq. (31)] assigned to CM atoms in Fig. 31 must be reduced to decrease the unrealistic large electrostatic potential applied to the QM electrons. Such reduction of the CM point charges is rather arbitrary and has no criterion for convergence. Since the linkatoms can terminate the dangling bonds very well and approximately restore the correct charge distribution at the QM/CM boundary (see Fig. 33), the additional charge polarization induced by CM point charges is not necessary. Thus, we did not use Eq. (2106) to include the electrostatic interaction between QM electrons and CM atoms. S (a) H H(link) O Si (b) Figure 33. The structures of the training system with isodensity surfaces, obtained from (a) allquantum calculation, and (b) combined QM/CM calculation. Fig. 34 shows the total energy change due to the change of the Si(l)O(b) or Si(2) O(b) bond length. When the SiO bond in the quantum region varies, the energy curve obtained from the combined QM/CM calculations is very close to that from the all quantum calculations. On the other hand, if the SiO bond in the CM region is changed this energy curve is close to that from the allclassical calculations. This is an indication of a transparent quantumclassical interface. A more rigorous test of the interface is shown in Figures 35 and 36. Here, we show the forces acting on several atoms near the QM/CM boundary as functions of the Si(1)O(b) or Si(2)O(b) bond length. Again, a consistent pattern can be recognized from these figures. When the Si(1)O(b) bond length varies, the values of force acting on the Si(1) and O(b) atoms, obtained from the combined QM/CM calculations, are both very close those obtained from allquantum calculations. When the Si(2)O(b) bond length varies, the values of force acting on the Si(2) and O(b) atoms, obtained from the combined QM/CM method, are nearly identical to those from the allclassical calculations. Based on these tests (Figures 34, 35 and 36), one can see that the interface is transparent from both the quantum side and the classical side. To the atoms in the CM region, the atoms in the QM region behave like classical atoms. This can be anticipated from the treatment of the QM/CM interactions (see Eq. 299 and 2100). On the other hand, to the quantum atoms, which consists of ions and electrons, the classical atoms plus the linkatoms imitate an environment that is similar enough to the one provided by all quantum atoms. The difference among the three curves in any of the figures 3436 arises from the disagreement between the quantum and classical treatments of the interactions. The Si(1)O(b)Si(2) bond angle in Si207H6 cluster (Fig. 32) is 1800. However, the average SiOSi bond angle in amorphous silica is around 1460. Thus, it is necessary to further test the validity of the linkatom method in other silica clusters with different chemical environments. 0.05 ...  QM 0.04 QM/CM  CM / 0.03 0.02 ,' " (a) 0.06 0.05  QM QMICM 0.04 0.03  0.02  0.01  0.4 0.2 0 0.2 0.4 0.6 r (bohr) (b) Figure 34. The total energy change due to the change of (a) the Si(l)O(b) bond length or (b) the Si(2)O(b) bond length. The longdash red line, solid black line, and shortdash red line represent the results from the QM, combined QM/CM, and CM calculations, respectively. For each of these three methods, the energy of the fullyrelaxed cluster is shifted to zero. 0.3 SQM 0 \ QM/CM S 0.2  CM S 0.1 ",, 0  0.1 0.4 0.2 0 0.2 0.4 0.6 r (bohr) (a) QM 6. QM/CM  CM 4 0.1 / t l 0.3 0.4 0.2 0 0.2 0.4 0.6 r (bohr) (b) Figure 35. The force on (a) Si(1l) and (b) O(b), when the Si(l)O(b) bond length varies. The longdash red line, solid black line, and shortdash red line represent the results from the QM, combined QM/CM, and CM calculations, respectively. 0.2 0.2 0 r 0 0.2 r (bohr) 0.2 (bohr) (b) Figure 36. The force on (a) Si(1) and (b) O(b), when the Si(1)O(b) bond length varies. The longdash red line, solid black line, and shortdash red line represent the results from the QM, combined QM/CM, and CM calculations, respectively.  QM QM/CM  CM / , 0.4 0.1 0.1 0.3 0.4 0.2 0.1 0 0.1 0.2 0.4 Figure 37. The structure of the Si309H6 cluster. The gray, dark gray, and light gray spheres represent the Si, 0, and H atoms, respectively. To test the transferability of the QM/CM interface, the linkatom method is applied to a threemembered silica ring (Si309H6, see Fig. 37). The comparison of the optimized structures obtained from the QM, combined QM/CM, and CM methods again demonstrates the transparency of the interface (see Fig. 38), i.e., the structure parameters, generated by the combined QM/CM method, are in good agreement with the results from the allquantum calculations in the QM region, as well as with the results from the allclassical calculations in the CM region. The tests in small silica clusters demonstrate that the linkatom method can couple the QM and CM regions very well for silicon dioxide systems. The reason for this success comes mainly from the local nature of the SiO bonding and the proper construction of the Hamiltonian based on that property. QPM H H I I 0 0 S134.00 .10 3.103 (a) O Si H / H Si.o 0 H H H H I I Si 3.113 O O 0/ 10431 ,o / 137 3.19 H(link) H(link) 3.139 o.i, 99A4 o Si 0 0 1425 90 0 03.129 0 QM CM 0 0 Si 0 9999 140.01 14 3.129 1 SiSi O / o \ Figure 38. The structures of the cluster shown in Fig. 37, obtained from (a) allquantum calculation, (b) combined QM/CM calculation using the linkatom method, and (c) allclassical calculation (bottom panel). In (a) and (c), the lengths of the SiO bonds on the plane of the ring are identical and this ring has a three fold rotational symmetry. QM/CM (b) CM (c) * 0 Test of the PseudoAtom Method As already discussed in this method, the dangling bond at the QM/CM boundary is "absorbed" into a pseudoatom 0*, which has seven valence electrons, a nuclear charge of seven, and an effective core potential. The simplest scheme is to replace the oxygen atom at the QM/CM boundary by a fluorine atom. The interface is again tested on Si207H6 and Si309H6 clusters, as already done for the linkatom method. The optimized structures of these two clusters, obtained from the combined QM/CM method, are shown in Fig. 39. For Si207H6, the interface works very well, as compared with Fig. 32. The structure in the QM region is in good agreement with that from allquantum calculations, and the structure in the CM region is very similar to that from allclassical calculations. Thus, the pseudooxygenatom properly terminates the dangling bond on the boundary. However, a discrepancy appears for the Si309H6 cluster, as seen by comparison of Fig. 39 with Fig. 38. The difference in the optimized structures is mainly from the bond angle. For example, the SiO*Si bond angle obtained from the combined QM/CM method is 149.1 (see the bottom panel of Fig. 39), which is much bigger than the corresponding bond angle (134.00) from the allquantum calculation (see the upper panel of Fig. 38). Such disagreement arises because no bond angle dependence is built into the current pseudoatom method. Thus, the electrons near the QM/CM boundary do not feel any difference between structures with different SiO Si bond angle. This explains why such an interface works well for the Si207H6 cluster, in which the SiOSi bond angle is 1800 at the QM/CM boundary, but fails for the Si309H6 cluster. Varying the nuclear charges of O* atoms improves the SiO* bond lengths, but fails to improve the SiO*Si bond angles for the same reason. Thus, a SiO*Si bond angle dependent O* effective core potential must be designed if the pseudoatom method is to be used. The construction of a transferable O* effective potential in a dynamically evolving SiO2 system is rather complicated. Fortunately, such bondangle dependence is naturally included in the linkatom method if the linkatoms are constrained on the lines of the Si O bond on the QM/CM boundary. Based on the preceding considerations, we have used the linkatom method, which is very simple but effective as well, in the studies of the hydrolytic weakening of the amorphous silica. IQM/CM QM H \3074 x 3.045 3.145 H 0 Si(1) (b) H / 10648 IQM/CM CM 3062  Si(2) 109 81 H H I I 0 0 Si 3064 0 100.893 0 /13.62 ' \ 0 0 "10 3': 124 0 QM CM Figure 39. Cluster structures obtained from combined QM/CM calculations using the pseudoatom method. (a) Cluster shown in Fig. 31. (b) Cluster shown in Fig. 37. (a) (b) Hydrolysis of a TwoMembered Silica Ring on the Amorphous Silica Surface Generation of the Amorphous Silica Surface The amorphous silica sample was prepared by cooling a preheated flcristobalite sample via a classical MD simulation that uses the BKS potential (Eq. 31) with a correction to suppress an unphysical attraction at short interatomic distance. A bulk sample of/fcristobalite, consisting of 3000 atoms with periodic boundary conditions, first was heated and equilibrated at 8000 K using the NVT ensemble. Then the system was annealed continuously down to 300 K following a schedule suggested by Huff et al (1999, Table 3, cycle 2VIII). Fig. 310 shows the MD results for the 00, SiO and SiSi pairdistribution functions in the amorphous silica bulk. The positions of the first peak in gsio (r), goo (r) and gsisi (r) give the mean SiO bond length, nearestneighbor 00 and SiSi distances to be 1.61 A, 2.60 A and 3.09 A, respectively. The corresponding experimental values from neutron diffraction data are 1.608 A, 2.626 A, and 3.077 A (Wright, 1994). The density, average SiOSi angle and average OSiO angle of the resulting amorphous silica sample are 2.25 g/cm3, 146.40 and 109.30, respectively, which are in good agreement with experimental data and previous simulations (Huff, et al., 1999; Wright, 1994). A vacuum gap then was inserted in the zdirection of the periodic cell to generate a surface, which was relaxed further at 300 K using the NPT ensemble. The resulting surface was covered by silica rings with various sizes. The top layer is shown in Fig. 3 11. Six and fivemembered rings are the dominant species. There are also some larger rings and some smaller ones such as four, three, and twomembered rings on the surface. In general, the small rings are more strained and rigid. Consequently, all two membered rings are found to be perpendicular to and sticking a bit out of the surface to reduce the local strain. In addition, there are also some onecoordinated O atoms hanging on the surface. These defect sites (dangling bonds and twomembered rings) can be annihilated easily by water attack. Fig. 312 shows the density profile of atoms as a function of the distance from the center of the slab. The location of the original surface (before the relaxation) is indicated in the figure. The relaxed surface becomes a little rough, resulting in a decrease of density near the edge of the surface. There is a relatively higherdensity region on the surface at around 15 A. For crystals, such a feature is usually due to lattice contraction on the surface. However, for the amorphous silica, this highdensity region is due to the enrichment of small silica rings (two, three and fourmembered rings) on the surface. The SiO bond lengths in small rings are usually longer than those in bigger rings because of the enhanced strain in small rings. Thus, the highdensity region on the surface does not correspond to shorter bond lengths in the same region, as shown in Fig. 313(a). Instead, the average SiO bond length in the highdensity region is even slightly longer than that for the entire surface slab. However, small rings are more compact as compared to large rings, giving rise to the relatively higher density. Note that the sharp decline of the average bond length at the top of the surface [Fig. 313(a)] is due to those SiO bonds with onecoordinated O atoms. The enrichment of small rings on the surface also is evidenced by the decreasing average SiOSi bond angle on the surface, as shown in Fig. 313(b). The sharp decline at the top of the surface in Fig. 313(b) is due to those twomembered rings, which have smallest SiOSi bond angles. The concentration of twomembered rings on the surface is found to be 0.42/nm2, which is close to the experimental estimate (0.20.4/nm2) (Grabbe et al., 1995) 56 6 1. 4 0 0 2 0 A*i*i  25 S20 15  o 5 10  0 5 0 .* * 4 L_ S2 0 2 4 6 8 10 12 14 r(A) Figure 310. Partial pairdistribution functions for amorphous silica at 300 K. Figure 311. One layer of the amorphous silica surface. The gray and dark gray spheres represent the Si and O atoms, respectively. 0 5 10 15 20 z (A) Figure 312. Number density of atoms as a function of the distance from the center of the slab. The dash line indicates the original surface before the relaxation. 1.65 1.60 1.55 160 Cj C 140 CO 120 n 100 (a) (b) 0 5 10 15 20 z A) Figure 313. (a) The average SiO bond length, and (b) the average SiOSi bond angle, as functions of the distance from the center of the slab. The slab is divided into many layers, each of which has thickness of 1.22 A. The bond lengths and angles are averaged in each of these layers. Hybrid QM/CM Approach As already discussed, the hybrid scheme divides the system into QM and CM regions. The QM region is treated by pseudopotential density functional calculations using the generalized gradient approximation (GGA) and planewave basis sets (Barnett and Landman, 1993). We used the TroullierMartins pseudopotential (Troullier and Martins, 1993) and the PerdewBurkeErnzerhof (PBE) exchangecorrelation functional (Perdew, Burke, and Ernzerhof, 1996). The cutoff energy for the expansion of planewave basis is 31 Hartrees. The interactions between QM and CM regions are described by classical interatomic potential functions, which are also used to describe interactions between atoms in the CM region. See the discussion in the preceding chapter. The linkatom approach was used to saturate the valence structure in the QM region (i.e. a hydrogen atom is constrained on the line of a SiO bond at the QM/CM boundary to terminate a dangling bond). The formalism is shown in Chapter 2. The test calculations of ground state structures and bond strength on small clusters show that the linkatom method works very well for silica systems (Du and Cheng, 2003a). Two types of interatomic potential functions were used in our calculations for different purposes. The BKS potential (Eq. 31) is known to generate good lattice parameters and elastic constants for quartz and other SiO2 polymorphs (van Beest, Kramer, and van Santen, 1990). It has also been used to generate reasonable amorphous structures. We expect the BKS potential to provide a reasonably good elastic field in amorphous silica. However, it greatly overestimates the binding energy. Watanabe et al developed a StillingerWeber type potential for silica (Watanabe, et al., 1999), in which parameters are fitted to reproduce the cohesive energies of SiO bonds as a function of the coordination number of the O atom and also the total and deformation energies of a number of silica clusters. In the Watanabe potential for Si, O systems, the total interaction energy is given by S= C f ,(i,j)+I j f3(i, j, k) (32) I J>1 i j>i k>j where s (=50 kcal/mol) is an energy parameter,f2 is the pairwise interaction between ith andjth atoms, andf3 is the threebody term as a function of the i,j, kth atoms. The Watanabe potential gives much better energetic than does the BKS potential. However, the Watanabe potential has not been shown to produce reasonable crystalline and amorphous structures. From this comparison of these two potentials, we adopted a mixed scheme, in which the BKS potential was used in structure optimization calculations and the Watanabe potential was used in the energy calculation of those optimized structures. Thus, in structure optimization calculations, the interaction energy among the classical atoms, U((Rs}), was given by U({R, }) = (,,, (rs,) (33) s s'>s and the interaction energy between QM and CM atoms, V({Rn},{Rs}), was given by V({R, ,{RJ)= II ,(r.) (34) n s The forces acting on QM and CM atoms follow from by Eq. 2104 and 2105. After the forces were minimized, the total energy was calculated with the Watanabe potential. U({Rs}) and V({Rn},{Rs}) are given by U({fRs}) = I I f (rs,,) + I II Ff, (r,, r,,,r,,,,,) (35) s s>s s s s'>s s">s' V({R},{R,})= f2 fZ(rs) (36) n s s'>S n n'>n s Simulation Setup A 12000atom amorphous silica slab, consisting of four 3000atom unit cells generated by the method described previously, was used as a model surface slab. Since this slab was big enough, the periodic boundary condition was removed. While 10000 atoms at the bottom and four sides were held fixed in position, 2000 atoms at the top center were relaxed to minimize the force on each atom using the BKS potential (Eq. 3 1). A QM region that contains 31 Si and O atoms was then chosen on the surface and further relaxed by minimizing the forces on QM nuclei (as given by Eq. 2104), while the positions of all the CM atoms were fixed to simplify the problem. [The division of the system is depicted in Fig. 314(a).] Finally, the Watanabe potential was used to re calculate the interaction energy between the QM and CM atoms. Since the positions of all the CM atoms are fixed when relaxing the QM region, the total energy of the system can be calculated as a sum of EQM and EQM/CM. The interactions between CM atoms were excluded from the total energy. The water molecules were included later in the QM region to calculate adsorption energies, which are determined according to the expression Ead.orphon = E(surface + water) n x E(water) E(surface) (37) where n is the number of water molecules. A negative Eadsorpton value corresponds to a stable structure. The interactions between water and CM atoms were neglected in our calculations. In order to check the error in energetic introduced by the less accurate classical potential, neglect of the interactions between water and CM atoms, and disrupted charge transfer on the QM/CM boundary, the adsorption energies were calculated with two larger QM regions of 74 and 209 surface atoms, respectively [Fig. 314(b) and (c)]. In the rest of this chapter, the systems with 31, 74 and 209 surface atoms in QM regions are referred as system I, II and III, respectively. The QM subsystem of system I is relaxed to its ground state. In the calculations of water adsorption energies for system II and III, the optimized structure of system I is used. At this stage, no structural relaxations of system II and III are considered. We have estimated the energetic error of system I in a fixed configuration. The artificial strain in system I introduced by fixing all the CM atoms can be reduced by relaxation of more surface atoms. Figure 314. A schematic picture of a silica surface (side view). Three upper panels (a), (b) and (c) are enlarged pictures of the QM regions in system I, II and III, respectively. The gray, dark gray, and light gray spheres represent the Si, O, and H atoms, respectively. Results Adsorption of One Water Molecule on the TwoMembered Ring The calculated adsorption configuration of H20 on the amorphous silica surface is depicted in Fig. 315 (Only the 31 QM atoms on the surface are shown). In the upper left picture in Fig. 315, the water molecule is clearly bound to the surface. The Si(1)O(1), Si(1)O(2) and Si(1)O(3) bond lengths are 1.84, 1.81 and 1.82 A, respectively, as shown in Table 31. The average SiO bond length in amorphous silica is 1.61 A. In system I, the energies of 31 surface atoms in QM region calculated based on the atomic configurations before and after water physisorption, and the energy of the physisorbed water molecule are calculated as Eb(surf), E(surf) and E(water), respectively. The local surface energy increase is estimated as E(surf) Eb(surf) = 1.834 eV. The bonding strength between physisorbed water molecule and the surface can be estimated as EY(surf +water) E(surf) E"(water) = 1.846 eV, where E"(surf+ water) is the energy of the 31 surface atoms with the adsorbed water molecule. While the water molecule is bound to the surface, the bonding strength of local SiO network is weakened, which is evidence of the hydrolytic weakening. Table 31. Several SiO bond lengths in combined QM/CM model and cluster model (cluster B). The units are in A. Si(l)O(l) Si(l)0(2) Si(l)0(3) Combined QM/CM 1.84 1.81 1.82 surface model Cluster model 2.56 1.70 1.71 To compare our hybrid QM/CM approach with the cluster model, we have also studied the adsorption of water on two silica clusters as shown in Fig. 316(a) and (b). One is Si206H4 (cluster A), which is a twomembered ring. The other one is a 31atom silica cluster terminated by H atoms (cluster B), which is the QM subsystem in system For cluster B, the Si atom does not make a bond with the O atom in H20 (Si(1)O(1) distanceis 2.56 A), while our combined QM/CM model predicts such a bond. Some SiO distances are compared with those from the combined QM/CM surface model in Table 31. H Oi/ Si, 9  11589.0 H 11590.0 s ) 03 C 11590.5 11591.0 11591.5 0 10 20 30 40 50 60 Step Figure 315. Snapshots and potential energy of QM atoms along the dissociation path of one H20 molecule on a silica surface. Three pictures correspond to physisorption, transition and chemisorption states. The gray, dark gray, and light gray spheres represent the Si, 0 and H atoms, respectively. The valence electron density contours from the cluster model (cluster B) and the combined QM/CM model are compared in Fig. 317, 318 and 319. All density contours are drawn on the plane crossing the Si(1), 0(1) and 0(2) atoms. As shown in these figures, the watersurface interaction is significantly stronger in the combined QM/CM surface model than in the cluster model. In the combined QM/CM model, the water molecule is strongly perturbed. Electrons transfer from the 0(1) atom to its adjacent Si(1)O(1) and 0(1)H bond, and also from the H20 and Si(1) to 0(2). The surface electrons are redistributed to accumulate in the direction toward the water molecule. However, charge transfer in the cluster model is much less significant. Mulliken charge analysis shows 0.3 e charge transfer from the water molecule to the surface in the combined QM/CM model, but only 0.06 e charge transfer in the cluster model. (a) (b) vi (c) Figure 316. Cluster models: (a) one H20 molecule bound to a Si206H4 cluster, (b) one H20 molecule adsorbed on a silica cluster (31 Si and O atoms with dangling bonds saturated by 14 H atoms) and (c) two H20 molecules adsorbed on the same silica cluster in (b). The gray, dark gray, and light gray spheres represent the Si, 0, and H atoms, respectively. 66 01 1 140 0.025 A 10 3Sil V ec 2 05 0 0.4010 02 0$\ 1 0.400 0.10 00.50 1 0 1 2 x (A) (b) Figure 317: Valence electron number density contours on the plane crossing Si 1, 1 and 02 atoms in the combined QM/CM surface model (a) and the cluster model (b). Each contour represents twice/half the density of the adjacent contour lines. The values shown in figures are in units of e/A3. >02 lines. The values shown in figures are in units of e/8+3 67 0 1 .0E040 .OE3 2 R0.002 02 .OE3 Sil 1 0 1 2 x (A) (a) 3 01 1 0 1 2 3 X(A) (b) Figure 318. Contour plot of (a) excess electron density An' and (b) deficiency electron density An in the watersurface system relative to superposed electron densities of the surface and H20 on the plane crossing Si(1), 0(1) and 0(2) atoms in the combined QM/CM surface model. (An = n(surface + H20) n(surface) n(H20) ) Each contour represents twice/half the density of the adjacent contour lines. The values shown in figures are in units of e/A3. ____'__ 10 1 x (A) (a) 2 3 x (A) (b) Figure 319. Contour plot of (a) excess electron density An' and (b) deficiency electron density An in the watercluster system relative to superposed electron densities of the cluster and H20 on the plane crossing Si(1), 0(1) and 0(2) atoms in the cluster model. (An = n(cluster + HO) n(cluster) n(HO)) Each contour represents twice/half the density of the adjacent contour lines. The values shown in figures are in units of e/A3. 01 i.OE Si1 02  The adsorption energies are shown in Table 32. Following the procedure used above for analyzing both the local surface energy increase and the surfacewater bonding strength, the cluster energy increases by 0.14 eV and the clusterwater bonding strength is 0.12 eV (The H atoms that saturate dangling bonds of the cluster are fixed in this calculation). Comparing the energetic obtained in the combined QM/CM and cluster models, the cluster model predicts much weaker water bonding strength and less hydrolytic weakening effect. The significant structural difference implies a much higher reaction barrier for water dissociation in the cluster model than in the combined QM/CM model. Walsh et al reported such a reaction barrier ranging from 0.7 to 1.1 eV based on a small cluster calculation (Walsh et al., 2000). Our calculations show that the upper bound of the water dissociation barrier is 0.4 eV. The complete reaction pathway is shown in Fig. 315. According to experimental observations, the twomembered rings on the silica surface are annihilated upon water adsorption at room temperature (Bunker, 1989a; Bunker, 1989b), which suggests a significant interaction between water and two membered rings. The hybrid QM/CM method that takes into account the elastic field clearly produces more reasonable results than cluster model does. Table 32. Adsorption energies of one and two H20 molecules in the cluster model with or without fixing the positions ofH atoms that saturate the dangling bonds of the cluster. Si206H4 31atom 31atom cluster 31atom cluster cluster (free) (constrained) (constrained with 2 H20) Eadsorption(eV) 0.5 0.12 0.05 0.23 The reaction energies for physorption and chemisorption of H20 molecules are included in Table 33. The errors in adsorption energies in system I come from two sources: (1) disrupted charge transfer on the QM/CM boundary and neglect of the interactions between the water and CM atoms, (2) the artificial strain due to fixing the positions of CM atoms. By increasing the size of the QM region we can correct the first error. To reduce the second error, we must relax more surface atoms, which we have not done thus far. The first error in adsorption energies in system I can be evaluated as the difference between QM contributions to adsorption energies in system I and III (The QM region in system II is not sufficiently larger than that in system I for that task. Some boundary atoms in system I are also on the QM/CM boundary of system II. Hence, system II cannot be used to evaluate errors in system I). The errors from source (1) for physisorption of one H20 molecule and chemisorption of one and two H20 molecules are 0.125, 0.196 and 0.369 eV, respectively. Note that with the Watanabe potential the contributions from the interactions between QM and CM atoms to the adsorption energies in system III are zero for the following reasons: 1) Only the 31 QM atoms in the reaction center change their positions upon water adsorption in the simulations and the distances between these 31 atoms and CM atoms exceed the cutoff distance of the Watanabe potential. 2) For those QM atoms that do not change positions upon water adsorption, their contributions to adsorption energies are canceled out when subtracting energies according to Eq. 37. Within the limit of the Watanabe potential, the adsorption energies are solely from the contribution of QM atoms in system III and further increasing the size of QM region without relaxation of more atoms will not change the results much. The relaxation of more atoms may further lower the adsorption energies. For comparison, we also included the energetic results using the BKS potential in Table 33. The water adsorption energies are substantially overestimated if using the BKS potential. Apparently the adsorption energies have not reached convergence as the size of the QM region increases. Although the Watanabe potential is better than the BKS potential in terms of the energy calculation, its accuracy still needs improvement. The relatively low quality of the classical potential causes the oscillation and slow convergence of the adsorption energies. Table 33: Adsorption energies of water physisorbed and chemisorbed on a two membered ring on the amorphous silica surface. All the units are in eV. CMd QM Eads EM M EadsQM EadsQM/CM Eads Eads Eads (Watanabe) (BKS) (Watanabe) (BKS) Physisorption 0.097 0.595 1.337 0.69 1.24 System I Chemisorption 0.863 0.254 4.937 1.12 5.80 (1 H20) Chemisorption 1.234 0.228 4.653 1.46 5.88 (2 H20) Physisorption 0.229 0.088 1.178 0.32 0.95 System II Chemisorption 0.748 0.816 4.209 1.56 4.96 (1 H20) Chemisorption 1.142 0.858 3.827 2.00 4.97 (2 H20) Physisorption 0.028 0 1.020 0.03 1.05 System III Chemisorption 1.059 0 2.498 1.06 3.56 (1 H20) Chemisorption 1.603 0 2.243 1.60 3.85 (2 H20) The reaction energy of onewater chemisorption is 1.06 eV in our calculation. In principle, further relaxation of more surface atoms should lower the reaction energy. The experimental data on water adsorption on amorphous silica surface are largely dependent on how the sample is prepared and the pretreatment temperature. The water adsorption energy of the dry surface obtained from MD simulations should be compared with experimental results with high pretreatment temperature (which means fewer hydroxyls on the surface). Bolis et al reported the water dissociative adsorption energy as about  100 kJ/mol (1.04 eV) for a silica sample pretreated at 1073 K. Note that on a dehydroxylated surface, water can only be dissociatively adsorbed to strained SiO bonds and remaining hydroxyls. Our result is in reasonable agreement with these experimental data. Adsorption of Two Water Molecules on the TwoMembered Ring Cooperative effects (more than one water molecule participates in the reaction) have been reported on the hydrolysis reaction of alumina clusters (Witbrodt, Hase, and Schlegel, 1998), the Si (001) surface (Akagi and Tsukada, 1999), and Si02(H20)n clusters (Cheng, Barnett, and Landman, 2002). Substantial reduction of the reaction barrier was observed when additional water molecules were included in the simulations (Witbrodt, Hase, and Schlegel, 1998; Akagi and Tsukada, 1999; Cheng, Barnett, and Landman, 2002). Garofalini also observed cooperative events on the amorphous silica surface in classical MD simulations (Garofalini, 1990). Encouraged by the previous simulation results, we have used the combined QM/CM method to study the reaction energy and pathway of two H20 molecules interacting with the amorphous silica surface. The complete reaction pathway is shown in Fig. 320 and the reaction energy is shown in Table 33. A very interesting reaction pathway (hydrogen relay dissociation) is observed. Water molecules spontaneously dissociate and undergo a double hydrogen atom transfer process that breaks the siloxane bridge. The charge state of H30 on top of the twomembered ring shown in the upper middle panel of the Fig. 320 is +0.7 obtained using Mulliken charge analysis. The final configuration (bottomright panel in Fig. 320) consists of one H20 molecule with two silanols on the surface. The same reaction pattern has also been observed on the (0001) surface of aA1203 and the Si (001) surface (Witbrodt, Hase, and Schlegel, 1998; Akagi and Tsukada, 1999). The second water molecule acts as a catalyst that assists migrating a hydrogen atom from the first water molecule to the surface via a hydrogen relay process. However, in the cluster model such a reaction is not observed. (The optimized structure is shown in Fig. 316c and the adsorption energy is shown in Table 31.) The waterdimer cooperative reaction may be also energetically favored for other strained SiO bonds on the amorphous silica surface. Strained structures like a two menbered ring are more likely to buckle and may protrude on the surface so that water molecules can have better contact with the SiO bonds and cooperative reactions may take place. In the simulations, we merely minimize the energy to the ground state. The thermodynamic motion of water molecules is more complicated. Garofalini observed chain reaction of bond rupture that took place over a few different reactive sites in MD simulations (Garofalini, 1990). This observation suggests that more complex cooperative events that involve more than two water molecules may take place at room temperature. To study such phenomena, the QM region must include a larger surface area, where a hydrogenbonding network can exist among several water molecules to assist longer proton migration path. If a defect is not too far away from other defects or silanols, such longer proton migration path is possible. Adsorption of Water Molecules on the Silanols After the siloxane bridge breaks due to water attack, the local strain is reduced and the twomembered ring evolves into two silanols. Without further applying strain, this silanolcovered structure is stable under the presence of additional water. In system III, the adsorption energies for one and two H20 molecules on two surface silanols are 0.54 and 0.73 eV, respectively, and no bond rupture occurs. The configuration of two H20 molecules on top of two silanols is shown in Fig. 321. Each H20 molecule is bound to one silanol via a hydrogen bond, which is the origin of the enhancement of the hydrophilic character near silanol sites. H20 03 H2H 03 12056.5 \H4V4 12057.0 HH O H3 ( 12057.5 02 P 12058.0 0 12058.5 Si1 12059.0 0 50 100 150 200 Step Figure 320. Snapshots and potential energy of the QM atoms along the water dissociation path ("hydrogen relay process") on a silica surface. The gray, dark gray, and light gray spheres represent the Si, O, and H atoms, respectively. On a dehydroxylated surface, isolated silanols still exist (Zhuravlev, 1993; Zhuravlev, 2000). While most of the area on the dehydroxylated surface is hydrophobic, the remaining silanols and surface defects, such as twomembered rings, act as centers of adsorption upon rehydroxylation of the surface. The edges of these nucleation centers tend to trap water molecules and promote their dissociation. Subsequent rehydroxylation gradually expands the hydroxylated region on the surface. Such hydroxylated islands have been observed in STM studies of the Si surface (Chander et al., 1992; Anderson and Kohler, 1993). Akagi and Tsukada (1999) reported that two water molecules can dissociate on the Si (001) surface via a similar "hydrogen relay process" observed on the silica surface in this work, and also further adsorption of a water dimer on the site near the partially modified surface sites (SiH and SiOH) can be energetically more favored than on the clean surface. We expect that similar phenomena can happen on a dehydroxylated silica surface. Figure 321. Two H20 molecules on top of two silanols. The gray, dark gray, and light gray spheres represent the Si, 0, and H atoms, respectively. For silica samples pretreated at a temperature below 4000, complete rehydroxylation can be achieved (Silanol concentration can be restored to the original L value). (Zhuravlev, 1993; Zhuravlev, 2000; D'Souza and Pantano, 1999), since there are still plenty of silanols on the surface to adsorb water and assist the migration of H atoms. On the other hand, if silica samples are pretreated at higher temperature, the OH groups still exist but are separated by large distances. The rehydroxylation process becomes very slow and only partial rehydroxylation takes place (Zhuravlev, 1993; Zhuravlev, 2000; D' Souza and Pantano, 1999). Zhuravlev reported that it took about 5 years to complete rehydroxylation of a silica surface that was pretreated at 9000 (Zhuravlev, 1993; Zhuravlev, 2000). Discussion The combined QM/CM method is a very promising tool for studies of amorphous materials. The accuracy of the method can be improved systematically by increasing the size of the QM region, while such a feature is lacking in the flcristobalite model. For amorphous silica, charge transfer is local and the electronic response is shortranged. A simple linkatom method can serve as a good QM/CM interface, which simplifies the implementation of such a method. A major drawback in the current method is the lack of a good classical force field. Most force fields for silica are created to generate good equilibrium structures. The energetic information is not very well integrated into the potential functions. In addition, classical forces usually are correct only around the equilibrium distances. Larger stretching of the bonds cannot be properly described by these force fields. The SiO bonds described by DFT are stiffer than those described by the BKS potential. Therefore, the structure optimization and molecular dynamical simulations that include both QM and CM atoms as dynamic atoms (no fixing of the atomic positions of CM atoms) may result in undesired effects. Progress has been made to encode first principles cluster and crystalline information in semiempirical potentials for silica to ensure the sufficient transfer of information across the QM/CM interface in a chemically realistic multiscale simulation (AlDerzi, A. R., Cory, M. G., Runge, K., and Trickey, S. B., 2003, submitted; Flocke, N., Zhu, W. and Trickey, S. B., 2003, submitted, Zhu, W., Taylor, C. E., AlDerzi, A. R., Runge, K., Trickey, S. B., Zhu, T., Li, J., and Yip, S., 2003, submitted.). The potential energy functions have been fitted to not only structures, but also vibrational frequencies and forces obtained from aquartz and a set of clusters that represent the local environment of amorphous and crystalline silica. However, much effort is still needed to generate a good force field that gives accurate descriptions of structures, forces and energetic of silica systems. In addition to classical potentials, the small QM region we can currently handle also limits the applicability of the hybrid QM/CM methods in largescale dynamical problems, such as crack propagation, which need hundreds of atoms in the QM region. Previous multiscale simulations that permitted hundreds of atoms in the QM region were based on the tightbinding method (Broughton, et al., 1999), which is oversimplified and not able to describe bond breaking properly. A new attempt has been made to re parameterize the semiempirical Hamiltonian based on the input from highlevel quantum mechanical method, coupledcluster (CC) theory (Talor, et al., 2003). The new QM method (socalled "transfer Hamiltonian") is able to treat 5001000 atoms self consistently and retains the accuracy of CC method. In terms of future studies of the chemomechanical phenomena, the issues just summoned must be investigated carefully. To develop an accurate and flexible multi scale scheme in materials simulations, it remains for us to study a number of issues including more efficient quantum simulation, better embedding potentials, and higher quality QM/CM interfaces that allow sufficient and continuous chemical and mechanical information transfer across the interface. Conclusions In this work, we demonstrated a successful model to couple a QM region and a CM region for silicon dioxide systems. This combined QM/CM method has been applied to the study of hydrolysis of a twomembered ring on the amorphous silica surface. With this model, we have described the longrange response of an extended amorphous surface to water attack that neither the cluster nor the crystal/cristobalite models can treat. The simulations have revealed structure and energetic for both physisorption and chemisorption processes. The findings of cooperative proton transfer and a series of complex intermediate steps in the twowater case have significantly extended and advanced our understanding of the nature of the reaction. Our work focuses on understanding of the surface reactivity, the underlying mechanism of the hydrolytic weakening of silica, stress corrosion and crack growth in silica. It provides insight into the longterm performance of silica based material and is fundamentally important in silica fabrication and durability, polymer adhesion, cell biology, and industrial processes that utilize catalytic reactions. CHAPTER 4 ENERGETIC AND ELECTRONIC STRUCTURE OF THE CARBON NANOTUBE COMPLEXES Since the discovery of the carbon nanotube in 1991 (lijima, 1991), carbon nanotube research has become an interdisciplinary field that attracts many researchers worldwide. Carbon nanotubes have unusual structural and electronic properties (Saito et al., 1998), suggesting a wide variety of technological applications (Dresselhaus et al, 1996), including field effect transistors (FETs) (Tans et al., 1998; Martel et al., 1998; Derycke et al., 2002; Heinze et al., 2002; Appenzeller et al., 2002), hydrogen storage (Dillon et al., 1997; Liu etal., 1999), etc. The diameter of a carbon nanotube is of nanometer size while its length can be a few microns. The electronic structure of a singlewall carbon nanotube is either metallic or semiconducting, depending on its diameter and chirality (defined below). The energy gap of semiconducting nanotubes can be varied continuously from 1 eV to 0 eV, by varying the nanotube diameter. Another unique property of carbon nanotubes is their stiffness, corresponding to the upper limit of the best carbon fibers, which are commonly used as a strong lightweight material. Because of its very small size and special properties, the carbon nanotube has made great impact on semiconductor physics, nanoscience, and technology. In this chapter, I will first define the structure of the carbon nanotubes, and then discuss two important problems in the carbon nanotube research. One is the electronic structure of carbon nanopeapods (a carbon nanotube encapsulating a fullerene or metallofullerene chain) (Du and Cheng, 2003b). The other is the adsorption of bromine on the carbon nanotube, which is related to an important technical problem, i. e., the separation of the metallic and semiconducting carbon nanotubes (Chen, et al., 2003). Figure 41. A graphene sheet. The chiral vector Ch and the translational vector T are also shown. A (4, 2) nanotube can be constructed by rolling the graphene sheet into a cylinder such that the starting and ending points of the vector Ch are connected. The figure corresponds to Ch = (4,2), T = (4, 5), and dR =2. Structure of a SingleWall Carbon Nanotube A singlewall carbon nanotube (SWCNT) can be described as a graphene sheet rolled into a cylindrical shape so that the structure is onedimensional with axial symmetry, and in general exhibiting a spiral conformation, called chirality. The chirality can be defined by a chiral vector, Ch, given by Ch = na + ma (n,m) (41) where a, and a2 are the real space unit vectors of a graphene sheet, and n and m are integers (Fig. 41). A nanotube can be constructed by rolling a graphene sheet into a cylinder. Fig. 41 shows the construction of a (4, 2) nanotube. The circumferential length, L, of the carbon nanotube is just the length of the chiral vector: L = Ch = an + 2 + nm (42) where a is the length of the lattice vector a, or a2, and a = 1.42 A x/3 = 2.46 A. The translational vector T can be expressed as T= t1a, +t2a2 ( t2) (43) where t, and t2 are both integers. Since Ch*T = 0, the expressions for t1 and t2 are given by n+2m 2n+m tl= t2 = (44) dR dR where dR is the greatest common divisor of (n + 2m) and (2n + m) . Electronic Structure of the Carbon Nanopeapods Introduction Since the discoveries of carbon fullerenes in 1985 (Kroto et al., 1985) and nanotubes in 1991 (Iijima, 1991), these two carbon allotropes have attracted considerable attention in the scientific community because of their unique structures and novel properties. A carbon nanotube can be metallic or semiconducting, depending on its chiral vector (n, m). Chemical doping can further modify the electrical properties of carbon nanotube. For instance, potassium or bromine doping enhances conductivity (Lee et al., 1997). Alkali metal and alkalineearth intercalated fullerenes have also been studied extensively because of their complex phase diagram and superconductivity at temperatures surpassed only by the highTo cuprates (Yildirim et al., 2000; Gunnarsson, 1997). The empty spaces inside carbon fullerenes and nanotubes provide the possibility of modification of their properties by inserting atoms or molecules in these hollow spaces. A large number of species have been observed experimentally to sit stably inside fullerene cages (Bethune et al., 1993; Stevenson et al., 1999; Shinohara, 2000) or nanotubes (Fan et al., 2000; Mukhopadhyay et al., 2002; Shimoda et al., 2002). Theoretical calculations have revealed the details of energetic and electronic structure of fullerene and nanotube endohedral complexes (Dunlap et al., 1992; Wang et al., 1993; Cioslowski, 1995; Li and Tomlnek, 1994; Tominek and Li, 1995; Kobayashi and Nagase, 1996; Kobayashi and Nagase, 1998; Miyamoto et al., 1995; Miyake and Saito, 2002). Recently, Smith and co workers have successfully encapsulated C60 molecules in singlewall carbon nanotubes (Smith etal., 1998; Burteaux et al., 1999; Smith etal., 1999; Smith et al., 2000) and achieved highyield synthesis of such "peapods" (2000). More complex metallofullerene peapods ((Gd@C82)n@SWNTs, (La2@Cso)n@SWNTs, etc)1 have also been identified in transmission electron microscopy (TEM) images (Hirahara et al., 2000; Smith et al., 2000; Suenaga et al., 2000). 1 The symbol @ is used to indicate that species listed to the left of the @ symbol are encapsulated inside the species listed to the right of the symbol @. Encapsulation of fullerenes and metallofullerenes inside nanotubes enables the development of a new class of hybrid materials, which exhibits many interesting properties. The electrical resistivity, thermopower and thermal conductivity of peapods are different from those of empty SWNTs (Hirahara et al., 2000; Vavro et al., 2002; Pichler et al., 2001). Scanning tunneling microscopy (STM) studies show the spatially varied local density of states (LDOS) in C60@SWNT (Hornbaker et al., 2002) and the spatial modulation of the energy bandgap in (Gd@C82)n@SWNT (Lee et al., 2002). A temperatureinduced change from p to n conduction in (Dy@C82)n@SWNT (Chiu et al., 2001) and ambipolar fieldeffect transistor behavior of (Gd@C82)n@SWNT have also been reported (Shimada et al., 2002). An encapsulated C60 chain may also show superconductivity upon alkaline doping, in analogy with fullerene intercalation compounds (Yildirim et al., 2000; Gunnarsson, 1997). In order to fully understand how the different encapsulants affect physical properties of their host nanotubes, substantial theoretical studies of the electronic structures of peapods are highly desirable. Next, I report our work on electronic structures of several semiconducting and metallic nanopeapods. We manipulate the C60 induced impurity states inside the bandgap of the host semiconducting nanotubes by controlling the distance between the tube and C60 and the type of the encaged metal atom inside C60. Hybridization of states and the effects of density and orientation of encapsulated C60 molecules have been analyzed. In addition, I discuss the possibility of superconductivity of potassiumdoped peapods. Method and Simulation Setup The electronic structure calculations were performed using density functional theory (DFT) with the local density approximation (LDA). The electronion interactions were described by ultrasoft pesudopotentials (Vanderbilt, 1990). The valence wave functions were expanded in a planewave basis set with an energy cutoff of 286 eV, which was tested and found to give good energy convergence. The super cell was chosen such that the distance between adjacent nanotube walls is longer than 6.7 A. We use two k points in the irreducible Brillouin zone (BZ). All calculations were performed using the Vienna abinitio simulation package (VASP) (Kresse and Jurthmueller, 1996a and 1996b). Structures and Energetics As the model systems for our theoretical investigation, we choose (16, 0), (17, 0), (19, 0), (10, 10) nanotubes and their corresponding fullerene or metallofullerene encapsulated peapods. A commensurability condition is imposed between the periodicity of the nanotube and that of the C60 chain. The optimized lattice parameter c is 9.79 A for the (10, 10) nanotube, and 12.68 A for all zigzag tubes. In order to study (17, 0) peapods with different C60 intermolecular distance, we use two different lattice lengths (triple and quintuple periodicity of a zigzag nanotube) for the (17, 0) peapod. These two different unit cells contain one and two C60 molecules with optimized lattice parameters of 12.68 A and 21.14 A, respectively. The (17, 0) peapods with smaller and larger lattice parameters are labeled as (17, 0)a and (17, 0)b peapods in the rest of the paper. Fig. 42 shows the structures of the (10, 10), (17, 0)a and (17, 0)b peapods. The cohesive energy per C60 molecule is defined as E (peapod) E (nanotube) E(Nx C60 ) AE = (45) N The encapsulations of C60 molecules for (10, 10), (16, 0), (17, 0)a, (17, 0)b, and (19, 0) peapods are all exothermic, as shown in Table 41. The C60 intermolecular distances for (16, 0), (17, 0)a, and (19, 0) are same. (a) C60@ (10, 10) (c) C60@ (17, 0)b Figure 42. The structures of (a) (10, 10), (b) (17, 0)a, and (c) (17, 0)b nanopeapods. Table 41. The cohesive energy per C60 molecule for (10, 10), (16, 0), (17, 0)a, (17, 0)b, and (19, 0) peapod. (10, 10) (16, 0) (17, 0)a (17, 0)b (19, 0) AE(eV) 1.73 0.78 1.71 1.71 0.84 Manipulation of FullereneInduced Impurity States in M@C6o@SWNTs In an isolated C60 molecule, there is a threefold degenerate lowest unoccupied tlu state. Fig. 43 shows the electronic band structures of the C60@(17, 0) peapod. The band structure of the (10, 10) peapod obtained in this work is similar to a previous result (Okada et al., 2001). From Fig 43(b) and (c), the (17, 0) peapod remains a (b) C60@ (17, 0)a semiconductor. In the (17, 0)' peapod, the flat energy band inside the band gap is derived from the tlu state of C60. There is nearly no energy dispersion for the tiuderived band because of the long centertocenter distance dc between two C60 molecules, while there is a very small energy dispersion in the (17, 0)b peapod. When dc is large, the C60 chain is broken and not conductive so that the tiuderived states behave as impurity states [see Fig. 43(b)] and carriers are distributed on the nanotube wall. When d, decreases, these impurity states become impurity bands [see Fig. 43(c)] and conductive, resulting in most holecarriers distributed on the nanotube wall and most electroncarriers on the C60 chain at a temperature not high enough to make the intrinsic carrier concentration dominant. (a) (17, O)a (b) C60@(17, O)a (c) C60@(17, 0 S .0 *.* > ..* .* 0.5 E g 0.0 ** 0.5 *, * 1.0 F X F X 0.44 0.42 **, *. 0.40 ** ,* 0.38 " e.*. 0.36 r X Figure 43. Energy band structures of empty and C60 encapsulated (17, 0) nanotube. (a) Empty nanotube, (b) C60@(17, 0)a, (c) C60@(17, 0)b, only C60 tiuderived impurity bands inside the band gap are shown. The origins of the energies are set at the top of the valence band (Ev). Note the size of the firstBrillouin zone for the (17, 0)b peapod is 3/5 of that for the (17, 0)a peapod because of its larger unit cell. The small tiu bandwidth in the (17, 0)b peapod is not only due to the relatively larger dc, but also due to the ordered orientation of the C60 molecules inside the nanotube. On rotating one of the two C60 molecules in the unit cell of the (17, 0)b peapod around the )b tube axis by 1800, we find a doubled tiu bandwidth and the reduction of the tlu DOS peak by about 20%, while the total energy essentially remains the same. If d is about 10 A, as found in experiment, and the random orientation of C60 molecules at room temperature is taken into account, we expect a bigger tiu bandwidth than that shown in Fig. 43(c). (a) C60@(16, 0) (b) C6o)a (17, O)a (c) C~o( i(19, 0) 2.05 2.0 1 2.0., ..0 S.0 .0 O i 1.0 ,..,:1' 1.0 .1.0 ...* 0 0 ***0 *******05 o0.5 0.5 0.5 0.0 0. 0.0 .. 0.5 *::. 0.5 0.5 1.0 ill1 .....0 1.0 ' F X F X r X Figure 44. Energy band structures of(a) C60@(16, 0)a, (b) C60@(17, O)a, (c) C60@(19, 0)a. The C60 intermolecular distances in these three peapods are same. The origins of the energies are set at the top of the valence band (E,). While dc and the orientation of C60 molecules largely determine the tlu band dispersion, the space between the C60 and the nanotube wall is a key factor for the energy levels of tiuderived states. Our calculations show that the energy levels of tiuderived states are the inside conduction band in C60@(16, 0), about 2/3 of the bandgap (0.6 eV) above the top of the valence band (E,) in C60@(17, 0), and about 1/3 of the bandgap (0.48 eV) above E, in C60@(19, 0) (see Fig. 44). The diameters of the (16, 0), (17, 0) and (19, 0) nanotubes are 12.48, 13.26 and 14.80 A, respectively. It has been reported in a recent paper that energy levels of tiuderived states are 2/5 of the bandgap (0.5 eV) above E, in C60@(14, 7) (Miyake and Saito, 2003). The diameter of the (14, 7) nanotube is 14.48 A. It is evident that the tiuderived states shift down toward the top of the valence 