<%BANNER%>

Theoretical Modeling and Design of Complex Materials


PAGE 1

THEORETICAL MODELING AND DESIGN OF COMPLEX MATERIALS By MAO-HUA 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

PAGE 2

Copyright 2003 by Mao-Hua Du

PAGE 3

This dissertation is dedicated to my late grandfather, Mei-Lin Du.

PAGE 4

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 Hai-Ping Cheng, for her guidance and encouragement throughout my graduate career. Her open-minded 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. Lin-Lin 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 parents-in-law, for their help and wonderful food I iv

PAGE 5

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. v

PAGE 6

TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES...........................................................................................................viii LIST OF FIGURES...........................................................................................................ix ABSTRACT......................................................................................................................xii CHAPTER 1 INTRODUCTION........................................................................................................1 Density Functional Theory...........................................................................................1 Complex Systems.........................................................................................................2 Chemo-Mechanical Phenomena in Amorphous Silica..........................................3 Nanopeapods: Materials by Design.......................................................................3 2 METHODS...................................................................................................................5 Total-Energy Pseudopotential Calculations.................................................................5 Born-Oppenheimer Approximation.......................................................................6 Density Functional Theory....................................................................................6 Local-Density Approximation...............................................................................9 Generalized Gradient Correction.........................................................................10 Pseudopotential Approximation..........................................................................10 The Hellmann-Feyman Theorem........................................................................13 Molecular Dynamics...........................................................................................14 Density Functional Theory for Periodic Systems.......................................................14 Blochs Theorem.................................................................................................14 k-point Sampling.................................................................................................15 Momentum-Space Formalism.............................................................................15 Dual-Space Formalism........................................................................................20 Density Functional Theory for Finite Systems: A Dual-Space Formulism................21 Combining the Density Functional Theory with Classical Potential Functions: A Multiscale Approach..............................................................................................27 Introduction.........................................................................................................27 Partitioning of the System...................................................................................31 Total Energy........................................................................................................32 The QM/CM interface.........................................................................................35 vi

PAGE 7

The link-atom method..................................................................................35 The pseudo-atom method.............................................................................36 3 MULTISCALE SIMULATIONS OF AMORPHOUS SILICA.................................37 Introduction.................................................................................................................37 Test of the Link-Atom Method...................................................................................42 Test of the Pseudo-Atom Method...............................................................................52 Hydrolysis of a Two-Membered Silica Ring on the Amorphous Silica Surface........54 Generation of the Amorphous Silica Surface......................................................54 Hybrid QM/CM Approach..................................................................................59 Simulation Set-up................................................................................................61 Results.........................................................................................................................63 Adsorption of One Water Molecule on the Two-Membered Ring......................63 Adsorption of Two Water Molecules on the Two-Membered Ring...................72 Adsorption of Water Molecules on the Silanols..................................................73 Discussion...................................................................................................................76 Conclusions.................................................................................................................78 4 ENERGETICS AND ELECTRONIC STRUCTURE OF THE CARBON NANOTUBE COMPLEXES......................................................................................79 Structure of a Single-Wall Carbon Nanotube.............................................................80 Electronic Structure of the Carbon Nanopeapods......................................................81 Introduction.........................................................................................................81 Method and Simulation Setup.............................................................................83 Structures and Energetics....................................................................................84 Manipulation of Fullerene-Induced Impurity States in M@C 60 @SWNTs.........85 Electronic Structure of (K x C 60 ) n @SWNTs..........................................................92 Bromine Doping of Single-Wall Carbon Nanotubes: Separating Metallic from Semiconducting Nanotubes...................................................................................95 Conclusions...............................................................................................................100 5 CONCLUSIONS......................................................................................................102 APPENDIX A TOTAL ENERGY OF A PERIODIC SYSTEM IN MOMEMTUM-SPACE REPRESENTATION...............................................................................................104 B POPULATION ANALYSIS BASED ON ELECTROSTATIC POTENTIALS.....112 LIST OF REFERENCES.................................................................................................115 BIOGRAPHICAL SKETCH...........................................................................................123 vii

PAGE 8

LIST OF TABLES Table page 3-1. Several Si-O bond lengths in combined QM/CM model and cluster model (cluster B). ............................................................................................................................63 3-2. The adsorption energies of one and two H 2 O molecules in the cluster model with or without fixing the positions of H atoms that saturate the dangling bonds of the cluster.......................................................................................................................69 3-3. The adsorption energies of water physisorbed and chemisorbed on amorphous silica surface. ....................................................................................................................71 4-1. The cohesive energy per C 60 molecule for (10, 10), (16, 0), (17, 0) a (17, 0) b and (19, 0) peapod...........................................................................................................85 viii

PAGE 9

LIST OF FIGURES Figure page 2-1. The division of a system into QM and CM regions...................................................31 3-1. The structure of the Si 2 O 7 H 6 cluster. .........................................................................42 3-2. The structures of the cluster shown in Fig. 3-1, obtained from (a) all-quantum calculation, (b) combined QM/CM calculation using the link-atom method, and (c) all-classical calculation............................................................................................43 3-3. The structures of the training system with iso-density surfaces, obtained from (a) all-quantum calculation, and (b) combined QM/CM calculation..................................45 3-4. 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 3-5. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. .........48 3-6. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. .........49 3-7. The structure of the Si 3 O 9 H 6 cluster. .........................................................................50 3-8. The structures of the cluster shown in Fig. 3-7, obtained from (a) all-quantum calculation, (b) combined QM/CM calculation using the link-atom method, and (c) all-classical calculation (bottom panel). ..................................................................51 3-9. Cluster structures obtained from combined QM/CM calculations using the pseudo-atom method. ...........................................................................................................53 3-10. Partial pair-distribution functions for amorphous silica at 300 K............................56 3-11. One layer of the amorphous silica surface. ..............................................................57 3-12. Number density of atoms as a function of the distance from the center of the slab. 58 3-13. (a) The average Si-O bond length, and (b) the average Si-O-Si bond angle, as functions of the distance from the center of the slab. ..............................................58 3-14. A schematic picture of a silica surface (side view). ................................................62 ix

PAGE 10

3-15. Snapshots and potential energy of QM atoms along the dissociation path of one H 2 O molecule on a silica surface. ...........................................................................64 3-16. Cluster models..........................................................................................................65 3-17: Valence electron number density contours on the plane crossing Si1, O1 and O2 atoms in the combined QM/CM surface model (a) and the cluster model (b).........66 3-18. Contour plot of (a) excess electron density and (b) deficiency electron density in the water-surface system relative to superposed electron densities of the surface and H 2 O on the plane crossing Si(1), O(1) and O(2) atoms in the combined QM/CM surface model...........................................................................................................67 3-19. Contour plot of (a) excess electron density and (b) deficiency electron density in the water-cluster system relative to superposed electron densities of the cluster and H 2 O on the plane crossing Si(1), O(1) and O(2) atoms in the cluster model...........68 3-20. Snapshots and potential energy of the QM atoms along the water dissociation path (hydrogen relay process) on a silica surface. .......................................................74 3-21. Two H 2 O molecules on top of two silanols. ...........................................................75 4-1. A graphene sheet. The chiral vector and the translational vector are also shown.....80 4-2. The structures of (a) (10, 10), (b) (17, 0) a and (c) (17, 0) b nanopeapods................85 4-3. Energy band structures of empty and C 60 encapsulated (17, 0) nanotube. ...............86 4-4. Energy band structures of (a) C 60 @(16, 0) a (b) C 60 @(17, 0) a (c) C 60 @(19, 0) a ....87 4-5. (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 t 1u -derived energy band are shown for (c) a (17, 0) a semiconducting peapod and (d) a (10, 10) metallic peapod. .............................................................89 4-6. Energy band structures of (a) K@C 60 @(17, 0) a (b) Ca@C 60 @(17, 0) a (c) Y@C 60 @(17, 0) a (spin-up), (d) Y@C 60 @(17, 0) a (spin-down), and (e) K@C 60 @(16, 0). ......................................................................................................90 4-7. The geometry of K 3 C 60 @(17, 0) b ............................................................................93 4-8. The density of states of (a) C 60 @(10, 10), (b) K 1 C 60 @(10, 10), (c) K 2 C 60 @(10, 10), (d) K 3 C 60 @(10, 10), (e) K 4 C 60 @(10, 10), and (f) K 5 C 60 @(10, 10). .......................94 4-9. The density of states of (a) (C 60 ) 2 @(17, 0) b (b) (K 1 C 60 ) 2 @(17, 0) b (c) (K 2 C 60 ) 2 @(17, 0) b (d) (K 3 C 60 ) 2 @(17, 0) b .............................................................95 4-10. 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 x

PAGE 11

4-11. 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 4-12. 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 xi

PAGE 12

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 Mao-Hua Du December 2003 Chair: Hai-Ping 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 two-membered 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 xii

PAGE 13

thus control the type of the majority carrier (p or n) and the carrier density of a semiconducting peapod. The possibility of superconductivity of potassium-doped peapods is also discussed. xiii

PAGE 14

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 mid-1960s (Hohenberg and Kohn, 1964; Kohn and Sham, 1965). In that work, they proved that the ground-state many-electron problem can be replaced by an exactly equivalent set of self-consistent one-electron equations (Kohn-Sham equations); and that the many-body part of the problem can be included by an effective exchange-correlation potential. Therefore, the complex many-body problem is essentially reduced to the 1

PAGE 15

2 problem of a non-interacting inhomogeneous electron gas in an effective potential. If one knows all of the appropriate electronic interactions, for example, the electron-electron electrostatic interaction, the exchange-correlation, the ion-ion interaction, and the electron-ion 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 cost-effective first-principles method and has reproduced a variety of ground-state 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

PAGE 16

3 and fullerene-based hybrid materials, in which simple structures with different properties are self-assembled to form complicated structures with enhanced properties. Chemo-Mechanical Phenomena in Amorphous Silica The framework of solid state physics was built based on Blochs 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 so-called chemo-mechanical 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 chemo-mechanical phenomena are still in a primitive state because the chemical processes at nanoscale and the mechanical properties at mesoor macroscopic-scale 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 link-atom method (Du and Cheng, 2003a). This multi-scale 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, zero-dimensional dots or nanocrystals, one-dimensional wires, and two-dimensional films often have unusual properties distinctly different from those of the

PAGE 17

4 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 Feynmans 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 self-assembly 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 fullerene-induced states inside the band gap of the host-semiconducting 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 potassium-doped peapods.

PAGE 18

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 total-energy calculation, and this literature is still growing. A thorough review of this field is certainly beyond the scope of this work and the authors knowledge. Instead, a short introduction to some of the important ingredients in the total-energy calculations is given. In this chapter, I first introduce the Born-Oppenheimer approximation, density functional theory, the pseudopotential method, the Hellmann-Feynman theorem, and molecular dynamics. An excellent review of the total-energy 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; Barnett and Landman, 1993). Only the plane-wave based implementation is discussed. Finally, I describe the hybrid method that combines DFT with classical potential functions (Du and Cheng, 2003a). Total-Energy Pseudopotential Calculations To perform accurate and efficient total-energy calculations, a series of simplifications and approximations is needed. These include the Born-Oppenheimer 5

PAGE 19

6 approximation to separate the electronic and nuclear wave-functions, density functional theory to model the electron-electron interaction, and the pseodopotential method in conjunction with a plane-wave basis set to model the electron-ion interaction. Born-Oppenheimer 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 many-body wave-functions. This so-called Born-Oppenheimer 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 many-body 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 total-energy functional of the electron density yields the ground-state energy of the system, and the corresponding electron density is the exact ground-state electron density. Kohn and Sham (1965) then showed that the many-electron problem can be replaced by an exactly equivalent set of self-consistent one-electron equations (Kohn-Sham equations). To simplify the notation, I will present the total energy functional and Kohn-Sham equations for the spin-unpolarized case. They can be easily extended to the spin-polarized 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

PAGE 20

7 (2-1) [{},()]({})[()]totalIionIelecEEERrRr Here, is the Coulomb energy associated with the interactions among the ions at positions{, given by ({})ionIER}IR IJionIJIJZZERR (2-2) where and are the charges of the Iand Jth ions, and Hartree atomic units are used. is the ground-state energy of the electrons. The ground-state electronic energy is given by IZelecE IZ)] [(r (2-3) eleceeIeeETEE where T is the kinetic energy, is the electron-ion interaction energy, and is the electron-electron interaction energy. These terms are given by e eIE eeE 21()2eiiiiTfdr 3r (2-4) (2-5) 3()()eIionEVdrr 331()(')'[(2'eeXCEdrdrErrrrr )] (2-6) 2()()iiifrr (2-7) where V is the static total electron-ion potential, is the electron density, is the electron occupation number, and is the exchange-correlation energy functional. ion ()r if [()]XCEr The set of one-electron states that minimize the energy functional is given by the solutions to the Kohn-Sham equations (Kohn and Sham, 1965): i elecE

PAGE 21

8 21()()()()()2ionHXCiiiVVVrrrr r (2-8) where V is the Hartree potential of the electrons given by ()Hr 3(')'HVrrr (2-9) dr The exchange-correlation potential, V, is given by the functional derivative ()XCr ()()()XCXCEVrrr (2-10) The Kohn-Sham equations represent a mapping of an interacting inhomogeneous electron gas in the presence of an external potential V onto a non-interacting inhomogeneous electron gas moving in an effective potential If the exchange-correlation functional were known exactly, the Kohn-Sham equations would give the exact ground-stateand. ()r ()()()()effHXCVVVVrrr r elecE The Kohn-Sham equations are nonlinear and must be solved self-consistently. Once the eigenvalues and the wave functions are obtained, the total electronic energy can be determined from Eq. (2-3) to (2-7), or alternatively from the formula i i elecE 3331()(')'()()()2'NeleciiXCXCiEfdrdrEVdrrrrrrrr (2-11) where N is the number of occupied Kohn-Sham orbitals. Here 21()2NNiiiieffiiiffVr (2-12) The sum of the single-particle Kohn-Sham eigenvalues does not give the total energy, because it overcounts the electron-electron interaction and the exchange-correlation energy.

PAGE 22

9 Local-Density Approximation The Kohn-Sham equations [Eq. (2-8)] are exact provided that an exact exchange-correlation functional is known. However, the exact is usually not known, hence an approximation must be used. The simplest approximation is the local-density approximation (LDA) proposed by Kohn and Sham (Kohn and Sham, 1965). XCE XCE In LDA, the exchange-correlation functional can be written as (2-13) 3()()XCXCErr dr Thus, ()()()()()XCXCXCEVrrrr rr (2-14) The local density approximation assumes that the exchange-correlation energy functional is purely local. The exchange-correlation energy per electron is equal to the exchange-correlation energy per electron in a homogeneous electron-gas that has the same density as the electron gas at r. Thus, XC (2-15) hom()()()LDALDAXCXCXCrr A commonly used parameterization for is the Vosko-Wilk-Nusair parameterization (Vosko et al., 1980, 1982) of the Ceperly-Alder 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 exchange-correlation energy. LDAXC Considering the inexact nature of the local density approximation, it is remarkable that calculations using LDA have been so successful. Generally, total-energy differences between related structures can be believed to within a few percent and structural

PAGE 23

10 parameters to at least within a tenth of an 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 exchange-correlation 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 for. In it 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. XCE XCE Pseudopotential Approximation The total-energy density functional method is often implemented in a plane-wave basis set, which is discussed later. Here, I introduce the pseudopotential approximation that greatly reduces the computational cost of the plane-wave based DFT method. A large number of plane-waves are needed to expand tightly bound core orbitals and rapidly oscillating valence wave-functions 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

PAGE 24

11 pseudopotential that acts on pseudo-wave-functions rather than the true valence wavefunctions. For a central field atom, the radial Kohn-Sham equation is given by 2221(1)()()()22effnlnlnldllVrrRrrRdrr r (2-16) The pseudo-wave-functions are obtained from an all-electron atomic calculation that self-consistently solves Eq. (2-16). Once the pseudo-wave-function is obtained, the screened pseudopotential is given by the inversion of Eq. (2-16), 2,22(1)1()22()PPscrPPlllPPllldVrrRrdr rRr (2-17) where is the radial pseudo-wave-function. The ionic pseudopotential can be obtained by subtracting the Hartree V and exchange-correlation V potentials calculated from the valence pseudo-wave-functions from the screened pseudopotential ()PPlRr()Vr ()PPHr ()PPXCr ,PPscrl (2-18) ,()()()()PPPPscrPPPPllHXCVrVrVrVr The pseodopotential needs to meet a set of requirements. The standard requirement is norm-conservation, which means that the pseudo-wave-functions are identical to the real wave-functions outside a core region, and that they generate same electron density inside the core region. A widely used norm-conserving pseudopotential is the Troullier-Martins pseudopotential (Troullier and Martins, 1991a). 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

PAGE 25

12 wave-function. The pseudopotential operator can be decomposed into local and nonlocal parts, (2-19) ,,()()()PPPPlcPPnlcllVrVrVrP lrr,r where V is the local potential, and the semi-local potential V is given by ,()PPlc ,()PPnlcl (2-20) ,()()()PPnlcPPPPloclllVrVrV Here, is the operator that projects out the lth angular momentum component from the wave-function. The nonlocal potential is short-ranged, since at large r, the ionic potential is independent of angular momentum and approaches the value of lP Zr where Z is the core charge. In principle, the local potential can be arbitrarily chosen. However, since the summation in Eq. (2-19) 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. (2-20), ,,,,()()()()()()()()PPnlcPPPPPPnlcllllKBnlclPPPPnlcPPlllVrrrVrVrrVrr (2-21) where V is given by Eq. (2-20), and is the atomic pseudo-wave-function for the angular momentum l. This form of the non-local potential substantially reduces the number of integrals that need to be evaluated. ,()PPnlcl r ()PPlr The norm-conserving condition requires that the total pseudocharge inside the core region match that of the all-electron wavefunction. Thus for those highly localized valence orbitals (e.g., 2p for first-row atoms and 3d for transition-metal atoms), a large

PAGE 26

13 number of plane-waves are still needed. Vanderbilt relaxed the norm-conserving rule and constructed ultrasoft pseudopotentials with a lower cutoff energy for plane-wave 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 first-row elements by a factor of 2 to 4 (Kresse and Furthmuller, 1996a). The Hellmann-Feyman Theorem Once the total energy of the system is obtained via Eq. (2-1), the force on ion I,f, is given by I iIiiiidEEEEdIIIfRRR iIR (2-22) If are eigenstates of the Hamiltonian, and given that i iE is just the last two terms in (2-22) can be rewritten as iH 0iiiiiiiiiiiiiiiiiHHIIIIiRRRRR i (2-23) Thus, IEIfR (2-24) This is the Hellmann-Feynman theorem (Hellmann, 1937; Feynman, 1939), which holds for any derivative of the total energy.

PAGE 27

14 Molecular Dynamics The motion of ions on the Born-Oppenheimer potential energy surface is governed by the forces obtained by Eq. (2-18) via integration of the Newtonian equations of motion (2-25) ImIIRf The integration can be performed using the Verlet algorithm (Allen and Tildesley, 1987). In contrasted to empirical potentials, the forces determined by ab-initio 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 Blochs theorem in order to perform a realistic total-energy calculation. Blochs Theorem A wave-function in a periodic system can be written in the following form (Ashcroft and Mermin, 1976) (2-26) ,()exp()iriukkrr ,ikRGr Here, u is a periodic function that satisfies the condition ,()ikr (2-27) ,,()()iiuunkkrr where R is any lattice vector and n is an integer. Such a function can be expanded in a discrete plane-wave basis set with wave vectors as reciprocal lattice vectors of the crystal: (2-28) ,,()exp()iiucikGGr

PAGE 28

15 Therefore the electronic wave-function can be written as 2,,1,2()exp()cutiiEcikkGGkGr kGr (2-29) The summation is over all reciprocal lattice vectors up to a cutoff energy. But in principle. cutE cutE k-point 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 k-point. Therefore, the Bloch theorem transforms the problem with an infinite number of electronic wave-functions to one with a finite number of electronic wave-functions 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 Blochs theorem, k-point 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. Momentum-Space Formalism In this section, I show the derivation of the momentum-space 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:

PAGE 29

16 (2-30) totaleeIHXCionETEEEE where 321()()2eiiiiTfdrr rr (2-31) (2-32) 3,,()()()PPeIiiiisEfdrVsRrrRt 31()()2HEdrVrr Hdr (2-33) (2-34) 3XCXCErr '',,'2ssionssZZNERss'ttR (2-35) 2()()iiifr rlrP (2-36) The symbols and are, respectively, the electron occupation number and the valence pseudo-wave-function for the one-electron state i. N is the number of unit cells in the crystal, and R is any lattice vector. V is the Hartree potential, given by Eq. (2-9). The symbols and are, respectively, the core charge and the basis vector for atom s in the unit cell. The prime in Eq. (2-34) means that the s=s term is excluded in the summation if the sand sth ions are in the same unit cell. V is the angular-momentum dependent pseudopotential and is given by ifsZ ()irst ()Hr PP (2-37) ,,()()()()PPPPPPlcPPnlclllllVVrPVrVr The momemtum-space representation of each energy term per unit cell in Eq. (2-30) is given by

PAGE 30

17 22,,12eiiiTfckkkGkGkG (2-38) 220()42cHEGGG (2-39) (2-40) XCcXCEGG G ,0('),,,,',,,',,,'3,1siPPlceIcssiPPnlcciiislilssPPlcesscsEeVfcceVZQdrVsGtGGGtkkkGkGkGkGkGGGGrrt (2-41) ''',''2'222012241expcos4{}ssionssssssssssccerfcEZZRGttRttRGGttG (2-42) where 31cicedGrGr r (2-43) ('),,,','1iiiiicfcceGGrkkkGkGkGGr (2-44) ,3,1()PPlcPPlcisscVdrVGrG er (2-45) 31iXCXCcedGrGr r (2-46) ,,2,,,',04(21)(cos)'PPnlcPPnlcslllsllcVlPjrVrjGkGkGkG rrdr (2-47) and

PAGE 31

18 ()(cos'kGkGkGkG ') (2-48) Here, is the volume of the unit cell, is the normalized weight for each k point, is the total charge of electrons in the unit cell, and 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. c k eQ lj lP If the Kleinman-Bylander type of nonlocal pseudopotential is used, Eq. (2-47) becomes ,2,,,',0,204(21)(cos)()()()'()PPPPnlcPslsllllPsPPnlcPsclllPPnlcPsllllPVjrRrVrRrjrVrRrrdrGkGkGkG VrRrrdr (2-49) where is the radial pseudo-wave-function. If there are n different G vectors, for each l and each k point, one reduces the required number of integrals from n(n+1)/2 to n by using the Kleinman-Bylander nonlocal pseudopotential. ()PslRr The integrals in Eq. (2-49) 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 short-ranged and thus is only nonzero within a sphere around each ion, the number of operations necessary to evaluate the integral in real-space does not increase with the system size. Therefore, it is more efficient to evaluate the integrals in Eq. (2-49) in real space for large systems. In order to do so, the wave-function for each angular momentum must be represented in real space. However, the discrete real-space Fourier grid introduces errors when projecting angular momentum components of the wave-functions

PAGE 32

19 because of the non-uniqueness of the spherical harmonic projection in a discrete real-space grid (Payne et al., 1992). To reduce these errors, King-Smith, Payne and Lin (1991) proposed a procedure for improving the accuracy of the real-space calculation of the nonlocal contribution to the total energy. Tests of this method showed good agreement between realand reciprocal-space calculations (King-Smith, Payne and Lin, 1991). The Kohn-Sham equations in momentum-space representation can be derived variationally from the expression for the total energy in Eq. (2-30), 2',,,',',,'12iiiiVccGGkGGkGkkGGkG (2-50) where (2-51) ('),,,,,',,,'(')(')(')iPPlcPPnlciHXCsslslVVVeVVGGkGGkGkGGGGGGG Here, V is given by Eq. (A.10). HG In practice, the divergent terms V and V are set equal to zero when solving Eq. (2-50). (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. (2-41)] should be added to the total energy. Once Eq. (2-50) is solved, the total energy can be obtained, 0H ,0PPlcs ,,3,1totaliiHXCcXCionisPPlcesscsEfEEVZQdrVkkkkGGGrrt E (2-52)

PAGE 33

20 where and are given by Eq. (2-39), (2-40) and (2-42), respectively. This expression is simpler to evaluate than Eq. (2-30), (2-38)-(2-42), since the double summation over reciprocal lattice vectors in Eq. (2-41) is absorbed in the simple summation of the eigenvalues of the occupied states in Eq. (2-52). HE XCE ionE The Hellmann-Feynman force on each ion is given by (2-53) ,ssesFFF ,I where (2-54) ,,0('),,,,',,,',,,'()'seseIiPPlcciPPnlcciiislilsEieVifcceVsGrGGGtkkkGkGkGkGkGGFGGGGG and 22',2''22'0''32''4expsin42{}sssIsionsssssscssssssssEZZerfceGttRRFGGGttGttRttRttRttR (2-55) ,seF is the force contribution from the valence electrons, and is the force contribution from the ions. ,sIF Dual-Space Formalism The procedure of evaluating the kinetic energy in momentum space and the potential energy in real space can be thought of as a dual-space formalism for pseudopotential calculations (Martin and Cohen, 1988; Troullier and Martin, 1991b; Wentzcovitch and Martin, 1991) in contrast with momentum-space formalism (Ihm, Zunger and Cohen, 1979) in which all the calculations are carried on in momentum

PAGE 34

21 space. The dual-space 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 dual-space formalism for the pseudopotential calculations of finite systems. Density Functional Theory for Finite Systems: A Dual-Space Formulism For the study of finite systems, the method described in Section (2-2) 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 point charges decreases as 1L and the corresponding value for dipole interaction is 31L where 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 plane-wave 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 multipole moments efficiently, though they may become prohibitive for a supercell method. 3L Following the notation of Barnett and Landman (1993), the eigenfunctions of the Kohn-Sham equations are denoted where is the spin index. The corresponding density of spin is given by j

PAGE 35

22 2()()jjjfr re (2-56) where are the occupation numbers, with and is the number of valence electrons. jf ,jjfN eN In a finite system, the total energy is expressed as usual by Eq. (2-30). Each energy term in Eq. (2-30) is given by 2,12ejjTfjj (2-57) ,PPeIjsjsEfjVj (2-58) 31()()2HEdrVrr Hdr (2-59) (2-60) 3XCXCErr ''ssionssZZEss'tt (2-61) where V is the Hartree potential, given by Eq. (2-9), V is the pseudopotential, which can be decomposed into local and nonlocal parts, H PPs (2-62) ,,(,')()(,')PPPPlcPPnlcsssVVVrrrrr Thus, (2-63) lcnlceIeIeIEEE where 3,lcPPlceIsssEdrVrr t (2-64)

PAGE 36

23 ,,,,,,2,,,,,,023,,,,,()()()()()PPnlcPPnlcslssllmsllmslsnlceIjPPnlcjslmslsljslslmsjjslmjVRrYrRrYrVjEfdrrRrVrfFdrKrtrtrtr (2-65) ,2,,,01()slPPnlcslslFdrrRrVr (2-66) (2-67) ,,,,PPnlcslmslsllmKVrRrYr r In Eq. (2-65), the Kleinman-Bylander non-local pseudopotential is used (Kleinman and Bylander, 1982). The Kohn-Sham Hamiltonian is given by 212eIHXCHVV V (2-68) where XCXCXCEV (2-69) ,,,,,,',''lcnlceIeIeIPPlcssslslmsslmssslmVVVVFKKrrrrrrtrtrt (2-70) and V is the Hartree energy, given by Eq. (2-9). H Barnett and Landman (1993) suggested a dual-space 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 jHk is also evaluated in the dual-space,

PAGE 37

24 2323333233,3,1,'212''','12'''eIHXCeIHXCPPlcPPnlcHXCeIeIjHkjVVVkdkjkdkdrdrjVVVkdkjkdkjdrVVVkdrVkrrrrkkkkrrrrrrrrkkkkrrrrrrrrrr z (2-71) In practice, the integrals over real and momentum spaces are replaced by summations over the discrete real and momentum-space grids. The size of the real-space cell is defined as follows (2-72) ,,;0,0,0xyxyzxLxLxLr The wave-function ()jjrr is expanded in a plane-wave basis within the cell, i.e., max1()()ijje g rggr g (2-73) for r in the cell, and (2-74) ()0jr for r outside the cell. In Eq. (2-73), is the volume of the cell, and is the cutoff momentum in momentum space. The momentum-space grid (g-space) is defined by xyzLLL maxg 2,,yxzxyzkkkLLLg (2-75) where are integers satisfying ,,xyzk

PAGE 38

25 22Nk N (2-76) and maxyxzxyzNNNLLL gr (2-77) For a finite system, the wave-function can be chosen to be real in real-space, i.e., (2-78) ()()jjr and thus, ()()jj g g (2-79) The density can also be expanded in a Fourier series but requires a doubled momentum-cutoff with respect to the wave-functions, i.e., r (')'1'1irjjjifeDeggggGrGrgG g (2-80) The G-space grid in Eq. (2-80) is defined as 2,,yxzxyzmmmLLLG (2-81) where are integers satisfying ,,xyzm 222MMmandMN (2-82) Now, a real-space S-grid can be defined as a dual of the G-space grid. The density and the Kohn-Sham potential can be evaluated on this S-space grid, i.e.,

PAGE 39

26 ,,yyxxzzxyznLnLnLMMMS (2-83) where are integers satisfying ,,xyzn 01 (2-84) nM The wave-function can be evaluated on the S-grid as 1ijjjeMM g SgSSg (2-85) where The density on the S-grid is given by xyMMMM z 2jjjDfS S (2-86) Now the different energy terms in Eq. (2-59), (2-60), (2-64) and (2-65) can be directly evaluated on the S-grid as follows 12HESS HVSS (2-87) (2-88) XCXCESS ,lcPPlceIsssEVSSS t (2-89) 2,,,,,nlceIjslslmsjjslmEfFKSStS (2-90) The kinetic energy is evaluated in momentum space. The Hamiltonian matrix element in Eq. (2-71) is written as

PAGE 40

27 2,,,,,,,'121''{}jjPPlcjHXCslssislslmsslmsjslmSjHkgVVVMFKKeggSgSgg j g SSStStStS S (2-91) The force on each ion is given by '''sssseIsssZZEFrr (2-92) 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

PAGE 41

28 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 mesoand 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 first-order 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 bond-breaking and bond-forming. 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

PAGE 42

29 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 length-scales 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 all-quantum 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 long-range 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

PAGE 43

30 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 multi-scale method, which couples finite element, molecular dynamics, and semi-empirical tight-binding methods, to study crack propagation in silicon. However, the tight-binding method is over-simplified and not able to describe bond breaking properly. Sauer and Sierka (2000) have developed a combined quantum mechanics-interatomic 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 multi-scale simulation techniques, we have developed a combined quantum-mechanical and classical-mechanical (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-quantum-mechanical 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.

PAGE 44

31 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 2-1. 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. Figure 2-1. The division of a system into QM and CM regions.

PAGE 45

32 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., (2-93) /,QMnCMsQMCMnsEEEErrrrr Here, the QM nuclei, CM atoms are labeled by indices n and s, respectively. In Eq. (2-93), only is explicitly dependent on the electron density (r), which is obtained by solving the Kohn-Sham equations self-consistently QME 21()[()]()()2'()ixceIiEVdr'rrr'rrrr ir (2-94) 2()()iiifr r (2-95) where V is the ionic pseudopotential, E eIr xc is the exchange-correlation energy, and are the electron occupation number. The different energy components can be computed as following: if 22'''11()()221()()[()]2nnQMnniiinnnninniieIixcinZZEmdffVdddERrRRrr'rrrrrr'rrr' rr (2-96) 21({})2CMssssEmURR (2-97) (2-98) /({},{})QMCMnsEVRR where is the core charge of the n-th ion. U({R nZ s }) and V({R n },{R s }) are described by sums of classical interatomic potential functions in pairwise or three-body form.

PAGE 46

33 The particles inside the QM region are treated as ions and electrons when computing However, for computing, 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 ensures the continuity of the mechanical properties across the QM/CM boundary. If simple pair-wise classical interatomic potential functions are used, U({R QME /QMCME /QMCME ij s }) and V({R n },{R s }) are expressed as (2-99) '''({})()sssssUR sssrr3ssssrr')M (2-100) ({},{})nsnsnsnsVRR One can also use three-body potential functions having the form (2-101) 2(,)(,,)ijiijikjfijfijk where is an energy parameter, f 2 is the pairwise interaction between the i-th and j-th atoms, and f 3 is the three-body term as a function of the i, j, k-th atoms. In this case, U({R s }) and V({R n },{R s }) are expressed as (2-102) 2'3'''''''''''({})()(,,)sssssssssssssUfrfrR (2-103) 23''3'''({},{})()(,,)(,,)nsnsnsnsnsssnnnsnsnsssnnnsVfrfrrrfrrrRR The forces acting on the QM and CM atoms are respectively given by (2-104) /(nnQMQMCEEF

PAGE 47

34 and (2-105) /(ssCMQMCEEF )M In many other hybrid QM/CM methods, the QM-CM 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., 3/,selectrostaticnsQMCMsnssnqZqEdrrrRRR s (2-106) where is the partial charge of s-th classical atom, and is the n-th quantum ion. In these methods, includes the electrostatic interactions between the QM electrons and CM point charges and is explicitly dependent on the electron density 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, Newtons 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 QM-CM and CM-CM 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. sq nZ /QMCME r

PAGE 48

35 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 link-atom and the pseudo-atom methods. The link-atom 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 link-atoms generally are hydrogen atoms. They are attached to the broken covalent bonds on the QM/CM boundary to saturate the valence. Since the electronegativity of H 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 H-termination for clusters of silica and zeolites. In our method, the link-atoms 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, (2-107) /,,,QMnLCMsQMCMnsEEEErrrrrr where L stands for the link-atom. If the atoms and we define the constraint such that the link-atom is at a constant distance from along the bond. The bond nAQM sBCMnA LC nA sB nA LC

PAGE 49

36 length is optimized to minimize the force on C in some relevant small cluster calculations. L The pseudo-atom method The alternative to using hydrogen atoms to terminate the dangling bonds at the boundary is to use pseudo-atoms to absorb these bonds. Taking the silicon dioxide system as an example, a pseudo-oxygen-atom 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 pseudo-atom 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 pseudo-atom method, the dangling bonds are absorbed into the pseudo-atoms. However, the orientation of the dangling bond is not taken into account if simply using a real atom as the pseudo-atom. If one chooses to re-parameterize the pseudopotential, the bond-angle dependence is still difficult to include in the pseudopotential. Thus, it is one advantage of the link-atom method since such bond-angle dependence is naturally included if the link-atoms are constrained to lie on the line of the bonds that connect the QM and CM parts.

PAGE 50

CHAPTER 3 MULTISCALE SIMULATIONS OF AMORPHOUS SILICA Introduction Silica is the most abundant chemical species in the earths 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 semiconductor-ceramic research. Amorphous silica is abundant in the mineral and rock. The study of water-rock 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 water-rock 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 water-silica interaction (Zhuravlev, 1993; Zhuravlev, 2000; Bunker, et al., 1989a, Bunker, et al. 1989b, Chuang 37

PAGE 51

38 and Maciel, 1996; Chuang and Maciel, 1997; Bolis, et al., 1991; Chiang, Zegarski, and Dubois, 1993; DSouza 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 chemo-mechanical 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, Si-O bonds are known to be inert to chemicals (H 2 O, NH 3 etc) that can cause stress corrosion in silica glass. However, strained Si-O 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 (Si-O-H) with coverage of about 4.9 silanols/nm 2 (Zhuravlev, 1993; Zhuravlev, 2000). (The number of Si atoms on the surface is about 8 Si/nm 2 .) H 2 O 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 pre-treatment at temperature higher than 400C) 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; DSouza and Pantano, 1999). On the dehydroxylated silica surface,

PAGE 52

39 some highly strained defects consisting of a two-membered ring can be identified by IR peaks at 888 and 908 cm -1 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 two-membered ring structure creates a local reactive site that is vulnerable to chemisorption of gas-phase species. The reaction between water and strained Si-O bonds is an important aspect of surface hydroxylation. Reliable computational models for water-silica 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 co-workers (2000) studied the activation barriers of water dissociative reactions with a number of clusters that represent defect sites on the amorphous surface. For the SiO 2

PAGE 53

40 system, the charge transfer is local and the electronic response is short-ranged. 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 co-workers (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 -cristobalite surface to model a hydroxylated amorphous surface (Ceresoli, et al., 2000; Iarori, 2001). This approach was adopted because the density and refractive index of -cristobalite 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 co-workers (2001) have used the -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-, six-membered rings) but the -cristobalite 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 long-range

PAGE 54

41 order and hence long-range surface response in the -cristobalite 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 -cristobalite 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 multi-scale schemes usually involve local chemical processes in systems that do not have long-range order. Apparently, multi-scale 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 multi-scale 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 chemo-mechanical phenomena in amorphous silica, in particular, stress-induced 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.

PAGE 55

42 In the following part of this chapter, I will first show the results of test of the QM/CM interface (the link-atom and the pseudo-atom methods) on small silica clusters (Du and Cheng, 2003a), and then apply the combined QM/CM method on the studies of hydrolysis of the two-membered 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 Link-Atom Method A small silica cluster (Si 2 O 7 H 6 see Fig. 3-1) 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. Figure 3-1. The structure of the Si 2 O 7 H 6 cluster. The gray, dark gray, and light gray spheres represent the Si, O, 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. 3-2). The link-atom (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

PAGE 56

43 systems, the link-atoms will be constrained on the line of O()-Si() bonds with a fixed O()-H(link) bond length of 1.82 bohr. The Beest-Kramer-Santen (BKS) potential is chosen to describe the classical interatomic interactions (van Beest, Kramer, van Santen, 1990). For atoms i and j, the BKS potential has the form of QM CM QM 6exp()ijijijijijijijijqqcAbrrr (3-1) where r ij is the distance between the vi and j-th atoms. (a) (b) (c) Figure 3-2. The structures of the cluster shown in Fig. 3-1, obtained from (a) all-quantum calculation, (b) combined QM/CM calculation using the link-atom method, and (c) all-classical calculation.

PAGE 57

44 The optimized structures obtained from the all-quantum calculation, combined QM/CM calculation, and all-classical calculation, are shown in Fig. 3-2. 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. 3-1, indicate that the QM part of the hybrid system has a structure very similar to that in the all-quantum calculation, and the corresponding CM part has a structure nearly identical to that in the all-classical 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. 3-3. If including the electrostatic interaction between electrons in the QM region and atoms in the CM region [see Eq. (2-106)], the partial charges [ and in Eq. (3-1)] assigned to CM atoms in Fig. 3-1 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 link-atoms can terminate the dangling bonds very well and approximately restore the correct charge distribution at the QM/CM boundary (see Fig. 3-3), the additional charge polarization induced by CM point charges is not necessary. Thus, we did not use Eq. (2-106) to include the electrostatic interaction between QM electrons and CM atoms. iq jq

PAGE 58

45 Figure 3-3. The structures of the training system with iso-density surfaces, obtained from (a) all-quantum calculation, and (b) combined QM/CM calculation. Fig. 3-4 shows the total energy change due to the change of the Si(1)-O(b) or Si(2)-O(b) bond length. When the Si-O 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 Si-O bond in the CM region is changed this energy curve is close to that from the all-classical calculations. This is an indication of a transparent quantum-classical interface.

PAGE 59

46 A more rigorous test of the interface is shown in Figures 3-5 and 3-6. 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 all-quantum 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 all-classical calculations. Based on these tests (Figures 3-4, 3-5 and 3-6), 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. 2-99 and 2-100). On the other hand, to the quantum atoms, which consists of ions and electrons, the classical atoms plus the link-atoms 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 3-4~3-6 arises from the disagreement between the quantum and classical treatments of the interactions. The Si(1)-O(b)-Si(2) bond angle in Si 2 O 7 H 6 cluster (Fig. 3-2) is 180. However, the average Si-O-Si bond angle in amorphous silica is around 146. Thus, it is necessary to further test the validity of the link-atom method in other silica clusters with different chemical environments.

PAGE 60

47 (a) (b) Figure 3-4. 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. The long-dash red line, solid black line, and short-dash 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 fully-relaxed cluster is shifted to zero.

PAGE 61

48 (a) (b) Figure 3-5. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. The long-dash red line, solid black line, and short-dash red line represent the results from the QM, combined QM/CM, and CM calculations, respectively.

PAGE 62

49 (a) (b) Figure 3-6. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. The long-dash red line, solid black line, and short-dash red line represent the results from the QM, combined QM/CM, and CM calculations, respectively.

PAGE 63

50 Figure 3-7. The structure of the Si 3 O 9 H 6 cluster. The gray, dark gray, and light gray spheres represent the Si, O, and H atoms, respectively. To test the transferability of the QM/CM interface, the link-atom method is applied to a three-membered silica ring (Si 3 O 9 H 6 see Fig. 3-7). 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. 3-8), i.e., the structure parameters, generated by the combined QM/CM method, are in good agreement with the results from the all-quantum calculations in the QM region, as well as with the results from the all-classical calculations in the CM region. The tests in small silica clusters demonstrate that the link-atom 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 Si-O bonding and the proper construction of the Hamiltonian based on that property.

PAGE 64

51 (a) (b) (c) Figure 3-8. The structures of the cluster shown in Fig. 3-7, obtained from (a) all-quantum calculation, (b) combined QM/CM calculation using the link-atom method, and (c) all-classical calculation (bottom panel). In (a) and (c), the lengths of the Si-O bonds on the plane of the ring are identical and this ring has a three-fold rotational symmetry.

PAGE 65

52 Test of the Pseudo-Atom Method As already discussed in this method, the dangling bond at the QM/CM boundary is absorbed into a pseudo-atom O*, 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 Si 2 O 7 H 6 and Si 3 O 9 H 6 clusters, as already done for the link-atom method. The optimized structures of these two clusters, obtained from the combined QM/CM method, are shown in Fig. 3-9. For Si 2 O 7 H 6 the interface works very well, as compared with Fig. 3-2. The structure in the QM region is in good agreement with that from all-quantum calculations, and the structure in the CM region is very similar to that from all-classical calculations. Thus, the pseudo-oxygen-atom properly terminates the dangling bond on the boundary. However, a discrepancy appears for the Si 3 O 9 H 6 cluster, as seen by comparison of Fig. 3-9 with Fig. 3-8. The difference in the optimized structures is mainly from the bond angle. For example, the Si-O*-Si bond angle obtained from the combined QM/CM method is 149.1 (see the bottom panel of Fig. 3-9), which is much bigger than the corresponding bond angle (134.0) from the all-quantum calculation (see the upper panel of Fig. 3-8). Such disagreement arises because no bond-angle dependence is built into the current pseudo-atom method. Thus, the electrons near the QM/CM boundary do not feel any difference between structures with different Si-O-Si bond angle. This explains why such an interface works well for the Si 2 O 7 H 6 cluster, in which the Si-O-Si bond angle is 180 at the QM/CM boundary, but fails for the Si 3 O 9 H 6 cluster. Varying the nuclear charges of O* atoms improves the Si-O* bond lengths, but fails to improve the Si-O*-Si bond angles for the same reason. Thus, a Si-O*-Si bond

PAGE 66

53 angle dependent O* effective core potential must be designed if the pseudo-atom method is to be used. The construction of a transferable O* effective potential in a dynamically evolving SiO 2 system is rather complicated. Fortunately, such bond-angle dependence is naturally included in the link-atom method if the link-atoms are constrained on the lines of the Si-O bond on the QM/CM boundary. Based on the preceding considerations, we have used the link-atom method, which is very simple but effective as well, in the studies of the hydrolytic weakening of the amorphous silica. ( a ) ( b ) Figure 3-9. Cluster structures obtained from combined QM/CM calculations using the pseudo-atom method. (a) Cluster shown in Fig. 3-1. (b) Cluster shown in Fig. 3-7.

PAGE 67

54 Hydrolysis of a Two-Membered Silica Ring on the Amorphous Silica Surface Generation of the Amorphous Silica Surface The amorphous silica sample was prepared by cooling a preheated -cristobalite sample via a classical MD simulation that uses the BKS potential (Eq. 3-1) with a correction to suppress an unphysical attraction at short interatomic distance. A bulk sample of -cristobalite, 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 2-VIII). Fig. 3-10 shows the MD results for the O-O, Si-O and Si-Si pair-distribution functions in the amorphous silica bulk. The positions of the first peak in g Si-O (r), g O-O (r) and g Si-Si (r) give the mean Si-O bond length, nearest-neighbor O-O and Si-Si distances to be 1.61 2.60 and 3.09 respectively. The corresponding experimental values from neutron-diffraction data are 1.608 2.626 and 3.077 (Wright, 1994). The density, average Si-O-Si angle and average O-Si-O angle of the resulting amorphous silica sample are 2.25 g/cm3, 146.4 and 109.3, 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 z-direction 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. Sixand five-membered rings are the dominant species. There are also some larger rings and some smaller ones such as four-, three-, and two-membered 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

PAGE 68

55 reduce the local strain. In addition, there are also some one-coordinated O atoms hanging on the surface. These defect sites (dangling bonds and two-membered rings) can be annihilated easily by water attack. Fig. 3-12 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 higher-density region on the surface at around 15 For crystals, such a feature is usually due to lattice contraction on the surface. However, for the amorphous silica, this high-density region is due to the enrichment of small silica rings (two-, threeand four-membered rings) on the surface. The Si-O bond lengths in small rings are usually longer than those in bigger rings because of the enhanced strain in small rings. Thus, the high-density region on the surface does not correspond to shorter bond lengths in the same region, as shown in Fig. 3-13(a). Instead, the average Si-O bond length in the high-density 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. 3-13(a)] is due to those Si-O bonds with one-coordinated O atoms. The enrichment of small rings on the surface also is evidenced by the decreasing average Si-O-Si bond angle on the surface, as shown in Fig. 3-13(b). The sharp decline at the top of the surface in Fig. 3-13(b) is due to those two-membered rings, which have smallest Si-O-Si bond angles. The concentration of two-membered rings on the surface is found to be 0.42/nm 2 which is close to the experimental estimate (0.2-0.4/nm 2 ) (Grabbe et al., 1995)

PAGE 69

56 Figure 3-10. Partial pair-distribution functions for amorphous silica at 300 K.

PAGE 70

57 Figure 3-11. One layer of the amorphous silica surface. The gray and dark gray spheres represent the Si and O atoms, respectively.

PAGE 71

58 Figure 3-12. 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. ( a ) ( b ) Figure 3-13. (a) The average Si-O bond length, and (b) the average Si-O-Si 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.

PAGE 72

59 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 Troullier-Martins pseudopotential (Troullier and Martins, 1993) and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation 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 link-atom approach was used to saturate the valence structure in the QM region (i.e. a hydrogen atom is constrained on the line of a Si-O 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 link-atom 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. 3-1) is known to generate good lattice parameters and elastic constants for quartz and other SiO 2 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 Stillinger-Weber type potential for silica (Watanabe, et al., 1999), in which parameters are fitted to reproduce the cohesive energies of Si-O bonds as a function of the coordination number of the O atom and also the total and

PAGE 73

60 deformation energies of a number of silica clusters. In the Watanabe potential for Si, O systems, the total interaction energy is given by (3-2) 3(,)(,,)sijiijikjfijfijk where (=50 kcal/mol) is an energy parameter, f 2 is the pairwise interaction between i-th and j-th atoms, and f 3 is the three-body term as a function of the i, j, k-th atoms. The Watanabe potential gives much better energetics 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({R s }), was given by U (3-3) '''({})()sssssR sssrrsssrr' and the interaction energy between QM and CM atoms, V({R n },{R s }), was given by V (3-4) ({},{})nsnsnsnsRR The forces acting on QM and CM atoms follow from by Eq. 2-104 and 2-105. After the forces were minimized, the total energy was calculated with the Watanabe potential. U({R s }) and V({R n },{R s }) are given by Uf (3-5) 2'3'''''''''''({})()(,,)ssssssssssssssrfrR (3-6) 23''3'''({},{})()(,,)(,,)nsnsnsnsnsssnnnsnsnsssnnnsVfrfrrrfrrrRR

PAGE 74

61 Simulation Set-up A 12000-atom amorphous silica slab, consisting of four 3000-atom 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. 2-104), while the positions of all the CM atoms were fixed to simplify the problem. [The division of the system is depicted in Fig. 3-14(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 E QM and E QM/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 (3-7) ()()()adsorptionEEsurfacewaternEwaterEsurface where n is the number of water molecules. A negative E adsorption 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 energetics 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. 3-14(b) and (c)]. In the

PAGE 75

62 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 3-14. 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.

PAGE 76

63 Results Adsorption of One Water Molecule on the Two-Membered Ring The calculated adsorption configuration of H 2 O on the amorphous silica surface is depicted in Fig. 3-15 (Only the 31 QM atoms on the surface are shown). In the upper left picture in Fig. 3-15, 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 respectively, as shown in Table 3-1. The average Si-O bond length in amorphous silica is 1.61 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 E b (surf), E a (surf) and E a (water), respectively. The local surface energy increase is estimated as E a (surf) E b (surf) = 1.834 eV. The bonding strength between physisorbed water molecule and the surface can be estimated as E a (surf + water) E a (surf) E a (water) = -1.846 eV, where E a (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 Si-O network is weakened, which is evidence of the hydrolytic weakening. Table 3-1. Several Si-O bond lengths in combined QM/CM model and cluster model (cluster B). The units are in Si(1)-O(1) Si(1)-O(2) Si(1)-O(3) Combined QM/CM surface model 1.84 1.81 1.82 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. 3-16(a) and (b). One is Si 2 O 6 H 4 (cluster A), which is a two-membered ring. The other one is a 31-atom silica cluster terminated by H atoms (cluster B), which is the QM subsystem in system

PAGE 77

64 For cluster B, the Si atom does not make a bond with the O atom in H 2 O (Si(1)-O(1) distanceis 2.56 ), while our combined QM/CM model predicts such a bond. Some Si-O distances are compared with those from the combined QM/CM surface model in Table 3-1. Figure 3-15. Snapshots and potential energy of QM atoms along the dissociation path of one H 2 O 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, O 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. 3-17, 3-18 and 3-19. All density contours are drawn on the plane crossing the Si(1), O(1) and O(2) atoms. As shown in these figures, the water-surface interaction is significantly stronger in the combined QM/CM surface model than in the cluster model. In the combined QM/CM model, the water

PAGE 78

65 molecule is strongly perturbed. Electrons transfer from the O(1) atom to its adjacent Si(1)-O(1) and O(1)-H bond, and also from the H 2 O and Si(1) to O(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. Figure 3-16. Cluster models: (a) one H 2 O molecule bound to a Si 2 O 6 H 4 cluster, (b) one H 2 O molecule adsorbed on a silica cluster (31 Si and O atoms with dangling bonds saturated by 14 H atoms) and (c) two H 2 O molecules adsorbed on the same silica cluster in (b). The gray, dark gray, and light gray spheres represent the Si, O, and H atoms, respectively.

PAGE 79

66 (a) (b) Figure 3-17: Valence electron number density contours on the plane crossing Si1, O1 and O2 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/ 3

PAGE 80

67 (a) (b) Figure 3-18. Contour plot of (a) excess electron density and (b) deficiency electron density in the water-surface system relative to superposed electron densities of the surface and H n n(surface 2 O on the plane crossing Si(1), O(1) and O(2) atoms in the combined QM/CM surface model. () Each contour represents twice/half the density of the adjacent contour lines. The values shown in figures are in units of e/ ) ()()22OHnsurfacenOHnn 3

PAGE 81

68 (a) (b) Figure 3-19. Contour plot of (a) excess electron density and (b) deficiency electron density in the water-cluster system relative to superposed electron densities of the cluster and H n()HO n 2 O on the plane crossing Si(1), O(1) and O(2) atoms in the cluster model. () Each contour represents twice/half the density of the adjacent contour lines. The values shown in figures are in units of e/ 22()nnclusternclusternHO () 3

PAGE 82

69 The adsorption energies are shown in Table 3-2. Following the procedure used above for analyzing both the local surface energy increase and the surface-water bonding strength, the cluster energy increases by 0.14 eV and the cluster-water bonding strength is -0.12 eV (The H atoms that saturate dangling bonds of the cluster are fixed in this calculation). Comparing the energetics 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. 3-15. According to experimental observations, the two-membered 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 3-2. Adsorption energies of one and two H 2 O molecules in the cluster model with or without fixing the positions of H atoms that saturate the dangling bonds of the cluster. Si 2 O 6 H 4 31-atom cluster (free) 31-atom cluster (constrained) 31-atom cluster (constrained with 2 H 2 O) E adsorption (eV) -0.5 -0.12 0.05 -0.23 The reaction energies for physorption and chemisorption of H 2 O molecules are included in Table 3-3. 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

PAGE 83

70 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 H 2 O molecule and chemisorption of one and two H 2 O 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 cut-off 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. 3-7. 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 3-3. The water adsorption energies are substantially overestimated if using the

PAGE 84

71 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 3-3: Adsorption energies of water physisorbed and chemisorbed on a two-membered ring on the amorphous silica surface. All the units are in eV. E ads QM E ads QM/CM (Watanabe) E ads QM/CM (BKS) E ads (Watanabe) E ads (BKS) Physisorption 0.097 0.595 -1.337 0.69 -1.24 Chemisorption (1 H 2 O) -0.863 -0.254 -4.937 -1.12 -5.80 System I Chemisorption (2 H 2 O) -1.234 -0.228 -4.653 -1.46 -5.88 Physisorption 0.229 0.088 -1.178 0.32 -0.95 Chemisorption (1 H 2 O) -0.748 -0.816 -4.209 -1.56 -4.96 System II Chemisorption (2 H 2 O) -1.142 -0.858 -3.827 -2.00 -4.97 Physisorption -0.028 0 -1.020 -0.03 -1.05 Chemisorption (1 H 2 O) -1.059 0 -2.498 -1.06 -3.56 System III Chemisorption (2 H 2 O) -1.603 0 -2.243 -1.60 -3.85 The reaction energy of one-water 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

PAGE 85

72 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 Si-O bonds and remaining hydroxyls. Our result is in reasonable agreement with these experimental data. Adsorption of Two Water Molecules on the Two-Membered 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 SiO 2 (H2O) 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 H 2 O molecules interacting with the amorphous silica surface. The complete reaction pathway is shown in Fig. 3-20 and the reaction energy is shown in Table 3-3. 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 H 3 O on top of the two-membered ring shown in the upper-middle panel of the Fig. 3-20 is +0.7 obtained using Mulliken charge analysis. The final configuration (bottom-right panel in Fig. 3-20) consists of one H 2 O molecule with two silanols on the surface. The same reaction pattern has also been observed on the (0001)

PAGE 86

73 surface of -Al 2 O 3 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. 3-16c and the adsorption energy is shown in Table 3-1.) The water-dimer cooperative reaction may be also energetically favored for other strained Si-O 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 Si-O 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 hydrogen-bonding 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 two-membered ring evolves into two silanols. Without further applying strain, this silanol-covered structure is stable under the presence of additional water. In system III, the adsorption energies for one and two H 2 O molecules on two surface silanols are -0.54 and -0.73 eV, respectively, and no bond rupture occurs. The configuration of two H 2 O

PAGE 87

74 molecules on top of two silanols is shown in Fig. 3-21. Each H 2 O molecule is bound to one silanol via a hydrogen bond, which is the origin of the enhancement of the hydrophilic character near silanol sites. Figure 3-20. 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 two-membered 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

PAGE 88

75 have been observed in STM studies of the Si surface (Chander et al., 1992; Anderson and Khler, 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 (Si-H and Si-OH) can be energetically more favored than on the clean surface. We expect that similar phenomena can happen on a dehydroxylated silica surface. Figure 3-21. Two H 2 O molecules on top of two silanols. The gray, dark gray, and light gray spheres represent the Si, O, and H atoms, respectively. For silica samples pretreated at a temperature below 400, complete rehydroxylation can be achieved (Silanol concentration can be restored to the original

PAGE 89

76 value). (Zhuravlev, 1993; Zhuravlev, 2000; DSouza 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; DSouza and Pantano, 1999). Zhuravlev reported that it took about 5 years to complete rehydroxylation of a silica surface that was pretreated at 900 (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 -cristobalite model. For amorphous silica, charge transfer is local and the electronic response is short-ranged. A simple link-atom 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 energetics 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 Si-O 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 semi-empirical potentials

PAGE 90

77 for silica to ensure the sufficient transfer of information across the QM/CM interface in a chemically realistic multi-scale simulation (Al-Derzi, 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., Al-Derzi, 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 -quartz 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 energetics 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 large-scale dynamical problems, such as crack propagation, which need hundreds of atoms in the QM region. Previous multi-scale simulations that permitted hundreds of atoms in the QM region were based on the tight-binding method (Broughton, et al., 1999), which is over-simplified and not able to describe bond breaking properly. A new attempt has been made to re-parameterize the semi-empirical Hamiltonian based on the input from high-level quantum mechanical method, coupled-cluster (CC) theory (Talor, et al., 2003). The new QM method (so-called transfer Hamiltonian) is able to treat 500-1000 atoms self-consistently and retains the accuracy of CC method. In terms of future studies of the chemo-mechanical 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

PAGE 91

78 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 two-membered ring on the amorphous silica surface. With this model, we have described the long-range 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 energetics for both physisorption and chemisorption processes. The findings of cooperative proton transfer and a series of complex intermediate steps in the two-water 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 long-term 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.

PAGE 92

CHAPTER 4 ENERGETICS AND ELECTRONIC STRUCTURE OF THE CARBON NANOTUBE COMPLEXES Since the discovery of the carbon nanotube in 1991 (Iijima, 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 et al., 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 single-wall 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 79

PAGE 93

80 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 4-1. A graphene sheet. The chiral vector and the translational vector 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 C are connected. The figure corresponds to C, T, and hC T h 4,2h 4,5 2Rd Structure of a Single-Wall Carbon Nanotube A single-wall carbon nanotube (SWCNT) can be described as a graphene sheet rolled into a cylindrical shape so that the structure is one-dimensional with axial symmetry, and in general exhibiting a spiral conformation, called chirality. The chirality can be defined by a chiral vector, given by hC

PAGE 94

81 (4-1) 12(,)hnmnmCaa where a and are the real space unit vectors of a graphene sheet, and n and m are integers (Fig. 4-1). A nanotube can be constructed by rolling a graphene sheet into a cylinder. Fig. 4-1 shows the construction of a (4, 2) nanotube. 1 2a The circumferential length, of the carbon nanotube is just the length of the chiral vector: L 22hLanmC nm (4-2) where is the length of the lattice vector a or a, and a 1 2 1.42a 32.46 The translational vector can be expressed as T (4-3) 112212,ttttTaa where t and are both integers. Since CT, the expressions for t and are given by 1 2t 0h 1 2t 1222,RRnmnmttdd ) (4-4) where is the greatest common divisor of and (2. Rd (2nm )nm 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

PAGE 95

82 nanotube. For instance, potassium or bromine doping enhances conductivity (Lee et al., 1997). Alkali metal and alkaline-earth intercalated fullerenes have also been studied extensively because of their complex phase diagram and superconductivity at temperatures surpassed only by the high-T c 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 energetics and electronic structure of fullerene and nanotube endohedral complexes (Dunlap et al., 1992; Wang et al., 1993; Cioslowski, 1995; Li and Tomnek, 1994; Tomnek 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 C 60 molecules in single-wall carbon nanotubes (Smith et al., 1998; Burteaux et al., 1999; Smith et al., 1999; Smith et al., 2000) and achieved high-yield synthesis of such peapods (2000). More complex metallofullerene peapods ((Gd@C 82 ) n @SWNTs, (La 2 @C 80 ) 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 @.

PAGE 96

83 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 C 60 @SWNT (Hornbaker et al., 2002) and the spatial modulation of the energy band-gap in (Gd@C 82 ) n @SWNT (Lee et al., 2002). A temperature-induced change from p to n conduction in (Dy@C 82 ) n @SWNT (Chiu et al., 2001) and ambipolar field-effect transistor behavior of (Gd@C 82 ) n @SWNT have also been reported (Shimada et al., 2002). An encapsulated C 60 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 C 60 induced impurity states inside the band-gap of the host semiconducting nanotubes by controlling the distance between the tube and C 60 and the type of the encaged metal atom inside C 60 Hybridization of states and the effects of density and orientation of encapsulated C 60 molecules have been analyzed. In addition, I discuss the possibility of superconductivity of potassium-doped peapods. Method and Simulation Setup The electronic structure calculations were performed using density functional theory (DFT) with the local density approximation (LDA). The electron-ion interactions were described by ultra-soft pesudopotentials (Vanderbilt, 1990). The valence wave

PAGE 97

84 functions were expanded in a plane-wave 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 We use two k points in the irreducible Brillouin zone (BZ). All calculations were performed using the Vienna ab-initio 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 C 60 chain. The optimized lattice parameter c is 9.79 for the (10, 10) nanotube, and 12.68 for all zigzag tubes. In order to study (17, 0) peapods with different C 60 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 C 60 molecules with optimized lattice parameters of 12.68 and 21.14 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. 4-2 shows the structures of the (10, 10), (17, 0) a and (17, 0) b peapods. The cohesive energy per C 60 molecule is defined as 60EpeapodEnanotubeENCEN (4-5) The encapsulations of C 60 molecules for (10, 10), (16, 0), (17, 0) a (17, 0) b and (19, 0) peapods are all exothermic, as shown in Table 4-1. The C 60 intermolecular distances for (16, 0), (17, 0) a and (19, 0) are same.

PAGE 98

85 (a) C 60 @ (10, 10) (b) C 60 @ (17, 0) a (c) C 60 @ (17, 0) b Figure 4-2. The structures of (a) (10, 10), (b) (17, 0) a and (c) (17, 0) b nanopeapods. Table 4-1. The cohesive energy per C 60 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) EeV -1.73 -0.78 -1.71 -1.71 -0.84 Manipulation of Fullerene-Induced Impurity States in M@C 60 @SWNTs In an isolated C 60 molecule, there is a three-fold degenerate lowest unoccupied t 1u state. Fig. 4-3 shows the electronic band structures of the C 60 @(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 4-3(b) and (c), the (17, 0) peapod remains a

PAGE 99

86 semoconductor. In the (17, 0) a peapod, the flat energy band inside the band gap is derived from the t 1u state of C 60 There is nearly no energy dispersion for the t 1u -derived band because of the long center-to-center distance d c between two C 60 molecules, while there is a very small energy dispersion in the (17, 0) b peapod. When d c is large, the C 60 chain is broken and not conductive so that the t 1u -derived states behave as impurity states [see Fig. 4-3(b)] and carriers are distributed on the nanotube wall. When d c decreases, these impurity states become impurity bands [see Fig. 4-3(c)] and conductive, resulting in most hole-carriers distributed on the nanotube wall and most electron-carriers on the C 60 chain at a temperature not high enough to make the intrinsic carrier concentration dominant. Figure 4-3. Energy band structures of empty and C 60 encapsulated (17, 0) nanotube. (a) Empty nanotube, (b) C 60 @(17, 0) a (c) C 60 @(17, 0) b only C 60 t 1u -derived impurity bands inside the band gap are shown. The origins of the energies are set at the top of the valence band (E v ). Note the size of the first-Brillouin 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 t 1u bandwidth in the (17, 0) b peapod is not only due to the relatively larger d c but also due to the ordered orientation of the C 60 molecules inside the nanotube. On rotating one of the two C 60 molecules in the unit cell of the (17, 0) b peapod around the

PAGE 100

87 tube axis by 180, we find a doubled t 1u bandwidth and the reduction of the t 1u DOS peak by about 20%, while the total energy essentially remains the same. If d c is about 10 as found in experiment, and the random orientation of C 60 molecules at room temperature is taken into account, we expect a bigger t 1u bandwidth than that shown in Fig. 4-3(c). Figure 4-4. Energy band structures of (a) C 60 @(16, 0) a (b) C 60 @(17, 0) a (c) C 60 @(19, 0) a The C 60 intermolecular distances in these three peapods are same. The origins of the energies are set at the top of the valence band (E v ). While d c and the orientation of C 60 molecules largely determine the t 1u band dispersion, the space between the C 60 and the nanotube wall is a key factor for the energy levels of t 1u -derived states. Our calculations show that the energy levels of t 1u -derived states are the inside conduction band in C 60 @(16, 0), about 2/3 of the band-gap (0.6 eV) above the top of the valence band (E v ) in C 60 @(17, 0) and about 1/3 of the band-gap (0.48 eV) above E v in C 60 @(19, 0) (see Fig. 4-4). The diameters of the (16, 0), (17, 0) and (19, 0) nanotubes are 12.48, 13.26 and 14.80 respectively. It has been reported in a recent paper that energy levels of t 1u -derived states are 2/5 of the band-gap (0.5 eV) above E v in C 60 @(14, 7) (Miyake and Saito, 2003). The diameter of the (14, 7) nanotube is 14.48 It is evident that the t 1u -derived states shift down toward the top of the valence

PAGE 101

88 band as the tube diameter increases from 12.48 to 14.80 A similar trend has also been found in metallic nanotubes (Okada et al., 2001). Fig. 4-5(a) shows the calculated LDOS of the nanotube wall along the axis of the (17, 0) a peapod. Although the total DOS is only projected on the nanotube, the feature from the C 60 t 1u state still appears in the LDOS of the nanotube wall and gradually disappears when the distance to the center of the C 60 molecule increases. This feature is the signature of the hybridization of the states on the tube and the buckyball. This result is consistent with a recent STM study on the semiconducting peapod (Hornbaker et al., 2002). Analysis of hybridized t 1u wavefunctions shows that they are distributed on both C 60 molecules and on the tube wall and spatially varied [see Fig. 4-5(c)-(d)]. Also, they are more localized around the C 60 molecule in the semiconducting peapod than in the metallic peapod, which suggests that it would be relatively more difficult for the STM to image the spatial variation of the LDOS in a metallic peapod. Encapsulated fullerenes or other species inside a semiconducting nanotube can be treated as impurities, whose chemical properties have an important influence on the electronic properties of the host semiconducting nanotube. In principle, a semiconducting peapod can be por n-type, depending on whether the encapsulants are electron acceptors or donors. C 60 has a high electron affinity and thus is an electron acceptor. It is clear from the band structure of C 60 @(17, 0) a [Fig. 4-3(b)] that C 60 t 1u -derived states are deep acceptor states. Encaging a metal atom inside a C 60 molecule can change the electron affinity of the molecule and result in a change of the impurity energy level inside the band gap of the a (17, 0) nanotube.

PAGE 102

89 Figure 4-5. (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 t 1u -derived energy band are shown for (c) a (17, 0) a semiconducting peapod and (d) a (10, 10) metallic peapod. Each contour represents twice/half of the density of the neighboring line. Fig. 4-6 (a)-(d) show the band structures of three metallofullerene peapods [K@C 60 @(17, 0) a Ca@C 60 @(17, 0) a Y@C 60 @(17, 0) a ]. (We performed the spin-polarized calculations for the Y@C 60 @(17, 0) a and found one electron spin per C 60 molecule.) The encaged K, Ca and Y atoms transfer about one, two and three electrons to the C 60 cage, respectively, as estimated by counting the electron occupation numbers of t 1u -derived states. The energy levels of t 1u -derived states in K@C 60 @(17, 0) a [Fig. 4-6(a)] are about 0.1 eV lower than those in C 60 @(17, 0) a [Fig. 4-3(b)]. As the charge states of the encaged ions change from +1, +2 to +3 (the whole system is neutral), three t 1u -derived energy levels downshift toward the valence band, which implies easier electron

PAGE 103

90 transfer from nanotube to encapsulated C 60 s. The experimentally measured electron affinity of C 60 is 2.667 eV (Brink et al., 1995). Our calculated electron affinities for C 60 and K@C 60 are 2.69 and 2.87 eV, respectively. 2 It has been found experimentally that metallofullerenes usually have higher electron affinity than empty C 60 (Boltalina et al., 1997). Therefore, K@C 60 Ca@C 60 and Y@C 60 are better electron acceptors than empty C 60 which results in a downshift of C 60 t 1u -derived states toward the valence band upon encaging these metal atoms. Figure 4-6. Energy band structures of (a) K@C 60 @(17, 0) a (b) Ca@C 60 @(17, 0) a (c) Y@C 60 @(17, 0) a (spin-up), (d) Y@C 60 @(17, 0) a (spin-down), and (e) K@C 60 @(16, 0). The origins of the energies are set at the top of the valence band (E v ). The Fermi level E F is indicated in (a) and (e), and equal to E v in other figures. In (d), one of the t 1u derived states is a shallow acceptor state slightly above E v As already discussed, decreasing the tube-fullerene distance can raise the energy levels of the t 1u -derived states. Fig. 4-6(e) shows the band structure of K@C 60 @(16, 0). The partially occupied t 1u states are below and near the bottom of the conduction band (E c ). According to the band structures shown in Fig. 4-6(c)-(e), Y@C 60 @(17, 0) is a p-type semiconductor and K@C 60 @(16, 0) is a n-type semiconductor. 2 We have used Troullier-Martin pseudopotential and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation potential for electron affinity calculations. For details of the method, see Barnett and Landman, Phys. Rev. B 48, 2081 (1993)

PAGE 104

91 Apparently, the tube diameter and the type of the encaged metal atom are two key factors that determine the energy level and the electron occupation number of impurity states inside the band-gap of the host semiconducting nanotube. Next, consider the consequence of encapsulating other types of metallofullerenes inside a semiconducting nanotube. There are many types of endohedral metallofullerenes that have been synthesized in larger quantity and with higher purification than M@C 60 They have different chemical and electronic properties, depending on the fullerene size, the type and the number of metal atoms encaged (Shinohara, 2000; Cioslowski, 1995). For example, C 82 has a higher first electron affinity than C 60 and its lowest unoccupied molecular obital (LUMO) is nondegenerate in contrast to C 60 If M@C 82 has a charge state of M 3+ C 82 3, the energy level of the singly occupied LUMO+1 will determine whether the peapod is a por n-type semiconductor. The space between C 82 and the nanotube wall is expected to play a key role of determining the energy levels of LUMO and LUMO+1 of C 82 In general, the tube-fullerene distance and the chemical properties (electron affinity and charge state) of the encapsulated metallofullerenes determine the energy level and the electron occupation number of the fullerene-induced impurity states inside the band-gap. Manipulation of these impurity states enables us to find a series of semiconducting peapods with Fermi levels ranging from the valence band to the conduction band, and consequently to control the type of majority carrier (p or n) and carrier density of semiconducting peapods. The family of endohedral metallofullerenes (M x n+ @C 2y n) behaves like a superatom with tunable size and chemical properties. Such a large family of metallofullerenes provides us an opportunity to fine-tune the electronic properties of semiconducting nanotubes to desired conditions.

PAGE 105

92 Electronic Structure of (K x C 60 ) n @SWNTs Potassium atoms also can be encapsulated in the interstitial space between C 60 and nanotube wall (Fig. 4-7). Such encapsulation is exothermic if the K-K distance is not too small. For the (10, 10) peapod, the cohesive energies for adding one additional potassium atom up to five K atoms are .92, -3.02, -3.05, -2.78, and .57 eV, respectively. For the (17, 0) b peapod, the cohesive energies for adding one additional pair of potassium atoms, up to three pairs, are .56, -5.83, and .48 eV. In analogy with the potassium-intercalated C 60 crystal, a doped C 60 chain may show superconductivity if the narrow t 1u band is nearly half filled. With more encapsulated potassium atoms, more electrons are transferred to both C 60 molecules and the nanotube, raising the Fermi level and the t 1u derived bands, as shown in Figures 4-8 and 4-9, while the bands associated with the nanotube remain nearly unchanged. In K x C 60 @(17, 0) peapods, t 1u -derived bands and E F are both shifted into the conduction band. Fig. 4-8(d) and 4-9(d) show the DOS [N(E)] of K 3 C 60 @(10, 10) and (K 3 C 60 ) 2 @(17, 0) b The N(E F ) per C 60 in either the (10, 10) or the (17, 0) b peapod is high and comparable to that of K 3 C 60 crystal in the superconducting phase (about 7.2 states/eV-spin) (Gunnarsson, 1997). (Note that there are two C 60 molecules per unit cell in the (17, 0) b peapod.) The N(E F ) per C 60 in the (17, 0) b peapod is higher than that in the (10, 10) peapod due to the relatively larger C 60 intermolecular distance in the (17, 0) b peapod. Thus, the critical temperature (T c ) for alkali-doped peapods may also be comparable or even higher than that for the alkali-intercalated C 60 crystal. However, both lattice distortion of the C 60 chain due to the weak C 60 -C 60 interaction and incommensurate charge filling of the C 60 molecules due to the alkali atom vacancies are generally expected. How these problems will affect the properties of peapods is still an open question. Future studies on the

PAGE 106

93 structural and phase transformation in these dimensionally constrained systems would be very interesting. (a) (b) Figure 4-7. The geometry of K 3 C 60 @(17, 0) b (a) Side view. (b) Top view. The light and dark gray spheres represent potassium and carbon atoms, respectively.

PAGE 107

94 (a) (b) (c) (d) (e) (f) Figure 4-8. The density of states of (a) C 60 @(10, 10), (b) K 1 C 60 @(10, 10), (c) K 2 C 60 @(10, 10), (d) K 3 C 60 @(10, 10), (e) K 4 C 60 @(10, 10), and (f) K 5 C 60 @(10, 10). The Fermi energies are all shifted to zero.

PAGE 108

95 (a) (b) (c) (d) Figure 4-9. The density of states of (a) (C 60 ) 2 @(17, 0) b (b) (K 1 C 60 ) 2 @(17, 0) b (c) (K 2 C 60 ) 2 @(17, 0) b (d) (K 3 C 60 ) 2 @(17, 0) b The fermi energies are all shifted to zero. Bromine Doping of Single-Wall Carbon Nanotubes: Separating Metallic from Semiconducting Nanotubes The carbon nanotube is a promising material for the future nano-electronics, in which millions of nanotubes can be self-assembled into functional circuits. Such application requires the availability of pure semiconducting carbon nanotubes. However, all known methods of nanotube synthesis (carbon arc, laser vaporization, CVD) produce a mixture of both metallic and semiconducting nanotubes. Thus, the development of

PAGE 109

96 separation techniques is of great technological importance to the progress of nano-electronics. Recently, Chen et al. (2003) have devised a method to achieve such separation. In it, surfactant-stabilized nanotubes are doped with a certain amount of bromine in an aqueous environment. If the binding of bromine with one type of nanotube is stronger than that with the other type, centrifugation of such a solution may result in the abundance of one type of nanotube in either the sediment or the supernatant. In order to quantify the binding strength of bromine with either metallic or semiconducting nanotubes, we have calculated the binding energy of a bromine atom or ion with (10, 0), (9, 0) and (6, 6) carbon nanotubes, using density functional theory (DFT) with both the local density approximation (LDA) and the generalized gradient approximation (GGA). The electron-ion interaction is described by the ultra-soft pesudopotential (Vanderbilt, 1990). The valence wave functions are expanded in plane-wave basis set with an energy cutoff of 286 eV, which has been tested to give good energy convergence. Integration over the irreducible Brillouin zone (BZ) is carried out using two k points. All calculations were performed using the Vienna ab-initio simulation package (VASP) (Kresse and Jurthmueller, 1996a and 1996b). We use the supercell approach to model the Br-nanotube system. The nearest Br-Br distances are 17.040 for (10, 0) and (9, 0) nanotubes, and 17.217 for (6, 6) nanotube. Among top, hollow and bridge sites 3 on the nanotube surface, the Br atom is found to be most stable on the top site. The Br binding energy on the nanotube is given by (4-6) ()()EEnanotubeBrEnanotubeEBr () 3 The top, hollow, and bridge sites refer to sites on top of a carbon atom, the center of a carbon hexagon, and the middle point of a C-C bond, respectively.

PAGE 110

97 The binding energies of a Br atom at the top site on (10, 0), (9, 0) and (6, 6) nanotubes are shown in Table 4-2. The calculations using the local density approximation (LDA) and the generalized gradient approximation (GGA) show a consistent trend. The bromine binds more tightly to metallic nanotubes than to semiconducting nanotubes. Figs. 4-10~ 4-12 show the band structure of (10, 0), (9, 0) and (6, 6) nanotubes before and after bromine adsorption. Clearly, the electron depletion in the valence band of the (9, 0) and (6, 6) nanotubes is much more than that for the (10, 0) nanotube, indicating more electron transfer in metallic nanotubes than in semiconducting ones. However, the charge state of Br in the aqueous solution is unlikely to be zero. We have also calculated binding energies for several negatively charged Br ions with nanotubes. The bromine ions Br and Br 0.5bind with neither metallic nor semiconducting nanotubes. For Br 0.2, the binding energies are .16, -0.41 and .31 eV for (10, 0), (9, 0) and (6, 6) nanotubes, respectively. The binding of a Br ion is stronger for metallic nanotubes than for semiconducting nanotubes, consistent with the behavior of the Br atom. Table 4-2. The cohesive energy per C 60 molecule for (10, 10), (16, 0), (17, 0) a (17, 0) b and (19, 0) peapod. (10, 0) (9, 0) (6, 6) Diameter () 7.8 7.1 8.1 EeV -LDA -0.96 -1.23 -1.17 EeV -GGA -0.70 -0.90 -0.88 Although the adsorption of bromine on the carbon nanotube in the aqueous environment is a complicated dynamical process, we expect that the Br-nanotube binding is stronger and bromine coverage on the surface is higher for metallic nanotubes than for semiconducting nanotubes based on the simple calculations described above and the following arguments.

PAGE 111

98 Figure 4-10. 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. The Fermi energies are all shifted to zero. Figure 4-11. 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. The Fermi energies are all shifted to zero.

PAGE 112

99 Figure 4-12. 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. The Fermi energies are all shifted to zero. First, define as the fractional surface coverage / (4-7) tSS where S is the Br concentration on the nanotube surface, and S t is the total concentration of available adsorption sites. The kinetic equation for Br adsorption and desorption on the nanotube surface can be expressed as / (4-8) addesddtRR where and are adsorption and desorption rates, respectively. In the solution, Br ions are bound to water or surfactant molecules via electrostatic interaction. When a Br ion is near the nanotube surface, it may undergo a dissociative adsorption process, in which the Br ion transfers charge to surrounding solvent molecules while it takes charge from the nanotube surface, resulting in the dissociation of the Br ion from the solvent adR desR

PAGE 113

100 molecules and adsorption of the Br ion on the nanotube surface. At low surface coverage of Br, the expression for the rate equation can be approximately given by (4-9) /(1)exp(/)exp(/)AddtAEkTBEkT DE Here and are the activation energies for Br adsorption and desorption, respectively. where is the concentration of unoccupied adsorption sites with no bromine ion adjacent. A and B are functions that do not depend on and are also independent of at low surface coverage. can be approximately expressed as where is the heat of adsorption (Q). Solving Eq (1) for the steady state () yields AE DEEA /tSSQ/dtd S AE DE DE Q 0 11exp(/)CQk T (4-10) whereC. The relatively higher heat of adsorption of bromine on metallic nanotubes results in higher surface coverage of Br on metallic nanotubes than on semiconducting nanotubes. AB/ Consistent with these calculations are the experimental results, which show a repeated enhancement in the semiconducting nanotube concentration (relative to the metallic concentration) in the Br exposed supernatant with respect to the sediment. Conclusions This chapter reported two important problems in the carbon nanotube research, i.e., the electronic structure of the carbon nanopeapods and bromine doping of the carbon nanotubes, which is an important technique in separation of the metallic and semiconducting nanotubes.

PAGE 114

101 Carbon nanotubes have unusual structural and electronic properties. These properties can be further modified to desired conditions via doping. Encapsulation of fullerenes or metallofullerenes has been shown to dope the host semiconducting nanotubes by inducing impurity states inside the energy band-gap. We have systematically investigated the effects of two key factors, the tube diameter and the type of the encaged metal atom inside C 60 on the energy level and the electron occupation number of the C 60 -induced 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 find that the potassium doping of the peapod significantly increases DOS at Fermi level and makes the K-doped peapod a candidate superconductor. These findings give insight to the stable doping of the semiconducting nanotubes, and the control of the electronic properties of complex nano-scale materials. Another problem of interest is the separation of the metallic from semiconducting carbon nanotubes via bromine doping. Our calculations confirm the experimental results, predicting a stronger binding strength of bromine with metallic nanotubes than that for semiconducting nanotubes. This finding is important to the understanding of the separation mechanism and the design of the separation procedure.

PAGE 115

CHAPTER 5 CONCLUSIONS Theoretical modeling and computer simulation have emerged as powerful tools for studying and predicting material properties. This dissertation gives novel examples of devising computational methods for studying complex material phenomena and even designing new materials with enhanced properties. Such abilities of the theoretical modeling and design of the complex materials owe much to the development of the density functional theory, which raises the role of theory in chemistry, surface science, and materials physics to a new stage. In this work, the use of density functional theory is extended to study chemo-mechanical phenomena with the help of classical mechanics such that the long-range elastic field is incorporated into the theoretical treatment of chemical reactions. Within this hybrid quantum-classical framework, the material phenomena involving in the interaction between different length-scales can be described and predicted. The applications of such a method enabled us to discover a new reaction pathway in the water-silica surface interaction, and thus greatly advanced our understanding of the mechanism of hydrolytic weakening. This approach is very promising in studying complex phenomena in amorphous solids. Density functional theory also is the most frequently used theory in nano-science and technology because of its ability to describe the electronic structure and interatomic forces in systems with hundreds and sometimes thousands of atoms with high accuracy. We have used DFT to study the electronic structure of carbon nanopeapods, and found 102

PAGE 116

103 ways to modify the electronic properties of peapods by changing the space between nanotube and encapsulated fullerenes, and by inserting metal atom(s) inside fullerene cage. This is the first systematic study of the effects of two key factors, the encapsulated metal atoms and the internal space, on the electronic properties of peapods. This is also an example of Materials by Design, in which the materials are designed by constructing fundamental building blocks such that a specific desired function could be achieved. Here, again, theory is proven to be extremely useful in studying and designing complex materials. Nowadays, theory in solid state physics and materials science is not only used to explain the phenomena observed in experiments, but also used to predict phenomena for which experiments are difficult to perform or not even feasible. For example, the reactivity of a particular site on an amorphous surface can be studied by theory while experiments usually can measure average values over all surface sites; the structural, electronic, and transport properties of the nano-scale materials can be studied by theory while the synthesis, fabrication, detection and measurement of these materials are extremely difficult. This dissertation not only reports important method development and applications in studying complex materials, but also demonstrates the important role of theory in the development of materials science and nanotechnology.

PAGE 117

APPENDIX A TOTAL ENERGY OF A PERIODIC SYSTEM IN MOMENTUM-SPACE REPRESENTATION Here we will show how to write each term in Eq. (2-30) in the momentum-space representation. A periodic function can be Fourier transformed between real and momentum-space using following equations: (A-1) ()()expFFiGrGG r 31()()expccFFiGrG drr (A-2) where is the volume of the real-space unit cell. c As discussed in Chapter 2, the Bloch wavefunction can be written as 2,,1,2()exp()cutiiEcikkGGkGr kGr and its normalization condition is (A-3) 3,,1ciidrkkrr Since '3'1cicedrGGrGG ,,1iiccckGkGG (A-4) Denote by the normalized weight of each k point (), such that k 1kk (A-5) 3,ciirdrfkkk Then the electron density in both real and momentum spaces are, respectively, given by 104

PAGE 118

105 ('),,,','1iiiiicfcceGGrkkkGkGkGGr (A-6) 31cicedGrGr r (A-7) The kinetic energy per unit cell is straightforward to write as follows 22,,12eiiiTfckkkGkGkG (A-8) Using the identity 324iedrrGGr (A-9) the Hartree potential ( 3'''Hrrrr V) in momentum-space representation is given by dr 321()()ciHcVVedrGGrGGr 4() (A-10) Thus, the Hartree energy per unit cell can be written as follows 232()1()()422cHHEVdrGGGrr (A-11) Apparently, is divergent when G. We will discuss how to treat this divergent term later. The electron-ion interaction energy per unit cell is given by HE 0 3,,,3()(,,,,',,,'1()()()1()PPeIiislsliilsiPPiiislsliilsEfdrVPNfdrceVPceNRkGrkGrkkkGkGkGRGrrRtrrRt ') where is the number of unit cells in the crystal, is the lattice vector, and t is the basis vector for atom s in the unit cell. Change to, becomes N R s r srRt eIE

PAGE 119

106 ()()(')()3,,,,',,,,'(')3(')(),,,',,,'1()11()iiPPeIiislliilsiiiPPciiislscEfdrceVPceNfdrcceeeVPeNssskGrRtkGrRtkkkGkGkRGGGGtGGRkGrkGrkkkGkGkGGRrr,('),,,',,,',,,'iliPPciiislilsfcceVsGGtkkkGkGkGkGkGG (')i where ()()3,,,',1()PPiPPislsllcVeVPekGrkGrkGkGr dr (A-12) To evaluate the matrix element'V, we need to decompose the plane-waves into angular-momentum channels. Using the identity ,,,PPslkGkG 4()()(21lllmmlPkrYrYkl )lm() (A-13) (21)()()4()()illllllllmlmlmlelijkrPkrijkrYrYkkr where and are spherical Bessel functions and Legendre polynomials, respectively. Thus, lj lP 3',,,''''''',22'','2''014'4()(4)'sin()lPPlslllmlmlmlclPPlslllmlmmlllPPlslllclmlmVdrijrYrYVrijrYrYidrrijrVrjrddYrYrkGkGkGkGkGkGkGkGkGkG ''0'lllmlmmlmlYYkGkGkGkG

PAGE 120

107 Using the identity and 2''''00sin()lmlmllmmddYrYr 21cos4lmllYYPkGkG llmmlkGkG where ()(cos'kGkGkGkG ') 2,,,',04(21)(cos)'PPPPslllsllcVlPjrVrjGkGkGkG rrdr (A-14) Decompose the pseudopotential into a pure local and a non-local part, so that the electron-ion interaction energy per unit cell can be rewritten as (A-15) lcnlceIeIeIEEE where ,3,,1slcPPlceIsssiPPlccssEVNeVRGtGrrRtGG dr (A-16) ,3,1()PPlcPPlcisscVdrVGrG erdr (A-17) and (A-18) ('),,,,',,,',,,'inlcPPnlceIciiislilsEfcceVsGGtkkkGkGkGkGkGG The exchange-correlation energy is given by (A-19) 33XCXCiXCcXCGEdreGrGrrGrGG Now we have derived the expressions forT, and However, and are individually divergent. is apparently divergent when e eIE HEHE XCE 0HEG 0eIEG

PAGE 121

108 0G The reason that is divergent is because its local potential is a pure Coulomb potential of the form 0eIEG Zr for large r, where Z is the core charge. On taking the Fourier transform, the local pseudopotential diverges as 2ZG [see Eq. (A-9)] whenG. We can decompose the local potential into a Coulomb term and a remaining part, 0 () ) ionG ,,,,,,()()()PPlcPPlcnonCoulombPPCoulombsssPPlcnonCoulombsssVVVZVrrrrt (A-20) r The Coulomb contribution to the electron-ion interaction energy per unit cell in the momentum space representation is given by ,33300(1()''1''cCoulombPPCoulombeIcsssscssescsEVZdrdrZQdrGGrrtrt 0)G (A-21) where Q is the total charge of the electrons in one unit cell. The G term of the Hartree energy is given by e 0 33333(0(0)(0)211('')'2'11('')'''2'''cHHcecEGVdrdrdrQdrdrGGrrrrrrr (A-22) ''' In analogy with Eq. (A-21), the G term of is given by 0 ionE '''1(0)sIssscsZEQtt (A-23) s

PAGE 122

109 where Q is the total charge of the ions. Thus, the sum of equations (A-21), (A-22) and (A-23) is I 333'''(0)(0)(0)11('')'''''2'''CoulombeIHionsseIsscssEEEZZQdrdrdrQGGGrrtrrtt sssI (A-24) If the net charge of the unit cell is zero (), the Eq. (A-24) becomes eQQ 333'''(0)(0)(0)11('')'''''2'''0CoulombeIHionssesscssCoulombeaveEEEZZQdrdrdrQVGGGrrtrrtt sss (A-25) where V is the cell average of the electrostatic Coulomb potential, which is set to be zero. Thus, for a charge-neutral system, the divergent Coulomb contributions to the total energy from, and cancel exactly. Coulombave 0G eIE HE ionE The ion-ion interaction energy is given by the Ewald summation (Ewald 1921), which consists of a real-space and a reciprocal-space summation. The term is divided between these two summations. Therefore, when calculating the Coulomb energy of the ionic system, the term in the reciprocal-space summation should be omitted, and two terms should be added to remove the remaining part of the average electrostatic potential contribution in the lattice energy. The correct ion-ion interaction energy is given by 0G 0G ''',''2'222012241expcos4{}ssionssssssssssccerfcEZZRGttRttRGGttG (A-26)

PAGE 123

110 where is a parameter controlling the convergence of the Ewald summation, and is the complementary error function, given by (erfcx ) 22expxerfcxtdt (A-27) In summary, the final expression for each energy term per unit cell in Eq. (2-30) is given by 22,,12eiiiTfckkkGkGkG (A-28) 220()42cHEGGG (A-29) (A-30) XCcXCEGG G ,0('),,,,',,,',,,'3,1siPPlceIcssiPPnlcciiislilssPPlcesscsEeVfcceVZQdrVsGtGGGtkkkGkGkGkGkGGGGrrt (A-31) ''',''2'222012241expcos4{}ssionssssssssssccerfcEZZRGttRttRGGttG (A-32) where 31cicedGrGr r (A-33) (A-34) ('),,,','iiiiifcceGGrkkkGkGkGGr

PAGE 124

111 ,3,1()PPlcPPlcisscVdrVGrG er (A-35) ,,2,,,',04(21)(cos)'PPnlcPPnlcslllsllcVlPjrVrjGkGkGkG rrdr (A-36) and ()(cos'kGkGkGkG ') (A-37) The last term in the Eq. A-31 is the G component of the first term excluding the Coulomb energy. 0

PAGE 125

APPENDIX B POPULATION ANALYSIS BASED ON ELECTROSTATIC POTENTIALS Population analysis is often used to determine the number of electrons associated with a particular atom or orbital. Some popular methods include Mulliken analysis, natural population analysis, topological electron density analysis, etc (see a review by Bachrach, 1994). Here, we discuss briefly a different approach that is based on electrostatic potentials (Chirlian and Francl, 1987). The value of the electrostatic potential in real space can be obtained easily in any quantum mechanical calculation. An alternative method for obtaining the atomic charge is to fit the quantum mechanical potential to a potential generated by point charges centered on the atomic nuclei. This potential V is given by M 1NjMjjqVrrR (B-1) where is the number of atoms in the system, and is the charge associated with jth atom. Although it is possible to fit the quantum mechanical potential to multipoles of the atomic charge (William, 1988, 1991), we will confine our discussion to the monopole expansion alone. N jq To fit the quantum mechanical potential accurately, one must select a large number of grid points, on each of which the quantum mechanical potential is fitted. Since only the monopole expansion is used, these points must be outside the volume within the van der Waals radii of the molecules of interest. 112

PAGE 126

113 Now, we define the function y, (B-2) 2MkiMiyqVVrr ij where M is the total number of selected points, and V is the quantum mechanical potential on point i. A fit can be obtained if one minimizes y. Using a Lagrangian multiplier method, the Lagrangian function z is defined as ir (B-3) jjzqyqgq where is the Lagrangian multiplier, and is the constraint such that the total molecular charge is reproduced, i. e., g (B-4) 1NjjjgqqQ where Q is the net charge of the system. The best fit is obtained by solving the equations 0z (B-5) and 0kzq (B-6) Using Eq. (B-1) to (B-4), we obtain 211MNNjjijijjijqzqVqQrrR (B-7) (N+1) linear equations can be obtained from Equations B-5 and B-6,

PAGE 127

114 (B-8) 1122,0NMjijjiNjjqBkAkkNqQ 1 where MiiikVAkrrR (B-9) and 1|||ijijikBkrRrR | (B-10) The atomic charges can be obtained by solving these equations using standard library routines. 12,,,Nqqq

PAGE 128

LIST OF REFERENCES Akagi, K. and Tsukada, M., 1999, Surf. Sci. 438, 9. Allen, M. P., and Tildesley, D. J., 1987, Computer Simulation of Liquids (Oxford) Andersohn, L. and Khler, U., 1993, Surf. Sci. 284, 77. Antes, I. and Thiel, W., 1998, in Combined Quantum Mechanical and Molecular Mechanical Methods, edited by Gao, J., and Thompson, M. A., (American Chemical Society, Washington DC). Antes, I. and Thiel, W., 1999, J. Phys. Chem. A 103, 9290. Ashcroft, N. W., and Mermin, N. D., 1976, Solid State Physics (Saunders College Publishing, Fort Worth). Bachrach, S. M., 1994, in Reviews in Computational Chemistry, Vol. 5, edited by Lipkowitz, K. B., and Donald, D. B., (VCH Publisher, Inc., New York, NY) Barnett, R. N., and Landman, U., 1993, Phys. Rev. B 48, 2081. Becke, A. D., 1992, J. Chem. Phys. 96, 2155. Bethune, D. S., Johnson, R. D., Salem J. R., Devries, M. S., Yannoni C. S., 1993, Nature 366, 123. Bolis, V., Fubini, B., Marchese, L., Martra, G., and Costa, D., 1991, J. Chem. Soc. Faraday Trans. 87, 497. Boltalina, O. V., Ioffe, I. N., Sorokin, I. D., and Sidorov, L. N., 1997, J. Phys. Chem. A 101, 9561. Brink, C., Anderson, L. H., Hvelplund, P., Mathur, D., Voldstad, J. D., 1995, Chem. Phys. Lett. 233, 52. Broughton, J. Q., Abraham, F. F., Bernstein, N., and Kaxiras, E., 1999, Phys. Rev. B. 60, 2391. Bunker, B. C., Haaland, D. M., Ward, K. J., Michalske, T. A., Smith, W. L., Binkley, J. S., Melius, C. F., and Balfe, C. A., 1989a Surf. Sci. 210, 406. 115

PAGE 129

116 Bunker, B. C., Haaland, D. M., Michalske, T. A., and Smith, W. L., 1989b, Surf. Sci. 222, 95. Burteaux, B., Claye, A., Smith, B. W., Monthioux, M., Luzzi, D. E., Fischer, J. E., 1999, Chem. Phys. Letts 310, 21. Ceperley, D. M., and Alder, B. J., 1980, Phys. Rev. Lett. 45, 566. Ceresoli, D., Bernasconi, M., Iarlori, S., Parrinello, M., and Tosatti, E., 2000, Phys. Rev. Lett. 84, 3887. Chaban, G. M. and Gerber, R. B., 2001, J. Chem. Phys. 115, 1340. Chadi, D. J., and Cohen, M. L., 1973, Phys. Rev. B. 8, 5747. Chander, M., Li, Y. Z., Patrin, J. C., and Weaver, J. H., 1992, Phys. Rev. B 48, 2493. Chen, Z., Du, X., Du, M.-H., Ranchen, C. D., Cheng, H.-P., and Rinzler, A. G., 2003, Nano Letters, 3, 1245. Cheng, H.-P., Barnett, R. N., and Landman, U., 2002, J. Chem. Phys. 116, 9300. Chiang, C-M., Zegarski, B. R., and Dubois, L. H., 1993, J. Phys. Chem. 97, 6948. Chirlian, L. E., and Francl, M. M., 1987, J. Comput. Chem. 8, 894. Chiu, P. W., Gu G., Kim G. T., Philipp G., Roth S., Yang S. F., Yang S., 2001, Appl. Phys. Lett. 79, 3845. Chuang, I-S., and Maciel, G. E., 1996, J. Am. Chem. Soc. 118, 401. Chuang, I-S., and Maciel, G. E. 1997, J. Phys. Chem. B 101, 3052. Cioslowski, J., 1995, Electronic Structure calculations on Fullerenes and Their Derivatives (Oxford University Press, Inc., New York) Cui, Q., and Karplus, M., 2002, J. Phys. Chem. B 106, 1768. Dillon, A. C., Jones, K. M., Bekkadahl, T. A., Kiang, C. H., Bethune, D. S., and Heben, M. J., 1997, Nature (London) 386, 377. Dreizler, R. M., and Gross, E. K. U., 1990, Density Functional Theory (Springer-Verlag, Berlin) Dresselhaus, M. S., Dresselhaus, G., and Eklund, T. A., 1996, Science of Fullerenes and Carbon Nanotubes (Academic Press, San Diego) DSouza, A. S., and Pantano, Carlo G., 1999, J. Am. Ceram. Soc. 82, 1289.

PAGE 130

117 Du, M.-H., and Cheng, H.-P., 2003a, Int. J. Quant. Chem. 93, 1. Du, M.-H., and Cheng, H.-P., 2003b, Phys. Rev. B, 68, 113402. Du. M.-H., Kolchin, A., and Cheng, H.-P., 2003, J. Chem. Phys., 119, 6418. Dunlap, B. I., Ballester, J. L., and Schmidt, P. P., 1992, J. Phys. Chem, 96, 9781. Eichinger, M., Tavan, P., Hutter, J. and Parrinello M., 1999, J. Chem. Phys. 110, 10452. Ewald, P. P., 1921, Ann. Phys. (Leipzig) 64, 253. Fan, X., Dickey E. C., Eklund P. C., Williams K. A., Grigorian L., Buczko R., Pantelides S. T., Pennycook S. J., 2000, Phys. Rev. Lett. 85, 4621. Fermi, E., 1927, Rend. Accad. Lincei 6, 602. Fermi, E., 1928a, Z. Phys. 48, 73. [A translation into English may be found in March (1975).] Fermi, E., 1928b, Rend. Accad. Lincei 7, 342. Feynman, R. P., 1939, Phys. Rev. 56, 340. Field, M. J., Bash, P. A., and Karplus, M., 1990, J. Comput. Chem., 11, 700. Fubini, B., 1998, The surface properties of silica, edited by Legrand, A. P., (John Wiley & Sons Ltd, Chichester, England) Garofalini, S. H., 1990, J. Non-Cryst. Solids 120,1. Grabbe, A., Michalske, T. A., and Smith, W. L., 1995, J. Phys. Chem. 99, 4648. Griggs, D., 1967, Geophysical Journal of Royal Astronomical Society 14, 19. Gunnarsson, O., 1997, Rev. Mod. Phys. 69, 575. Hall R. J., Hindle S. A., Burton N. A. and Hillier. I. H., 2000, J. Comput. Chem. 21, 1433. Hamann, D. R., 1996, Phys. Rev. Lett. 76, 660. Hammer, B., Jacobsen K. W., and Norskov, J. K., 1993, Phys. Rev. Lett. 70, 3971. Hammer, B. and Scheffler, M., 1995, Phys. Rev. Lett. 74, 3487. Heine, V., and Cohen, M. L., 1970, Solid State Physics Vol. 24, p. 37. Hellmann, H., 1937, Einfuhrung in die Quantumchemie (Deuticke, Leipzig)

PAGE 131

118 Hirahara, K., Suenaga, K., Bandow, S., Kato, H., Okazaki, T., Shinohara, H., and Iijima, S., 2000, Phys. Rev. Lett. 85, 5384. Hochella, M. F. Jr., and White, A. F., 1990, in Review in Mineralogy and Geochemistry, Vol. 23, edited by Ribbe, P. H. (BookCrafters, Inc., Chelsea, Michigan) Hohenberg, H., and Kohn, W., 1964, Phys. Rev. 136, 864B. Hornbaker, D. J., Kahng, S.J., Misra, S., Smith, B. W., Johnson, A. T., Mele, E. J., Luzzi, D. E., and Yazdani, A., 2002, Science 295, 828. Huff. N. T., Demiralp, E., Cagin, T., and Goddard, W. A., III, J. Non-Cryst. Solids 253, 133 (1999) Iarori, S., Ceresoli, D., Bernasconi, M., Donadio, D., and Parrinello, M., 2001, J. Phys. Chem. B 105, 8007. Ihm, J., Zunger A., and Cohen, M. L., 1979, J. Phys. C 12, 4409. Iijima, S., 1991, Nature 354, 56. Jones, R. O., and Gunnarsson, O., 1989, Rev. Mod. Phys. 61, 689. Kairys, V. and Jensen, J. H., 2000, J. Phys. Chem. A 104, 6656. King-Smith, R. D., Payne, M. C., and Lin, J. S., 1991, Phys. Rev. B 44, 13063 Kleinman, L., and Bylander, D. M., 1982, Phys. Rev. Lett. 4, 1425. Kobayashi, K., and Nagase, S., 1996, Chem. Phys. Lett. 262, 227. Kobayashi, K., and Nagase, S., 1998, Chem. Phys. Lett. 282, 325. Kohn, W., and Sham, L. J., 1965, Phys. Rev. 140, 1133A. Kresse, G., and Furthmueller, 1996a, J. Comput. Mat. Sci. 6, 15. Kresse, G. and Furthmueller, 1996b, Phys. Rev. B 54, 11169. Kresse, G. and Hafner, J., 1994, J. Phys.: Condens. Matter 6, 8245. Kroto, H. W., Heath, J. R., OBrien, S. C., Curl, R. F. and Smalley, R. E., 1985, Nature 318, 162. Laasonen, K., Car, R., Lee, C., and Vanderbilt, D., 1991, Phys. Rev. B 43, 6796 Lee, J., Kim H., Kahng S. J., Kim G., Son Y. W., Ihm J., Kato H., Wang Z. W., Okazaki T., Shinohara H., Kuk Y., 2002, Nature, 415, 1005.

PAGE 132

119 Lee, R. S., Kim, H. J., Fischer, J. E., Thess, A., and Smalley, R. E., 1997, Nature 388, 17. Li, Y. S., and Tomnek, D., 1994, Chem. Phys. Lett. 221, 453. Liu, C., Fan, Y. Y., Liu, M., Cong, H. T., Cheng, H. M., and Dresselhaus, M. S., 1999, Science 286, 1127. Makov, G., and Payne, M. C., 1995, Phys. Rev. B 51, 4014. March, N. H., 1975, Self-Consistent Fields in Atoms (Oxford: Pergamon). Martel, R., Schmidt, T., Shea, H. R., Hertel, T., and Avouris, Ph., 1998, Appl. Phys. Lett. 73, 2447. Martins, J. L., and Cohen, M. L., 1988, Phys. Rev. B 37, 6134. Michalske, T. A., and Bunker, B. C., 1984, J. Appl., Phys. 56, 2687. Miyake, T., and Saito, S., 2002, Phys. Rev. B, 65, 165419. Miyake, T., and Saito, S., 2003, Solid State Commun., 125, 201. Miyamoto, Y. Rubio, A., Blase, X., Cohen, M. L., Louie, S. G., 1995, Phys. Rev. Lett. 74, 2993. Monkhorst, H. J., and Pack, J. D., 1976, Phys. Rev. B 13, 5188. Mukhopadhyay, I., Kawasaki S., Okino F., Govindaraj A., Rao C. N. R., Touhara H., 2002, Physica B 323, 130. Nicoll, R. A., Hindle, S. A., MacKenzie G., Hillier, I. H., and Buton, N. A., 2001, Theor. Chem. Acc. 106, 105. Ogata, S., Lidorikis, E., Shimojo, F., Nakano, A., Vashishta, P., and Kalia, R., 2001, Comput. Phys. Commun. 138, 143. Okada, S., Saito, S., and Oshiyama, A., 2001, Phys. Rev. Lett. 86, 3835. Parr, R. G., and Yang, W., 1989, Density Functional Theory of Atoms and Molecules (Oxford, New York) Payne, M. C., Teter, M. P., Allan, D. C., Arias, T. A., and Joannopoulos, J. D., 1992, Rev. Mod. Phys. 64, 1045. Penev, E., Kratzer, P., and Scheffler, M., 1999, J. Chem. Phys. 110, 3986. Perdew, J. P., Burke, K., and Ernzerhof, M., 1996, Phys. Rev. Lett. 77, 3865.

PAGE 133

120 Perdew, J. P., Chevary, J. A., Vosko, S. H., Jackson, K. A., Pederson, M. R., Singh, D. J., and Fiolhais, C., 1992, Phys. Rev. B 46, 6671. Philipsen, P. H. T., te Velde, G., and Barerends, E. J., 1994, Chem. Phys. Lett. 226, 583. Phillips, J. C., 1958, Phys. Rev. 112, 685. Pichler, T., Kuzmany, H., Kataura, H., and Achiba, Y., 2001, Phys. Rev. Lett. 87, 267401. Proynov, E. I., Ruiz, E., Vela, A., and Salahub, D. R., 1995, Int. J. Quantum Chem. S29, 61. Reuter, N., Dejaegere, A., Maigret, B., and Karplus, M., 2000, J. Phys. Chem. A 104, 1720. Saito, R., Dresselhaus, G., and Dresselhaus, M. S., 1998, Physical Properties of Carbon Nanotubes, (Imperial College Press, London) Sauer, J., and Sierka, M., 2000, J. Comput. Chem. 21, 1470. Shimada, T., Okazaki T., Taniguchi R., Sugai T., Shinohara H., Suenaga K., Ohno Y., Mizuno S., Kishimoto S., Mizutani T., 2002, Appl. Phys. Lett. 81, 4067. Shimoda, H., Gao B., Tang X. P., Kleinhammes A., Fleming L., Wu Y., Zhou O., 2002, Physica B 323, 133. Shinohara, H., 2000, Rep. Prog. Phys. 63, 843. Slater, J. C., 1951, Phys. Rev. 81, 385. Smith, B. W., and Luzzi, D. E., 2000, Chem. Phys. Lett. 321, 169. Smith, B. W., Luzzi, D. E., and Achiba, Y., 2000, Chem. Phys. Lett. 331, 137. Smith, B. W., Monthioux, M., and Luzzi, D. E., 1998, Nature 396, 323. Smith, B. W., Monthioux, M., and Luzzi, D. E., 1999, Chem. Phys. Lett 315, 31. Smith, B. W., Russo, R. M., Chikkannanavar, S. B., and Luzzi, D. E., 2002, J. Appl. Phys. 91, 9333. Stevenson S., Rice G, Glass T, Harich K, Cromer F, Jordan M.R., Craft J, Hadju E, Bible R, Olmstead M. M., Maitra K, Fisher A. J., Balch A. L., Dorn H.C., 1999, Nature 401, 55. Suenaga, K., Tenc, M., Mory, C., Colliex, C., Kato, H., Okazaki, T., Shinohara, H., Hirahara, K., Bandow, S., and Iijima, S., 2000, Science 290, 2280.

PAGE 134

121 Tan, S., Verschueren, A., and Dekker, C., 1998, Nature (London) 393, 49. Taylor, C. E., Cory, M., Bartlett, R. J., and Thiel, W., 2003, Comput. Mat. Sci. 27, 204 (2003). Thomas, L. H., 1927, Proc. Camb. Phil. Soc. 23, 542. [Reprinted in March (1975).] Titmuss, S. J., Cummins, P. L., Bliznyuk, A. A., and Rendell, A. P., and Gready, J. E., 2000, Chem. Phys. Lett. 320, 169. Tomnek, D., and Li, Y. S., 1995, Chem. Phys. Lett. 243, 42. Troullier, N., and Martins, J. L., 1991a, Phys. Rev. B 43, 1993. Troullier, N., and Martins, J. L., 1991b, Phys. Rev. B 43, 8861. Vanderbilt, D., 1990, Phys. Rev. B 41, 7892. Vavro, J., Llaguno M. C., Satishkumar B. C., Luzzi D. E., Fischer J. E., 2002, Appl. Phys. Lett. 80, 1450. Vosko, S. H., Wilk, L., and Nusair, M., 1980, Can. J. Phys. 58, 1200. Vosko, S. H., and Wilk, L., 1982, J. Phys. C 15, 2139. van Beest, B. W. H., Kramer, G. J., van Santen, R. A., 1990, Phys. Rev. Lett. 64, 1955. Walsh, T. R., Wilson, M., Sutton, A. P., 2000, J. Chem. Phys. 113, 9191. Wang, Y., Tomnek, D., and Ruoff, R. S., 1993, Chem. Phys. Lett. 208, 79. Watanabe, T., Fujiwara, H., Noguchi, H., Hoshino, T., and Ohdomari, I., 1999, Jpn. J. Appl. Phys. Part 2, 38, L366. Wentzcovitch, R. M., and Martins, J. L., 1991, Solid State Commun. 78, 831. William, D. E., 1988, J. Comput. Chem. 9, 745. William, D. E., 1991, in Reviews in Computational Chemistry, Vol. 2, edited by Lipkowitz, K. B., and Boyd, D. B., (VCH Publisher, Inc., New York, NY) Witbrodt, J. M., Hase, W. L., and Schlegel, H. B., 1998, J. Phys. Chem. B 102, 6539. Wright, A. C., 1994, J. Non-Cryst. Solids 179, 84. Yildirim, T., Zhou, O., and Fischer, J. E., 2000, in The Physics of Fullerene-Based and Fullerene-Related Materials, edited by Andreoni, W. (Kluwer Academic Publishers, The Netherlands)

PAGE 135

122 Yin, M. T., and Cohen, M. L., 1982a, Phys. Rev. B 25, 7403. Yin, M. T., and Cohen, M. L., 1982b, Phys. Rev. B 26, 3259. Zhang, Y., Lee, T.-S., and Yang, W., 1999, J. Chem. Phys. 110, 46. Zhuravlev, L.T., 1993, Colloids and Surfaces A 74, 71. Zhuravlev, L.T., 2000, Colloids and Surfaces A 173, 1.

PAGE 136

BIOGRAPHICAL SKETCH Mao-Hua Du was born on February 3, 1976 and grew up in Jiading, a small town near the city of Shanghai in China. Mao-Hua enjoyed his time in Jiading, which was a small and quiet town in the 1980s, surrounded by the rural area near the Yangzi River and in the vicinity of the great city of Shanghai. He went to schools in Jiading and developed interests in history and geography. However, Mao-Hua decided to keep his favorite, history, as a hobby, instead of a career. He went to Shanghai to attend Fudan University, majoring in Physics. Mao-Hua graduated in the summer of 1998 with a B. S. degree in Physics. Mao-Hua left his homeland and came to the United States to pursue his doctoral degree. He has been working with Professor Hai-Ping Cheng at the University of Florida for the last five years. During his stay here, he visited France for five weeks in the summer of 2000 to attend the Les Houches summer school. He gained much experience in programming for numerical computation and molecular visualization at UF. He also served one year as the Vice President of the Friendship Association of Chinese Students and Scholars at the University of Florida. 123


Permanent Link: http://ufdc.ufl.edu/UFE0001540/00001

Material Information

Title: Theoretical Modeling and Design of Complex Materials
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0001540:00001

Permanent Link: http://ufdc.ufl.edu/UFE0001540/00001

Material Information

Title: Theoretical Modeling and Design of Complex Materials
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0001540:00001


This item has the following downloads:


Full Text












THEORETICAL MODELING AND DESIGN OF COMPLEX MATERIALS


By

MAO-HUA 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

Mao-Hua Du

































This dissertation is dedicated to my late grandfather, Mei-Lin 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 Hai-Ping Cheng, for her guidance and encouragement throughout my

graduate career. Her open-minded 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. Lin-Lin 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 parents-in-law, 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
Chemo-Mechanical Phenomena in Amorphous Silica.......................................3
Nanopeapods: M materials by Design.......................... ........... ............... 3

2 M E T H O D S ............................................................................. 5

Total-Energy Pseudopotential Calculations ...................................... ............... 5
Born-Oppenheimer Approximation...................... ... ......................... 6
Density Functional Theory ...................... ......... ............................. 6
Local-Density Approximation..................... ...... ............................ 9
G eneralized G radient Correction...................................... ........ ............... 10
Pseudopotential Approxim ation ....................................................................... 10
The H ellm ann-Feym 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 ual-Space Form alism .................... ........................ .. ............... ... 20
Density Functional Theory for Finite Systems: A Dual-Space 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 link-atom m ethod ....................................................... .....................35
T he pseudo-atom m ethod ............................................................................36

3 MULTISCALE SIMULATIONS OF AMORPHOUS SILICA.............................37

In tro d u ctio n ........................................................................................................... 3 7
Test of the Link-Atom Method ........... .. ......... ................... 42
Test of the Pseudo-Atom M ethod.........................................................................52
Hydrolysis of a Two-Membered Silica Ring on the Amorphous Silica Surface........54
Generation of the Amorphous Silica Surface........................... ............... 54
Hybrid QM /CM Approach ............................................................................ 59
Sim ulation Set-up .......................... ............ ............... .... ....... 61
R results .................... .... ..... ...... ....... ........ ... ..................... .......... ...... 63
Adsorption of One Water Molecule on the Two-Membered Ring....................63
Adsorption of Two Water Molecules on the Two-Membered 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 Single-W 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 Fullerene-Induced Impurity States in M@C60@SWNTs .........85
Electronic Structure of (KxC60)n@ SW NTs ........................................................ 92
Bromine Doping of Single-Wall 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 MOMEMTUM-SPACE
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

3-1. Several Si-O bond lengths in combined QM/CM model and cluster model (cluster
B ) ....................................... ......................................... 6 3

3-2. 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

3-3. The adsorption energies of water physisorbed and chemisorbed on amorphous silica
surface. .............................................................................7 1

4-1. 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

2-1. The division of a system into QM and CM regions. ................ .............. 31

3-1. The structure of the Si207H6 cluster. ........................................................................42

3-2. The structures of the cluster shown in Fig. 3-1, obtained from (a) all-quantum
calculation, (b) combined QM/CM calculation using the link-atom method, and (c)
all-classical calculation. ......................... .................. .................. ........... 43

3-3. The structures of the training system with iso-density surfaces, obtained from (a) all-
quantum calculation, and (b) combined QM/CM calculation........................ 45

3-4. 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

3-5. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. .........48

3-6. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies. .........49

3-7. The structure of the Si309H6 cluster. ................................ .. ......................... 50

3-8. The structures of the cluster shown in Fig. 3-7, obtained from (a) all-quantum
calculation, (b) combined QM/CM calculation using the link-atom method, and (c)
all-classical calculation (bottom panel). ................ .......................... ............... 51

3-9. Cluster structures obtained from combined QM/CM calculations using the pseudo-
atom m ethod ...................................................... ................. 53

3-10. Partial pair-distribution functions for amorphous silica at 300 K ..........................56

3-11. One layer of the amorphous silica surface. ................................... ............... 57

3-12. Number density of atoms as a function of the distance from the center of the slab. 58

3-13. (a) The average Si-O bond length, and (b) the average Si-O-Si bond angle, as
functions of the distance from the center of the slab. ............................................58

3-14. A schematic picture of a silica surface (side view). .............................................. 62









3-15. Snapshots and potential energy of QM atoms along the dissociation path of one
H20 m olecule on a silica surface. ............................... .................................. 64

3-16. C lu ster m models .........................................................................65

3-17: 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

3-18. Contour plot of(a) excess electron density and (b) deficiency electron density in
the water-surface 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

3-19. Contour plot of (a) excess electron density and (b) deficiency electron density in
the water-cluster 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

3-20. Snapshots and potential energy of the QM atoms along the water dissociation path
("hydrogen relay process") on a silica surface. ............................................. 74

3-21. Two H20 molecules on top of two silanols. ...................................................75

4-1. A graphene sheet. The chiral vector and the translational vector are also shown.....80

4-2. The structures of (a) (10, 10), (b) (17, 0)a, and (c) (17, 0)b nanopeapods ...............85

4-3. Energy band structures of empty and C60 encapsulated (17, 0) nanotube. ...............86

4-4. Energy band structures of(a) C60@(16, 0)a, (b) C60@(17, 0)a, (c) C60@(19, 0)a. ....87

4-5. (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 tiu-derived energy band are shown for (c) a (17, 0)a semiconducting
peapod and (d) a (10, 10) m etallic peapod. ........................................ ..................89

4-6. Energy band structures of(a) K@C60@(17, 0), (b) Ca@C60@(17, 0), (c)
Y@C60@(17, O)a (spin-up), (d) Y@C60@(17, 0)a (spin-down), and (e)
K @ C 6o@ (16, 0). ......................................................................90

4-7. The geometry ofK3C6o@(17, )b. ................ ............................................. 93

4-8. 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

4-9. 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

4-10. 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









4-11. 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

4-12. 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

Mao-Hua Du

December 2003

Chair: Hai-Ping 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 two-membered 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 potassium-doped

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

mid-1960s (Hohenberg and Kohn, 1964; Kohn and Sham, 1965). In that work, they

proved that the ground-state many-electron problem can be replaced by an exactly

equivalent set of self-consistent one-electron equations (Kohn-Sham equations); and that

the many-body part of the problem can be included by an effective exchange-correlation

potential. Therefore, the complex many-body problem is essentially reduced to the









problem of a non-interacting inhomogeneous electron gas in an effective potential. If one

knows all of the appropriate electronic interactions, for example, the electron-electron

electrostatic interaction, the exchange-correlation, the ion-ion interaction, and the

electron-ion 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 cost-effective first-principles method and has reproduced

a variety of ground-state 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 fullerene-based hybrid materials, in which simple structures with different properties

are self-assembled to form complicated structures with enhanced properties.

Chemo-Mechanical 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 so-called chemo-mechanical 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 chemo-mechanical phenomena are still in a primitive state

because the chemical processes at nanoscale and the mechanical properties at meso- or

macroscopic-scale 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 link-atom method (Du and Cheng, 2003a). This multi-scale

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, zero-dimensional dots or nanocrystals, one-dimensional 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 self-assembly 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 fullerene-induced states inside the band gap of the

host-semiconducting 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 potassium-doped 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 total-energy 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 total-energy calculations is given. In this chapter, I first

introduce the Born-Oppenheimer approximation, density functional theory, the

pseudopotential method, the Hellmann-Feynman theorem, and molecular dynamics. An

excellent review of the total-energy 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 plane-wave based implementation is discussed.

Finally, I describe the hybrid method that combines DFT with classical potential

functions (Du and Cheng, 2003a).

Total-Energy Pseudopotential Calculations

To perform accurate and efficient total-energy calculations, a series of

simplifications and approximations is needed. These include the Born-Oppenheimer









approximation to separate the electronic and nuclear wave-functions, density functional

theory to model the electron-electron interaction, and the pseodopotential method in

conjunction with a plane-wave basis set to model the electron-ion interaction.

Born-Oppenheimer 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 many-body wave-functions. This so-called Born-Oppenheimer

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 many-body 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 total-energy functional of the electron density yields the ground-state

energy of the system, and the corresponding electron density is the exact ground-state

electron density. Kohn and Sham (1965) then showed that the many-electron problem can

be replaced by an exactly equivalent set of self-consistent one-electron equations (Kohn-

Sham equations). To simplify the notation, I will present the total energy functional and

Kohn-Sham equations for the spin-unpolarized case. They can be easily extended to the

spin-polarized 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)] (2-1)

Here, ELo,({R, })is the Coulomb energy associated with the interactions among the ions

at positions {R}) given by


En JZ (2-2)
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 ground-state energy of the electrons. The ground-state electronic

energy is given by

Ez,, = TE + E + E (2-3)

where Te is the kinetic energy, EeL is the electron-ion interaction energy, and E,, is the

electron-electron interaction energy. These terms are given by


T=If _- V2 ), r (2-4)
,ZfJy-( 2 (2-4)


E = JV,, (r)p (r)d3r (2-5)


E = p(r)p(')d rd r' + E [(r)] (2-6)
2 11 r-r'

p (r) = f ,(r)2 (2-7)


where Vo is the static total electron-ion potential, p (r) is the electron density, f is the

electron occupation number, and E.x[p(r)] is the exchange-correlation energy

functional.

The set of one-electron states y, that minimize the energy functional Eeec is given

by the solutions to the Kohn-Sham equations (Kohn and Sham, 1965):










V2 + o,"(r) + VH(r) + (r) (r) = s (r) (2-8)


where VH (r) is the Hartree potential of the electrons given by


VH P r (2-9)
V r r-r'

The exchange-correlation potential, V (r), is given by the functional derivative


V (r) (r) (2-10)
6p(r)

The Kohn-Sham equations represent a mapping of an interacting inhomogeneous

electron gas in the presence of an external potential V(r) onto a non-interacting

inhomogeneous electron gas moving in an effective potential

Vf (r) = V(r) + VH (r) + Vxc (r). If the exchange-correlation functional were known

exactly, the Kohn-Sham equations would give the exact ground-state p and Eiec .

The Kohn-Sham equations are nonlinear and must be solved self-consistently. Once

the eigenvalues s, and the wave functions y, are obtained, the total electronic energy

Ee ec can be determined from Eq. (2-3) to (2-7), or alternatively from the formula


c = c P(r)p(r )d3rd3r'+Ex [p(r)]- ixc(r)p(r) (2-11)
r2 r-r'

where N is the number of occupied Kohn-Sham orbitals. Here

N N
S(2-12)


The sum of the single-particle Kohn-Sham eigenvalues does not give the total energy,

because it overcounts the electron-electron interaction and the exchange-correlation

energy.









Local-Density Approximation

The Kohn-Sham equations [Eq. (2-8)] 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 local-density

approximation (LDA) proposed by Kohn and Sham (Kohn and Sham, 1965).

In LDA, the exchange-correlation functional can be written as

Ec = c(r) p (r)dr (2-13)

Thus,

C Ec [p (r)] = [p(r)sx(r)]
6p(r) p (r)

The local density approximation assumes that the exchange-correlation energy functional

is purely local. The exchange-correlation energy per electron s, is equal to the

exchange-correlation energy per electron in a homogeneous electron-gas that has the

same density as the electron gas at r. Thus,

a'(r)- f'A[p(r)] =s [p(r)] (2-15)


A commonly used parameterization for afLA is the Vosko-Wilk-Nusair parameterization

(Vosko et al., 1980, 1982) of the Ceperly-Alder 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 exchange-correlation

energy.

Considering the inexact nature of the local density approximation, it is remarkable

that calculations using LDA have been so successful. Generally, total-energy 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 exchange-correlation 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 total-energy density functional method is often implemented in a plane-wave

basis set, which is discussed later. Here, I introduce the pseudopotential approximation

that greatly reduces the computational cost of the plane-wave based DFT method.

A large number of plane-waves are needed to expand tightly bound core orbitals

and rapidly oscillating valence wave-functions 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 pseudo-wave-functions rather than the true valence

wavefunctions.

For a central field atom, the radial Kohn-Sham equation is given by

S 2+ t+ [p(r)] rRn,(r)=+ s[rR ,(r) (2-16)
2 dr2 2r2

The pseudo-wave-functions are obtained from an all-electron atomic calculation that self-

consistently solves Eq. (2-16). Once the pseudo-wave-function is obtained, the screened

pseudopotential is given by the inversion of Eq. (2-16),

PPcr / 1(1+ 1) 1 d2
--, [I +rRP (r)] (2-17)
2r2 2rRPP(r) dr2

where RPP(r) is the radial pseudo-wave-function. The ionic pseudopotential can be

obtained by subtracting the Hartree VPP(r) and exchange-correlation VP (r) potentials

calculated from the valence pseudo-wave-functions from the screened pseudopotential

PPsc (r).

Vr (r) = V PPs (- (r) VV4 (r) (2-18)

The pseodopotential needs to meet a set of requirements. The standard requirement

is "norm-conservation", which means that the pseudo-wave-functions are identical to the

real wave-functions outside a core region, and that they generate same electron density

inside the core region. A widely used norm-conserving 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









wave-function. The pseudopotential operator can be decomposed into local and nonlocal

parts,

V~p (r) = VPPc(r) + VfPPnc (r)P (2-19)


where VPPl'(r) is the local potential, and the semi-local potential JPPn"ic(r) is given by

SPP'nic(r) = VPP (r)- VPPloc (r) (2-20)

Here, P, is the operator that projects out the /th angular momentum component from the

wave-function. The nonlocal potential is short-ranged, 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. (2-19) 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. (2-20),

KC .f/ lc(r)(DP (r)) (PP (r)V "c(r)
V, Kc(r) =I (2-21)
(O ((r) ) OPP(/)c ( (r)


where PP'"f"(r) is given by Eq. (2-20), and OPP(r) is the atomic pseudo-wave-function

for the angular momentum 1. This form of the non-local potential substantially reduces

the number of integrals that need to be evaluated.

The norm-conserving condition requires that the total pseudocharge inside the core

region match that of the all-electron wavefunction. Thus for those highly localized

valence orbitals (e.g., 2p for first-row atoms and 3d for transition-metal atoms), a large









number of plane-waves are still needed. Vanderbilt relaxed the norm-conserving rule and

constructed ultrasoftt" pseudopotentials with a lower cutoff energy for plane-wave

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 first-row elements by a factor of 2 to 4

(Kresse and Furthmuller, 1996a).

The Hellmann-Feyman Theorem

Once the total energy of the system is obtained via Eq. (2-1), the force on ion I, f,

is given by

dE OE OE Oy OE 8Oy
df, C E C_ E (2-22)
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 (2-22) can be rewritten as


R,Hy, + Y R,
(RI / ORI/








f, (2-24)
aRI

This is the Hellmann-Feynman theorem (Hellmann, 1937; Feynman, 1939), which holds

for any derivative of the total energy.









Molecular Dynamics

The motion of ions on the Born-Oppenheimer potential energy surface is governed

by the forces obtained by Eq. (2-18) via integration of the Newtonian equations of motion


m, RI = fI (2-25)

The integration can be performed using the Verlet algorithm (Allen and Tildesley,

1987). In contrasted to empirical potentials, the forces determined by ab-initio 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 total-energy calculation.

Bloch's Theorem

A wave-function in a periodic system can be written in the following form

(Ashcroft and Mermin, 1976)

VY,k (r) = exp(ik.r)u,k (r) (2-26)

Here, u1,k (r) is a periodic function that satisfies the condition


uk (r) = uk (r + nR) (2-27)

where R is any lattice vector and n is an integer. Such a function can be expanded in a

discrete plane-wave basis set with wave vectors as reciprocal lattice vectors of the

crystal:

k (r)= cG exp(iG.r) (2-28)
G









Therefore the electronic wave-function can be written as

,,k (r) ,kG Iexp [i(k + G).r] (2-29)
G, k+G2 Ecu


The summation is over all reciprocal lattice vectors up to a cutoff energy Ek,. But in

principle Eut = o.

k-point 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 k-point. Therefore, the Bloch theorem transforms the problem with an

infinite number of electronic wave-functions to one with a finite number of electronic

wave-functions 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, k-point 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.

Momentum-Space Formalism

In this section, I show the derivation of the momentum-space 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,, (2-30)

where


Te =:ZfJd3' (r)f -V2 v(r) (2-31)


E = f d'r~ (r)V (r-R-ts)p,(r) (2-32)
i,R,s


EH = d3r p(r)VH(r) (2-33)


Ex = p (r)sc (r)d3r (2-34)


E =NY' zS-S (2-35)
o 2 R,s,s' It, t, + R


p(r) = f k,(r)12 (2-36)


The symbols / and y, (r) are, respectively, the electron occupation number and the

valence pseudo-wave-function for the one-electron 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. (2-9).

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. (2-34) 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 angular-momentum

dependent pseudopotential and is given by

PP (r)= ~PP(r = V(PPlc()c+ PPPnlc (r) (2-37)
1 /

The momemtum-space representation of each energy term per unit cell in Eq. (2-

30) is given by









T -- IC kk-- Ck+G 2(k+G)2 (2-38)
2 k k G

E p (G)
EH 47. (2-39)
2G)G

Exc = Q2 p*(G)Fxc (G) (2-40)
G

Ee = cZp (G) e- t rPPlc (G)
GO s
+ W .* c Yr-(G'-G)o-tsPPn (2-41
+ ""cZ k ,kek+G Lk+G'Zl L s,l,k+G, (2-41)
k,,l G,G'

+-QZl d'r VsPPic(r) s
s ir-ts


1n Ierfc( tt, -t,, + R) 2f
E, 2 Z tt, +R
(2-42)
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 (2-43)


P (r)= -ZZC kf, k Z Ck+GCk+G e(G' G)*r (2-44)
4 k i G,G'

VPPlc (G)= Jd3rVfPP (r)e-,G.r (2-45)


ec (G)=-' ex(r)e -G*rd3r (2-46)


V",. kG, (21+l)P(cosy)ji j(, (k+ GIr) (r) j,(k+G'r)rdr (2-47)









(k + G).(k + G')
cosy = (2-48)
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 Kleinman-Bylander type of nonlocal pseudopotential is used, Eq. (2-47)

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 (2-49)

x j^i (k+G'-r)(PPnc RIs 2dr


where RPs (r) is the radial pseudo-wave-function. 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 Kleinman-Bylander nonlocal pseudopotential.

The integrals in Eq. (2-49) 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

short-ranged and thus is only nonzero within a sphere around each ion, the number of

operations necessary to evaluate the integral in real-space does not increase with the

system size. Therefore, it is more efficient to evaluate the integrals in Eq. (2-49) in real

space for large systems. In order to do so, the wave-function for each angular momentum

must be represented in real space. However, the discrete real-space Fourier grid

introduces errors when projecting angular momentum components of the wave-functions









because of the non-uniqueness of the spherical harmonic projection in a discrete real-

space grid (Payne et al., 1992). To reduce these errors, King-Smith, Payne and Lin

(1991) proposed a procedure for improving the accuracy of the real-space calculation of

the nonlocal contribution to the total energy. Tests of this method showed good

agreement between real- and reciprocal-space calculations (King-Smith, Payne and Lin,

1991).

The Kohn-Sham equations in momentum-space representation can be derived

variationally from the expression for the total energy in Eq. (2-30),


ZI k+ G GG +k,G,G ,k+G =E,k C,k+G (2-50)
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. (2-50). (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. (2-41)] should be added to the

total energy. Once Eq. (2-50) 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) ) (2-52)
+ Q r Ve~o(r) -y









where EH, Ex and E,,, are given by Eq. (2-39), (2-40) and (2-42), respectively. This

expression is simpler to evaluate than Eq. (2-30), (2-38)-(2-42), since the double

summation over reciprocal lattice vectors in Eq. (2-41) is absorbed in the simple

summation of the eigenvalues of the occupied states in Eq. (2-52).

The Hellmann-Feynman force on each ion is given by

F Fs,+Fs, (2-53)

where
F, = -V'E5
= iQ p* (G) Ge-'G'rVPP'c (G) (2-54)
G#O
7i G-G )-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)] (2-55)
ss iC o GO 42)

Rt -t, )R r erfc( t(r tR t + R) 21-R 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.

Dual-Space Formalism

The procedure of evaluating the kinetic energy in momentum space and the

potential energy in real space can be thought of as a dual-space formalism for

pseudopotential calculations (Martin and Cohen, 1988; Troullier and Martin, 1991b;

Wentzcovitch and Martin, 1991) in contrast with momentum-space formalism (Ihm,

Zunger and Cohen, 1979) in which all the calculations are carried on in momentum









space. The dual-space 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 dual-space formalism for the pseudopotential calculations of finite

systems.

Density Functional Theory for Finite Systems: A Dual-Space Formulism

For the study of finite systems, the method described in Section (2-2) 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 plane-wave 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

Kohn-Sham equations are denoted xy,, where o is the spin index. The corresponding

density of spin o is given by









p, (r) = f (r)2 (2-56)


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. (2-30). Each energy

term in Eq. (2-30) is given by


=T f= J (Jy (-V2 jy) (2-57)


E ZZ = (Z J,(T S

EH = d3r p (r)VH(r) (2-59)


E, p (r);s (r)d3r (2-60)


Son = (2-61)
s' s' ts ts,

where VH is the Hartree potential, given by Eq. (2-9), VPP is the pseudopotential, which

can be decomposed into local and nonlocal parts,

JVPP = (r, r') = VPP" (r) + VJ," (r, r') (2-62)

Thus,
E1- = Ec + En (2-63)

where


ELc = d3r p(r)VPP"c (r t,) )


(2-64)









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) (2-65)


,j C s,l,m

1
F1, = (2-66)
fdr rRs, (r) V flC (r)

Ks,Rm (r)= V nc (r)R], (r) Ym (i) (2-67)

In Eq. (2-65), the Kleinman-Bylander non-local pseudopotential is used (Kleinman and

Bylander, 1982). The Kohn-Sham Hamiltonian is given by

1
H= V2 + +VH + VX (2-68)
2

where


V E a c(p)] (2-69)
6p ap

Vat (r,r') = Vic (r) + Vfc (r, r ')
= VPPc (r-t, )+ j ,tKs,, (r-ts)K,,m (r'-ts) (2-70)
s s l,m

and VH is the Hartree energy, given by Eq. (2-9).

Barnett and Landman (1993) suggested a dual-space 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


dual-space,









(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) (2-71)

= 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 momentum-space grids.

The size of the real-space cell is defined as follows

r (x,y,z); 0
The wave-function \ ,, (r) = (r jo) is expanded in a plane-wave basis within the cell,

i.e.,

1
y (r)= 4 (g)eg'r (2-73)
g9gmax

for r in the cell, and
vo (r) = 0 (2-74)

for r outside the cell.

In Eq. (2-73), Q is the volume of the cell, Q = LLyL and gmax is the cutoff

momentum in momentum space. The momentum-space grid (g-space) is defined by


g = 2 k,- k (2-75)
Lx Ly L,


are integers satisfying


where {k }a= ,y,










2 2

and

Nx N N ax (2-77)
Lx Ly Lz 7x

For a finite system, the wave-function can be chosen to be real in real-space, i.e.,

o v(r) = J (r) (2-78)

and thus,

J, (g) =j*T (-g) (2-79)

The density p, (r) can also be expanded in a Fourier series but requires a doubled

momentum-cutoff with respect to the wave-functions, i.e.,


p,,(r). 1 :I I^ (-0: (0 (g)g').r
g g (2-80)
1
SGD (G)eGr

The G-space grid in Eq. (2-80) is defined as


G=27r m my m (2-81)
[L L,)

where {m }a=,,y,z are integers satisfying


-M < m < and M = 2Na (2-82)
2 2

Now, a real-space S-grid can be defined as a dual of the G-space grid. The density

and the Kohn-Sham potential can be evaluated on this S-space grid, i.e.,










S= nx L nyL, nZ L
M 'My ',M
MMM A A
x y z/M


where n x,y
(. acL-x,v,


are integers satisfying


The wave-function can be evaluated on the S-grid as
The wave-function can be evaluated on the S-grid as


0 (S)=JM (S)


(2-84)


(2-85)


I (g)e'g.s
Jr g


where M = MXM The density on the S-grid is given by


D, (S)= If, (S)2 (2-86)


Now the different energy terms in Eq. (2-59), (2-60), (2-64) and (2-65) can be

directly evaluated on the S-grid as follows

1
EH (S)VH (S) (2-87)


Exc = p (S)Ec (S) (2-88)
S
s

Ec = p (S) I VPP,c (S t, |) (2-89)
S

2
Ecd = f Fs,, K, (S ts)yk, (S) (2-90)
j,0 s,i,m S

The kinetic energy is evaluated in momentum space. The Hamiltonian matrix element in


Eq. (2-71) is written as


(2-83)









1 2
(jo H ko) = 2 2 (-g) (g)


+ (-g) V, (S) + V, (S) + IP (S t) j (S)(2-91)


+ 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 -- (2-92)


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 first-order 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 bond-breaking and bond-forming. 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 length-scales 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

all-quantum 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 long-range

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 multi-scale method, which couples finite

element, molecular dynamics, and semi-empirical tight-binding methods, to study crack

propagation in silicon. However, the tight-binding method is over-simplified and not able

to describe bond breaking properly. Sauer and Sierka (2000) have developed a combined

quantum mechanics-interatomic 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 multi-scale simulation techniques, we have developed

a combined quantum-mechanical and classical-mechanical (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-

quantum-mechanical 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 2-1.

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 2-1. 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 ) (2-93)

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 Kohn-Sham equations self-consistently


V2 +V V,(r d' p(r' + 6E- ,(r) = s ,(r) (2-94)
L2 I r r' 6p(r)

p (r) => f ),(r)12 (2-95)

where V, (r) is the ionic pseudopotential, Ex is the exchange-correlation 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 }) (2-97)


EQmCM = V({Rn},{Rj}) (2-98)

where Z, is the core charge of the n-th ion. U({Rs}) and V({Rn},{Rs}) are described by

sums of classical interatomic potential functions in pairwise or three-body 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 pair-wise classical interatomic potential functions (, are used, U({Rs}) and


V({Rn}, {Rs}) are expressed as

U({R,})= ~-(SP,,) (2-99)
s s'>s

V({R },{R }) = I I ,D, (r, ) (2-100)
n s


One can also use three-body potential functions having the form

D = C f2 (i, j)+ jY f (i, j, k), (2-101)
1 J>1 i j>i k>j

where s is an energy parameter,f2 is the pairwise interaction between the i-th andj-th

atoms, andf3 is the three-body term as a function of the i,j, k-th atoms. In this case,

U({Rs}) and V({Rn}, {Rs}) are expressed as

U({RJ})= I-I F f:Ls23 +3(, ,,, ,,,) (2-102)
Ss'>s s s'>s s">s'

V({R },{R }J)= Y f2 (s)
(2-103)
+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) (2-104)









and

F, = -V,(EM + EQM M) (2-105)

In many other hybrid QM/CM methods, the QM-CM 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 (2-106)
s ir -R. n,s Rn -Rsl

where q, is the partial charge ofs-th classical atom, and Z, is the n-th 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 QM-CM and CM-CM 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 link-atom and the pseudo-atom methods.

The link-atom 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 link-atoms 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 H-termination for clusters of silica and zeolites. In our

method, the link-atoms 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,}) (2-107)

where L stands for the link-atom.

If the atoms A, e {QM} and B, e {CM}, we define the constraint such that the

link-atom 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 pseudo-atom method

The alternative to using hydrogen atoms to terminate the dangling bonds at the

boundary is to use pseudo-atoms to "absorb" these bonds. Taking the silicon dioxide

system as an example, a pseudo-oxygen-atom 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 pseudo-atom 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

pseudo-atom method, the dangling bonds are "absorbed" into the pseudo-atoms.

However, the orientation of the dangling bond is not taken into account if simply using a

real atom as the pseudo-atom. If one chooses to re-parameterize the pseudopotential, the

bond-angle dependence is still difficult to include in the pseudopotential. Thus, it is one

advantage of the link-atom method since such bond-angle dependence is naturally

included if the link-atoms 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

semiconductor-ceramic research.

Amorphous silica is abundant in the mineral and rock. The study of water-rock

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 water-rock 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 water-silica 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 chemo-mechanical 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, Si-O bonds are known to be inert to chemicals (H20, NH3,

etc) that can cause stress corrosion in silica glass. However, strained Si-O 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 (Si-O-H)

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 pre-treatment 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 two-membered ring can be identified by IR

peaks at 888 and 908 cm-1. 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 two-membered ring structure creates a local reactive site that is

vulnerable to chemisorption of gas-phase species.

The reaction between water and strained Si-O bonds is an important aspect of

surface hydroxylation. Reliable computational models for water-silica 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 co-workers (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 short-ranged. 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 co-workers (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 f-cristobalite surface to model a hydroxylated

amorphous surface (Ceresoli, et al., 2000; larori, 2001). This approach was adopted

because the density and refractive index of f-cristobalite 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 co-workers (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-, six-membered rings) but the f-cristobalite 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 long-range









order and hence long-range surface response in the fl-cristobalite 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 f-cristobalite 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 multi-scale schemes usually involve local

chemical processes in systems that do not have long-range order. Apparently, multi-scale

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 multi-scale 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 chemo-mechanical phenomena in amorphous silica, in

particular, stress-induced 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 link-atom and the pseudo-atom methods) on small silica clusters

(Du and Cheng, 2003a), and then apply the combined QM/CM method on the studies of

hydrolysis of the two-membered 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 Link-Atom Method

A small silica cluster (Si207H6, see Fig. 3-1) 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 3-1. 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. 3-2).

The link-atom (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 link-atoms 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 Beest-Kramer-Santen

(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


(3-1)


where r, is the distance between the vi andj-th 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 3-2. The structures of the cluster shown in Fig. 3-1, obtained from (a) all-quantum
calculation, (b) combined QM/CM calculation using the link-atom method,
and (c) all-classical calculation.









The optimized structures obtained from the all-quantum calculation, combined

QM/CM calculation, and all-classical calculation, are shown in Fig. 3-2. 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. 3-1, indicate that the QM part of the hybrid system has a

structure very similar to that in the all-quantum calculation, and the corresponding CM

part has a structure nearly identical to that in the all-classical 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. 3-3.

If including the electrostatic interaction between electrons in the QM region and

atoms in the CM region [see Eq. (2-106)], the partial charges [q, and q, in Eq. (3-1)]

assigned to CM atoms in Fig. 3-1 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 link-atoms can

terminate the dangling bonds very well and approximately restore the correct charge

distribution at the QM/CM boundary (see Fig. 3-3), the additional charge polarization

induced by CM point charges is not necessary. Thus, we did not use Eq. (2-106) to

include the electrostatic interaction between QM electrons and CM atoms.












S





(a)



H H(link) O




Si



(b)



Figure 3-3. The structures of the training system with iso-density surfaces, obtained from
(a) all-quantum calculation, and (b) combined QM/CM calculation.
Fig. 3-4 shows the total energy change due to the change of the Si(l)-O(b) or Si(2)-

O(b) bond length. When the Si-O 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 Si-O bond in the CM region is changed

this energy curve is close to that from the all-classical calculations. This is an indication

of a transparent quantum-classical interface.









A more rigorous test of the interface is shown in Figures 3-5 and 3-6. 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 all-quantum 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 all-classical

calculations.

Based on these tests (Figures 3-4, 3-5 and 3-6), 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. 2-99 and 2-100). On the other

hand, to the quantum atoms, which consists of ions and electrons, the classical atoms plus

the link-atoms 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 3-4-3-6

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. 3-2) is 1800. However, the

average Si-O-Si bond angle in amorphous silica is around 1460. Thus, it is necessary to

further test the validity of the link-atom 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 3-4. 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 long-dash red line, solid black line, and
short-dash 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 fully-relaxed cluster is shifted to zero.









0.3

S------QM
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 3-5. The force on (a) Si(1l) and (b) O(b), when the Si(l)-O(b) bond length varies.
The long-dash red line, solid black line, and short-dash 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 3-6. The force on (a) Si(1) and (b) O(b), when the Si(1)-O(b) bond length varies.
The long-dash red line, solid black line, and short-dash 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 3-7. 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 link-atom method is applied

to a three-membered silica ring (Si309H6, see Fig. 3-7). 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. 3-8), i.e., the structure

parameters, generated by the combined QM/CM method, are in good agreement with the

results from the all-quantum calculations in the QM region, as well as with the results

from the all-classical calculations in the CM region.

The tests in small silica clusters demonstrate that the link-atom 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 Si-O 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 3-8. The structures of the cluster shown in Fig. 3-7, obtained from (a) all-quantum
calculation, (b) combined QM/CM calculation using the link-atom method,
and (c) all-classical calculation (bottom panel). In (a) and (c), the lengths of
the Si-O 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 Pseudo-Atom Method

As already discussed in this method, the dangling bond at the QM/CM boundary is

"absorbed" into a pseudo-atom 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 link-atom method.

The optimized structures of these two clusters, obtained from the combined

QM/CM method, are shown in Fig. 3-9. For Si207H6, the interface works very well, as

compared with Fig. 3-2. The structure in the QM region is in good agreement with that

from all-quantum calculations, and the structure in the CM region is very similar to that

from all-classical calculations. Thus, the pseudo-oxygen-atom properly terminates the

dangling bond on the boundary. However, a discrepancy appears for the Si309H6 cluster,

as seen by comparison of Fig. 3-9 with Fig. 3-8. The difference in the optimized

structures is mainly from the bond angle. For example, the Si-O*-Si bond angle obtained

from the combined QM/CM method is 149.1 (see the bottom panel of Fig. 3-9), which is

much bigger than the corresponding bond angle (134.00) from the all-quantum

calculation (see the upper panel of Fig. 3-8). Such disagreement arises because no bond-

angle dependence is built into the current pseudo-atom method. Thus, the electrons near

the QM/CM boundary do not feel any difference between structures with different Si-O-

Si bond angle. This explains why such an interface works well for the Si207H6 cluster, in

which the Si-O-Si bond angle is 1800 at the QM/CM boundary, but fails for the Si309H6

cluster. Varying the nuclear charges of O* atoms improves the Si-O* bond lengths, but

fails to improve the Si-O*-Si bond angles for the same reason. Thus, a Si-O*-Si bond









angle dependent O* effective core potential must be designed if the pseudo-atom method

is to be used.

The construction of a transferable O* effective potential in a dynamically evolving

SiO2 system is rather complicated. Fortunately, such bond-angle dependence is naturally

included in the link-atom method if the link-atoms are constrained on the lines of the Si-

O bond on the QM/CM boundary. Based on the preceding considerations, we have used

the link-atom 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 3-9. Cluster structures obtained from combined QM/CM calculations using the
pseudo-atom method. (a) Cluster shown in Fig. 3-1. (b) Cluster shown in Fig.
3-7.


(a)


(b)









Hydrolysis of a Two-Membered Silica Ring on the Amorphous Silica Surface

Generation of the Amorphous Silica Surface

The amorphous silica sample was prepared by cooling a preheated fl-cristobalite

sample via a classical MD simulation that uses the BKS potential (Eq. 3-1) with a

correction to suppress an unphysical attraction at short interatomic distance. A bulk

sample of/f-cristobalite, 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 2-VIII).

Fig. 3-10 shows the MD results for the 0-0, Si-O and Si-Si pair-distribution functions in

the amorphous silica bulk. The positions of the first peak in gsi-o (r), go-o (r) and gsi-si (r)

give the mean Si-O bond length, nearest-neighbor 0-0 and Si-Si 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

Si-O-Si angle and average O-Si-O 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 z-direction 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 five-membered rings are the dominant species. There are also some larger

rings and some smaller ones such as four-, three-, and two-membered 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 one-coordinated O atoms hanging

on the surface. These defect sites (dangling bonds and two-membered rings) can be

annihilated easily by water attack.

Fig. 3-12 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 higher-density 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 high-density region is due to the

enrichment of small silica rings (two-, three- and four-membered rings) on the surface.

The Si-O bond lengths in small rings are usually longer than those in bigger rings

because of the enhanced strain in small rings. Thus, the high-density region on the

surface does not correspond to shorter bond lengths in the same region, as shown in Fig.

3-13(a). Instead, the average Si-O bond length in the high-density 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. 3-13(a)] is due to those

Si-O bonds with one-coordinated O atoms. The enrichment of small rings on the surface

also is evidenced by the decreasing average Si-O-Si bond angle on the surface, as shown

in Fig. 3-13(b). The sharp decline at the top of the surface in Fig. 3-13(b) is due to those

two-membered rings, which have smallest Si-O-Si bond angles. The concentration of

two-membered rings on the surface is found to be 0.42/nm2, which is close to the

experimental estimate (0.2-0.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 3-10. Partial pair-distribution functions for amorphous silica at 300 K.





















































Figure 3-11. 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 3-12. 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 3-13. (a) The average Si-O bond length, and (b) the average Si-O-Si 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 Troullier-Martins pseudopotential (Troullier and

Martins, 1993) and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation 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 link-atom approach was used to saturate the valence structure in the QM region

(i.e. a hydrogen atom is constrained on the line of a Si-O 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 link-atom

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. 3-1) 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 Stillinger-Weber type potential for silica (Watanabe, et

al., 1999), in which parameters are fitted to reproduce the cohesive energies of Si-O

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) (3-2)
I J>1 i j>i k>j

where s (=50 kcal/mol) is an energy parameter,f2 is the pairwise interaction between i-th

andj-th atoms, andf3 is the three-body term as a function of the i,j, k-th 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,) (3-3)
s s'>s

and the interaction energy between QM and CM atoms, V({Rn},{Rs}), was given by

V({R, ,{RJ)= II ,(r.) (3-4)
n s

The forces acting on QM and CM atoms follow from by Eq. 2-104 and 2-105.

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 I-I Ff, (r,, r,,,r,,,,,) (3-5)
s s>s s s s'>s s">s'


V({R},{R,})= f2 fZ(rs)
(3-6)

n s s'>S n n'>n s









Simulation Set-up

A 12000-atom amorphous silica slab, consisting of four 3000-atom 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. 2-104), while the

positions of all the CM atoms were fixed to simplify the problem. [The division of the

system is depicted in Fig. 3-14(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) (3-7)

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. 3-14(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 3-14. 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 Two-Membered Ring

The calculated adsorption configuration of H20 on the amorphous silica surface is

depicted in Fig. 3-15 (Only the 31 QM atoms on the surface are shown). In the upper left

picture in Fig. 3-15, 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 3-1. The average Si-O 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 Si-O network is weakened, which is evidence

of the hydrolytic weakening.

Table 3-1. Several Si-O 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. 3-16(a) and (b).

One is Si206H4 (cluster A), which is a two-membered ring. The other one is a 31-atom

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 Si-O

distances are compared with those from the combined QM/CM surface model in Table

3-1.


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 3-15. 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. 3-17, 3-18 and 3-19. All density contours

are drawn on the plane crossing the Si(1), 0(1) and 0(2) atoms. As shown in these

figures, the water-surface 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 3-16. 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 3-17: 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 .OE-3


2 R-0.002
02 .OE-3




















Sil
-1 0 1 2
x (A)

(a)




3-
01









-1 0 1 2 3
X(A)

(b)

Figure 3-18. Contour plot of (a) excess electron density An' and (b) deficiency electron
density An in the water-surface 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 3-19. Contour plot of (a) excess electron density An' and (b) deficiency electron
density An in the water-cluster 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 3-2. Following the procedure used

above for analyzing both the local surface energy increase and the surface-water bonding

strength, the cluster energy increases by 0.14 eV and the cluster-water 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. 3-15. According to experimental observations, the two-membered 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 3-2. 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 31-atom 31-atom cluster 31-atom 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 3-3. 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 cut-off 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. 3-7. 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 3-3. 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 3-3: 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 one-water 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 Si-O bonds

and remaining hydroxyls. Our result is in reasonable agreement with these experimental

data.

Adsorption of Two Water Molecules on the Two-Membered 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. 3-20 and the reaction energy is shown in Table 3-3. 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 two-membered ring shown in the upper-

middle panel of the Fig. 3-20 is +0.7 obtained using Mulliken charge analysis. The final

configuration (bottom-right panel in Fig. 3-20) consists of one H20 molecule with two

silanols on the surface. The same reaction pattern has also been observed on the (0001)









surface of a-A1203 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. 3-16c and the adsorption energy is shown in Table 3-1.)

The water-dimer cooperative reaction may be also energetically favored for other

strained Si-O 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 Si-O 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

hydrogen-bonding 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 two-membered ring evolves into two silanols. Without further applying strain, this

silanol-covered 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. 3-21. 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 3-20. 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 two-membered 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 (Si-H and Si-OH) can be energetically more favored

than on the clean surface. We expect that similar phenomena can happen on a

dehydroxylated silica surface.


Figure 3-21. 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 fl-cristobalite model. For

amorphous silica, charge transfer is local and the electronic response is short-ranged. A

simple link-atom 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 Si-O 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 semi-empirical potentials









for silica to ensure the sufficient transfer of information across the QM/CM interface in a

chemically realistic multi-scale simulation (Al-Derzi, 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., Al-Derzi, 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 a-quartz 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 large-scale dynamical

problems, such as crack propagation, which need hundreds of atoms in the QM region.

Previous multi-scale simulations that permitted hundreds of atoms in the QM region were

based on the tight-binding method (Broughton, et al., 1999), which is over-simplified and

not able to describe bond breaking properly. A new attempt has been made to re-

parameterize the semi-empirical Hamiltonian based on the input from high-level quantum

mechanical method, coupled-cluster (CC) theory (Talor, et al., 2003). The new QM

method (so-called "transfer Hamiltonian") is able to treat 500-1000 atoms self-

consistently and retains the accuracy of CC method.

In terms of future studies of the chemo-mechanical 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 two-membered ring on the amorphous silica surface. With

this model, we have described the long-range 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 two-water 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 long-term 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 single-wall 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 4-1. 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 Single-Wall Carbon Nanotube

A single-wall carbon nanotube (SWCNT) can be described as a graphene sheet

rolled into a cylindrical shape so that the structure is one-dimensional 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) (4-1)

where a, and a2 are the real space unit vectors of a graphene sheet, and n and m are

integers (Fig. 4-1). A nanotube can be constructed by rolling a graphene sheet into a

cylinder. Fig. 4-1 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 (4-2)


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) (4-3)

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 = (4-4)
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 alkaline-earth intercalated fullerenes have also been studied

extensively because of their complex phase diagram and superconductivity at

temperatures surpassed only by the high-To 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 single-wall carbon nanotubes

(Smith etal., 1998; Burteaux et al., 1999; Smith etal., 1999; Smith et al., 2000) and

achieved high-yield 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 band-gap in (Gd@C82)n@SWNT (Lee et al., 2002). A

temperature-induced change from p to n conduction in (Dy@C82)n@SWNT (Chiu et al.,

2001) and ambipolar field-effect 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 band-gap

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 potassium-doped peapods.

Method and Simulation Setup

The electronic structure calculations were performed using density functional

theory (DFT) with the local density approximation (LDA). The electron-ion interactions

were described by ultra-soft pesudopotentials (Vanderbilt, 1990). The valence wave









functions were expanded in a plane-wave 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 ab-initio 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. 4-2

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 = (4-5)
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 4-1. 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 4-2. The structures of (a) (10, 10), (b) (17, 0)a, and (c) (17, 0)b nanopeapods.
Table 4-1. 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 Fullerene-Induced Impurity States in M@C6o@SWNTs
In an isolated C60 molecule, there is a three-fold degenerate lowest unoccupied tlu

state. Fig. 4-3 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 4-3(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 tiu-derived band

because of the long center-to-center 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 tiu-derived states behave as impurity states [see

Fig. 4-3(b)] and carriers are distributed on the nanotube wall. When d, decreases, these

impurity states become impurity bands [see Fig. 4-3(c)] and conductive, resulting in most

hole-carriers distributed on the nanotube wall and most electron-carriers 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 4-3. 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 tiu-derived
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 first-Brillouin 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. 4-3(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 ill-1 .....0 -1.0 '
F X F X r X

Figure 4-4. 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 tiu-derived states. Our calculations show that the energy levels of tiu-derived

states are the inside conduction band in C60@(16, 0), about 2/3 of the band-gap (0.6 eV)

above the top of the valence band (E,) in C60@(17, 0), and about 1/3 of the band-gap

(0.48 eV) above E, in C60@(19, 0) (see Fig. 4-4). 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 tiu-derived states are 2/5 of the band-gap (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 tiu-derived states shift down toward the top of the valence