<%BANNER%>

Multi-Scale Modeling of Solids as a Composite of Quantum Mechanical (QM) and Classical Mechanical (CM) Domains


PAGE 1

MULTI-SCALE MODELING OF SOLIDS AS A COMPOSITE OF QUANTUM MECHANICAL (QM) AND CLASSICAL MECHANICAL (CM) DOMAINS By ADITI MALLIK A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005

PAGE 2

Copyright 2005 by Aditi Mallik

PAGE 3

This dissertation is dedicated to my parents.

PAGE 4

ACKNOWLEDGMENTS This dissertation owes its existence to my supervisory committee chair Professor James W. Dufty of University of Florida from whom I have learnt to perceive a research problem in proper context, translate rudimentary ideas into formal analysis, and check the logical consistency of the presentation. I cannot adequately thank him for his attention, emotional support, and inspiration that I have received ever since I joined the University of Florida, Physics Department in 2000. I would also like to thank my other supervisory committee members, Professors Samuel B. Trickey, Hai-Ping Cheng, John Mecholsky, Chris J. Stanton, and Keith Runge. Their valuable time and suggestions helped me to complete the work within a reasonable time frame. I am deeply indebted to Dr. Keith Runge for his unstinted interest in my research progress and especially his remarkable patience in answering all my nave questions on Quantum Chemistry. Thanks go to Dr. DeCarlos Taylor for helping me with the computer programming. The implication of this support will only be appreciated by those who frequently got lost in the programming-maze as I did. I am also grateful to Professor Hai-Ping Cheng for teaching me the parallel programming code on Density Functional Theory. Dr. Krishna Muralidharan is gratefully acknowledged for collaboration in his work on glass which developed into a part of a chapter in this dissertation. The Physics Department has been a wonderful place for my learning experience. It was a great experience to be a part of the Quantum Theory Project (QTP) group with people who have been good friends, sources of discussion and technical support. I wish to iv

PAGE 5

thank all former and current members of the QTP group especially Dr. Ling-Ling Wang, Dr. Moahua Du, Chao Cao, and Ying-Xao Wang for their help and mentorship. I am also grateful to my friend Aparna Baskaran who offered pertinent comments on a few chapters in this dissertation. I am deeply indebted to those who sustained me at a personal level. My heart-felt appreciation goes to my parents and all my family members back home in India. Without their constant emotional support, I simply could not have come this far. I wish to thank my husband, Samir, for making my stay in the US so enjoyable; not to mention the painstaking job of proof-reading this dissertation. This work was support by ITR grant 0325553. v

PAGE 6

TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iv LIST OF TABLES .............................................................................................................ix LIST OF FIGURES ...........................................................................................................xi GLOSSARY ....................................................................................................................xiv CHAPTER 1 INTRODUCTION........................................................................................................1 Definition of Multi-Scale Modeling.............................................................................1 Examples of Multi-Scale Modeling..............................................................................2 Methods for Different Length Scales...........................................................................5 Quantum Mechanical (QM)/Classical Mechanical (CM) Simulations........................6 Proposed Scheme..........................................................................................................7 Model System...............................................................................................................8 2 THE MATERIAL AND ITS FORMAL PARTITIONING.......................................11 Introduction.................................................................................................................11 The Idealized Quantum Solid.....................................................................................12 The Formal Partition and Composite Solid................................................................16 The Representative Classical Solid.....................................................................16 The Composite Quantum/Classical Solid............................................................18 Description of the Classical Subsystem......................................................................21 Description of the Quantum Subsystem.....................................................................21 Coulomb Effects of Environment........................................................................24 Electron Exchange with Environment.................................................................25 3 DESCRIPTION OF THE QUANTUM DOMAIN.....................................................28 Introduction.................................................................................................................28 Quantum Methods......................................................................................................29 Description of the QM Domain..................................................................................32 Different Termination Schemes..........................................................................32 Pseudo-atom Method...........................................................................................35 vi

PAGE 7

Calculation of Charge Densities..........................................................................37 Embedding Scheme.............................................................................................40 Other Options for Links.......................................................................................44 Conclusions.................................................................................................................47 4 CONSTRUCTION OF A NEW CLASSICAL POTENTIAL FOR THE CM DOMAIN....................................................................................................................49 Introduction.................................................................................................................49 Molecular Dynamics...................................................................................................50 Corrector-Predictor Algorithms...........................................................................50 Velocity Rescaling...............................................................................................52 Failure of Standard Potentials.....................................................................................52 Objective and Method for Constructing a New Potential...........................................54 Results.........................................................................................................................60 Finding Parameters for the Potential...................................................................60 Stress-Strain Curves............................................................................................64 Improving Stress-Strain Curves for Higher Strains.............................................67 Conclusions.................................................................................................................70 5 CONSTRUCTION OF THE COMPOSITE NANOROD..........................................72 Introduction.................................................................................................................72 Results for the Composite Rod...................................................................................73 Conclusions.................................................................................................................74 6 INDEPENDENCE OF MULTI-SCALE MODELING WITH RESPECT TO THE CHOICE OF UNDERLYING QM.............................................................................75 Introduction.................................................................................................................75 Review of DFT Method..............................................................................................76 Born-Oppenheimer Approximation.....................................................................76 Hohenberg-Kohn Theorems................................................................................76 Kohn-Sham Equations.........................................................................................77 Local Density Approximation.............................................................................78 Generalized Gradient Approximation.................................................................79 Troullier-Martin Pseudo-Potential.......................................................................80 Dual Space Formalism........................................................................................82 Computational Details................................................................................................82 Results and Discussion...............................................................................................83 Training of DFT-Pseudo-atoms...........................................................................83 Construction of Classical-DFT Potential for the CM Region.............................86 Conclusions.................................................................................................................89 vii

PAGE 8

7 APPLICATIONS OF THE MULTI-SCALE MODELING PROCEDURE TO OTHER SILICA SYSTEMS......................................................................................90 Introduction.................................................................................................................90 Systems Studied..........................................................................................................91 Silica Nanorings..................................................................................................91 Silica Nanoclusters..............................................................................................94 Notched Nanorod.................................................................................................95 Bulk Glass...........................................................................................................97 Results Obtained from New-DFT-Classical-Potential.............................................104 The 10-Membered Silica Nanoring...................................................................104 Bulk-Glass.........................................................................................................105 Conclusions...............................................................................................................106 8 COMPARISON OF THE DIFFERENT QUANTUM MECHANICAL METHODS...............................................................................................................107 Introduction...............................................................................................................107 The Nanorod.............................................................................................................108 Comparison of Structure....................................................................................108 Elastic Properties...............................................................................................109 Comparison of Different Types of DFT...................................................................110 Bulk Glass.................................................................................................................111 Discussion.................................................................................................................112 9 SUMMARY AND FUTURE RECOMMENDATIONS..........................................115 Conclusions...............................................................................................................115 Recommendations for Future Work.........................................................................118 APPENDIX A NEGLECT OF DIATOMIC DIFFERENTIAL OVERLAP (NDDO) MODEL......120 B GENETIC ALGORITHM........................................................................................124 C SCALING METHOD...............................................................................................128 D CHARGE TRANSFER POTENTIAL.....................................................................130 Charge-Transfer Potential for Old Transfer Hamiltonian Classical Potential #1.....131 Charge-Transfer Classical Potential for New-Density Functional Theory...............133 LIST OF REFERENCES.................................................................................................134 BIOGRAPHICAL SKETCH...........................................................................................140 viii

PAGE 9

LIST OF TABLES Table page 3-1. Comparison of charge densities of the ring with LA and ring with pseudo-atoms with that of the ring in the bulk calculated at planes parallel to the plane of the ring...........................................................................................................................39 3-2. Percentage charge difference with respect to bulk calculated with our method for various cases studied................................................................................................43 4-1. Different values of potential parameters obtained from GA which gave a stable rod.............................................................................................................................61 4-2. Parameters for New Potential in comparison to TTAM and BKS potential parameters................................................................................................................62 4-3. Comparison of structure of the rod obtained from the different potentials...............63 4-4. New parameters for classical potential-2 (NTH-2 Classical Potential).....................68 4-5. Structure obtained from NTH-2 Classical Potential in comparison to that of the TH nanorod..............................................................................................................69 6-1. Magnitude of force on the terminated Si atom in pyrosilicic acid with changing parameters of the dft-pseudo-atom...........................................................................84 6-2. Percentage charge density difference obtained with our scheme versus the charge density in the actual molecule for different cases studied in the DFT method........85 6-3. Parameters for New-Classical-DFT Potential in comparison to standard TTAM and BKS potential parameters..................................................................................87 6-4. Structure of the nanorod obtained using New-Classical-DFT, BKS potentials versus structure of DFT method...............................................................................87 7-1. Comparison of structure of the nanoring obtained from TH, NTH-2 Classical Potential and TTAM potential..................................................................................92 7-2. Structure comparison of 2-membered and 3-membered silica rings.........................94 7-3. Comparison of different properties of the glass obtained from the NTH-2 Classical Potential to those of the BKS glass.........................................................101 ix

PAGE 10

7-4. Comparison of structure of the nanoring obtained from DFT and New-DFTClassical Potential..................................................................................................105 7-6. Comparison of different properties of the glass obtained from New-DFT Classical Potential to those of BKS glass...............................................................105 8-1. Comparison of structure of the nanorod obtained from DFT and TH quantum methods..................................................................................................................108 8-2. Comparison of the 10-membered nanoring using two different types of DFT methods..................................................................................................................110 8-3. Comparison of structure of bulk glass using two different classical potentials......111 x

PAGE 11

LIST OF FIGURES Figure page 1-1.Use of different methods at different length scales in crack propagation....................2 1-2. Chemo-mechanical polishing of silica wafers.............................................................3 1-3. Structure of the nanorod. A) Front view. B) Top view. The red spheres are the Si atoms and the blue spheres are O atoms.....................................................................9 1-4. Partitioning of the nanorod into QM and CM domains. Localized (blue) electron cloud of the CM region shows the appropriateness of such partitioning.................10 3-1. Training of Pseudo-atom on smaller molecules, here Pyrosilicic acid......................36 3-2. Charge densities with different termination schemes. (A) Isolated ring (B) with links (C) with pseudo-atoms (D) in actual nanorod.................................................38 3-3. Comparison of magnitude of force on Si nuclei with different termination schemes....................................................................................................................40 3-4. Embedding Scheme: Approximating the rest (CM region) of the rod beyond the pseudo-atoms by dipoles..........................................................................................41 3-5. Effect of inclusion of the dipoles in the TH. (A) Improvement in the value of forces. (B) Reduction in the value of normalized difference of charge densities Open symbols, without dipoles; squares, with dipoles............................................41 3-6. Cross-sectional view. (A) Distorted ring. (B) Ring in equilibrium...........................42 3-7. Comparison of forces on Si nuclei with our scheme and LA to those in the actual molecule for various cases studied...........................................................................43 3-8. Partitioning of 3-membered silica ring......................................................................44 3-9. Deviated position of optimized LA...........................................................................45 3-10. Charge density with LA at optimized positions......................................................45 3-11. Position of LA beyond intermediate oxygen atoms................................................46 xi

PAGE 12

3-12. Comparison of force on Si nuclei with LA on oxygen to LA on silicon and our embedding scheme...................................................................................................47 4-1. Comparison of stress-strain curves for the standard pair potentials (BKS, TTAM) with that of Transfer Hamiltonian............................................................................54 4-2. Presence of local and global minima.........................................................................59 4-3. Comparison of Si-O pairwise interactions for the three potentials...........................63 4-4. Comparison of Youngs modulus at low strains for different classical potentials with that of TH.........................................................................................................64 4-5. Reduction of the wiggles with temperature...............................................................65 4-6. Stress-strain curves at lower strain rate of 2m/s and a temperature of 0.2K.............65 4-7. Histogram of number of modes vs. wave number (cm -1 )..........................................67 4-8. Behavior of NTH-2 Classical Potential near fracture................................................69 5-1. Comparison of stress-strain of the composite rod with those obtained from NTH-2 Classical Potential and TH....................................................................................74 6-1. Comparison of force on the Si nucleus with different termination schemes and dipoles for the DFT case..........................................................................................85 6-2. Comparison of adiabatic stress curves for DFT, new-DFT-classical and TTAM potential....................................................................................................................88 7-1. The 10-membered silica nanoring.............................................................................91 7-2. Stress-strain curve for the nanoring obtained from TH and NTH-2 Classical Potential....................................................................................................................93 7-3. The 2and 3-membered silica rings..........................................................................94 7-4. Silica Nanoclusters. A) Si 12 O 24 B) Si 18 O 36 C) Si 24 O 48 -Cage. D) Si 24 O 48 -Fullerene...................................................................................................................95 7-5. Front view of the 107 atom notched nanorod showing one defect (missing oxygen).....................................................................................................................96 7-6. Comparison of stress-strain curves for the notched rod obtained from different methods, showing the significant reduction in the yield stress with the presence of a single defect......................................................................................................97 7-7. Structure of glass obtained from the NTH-2 Classical Potential...............................98 xii

PAGE 13

78. Radial distribution function of the Si-O, O-O and Si-Si pairs in the glass.............100 7-9. Bond angle distribution of O-Si-O and Si-O-Si......................................................100 7-10. Stress-strain curves for NTH-2 and BKS glass.....................................................102 8-1. Comparison of elastic properties for DFT and TH nanorods..................................109 8-2. Comparison of elastic properties of the two glasses in the elastic region...............112 D-1. Charge transferred from Si to O as a function of distance.....................................131 D-2. Increase in the stiffness with introduction of the charge transfer term...................132 D-3. Increase in the charge of Si with strain...................................................................132 D-4. Increase in the stress with the charge transfer DFT potential.................................133 xiii

PAGE 14

GLOSSARY BAD Bond angle distribution CC Coupled Cluster CM Classical mechanics / mechanical DFT Density Functional Theory FE Finite Element GA Genetic Algorithm LA Link Atoms MD Molecular dynamics NDDO Neglect of Diatomic Differential Overlap NTH-2 New Transfer Hamiltonian Classical Potential # 2 QM Quantum mechanics / mechanical RDF Radial distribution function TB Tight Binding TH Transfer Hamiltonian TM Troullier-Martin Y Youngs modulus xiv

PAGE 15

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 MULTI-SCALE MODELING OF SOLIDS AS A COMPOSITE OF QUANTUM MECHANICAL (QM) AND CLASSICAL MECHANICAL (CM) DOMAINS By Aditi Mallik May 2005 Chair: James W. Dufty Major Department: Physics The simulation of non-equilibrium phenomena such as fracture in a macroscopic sample using ab initio quantum methods is not feasible because of the high computational time required. The alternative approach of a classical mechanical method does not have this limitation, but an accurate description cannot be made for important states far from equilibrium. A promising tool for studying such phenomena is the application of quantum techniques only to the small reactive regions and classical mechanics to describe the near-equilibrium bulk. Such methods are referred to here as quantum mechanical/classical mechanical (QM/CM) methods in multi-scale modeling. Our study addressed the problem of how to bridge these QM and CM domains by constructing an approximate semi-classical representation of the quantum (actual) solid as an equivalent composite system comprising a small QM domain plus the CM bulk region. Considering each domain as an open system, it is necessary to incorporate the relevant information about the state of the complimentary domain in a simple and xv

PAGE 16

accurate manner. For the QM domain, the CM environment is reduced in terms of trained pseudo-atoms and dipoles to account for the short-range chemical bonding and long-range Coulombic interactions, respectively; for the CM region a new classical pair potential has been constructed using quantum force data, such that this potential yields the same equilibrium structure and elastic properties as the QM region. This approach to multi-scale modeling has been applied to a SiO 2 nanorod and SiO 2 nanoring for critical testing. It is shown that the modeling reproduces the forces, charge densities, and elastic properties to within a few percent of the actual quantum values. Applications to hydrolytic weakening and fracture in bulk glass are discussed. Although our study was limited to silica-based materials, it should provide guidance for extension to other materials as well. xvi

PAGE 17

CHAPTER 1 INTRODUCTION Definition of Multi-Scale Modeling The field of multi-scale modeling has opened a new era to computational science for studying complex phenomena such as fracture, hydrolysis, enzymatic reactions, solute-solvent studies, hydrogen embrittlement, and many other chemo-mechanical processes in macroscopic samples. Often these phenomena require a very accurate description at one scale, while at another scale, one can apply a much less sophisticated method to get satisfactory results [Broughton et al., 1999]. In fact, it is necessary to have the less sophisticated description for the bulk sample because the higher accuracy methods are computationally too intensive to be applied to the bulk. These scales may be length or time or a combination of both. This scheme of combining different models at different scales to achieve a balance of accuracy, efficiency and realistic description is known as multi-scale modeling. It is accomplished by applying the very high accuracy method only in a small domain where it is needed the most and more approximate methods for the rest of the bulk where they are appropriate. The next example shows the need to combine different models. Crack propagation : For the case of crack propagation (Figure 1-1), one has to apply detailed quantum mechanics at the tip of the crack where the bonds are breaking between the atoms, causing a marked deformation of the electron cloud and charge transfer among ions. Classical mechanical methods cannot be applied in this region because classical Hamiltonians forgo any electronic treatment. But far from the crack tip, where the atoms 1

PAGE 18

2 are deformed less, atom interactions can be described well by classical potentials. It would be computationally infeasible to apply the detailed quantum mechanical methods to the whole sample; as the application of these methods are limited to only very small molecules (except for crystals). Thus, one must link different models for the different length scales; as no single theory is sufficient to describe the process of crack propagation in the bulk sample. Broughton et al., [1999] described multi-scale simulations in detail. Deformation in electron density, QM required Figure 1-1.Use of different methods at different length scales in crack propagation [Broughton et al., 1999]. Examples of Multi-Scale Modeling (i) Chemo-mechanical polishing (CMP): combines chemical and mechanical interactions to planarize silica wafer surfaces, using a slurry composed of solvents, wetting agents, other compounds and submicron-sized particles [Singh and Bajaj, 2002]. The chemical compounds interact with the material to form a chemically modified surface and the abrasives in the slurry mechanically interact with the chemically modified surface layers, resulting in material removal. Here interactions occur at both micro and nano-scale, so these two scales have to be treated generally with two different methods. At the microscale, the pad carrying the abrasive particles interacts with the surface of the

PAGE 19

3 wafer (Figure 2-2). At the nanoscale, the kinetics of the formation and removal of the thin surface layer controls CMP output parameters (such as removal rate, surface planarity, density of surface defects etc.). slurryChemically passivating layer Abrasive Figure 1-2. Chemo-mechanical polishing of silica wafers [Singh and Bajaj, 2002]. (ii) Tidal wave prediction : Clementi [1988] studied the interactions of several water molecules on an angstrom scale using quantum calculations. From this database, he found an empirical potential for molecular dynamics simulation, which was then used to determine the viscosity and density of water. Finally they used the density and viscosity of water in computational fluid dynamics to model tidal circulation in Buzzards Bay of Massachusetts. (iii) Atmospheric science: In this type of multi-scale modeling one has to combine different length as well as time scales. High accuracy computational methods are used to determine reaction barriers for simple chemical reactions which are then used in large rate equations coupled to grid computation to predict chemical meteorology.

PAGE 20

4 (iv) Embrittlement of pressure-vessels : Another multi-scale process in which the phenomenon occurs over different time scales is embrittlement of pressure-vessel steels of nuclear reactors caused by prolonged exposure to nuclear radiation [Odette et al., 2001]. In the preceeding examples of tidal waves and atmospheric science, the various scales are weakly coupled i.e., the behavior at one scale does not influence the phenomena at the others strongly. Rudd and Broughton [2000] described such models as serial multi-scale-modeling which means the computation of parameters at smaller scale for its use in more phenomenological models at a larger scale. The other type of multi-scale problem, in which the scales are strongly coupled is known as concurrent multi-scale modeling. In this type, all the different models applied on different scales must be nested with proper boundary conditions. For example in brittle fracture, the fracture takes place on a macroscopic scale with a propagation direction determined by local stresses in the material. At the atomic scale, the chemical bonds are taken far from the equilibrium. So the overall fracture mechanism depends on the instantaneous atomic bond lengths, at the time of strain application, yet the crack changes character as it propagates over macroscopic distances. Similarly, the CMP process can be conceptualized as the altered reactivity of a strained bond. This alteration occurs in combination with effects of local surface, neighboring atoms, and strain distribution. In the presence of a chemically reactive species, a stress that normally would produce a local elastic strain instead can lead to bond breaking. Our multi-scale modeling is more of the concurrent type, though some of the work can be considered as the serial type. Chapter 3 shows that it is necessary to

PAGE 21

5 include the information about other scales while studying the phenomena at one particular scale as in the concurrent method. However the parameterization of the quantum method (Transfer Hamiltonian) used is initially done on smaller molecules. The parameterization then is transferred to bigger molecules as in the serial method. Such distinctions are discussed more in chapter 3. Methods for Different Length Scales One can distinguish three different length scale levels in the context of multi-scale modeling of materials (1) The nanoscale: Here electrons play an important role so the phenomena are described using quantum mechanical methods like Density Functional theory (DFT) [Parr and Yang, 1989], post-Hartree-Fock methods [Stewart, 1990] and quantum Monte Carlo (QMC) [Foulkes et al., 2001]. (2) The atomic or microscale: Here phenomena are described in terms of classical interactions among atoms via Molecular Dynamics (MD) [Haile, 1992] or Monte Carlo (MC) simulations. (3) The macroscale: Here continuum fields such as density, temperature, velocity, stress-fields etc., are important players, and therefore finite-element (FE) methods [Hughes, 1987] are used to examine the large-scale properties of the material. We studied only one particular subclass of multi-scale modeling, known as hybrid quantum mechanical (QM)/classical mechanical (CM) problems, and their application to solid insulators, silica in particular. This scheme is sometimes referred to as QM/MM (molecular mechanics) simulation. In the literature, the classical MD region is further embedded in the FE continuum region [Broughton et al., 1999]. However this MD/FE coupling is unbalanced because of the differences in the dispersion relations and the

PAGE 22

6 underlying elastic properties between the MD and FE continuum region. The work on MD/FE coupling is beyond the scope of our study and is not discussed here. Quantum Mechanical (QM)/Classical Mechanical (CM) Simulations Over the past decade a large body of work has been done on this QM/CM topic. Warshel and Levitt [1976] laid out their first formalism on the QM/CM theory in their seminal paper. Aqvist and Warshel [1993] and Gao [1996] gave a complete overview of QM/CM simulation methods. Applications of QM/MM include study of (i) Biological systems like enzymes, DNAs and proteins [Amara et al., 1999; Monard and Merz, 1999; Hall et al., 2000; Titmuss et al., 2000; Houjou et al., 2001; Nicoll et al., 2001; Cui and Karplus, 2002; Dilabio et al., 2002; Gao and Truhlar, 2002; Murphy et al., 2002; Crespo et al., 2003; Lofere et al., 2003]. For studying enzymes, one applies the quantum method to the active sites in the enzyme molecules where ligands are binding, while the rest of the molecule that does not take part in the chemical reaction is treated classically. Gogonea [2002] reviewed applications of QM/CM methods to biological systems. (ii) Vibrational spectroscopy of complexes [Chabana and Gerber, 2001]. (iii) Electronic excitations [Thompson, 1996; Gao et al., 1996; Vries and Duijen, 1996]. (iv) Hydrolysis of silica [Jung et al., 2001; Du et al., 2003]. (v) Solute-solvent problems [Eichinger et al. 1999; Guo and Xia, 1992; Gao and Freindorf, 1997] in which the solute is described quantum mechanically and solvent that plays the role of an external electrostatic environment is treated classically. However if the QM/CM partitioning has to be done in a single molecule (as in our case) the partitioning becomes more complicated because one or more covalent bonds

PAGE 23

7 have to be cut between the QM and CM regions. Then a proper treatment has to be done to saturate the dangling bonds because these dangling bonds can give rise to incorrect charge densities in the QM domain (chapter 3) and to virtual energy bands in the energy gap [Sauer and Sierka, 2000], electronic states around the Fermi level which are localized at the peripheral region of the cluster containing dangling bonds that are otherwise not observed in bulk Si [Ogata and Belkada, 2004]. Proposed Scheme The challenge of the QM/CM simulation is to move from one length scale to another, as smoothly as possible. We have proposed a scheme for multi-scale modeling that would give a seamless coupling between the QM and CM regions. This is based on the following three criteria: (i) Modeling environmental effects for the quantum domain in an accurate and simple manner. In QM/CM partitioning methods for multi-scale modeling, one is often forced to introduce uncontrolled phenomenological effects of the environment (CM) in the QM domain. Chapter 3 shows that one should consider the precise effects of the CM region in the QM domain to get an accurate description of the QM domain when separated from its CM bulk. Here, accurate means that the quantities under study are indistinguishable whether determined from a length-scale coupled system or from a system comprising only the finest description scale. Our characterization is based on predicted quality of charge densities and forces (which implies correct geometry) in the QM domain [Mallik et al., 2004]. Previous research focused on obtaining only the correct geometry; so far, no studies have reported on getting the correct charge densities to the best of authors knowledge. The rationale behind selecting the charge density as a criterion to assess the quality of QM/CM embedding is explained in chapter 2. The

PAGE 24

8 details of implementations of embedding the QM domain in its classical environment are given in chapter 3 (ii) Developing a new classical potential for the classical region. Regarding the coupling to the CM region, most researchers consider some kind of handshaking region in which they average the force between two regions [Rudd and Broughton, 2000]. This technique of force matching by choosing the average of two forces is however very artificial. We avoided this obstacle by designing a new potential for the entire CM region, such that forces obtained from this potential are the same as those given by the QM domain [Mallik et al., 2005]. Chapter 4 proposes that each development of a multi-scaling algorithm should include the construction of its own classical potential. The potential is constructed using the information about the forces of the QM domain when strained as well as in equilibrium such that the structure and elastic properties predicted by this classical potential agree with those of the QM domain. The use of this new potential will prevent any mismatch of mechanical properties across the QM/CM boundary. (iii) Building a composite material by joining the QM and CM regions with suitable boundary conditions. The composite material is required to be indistinguishable from the one obtained from full quantum or classical calculation on the entire system for small strain elastic response. The accuracy of the proposed method is demonstrated in detail for a benchmark system: a silica nanorod. Model System The complete formalism has been illustrated in detail for a model system: a silica nanorod [Zhu et al., 2003]. This silica nanorod containing 108 nuclei has been chosen as

PAGE 25

9 our model system since a quantum treatment of the whole sample is possible to assess the quality and the physical basis for the proposed scheme and yet the system is big enough to yield bulk properties such as the Youngs modulus. The nanorod has the proper stoichiometric ratio of silicon to oxygen observed in real silica (1:2) and is considered a viable model for studying silica. Figure 1-3 shows front and top views of this nanorod. The nanorod consists of six Si 6 O 6 rings and each of the four planar ring shares, both above and below, a ring of oxygen atoms with adjacent six member rings. To alleviate any dangling bonds, the nanorod is terminated with capped rings whereby each silicon atom of the terminating ring is connected by bridging oxygen or two interstitial oxygens. The size of this nanorod can be readily adjusted by adding or removing Si 6 O 6 planar rings with corresponding O atoms. In the present study, one of the rings near the center of the rod is chosen to be the QM domain and rest of the rod is treated classically (Figure 1-4A). AB Figure 1-3. Structure of the nanorod. A) Front view. B) Top view. The red spheres are the Si atoms and the blue spheres are O atoms

PAGE 26

10 The localized nature of the valence electron charge density (the blue colored cloud in Figure 1-4) of the CM regions ensures the appropriateness of such a partitioning. The lighter shades represent higher electron densities and the black region in between shows that no overlap exists between the two clouds. QM CMCM Figure 1-4. Partitioning of the nanorod into QM and CM domains. Localized (blue) electron cloud of the CM region shows the appropriateness of such partitioning. The first, second and third criteria of multi-scale modeling are detailed in chapters 3, 4, and 5 respectively.

PAGE 27

CHAPTER 2 THE MATERIAL AND ITS FORMAL PARTITIONING Introduction In this chapter, we present the theoretical basis for multi-scale modeling. The objective is to provide a means to assess the approximations made more directly, and also to provide the basis for their extension to other systems (e.g., metals). Here, the problem of multi-scale modeling is introduced by stating what variables are most appropriate for control in the multi-scale modeling. First, the quantum description in terms of ions and electrons is given and the limitations for practical implementation are noted. Next, reduced self-consistent descriptions are given for the ions and electrons, each coupled to the other through their average charge density. In this form, the classical limit for the ions can be taken, allowing the use of MD simulation methods for their dynamics. The equation for the electrons is simplified exploiting the different time scales for electron and nuclear motion (the Born-Oppenheimer limit). The first description constitutes the idealized quantum solid for which the subsequent multi-scale modeling scheme is proposed. This idealized quantum solid is then represented as a composite of two domains, the larger bulk in which a classical representation is used and a smaller reactive domain in which the original quantum domain is retained. The objective of the modeling is then to construct a potential function that gives an accurate description in both the small and the classical domain. This has two components: the determination of a pair potential for the forces on ions in the classical domain and an accurate calculation of the charge density for the ions in the quantum domain. For the latter calculation, 11

PAGE 28

12 approximations are made in term of pseudo-atoms and dipoles to include the relevant effects of the classical part onto the quantum domain. The Idealized Quantum Solid At the fundamental level a simple atomic or molecular solid can be described in terms of ions with charge number for the corresponding atoms of species, and a set of electrons, with overall charge neutrality ( iN iZ i eN ieiiNZN ). To introduce the various levels of approximation we start with the density operator for the system as a whole. Properties of interest are given by the expectation value D A DATrAie, (2-1) where the trace is taken over all electron and ion degrees of freedom. The density operator obeys the Liouville-von Neumann equation 0],[DHiDt (2-2) The Hamiltonian operator H is comprised of the Hamiltonians and for the isolated systems of ions and electrons, respectively, and their interaction iH eH ieU ieeiUHHH (2-3) |'rr|)'r()r('rr eiieddU (2-4) Here )r(i and )r(e are the ion and electron charge density operators. This is the most fundamental level for a quantum description of the system. For a pure state, is a projection onto a state D which is governed by the corresponding Schrdinger equation. For a small system, this can be solved using high-quality quantum methods; however this

PAGE 29

13 is not a practical method for application to larger systems. Here only non-crystalline systems are considered in which one cannot apply periodic boundary conditions. In many cases the properties of interest are functions only of the ion degrees of freedom (e.g. structure), or they are purely electronic (e.g., optical),. Then a description is possible in terms of the reduced density operators for the ions and for the electrons, and, resulting from appropriate partial traces over all the electrons or ions, respectively iAA eAA iD eD eeeeiiiiADTrAADTrA, (2-5) DTrDDTrDieei (2-6) Their equations of motion follow directly from equation 2-2. 0)(1],[iiiiiiitUDDUDHiD (2-7) 0)(1],[eeeeeeetUDDUDHiD (2-8) where the potential energy operators coupling the nuclear and electronic degrees of freedom are )'r(~)r(|'rr|1'rreiiddU ,))()'r(()'r( ~ 1DTrDTreeee (2-9) )'r()r(~|'rr|1'rreieddU 1))()'r(()'r( ~ DTrDTriiii (2-10) This is similar to the microscopic ion electron coupling of 24 except that now the electron charge density operator )r( e is replaced by its conditional average, )'r( ~ e in the equation for and the ion charge density iD )r( i is replaced by its conditional average in the equation for The equations 2-7 and 2-8 are still exact but formal since eD

PAGE 30

14 these average charge densities are not determined by these equations alone. The simplest realistic approximation (mean field) is to neglect the direct correlations in the charge densities i and e i.e. make the replacement in equations 2-7 and 2-10 to get ieDDD i iiiDrTr)( e eeeDrTr)( (2-11) As a consequence, the potentials and become Hermitian and the average charge densities are now self-consistently determined by 2-7 and 2-8. Self-consistency is required since and are functionals of the average charge densities iU eU eD iD i and e respectively. The advantage of this reduced description is that the ions and electrons are described by separate equations, reducing the degree of difficulty of each. The large differences in electron and nuclear masses imply corresponding differences in time scales and thermal de Broglie wavelengths. Consequently, the equation for the ions admits a classical limit for the conditions of interest. The classical limit of equation 2-7 becomes 0}]),[{( ieiiitDUHD (2-12) where { } now denotes a Poisson bracket operation, and is a function of the ion positions and momenta iD )},P{},R({tDDiiii This equation now can be solved accurately and efficiently by molecular dynamics simulation methods, even for large systems. Implementation of equation 2-12 still requires calculation of the potential energy ][eiU which in turn requires determination of the electronic charge density from However, the general solution to 2-8 is a formidable problem: determination of the dynamics of electrons self-consistently in the presence of a changing ion charge density. eD

PAGE 31

15 Two simplifications are made to bring this problem under control. First, is the Born-Oppenheimer approximation and the second assumption is that only the lowest energy state contributes to at any given time. In principle, 2-12 is solved in time steps eD t At each time step the electronic charge density is computed for the ion configurations at that time step in order to recompute new forces for the next time step. The electronic charge density calculation is a ground state eigenvalue problem for the given ion configuration. While simple in principle, the number of time steps and number of electrons for each calculation easily can be unreasonably large for practical implementation. To summarize; the final description of this idealized solid consists of a set of point ions governed by the classical equation 2-12 for their probability distribution )},P{},R({tDDiiii 0}],[{ ieiiitDUHD (2-13) ),'()r(|'rr|1'rrtrddUeii )],r([)'r(),r(tDTrtieeee (2-14) The average electron density )r(e is determined from the ground state solution to 2-8. ,0],[eeeDUHi (2-15) ),'r(),r(|'rr|1'rreietddU ]),,r([)r(),r(ttDTrteiiii (2-16) The analysis proceeds stepwise. The classical equations 2-13 are solved analytically in discrete time steps for the atomic coordinates and momenta. At each step the electron problem 2-15 is solved for the electronic ground state using the ion configuration at the previous time step. From this ground state the electron charge density is determined. This gives the potential energy function and consequently the forces iU

PAGE 32

16 required to change the ion positions and momenta at the next time step. The process is repeated with a new electron charge density calculated at each time step using the new atomic configurations. All electron correlations are accounted for quantum mechanically in the eigenvalue problem; all ionic correlations are determined classically through direct solution of Newtons equations. All structural properties of interest can be calculated since the phase points for the ions are known at all times. The dynamics of the electrons is only coarse-grained as they are slaved to the time dependence of the ions. The foregoing treatment defines a quantum molecular dynamics representation of the idealized quantum solid. The resulting description allows an accurate classical treatment of the atomic structure while retaining relevant quantum chemistry for the interatomic forces due to electrons. If the scheme could be implemented in practice for bulk systems of interest over reasonable time intervals there would be little need for multi-scale modeling. With MD simulation the solution to 2-13 once has been provided is straightforward. So the problem has been reduced to a determination of the electron charge density. Unfortunately, the solution to 2-15 using realistic ab initio quantum methods for even a few hundred ions at each time step becomes prohibitively time intensive. iU The Formal Partition and Composite Solid The Representative Classical Solid In many cases of interest (e.g., equilibrium structure, thermodynamics) the computational limitations of ab initio quantum methods can be avoided through a purely classical representation of the solid. This approach entails an idealization that has many variants. It consists of the representation of the true potential energy function })R({iU

PAGE 33

17 by a suitably chosen function })R({icU Consequently, its form does not need to be computed at each time step and the speed and efficiency of classical molecular dynamics is not compromised. The problem with this approach lies in the choice for })R({icU In principle, an exact mapping for some fundamental property such as the free energy can be imposed ][][ccUFUF (2-17) where is the corresponding classical calculation. Generally, such exact methods can be inverted only perturbatively and lead to a sequence of effective many ion interactions involving increasingly more particles. A more practical method is to assume pairwise additivity for effective point atoms cF jiNNjiijicijVRU,|).RR(|21})({ (2-18) The exact determination of the pair potentials )|RR|(jiijV is now much more restricted as not all properties of interest have such a representation. Nevertheless pair properties such as the radial distribution functions )|r|( ijg might be used to determine the pair potentials. As the )|r|( ijg are unknown and difficult to calculate, the inversion is again difficult and not practical. Instead, experimental data often is used to fit a parameterized functional form chosen for the pair potentials. At this stage control over the approximation is lost and the method becomes phenomenological. This phenomenological approach has been and remains a valuable tool of materials sciences. However, in the context of multi-scale modeling it must be reconsidered carefully. Instead of comparing bulk properties of a solid and its classical representation,

PAGE 34

18 multi-scale modeling constructs classical and quantum models of subsystems and then requires fidelity at their interface. This is a severe test of the modeling assumptions in each subsystem. The approach proposed here confronts this issue directly in the construction of appropriate pair potentials for the problem considered. The Composite Quantum/Classical Solid The use of parameterized classical potentials })R({icU to avoid the time intensive quantum electron calculations cannot be justified for domains far from equilibrium, where an accurate solution to the quantum description for the electrons is required. It is presumed that there is some method for identifying such domains within the solid under given conditions for which quantum chemical effects should be treated in detail. The quantum solid then is represented as a composite of two domains, the larger bulk in which a classical representation is used and a smaller reactive domain in which the original quantum description is retained. The objective of the modeling described here is therefore to construct the potential function such that it gives an accurate description in both the reactive and non-reactive domains. This has two components, the determination of a pair potential for the forces on ions in the classical domain, and an accurate calculation of the charge density for the forces on ions in the quantum domain. iU The classical and quantum domains are defined by labeling all ions as classical or quantum, and spatial domains associated with the coordinates of each, denoted }R{ic and }R{iq respectively. The quantum domains are assumed small, to allow practical calculation of the electronic structure. In principle there could be several disconnected quantum domains, but to simplify the discussion we consider only one. It is assumed initially that the two sets are contiguous with a smooth interface and that

PAGE 35

19 diffusion or migration between them is not significant over the times of interest. In addition to the ions in the quantum domain, there are electrons where is determined by the condition that the quantum domain be charge neutral. The total average electron charge density then is decomposed as Dn Dn ),,r(),r())r()r()(,r(),r(tttteceqCQee (2-19) where )r( Q and )r( C are characteristic functions for and C. The boundaries of the quantum spatial domain Q are fixed by the choice of ions and Q QeqDtdrn),,r( (2-20) where Qencloses all }R{iq It is to be noted that the boundary of is not unique. To make the calculations simple, we chose the boundary of Q to be some smooth surface such as a plane, beyond which the charge density falls below a particular value, e.g. 10 Q -5 The complement of this domain is the classical domain This gives a corresponding decomposition of the total ion potential energy function for the ions (ion -ion Coulomb interactions) plus of 2-14 so that 2-13 becomes C iV iU ,0}),{ iiqiciitDVVKD (2-21) where is the ion kinetic energy and iK ',t))r(')r(('|rr|rd',t))r(')rr(')r(('|rr|'rd)r(rdVeqiQeciCiCic1121 (2-22) and ',t))r(')r(('|rr|rd',t))r(')rr(')r(('|rr|'rdrdVeciCeqiQiQiq11)r(21 (2-23)

PAGE 36

20 The first terms of the integrands on the right sides of 2-22 and 2-23 represent the interactions of the ions with a given subsystem in the presence of the average electronic charge density of that subsystem. The second terms represent the interactions of those ions with their complementary subsystem. Half of the ion-ion potential energy between the two subsystems has been associated with each potential in this decomposition so that the total force acting on the quantum domain by the classical domain is equal and opposite to that on the classical domain due to the quantum domain. The potential is due to ions in the classical domain. By definition these ions are in near equilibrium states and therefore this part of the potential should be represented well by classical pair potentials of the form 2-18, with appropriately chosen parameters. The potential is due to ions in the quantum domain. As this is the domain that can be far from equilibrium the average electronic charge densities must be calculated in detail from the quantum description 2-15. There are two parts to this charge density affecting the ions in the quantum domain, that due to the electrons of the quantum domain icV iqV Dn ),r(teq and that due to the surrounding classical domain ),r()r(tcei The equation governing the electrons of the quantum domain is coupled to this same classical domain average charge density. The scheme proposed here consists of modeling this charge density Dn ),r(tce as an accurate and practical representation for the environment of the quantum domain, allowing calculation of ),r(teq and therefore determining both and iqV .ciV

PAGE 37

21 Description of the Classical Subsystem The steps in constructing the pair potentials are the following. First the specific quantum method to be used in the multi-scale modeling is identified and applied to a large cluster or sample of the solid to be modeled. The forces on atoms for both equilibrium and near equilibrium states are then calculated quantum mechanically. Next, a simple functional form for the pair potentials is chosen and the parameters are adjusted for a best fit to the quantum force data. Chapter 4 details the construction of such pair potentials for the silica nanorod. This method of construction describes the modeling of the first term on the right side of equation 2-22 where both the ions and the average electron charge density refers to the classical domain, assumed in a near equilibrium state. The second term involves a coupling to ions and average charge density in the quantum domain ),(treq This charge density is computed for the forces in the quantum domain (next section) and hence the second term of 2-22 is known as well. In summary, the potential energy for the ions of the classical domain is determined from synthesized pair potentials among the ions and electrons of the classical domain, plus an interaction with the electrostatic potential of the ions and average electron charge density of the quantum domain. Description of the Quantum Subsystem The quantum domain is a subsystem of electrons localized about the designated ions defining that domain. Consider the reduced density operator for electrons defined by Dn Dn enNeneDTrDDeD)( (2-24)

PAGE 38

22 More specifically, this partial trace is defined in coordinate representation [Coleman and Yukalov, 2000] by eeDDeDDDNeNnnNnnenDddDr,..r,'r,..,'r||r,..r,r,...,rr..r'r,...'r||r,...,rDDDnn111n1)(1 (2-25) Clearly the full exchange symmetry among all electrons is preserved. However, this reduced density operator is not specific solely to the quantum domain defined above. For example, the diagonal elements DDnqnDr,..,r||r,...,r11 give the probability density to find electrons at the specified positions, and the latter can be chosen anywhere inside the system. Thus, only when the positions are restricted to the quantum domain does this reduced density operator represents the electrons of that domain. Similarly, if Dn Q )r()( Dne is the charge density operator for electrons its average is Dn )()()()()()()(DDDDDneneneeneeneDrTrDrTrr (2-26) where the trace in the second equality is over degrees of freedom. This average charge density represents the average contribution of electrons at any point Dn r Both DDnqnDr,..,r||r,..,r11 and )r()( Dne change with r since there is an absolute reference background set by the functional dependence on the ion charge density. Consequently, in all of the following discussion of this section DDnqnDr,..,r||r,..,r11 and )r()( Dne are considered only for positions within the chosen quantum domain. Accordingly this coordinate representation }r,...,r{1Dn defines an electron Hilbert space of functions defined over the volume of quantum domain In this context becomes the Dn Q )(DneD

PAGE 39

23 reduced density operator for the quantum domain and )r()( Dne its charge density operator qneDDD)( )r()r()( eqneD (2-27) The equation determining the reduced density operator for the quantum subsystem follows directly from this definition and equation 2-15 for eD .0]),[(qeqeDVKi (2-28) where is the kinetic energy for the electrons and eK Dn )),'r()'r((|'rr|1'r)),'r()'rr()'r((|'rr|1'r)r(r21tdtddVeqiCieqQeqQeq (2-29) The first term on the left side of 2-29 describes the isolated quantum subsystem, composed of the Coulomb interactions among the electrons and their coupling to the average charge density of the ions in the quantum domain. The second term is the interaction of these electrons with their environment, the total average charge density of the classical domain. There are two distinct types of contributions from this charge density of the classical domain. The first is associated with a subset of ions at the border of the QM/CM domains which describe chemical bonds in the full quantum solid. These ions locate regions where there is a highly localized electron charge density shared with the quantum domain, including both strong correlation and exchange effects, localized and non-uniform. The second type of contribution is associated with the atoms beyond next nearest neighbors for which the quantum correlation and exchange effects are much

PAGE 40

24 weaker, and the dominant effect is that of a polarized charge density due to the ions and electrons. So the total charge density for the classical region can be written as Rccieccrr)r()r()|Rr|()()()r( (2-30) where denotes the set of border ions in the environment for which a bond has been broken in identifying the quantum domain. Also, )|'|(Rr is a characteristic function specifying a domain centered on R such that it does not overlap neighboring ions. Its size is taken large enough to incorporate the bound electrons forming the atom for this nucleus in the quantum solid. The second term )(rc is the charge density for all remaining ions and electrons of the classical environment Coulomb Effects of Environment Since by definition, )'r( c does not include any of the chemical bonding with the quantum domain required for its valency, the electrostatic potential associated with it can be expected to have a regular multipole expansion 2d.r)(|'rr|1'rrrdcQ + (2-31) The leading monopole term is zero due to charge neutrality of the classical subsystems (the integral excludes the border ions). The dipole moment d for the entire environment is )'r('r'rd cd (2-32) These results are still quite formal but they provide the basis for the phenomenology proposed for this part of the modeling: replace all effects of the classical

PAGE 41

25 environment on the quantum system, exclusive of bonding, by an effective dipole representing the polarization of the medium by the quantum domain. The origin of this dipole in the above analysis shows that in general it will depend on the geometry and the state of both the classical and quantum domains (e.g., the dipole will change under conditions of strain). In the phenomenological application of this prescription the dipole must be supplied by some simpler means since the electronic contribution to )r( c is not known. Electron Exchange with Environment The contributions from )r( c in the regions where bonds have been cut require a more detailed treatment. Clearly, a necessary condition is that valence saturation should be restored in the quantum domain. However, that is not sufficient to assure the correct charge density or the correct forces within that domain. Instead, the charge density near the border ions responsible for bonding should induce a realistic charge density within the quantum domain. To see how this can be done first write the contribution from one such ion as )r()()|R'r|(|'rr|1'r pcrd (2-33) The potential )r( p represents the actual electrostatic and exchange effects of the ion of the border on the quantum domain. These effects are comprised largely of the ion at the site plus its closed shell electrons distorted by exchange and correlation effects of the site in the quantum domain with which it is bonding. Consequently, an appropriate pseudo-potential is introduced at each such border ion whose behavior in the direction of

PAGE 42

26 the quantum system is the same as that of the actual ion. These ions plus their pseudo-potentials will be called pseudo-atoms. To accomplish this, an appropriate candidate for the pseudo-atom is chosen based on valency of the particular pair of border ion and its neighbor in the quantum domain. Next, a small cluster containing that bonding pair is chosen for training the pseudoatom. The bond is then broken and the relevant member of the pair replaced by the pseudoatom. The training consists of parameterizing the pseudo potential to give the same forces as in the original cluster. For example, the pseudo-atom might consist of an ion plus its closed shell electrons chosen to satisfy the valency for the bond broken. The adjustable parameters refer to a characterization of the closed shell electron distribution. In this way, it is assured that the pseudo-atom not only gives the correct saturation of the dangling bond but also reproduces the forces within the cluster and hence gives a realistic representation of the charge distribution between the chosen pair. The primary assumption is that the bonding of interest is a local effect, so that the pseudo atom trained in the cluster can be used in the bulk solid with a similar accuracy. The training depends on the particular quantum method used to describe the solid, and is considered in more detail in chapters 3 and 6. In summary, the environmental effects on the quantum system are accounted for approximately by a dipole representing its polarization and a collection of pseudo-atoms located at the sites of ions where bonds have been cut QpieqeqQeqrddV)r(d.r))'r()'rr()'r((|'rr|1'r)r(r212 (2-34)

PAGE 43

27 This model now allows solution to 2-28 for the electron distribution in the quantum domain, including its coupling to the classical environment. Then )r( eq is calculated from 2-26 and 2-27. Finally, the desired potentials and of 2-22 and 2-23 are fully determined ciV iqV )),,'r()'r((|'rr|1'r)r('r21|)RR(|21,tddVVeqiQiCjjiNNiijicij (2-35) and QpeqiiQqirtddV)r(d.r)),'r()'rr()'r((|'rr|1'r)r(r212 (2-36) The solution to the classical ion motion 2-21 can then proceed via MD simulation. This completes the proposed scheme for multi-scale modeling of the idealized quantum solid. The verification of this modeling in the silica nanorod and related structures is considered in chapters 3 through 6.

PAGE 44

CHAPTER 3 DESCRIPTION OF THE QUANTUM DOMAIN Introduction This chapter illustrates the technique to meet the first criterion for the multi-scale modeling, proposed in chapter 1, which is to obtain the correct charge densities and forces in the QM domain when isolated from its CM bulk. This requires incorporating the information about the state of the CM environment into the QM domain. The CM region has two kinds of effects on the QM domain: a short-range electronic exchange interaction and long-range Coulombic interactions. Most previous studies of QM/CM simulations neglected the long-range interactions. However it will be shown here that both kinds of interactions must be taken into account for an accurate description of the QM domain. Also the effects of the CM environment should be modeled in a simple enough manner that it can be implemented in the quantum calculations without significantly increasing the computational time. A method for such modeling is shown in this chapter. In doing so, the chapter describes: (i) the quantum mechanical tool employed to describe the QM domain, (ii) an overview of various termination schemes used in the literature to treat the bond cutting region, (iii) the method proposed here for construction of a pseudo-atom to saturate dangling bonds and to account for short-range exchange interactions, and (iv) modeling the rest of the CM environment beyond the next neighbor interactions with lower order multipoles. Finally we compare the charge densities and forces in the QM domain obtained from our scheme of pseudo-atoms and dipoles with those obtained from the conventionally used link atoms (bond saturation with hydrogen atoms) and dipoles. 28

PAGE 45

29 The next section describes the different choices for the quantum mechanical method used to obtain the reference data against which our method of multi-scale modeling is to be judged. Quantum Methods At the quantum level, each possible choice for the ab initio method entails different approximations used to optimize the accuracy and speed of the calculations. Each approximation has its own advantages and limitations, but a general embedding scheme should be insensitive to these. To test how robust the proposed embedding scheme is, we chose two different methods for the underlying quantum mechanical description. The first method is the Transfer Hamiltonian (TH) method [Taylor et al., 2003] based upon the Neglect of Diatomic Differential Overlap (NDDO) approximation. The second method is the Born-Oppenheimer density functional theory (DFT) using the generalized gradient approximation (GGA) [Barnett and Landman, 1993]. Chapters 3 through 5 give the results obtained from the TH method; chapter 6 shows results of the DFT method. The next section briefly introduces the TH strategy as the details are published elsewhere [Taylor et al., 2003]. The DFT method is summarized in chapter 6. Most of the previous work on multi-scale modeling is based on tight binding (TB) models [Broughton et al. 1999; Rudd and Broughton, 2000] because of their ability to treat hundreds of atoms in the QM domain. However, the Hamiltonian in TB models, once formed, is not solved self-consistently. It is oversimplified, and in particular, it cannot account for charge transfer, which is essential to bond breaking. The underlying Hamiltonian in TB models is an extended Huckel Hamiltonian to which a short-ranged repulsion has been added.

PAGE 46

30 Transfer Hamiltonian (TH ): The goal of the TH strategy is to provide forces for realistic MD simulations of the quality of the coupled cluster singles and doubles (CCSD) method [Taylor et al., 2003] in computationally accessible times. Bartlett [1995] gave a comprehensive review of coupled cluster theory and its application to chemistry. The relevant aspects here are that CC theory is the present state-of-the-art method for high accuracy predictions on small molecules and that the method is computationally very costly. The TH is a low-rank, self-consistent, quantum mechanical single-particle operator whose matrix elements are given in terms of parameterizable functions. The functional form of the TH can be chosen to be of any semi-empirical functional form but here the NDDO (Neglect of diatomic differential overlap) form is chosen because of its better qualitative description over Intermediate Neglect of Differential overlap (INDO) Hamiltonians [Hsiao et al., 2001]. The NDDO approximation restricts charge distribution to basis functions on the same center, so this has at most two-atom interactions. Appendix A details the NDDO method. Because of the NDDO Hamiltonian form it is anticipated that the parameters are saturated for small cluster size and that these parameters can be then transferred to extended systems. An overview of the derivation of TH from coupled cluster theory is as follows. In coupled cluster (CC) theory, the wavefunction is written as oTe (3-1) where o is a single Slater determinant and T is an excitation operator that creates single, double, and so on up to n th excited states. T is given by bajijiabijaiiaibatatTTT,,,,21...21.... (3-2)

PAGE 47

31 where and are the standard creation/annihilation operators. In principle, the CC theory considers all possible electron excitations and correlations like post-Hartree Fock theories. However, for simplicity we retain only terms with single and double excitations. This is a standard working approximation in contemporary small molecule computation. The independent particle TH can be connected to many-body CC theory by inserting equation 3-1 into the Schrodinger equation, yielding ia ia ooTTEeHe (3-3) ooEH with (3-4) TTeHeH The energy and the forces given by the gradient of the energy obtained from equation 3-4 are exact up to the level of correlation included (here single and double excitations). The ground state eigenfunction of the effective Hamiltonian is a single Slater determinant, but it has electron-electron correlations corresponding to the truncation of 3-2 employed. To apply the TH to silica system for studying complex processes such as fracture, or hydrolysis where the role of electrons is important, Taylor et al. [2003] parameterized the TH for the following reactions. SiH 2SiH 2 6 3 Si(OH) 3 + SiO(OH) 3 (radical mode of dissociation) Si 2 O 7 H 6 + Si(OH) 3 + OSi(OH) 3 (ionic mode of dissociation) Since the TH models two different dissociation pathways for Si 2 O 7 H 6 electron state specificity is thus systematically introduced into TH. Therefore MD driven by TH forces is expected to predict fracture realistically because in fracture a variety of dissociations always occur. The TH thus is particularly trained to drive high-quality MD

PAGE 48

32 results while requiring only the computational intensity of semi-empirical methods which are orders of magnitude less demanding than CC calculations on the same systems. We will assume that this TH strategy serves as the framework within which we investigate our scheme for the embedding a QM region inside a CM region. In the current application, we have presumed that the forces that arise from the TH are the reference forces by which our proposed method for multi-scale modeling is to be judged. The next section shows the procedure for accurately representing the QM domain when separated from its CM environment. Description of the QM Domain Description of the QM domain requires incorporating the relevant effects of the bulk CM environment on the reactive QM region. There are two qualitatively different effects of the environment. The first is the short-range interaction at the QM/CM boundary, where the partitioning involves cutting covalent bonds between the two regions. The second is the long-range Coulomb interaction between the rest of the bulk and the QM region. As will be seen, both effects must be accounted for to describe the QM region accurately. A successful QM/CM treatment needs a strategy to treat: (i) the bond-cutting region, called termination schemes, and (ii) the long-range interactions through embedding. Different Termination Schemes Over the past decade there have been numerous proposals for different types of termination schemes. The description of all of these is beyond the scope of this work. Only a few of the more relevant ones are discussed in some detail. (i) Link Atom method : The predominant scheme used is the link atom (LA) method, first presented by Singh and Kollman [1986]. In this method, hydrogen atoms are

PAGE 49

33 added to the CM side of broken covalent bond to satisfy the valency of the QM system. There are many variations within the implementation of the LA method, for example, the double Link Atom method [Das et al., 2002], the Add-Remove Link Atom method [Swart, 2003] or the Scaled-Position-Link-Atom Method (SPLAM) [Eichinger et al., 1999]. Reuter et al. [2000] reported that in all cases the LA should be initially placed at approximately 0.97 from the QM atom. One can choose from two different options depending on whether the link interacts with the CM atoms or not. The type of links that do not interact with the CM, referred to as QQ Links give large errors in energy computation. The second type of link called HQ Links, when geometry optimizations are performed leads to incorrect bond length for the frontier QM/CM bond, which in turn can affect the rotational barrier around the frontier bond. So, one has to remove the interactions between the link and the classical group at the frontier. But then, to avoid the electric field problems mentioned above, it is necessary to prevent the frontier classical group from interacting with all QM atoms. This procedure is recommended only if the charge of the classical group at the frontier is small. In spite of all the above-mentioned difficulties, the LA method is still widely used in many types of QM/CM applications because of its simplicity in implementation. For the comparison with the method proposed here, we only studied consequences of the QQ types of LA in detail where the hydrogen atom is placed at a distance of 0.97 from the QM atom at the boundary (which is a Si atom in our case) and the interactions of the LA with the CM atoms are ignored. We also studied the effects of changing the position and alignment of LA. These results will be discussed in the last section.

PAGE 50

34 (ii) Frozen Orbital method : This approach also is known as the Local-Self-Consistent-Field (LCSF) method [Assfeld and Rivail, 1996]. In this method, the electron density in the cut bond is precalculated and then the electronic density along the frontier bond is represented by a strictly localized frozen atomic orbital, which has a preset geometry and electronic population. These frozen orbitals are not allowed to adapt during the QM/CM calculations. The frozen orbitals are defined by their hybridization coefficients and their electron population. These two parameters are determined by quantum chemical calculations on small model systems and are assumed to be transferable to a bigger system. This method avoids the use of any additional atoms but Reuter found it to be less robust compared to the LA method. It was seen that energy calculations are very sensitive to QM size and values of CM charges. He compared the results for proton affinities and deprotonation enthalpies for a number of molecules using LA and LSCF methods and found the results to be comparable. (iii) Generalized Hybrid Oribital method : Gao et al. proposed a refinement of this LSCF concept known as the Generalized Hybrid Orbital (GHO) Method [Gao et al., 1998; Pu et al., 2004] in which the electron density in the cut bond is allowed to readjust. In this method, the hybrid orbitals are divided into auxiliary and active orbitals, the latter of which are then optimized in the QM calculations. But both of these frozen orbital methods require substantial computational work and hence are difficult to implement. (iv) Connection Atom method : In this method [Antes and Thiel, 1999] a connection atom is developed to saturate a C-C bond, such that the connection atom mimics the effect of a methyl group. This connection atom interacts with the QM atom quantum mechanically and the interactions with the CM atoms are handled classically using a

PAGE 51

35 carbon force field. The parameters of the connection atom are determined using semi-empirical methods [Stewart, 1990] such as AM1, MNDO or PM3 designed to reproduce theoretical QM data for energies, geometry and net charges. About 30 different methyl hydrocarbons were used as reference molecules. The parameters adjusted are the orbital exponent s one-center one-electron energy one-center two-electron integral, resonance parameter ssU ssg ss and repulsion term The mathematical functional forms of these parameters can be found in any book on semi-empirical theory [Stewart, 1990] or papers by Dewar and Theil [1976]. There are numerous other termination schemes for the QM/CM boundary like the 'pseudobond' scheme [Zhang et al., 1999], IMOMM [Meseras and Morukuma, 1995; Humbel et al., 1996] and 'ONIOM' procedures [Svesson et al., 1996; Maseras et al., 1999; Dapprish et al., 1999], 'effective group potential (EGP) [Poteau, 2001]. However these methods are not discussed here. Pseudo-atom Method These are developed on the basis of the previously described TH method for saturating a bond terminating in Si (in silica systems). In keeping with the TH strategy, we have based the parameterization of the pseudo-atom by matching forces for the QM portion of the system. We intend that the short-range interactions, particularly the exchange interaction, will be taken into account using the pseudo-atom. The following NDDO parameters were modified to account for the near environment, based on a molecular model, pyrosilicic acid: one-center-one-electron integrals,; Coulomb integrals , ; exchange integral ; two-center-one-electron resonance integrals ssU spU ssg spg ppg sph s p [Dewar and Theil, 1976]. As already mentioned, because of the NDDO

PAGE 52

36 form it is anticipated that the parameters are saturated for small cluster size and that these parameters can be then transferred to extended systems. ModifiedF Figure 3-1. Training of Pseudo-atom on smaller molecules, here Pyrosilicic acid. Figure 3 -1 shows the parameterization scheme for the pseudo-atom on pyrosilicic acid. The part of the molecule within the dotted lines is replaced by a Fluorine (F) atom whose NDDO parameters are then adjusted to give the correct QM forces (which implies correct geometry as well) and charge density in other parts of the molecule (outside the dotted lines). The F atom is placed at the same position as the neighboring CM atom (O in this case) in the bond being cut and hence geometry optimization as in the LA method is not required. This method is similar to the previously described connection atom developed by Antes and Thiel [1999]. However the advantage of using pseudo-atoms is that those are trained to give the correct QM forces and charge densities for both the equilibrium Si-O bond being cut and for small stretches (up to 3-4% from the equilibrium) in the Si-O bondlength as well. This fact will allow us to use the pseudo-atoms even while studying dynamics in the system. The actual charge density from the TH method is calculated instead of the more commonly used Mulliken charge densities [Mulliken, 1955], which

PAGE 53

37 are known to have several common errors ( e.g., a Mulliken charge distribution always has equal apportioning of electrons between pairs of atoms, even if their electronegetavities are very different, which leads to quite unrealistic net atomic charge. In extreme cases some orbitals may contain a negative number of electrons and others more than two electrons). Also in all the other termination methods reported, the role of the more remote CM environment is not quite clear. Calculation of Charge Densities In the nanorod (1-4), the charge densities of the (i) ring with LA, and (ii) ring with pseudo-atoms are compared to the charge densities of the ring in bulk (i.e., the charge density for QM domain resulting from the benchmark QM calculation for the entire rod). The density is calculated on a planar grid placed parallel to the plane of ring (which is taken to be at z = 0). Visualization of the charge densities provides some insights into the descriptive differences between the pseudo-atom and LA. Figures 3-2(A), 3-2(B), 3-2(C) and 3-2(D) show the charge density of the isolated ring, the ring with LA, the ring with pseudo-atoms, and the ring in bulk respectively in the plane of the ring. The six red spots (high density) on the contour of Figure 3-2(A) correspond to the six oxygen nuclei of the ring. It is noted that the charge density found from superposition of single atom calculations of oxygen and silicon atoms is not same as that of the isolated ring. The six blue spots (low density) in the contour of Figure 3-2(B) correspond to the LAs. In the contour of Figure 3-2(C), an overlap between the spots indicates the bonding between the nuclei similar to the ring in the bulk [Figure 3-2(D)]. These plots also indicate that pseudo-atoms are a more realistic representation of the bulk compared to the LAs because they take into account the bonding of the ring.

PAGE 54

38 (A)(B) (D) (C) Figure 3-2. Charge densities with different termination schemes. (A) Isolated ring (B) with links (C) with pseudo-atoms (D) in actual nanorod. Table 3-1 shows the normalized square difference of the charge density for these two terminations with respect to the ring in bulk. Note that as the z coordinate gets larger we approach the pseudo-atom location where we would not expect the termination scheme to reproduce the exact result, hence we stop at z = 0.8 The right-most column of the Table 3-1 shows that, while the difference is least for the ring with pseudo

PAGE 55

39 atoms, still the discrepancy from the exact result is not less than 1% even in the best case and the relative agreement gets worse the further away from the ring we look. Hence, even though the pseudo-atom provides reasonable chemical bonding behavior, there is some significant effect from the rest of the system. Table 3-1. Comparison of charge densities of the ring with LA and ring with pseudo-atoms with that of the ring in the bulk calculated at planes parallel to the plane of the ring. Normalized difference Distance from plane of ring (in ) Density in ring in bulk (electons/Volume) Ring with Link-atoms Ring with Pseudo-atoms Z=0.2 942 0.05 0.02 Z=0.4 507 0.12 0.05 Z=0.6 246 0.32 0.19 Z=0.8 171 0.42 0.42 The forces on the nuclei of the ring with LA and pseudo-atoms were compared to those of the ring in bulk (the benchmark). Figure 3-3 shows the magnitude of the forces on one of the Si nuclei of the ring resulting from the LA or pseudo-atom termination schemes. Since the ring in the rod is in equilibrium, there should not be any forces as is evident from the graph. It is seen that neither the LA nor the pseudo-atom method alone is sufficient to represent the bulk, and additional contributions from the rest of the rod are necessary. These are the long-range Coulomb interactions for which the detailed quantum mechanics is not expected to be important.

PAGE 56

40 01002003004005006007008009001000force (in kcal/mol/A) Ring with linksRing with pseudo-atomsRing in rod Figure 3-3. Comparison of magnitude of force on Si nuclei with different termination schemes. Embedding Scheme To improve our embedding procedure, we describe the remainder of the environment as two dipoles for the top and bottom portions of the rod (Figure 3-4). The values of the dipole have been calculated using the TH-NDDO charge density for these two portions of the rod. These domains are taken to be charge neutral, but are polarized by the presence of the QM domain. The validity of the approximation of the rod by dipoles has been checked by finding the difference between the force on Si nuclei due to all the charges and that due to the dipole and it is found to be very small of (about 6.25 kcal/mol/ and is within our computational error bar which is about 50 kcal/mol/). It is found that even with dipoles placed externally above and below the ring, the force on Si nuclei is 364.7 kcal/mol/. Hence, it is clear that using pseudo-atoms or external dipoles alone will not lead to a satisfactory embedding.

PAGE 57

41 _+ QMRest of the rodRest of the rodPseudo-atoms CMCM QM Long range interactionsmodeled by dipolesShort range by Pseudo-atomsRest of the rod + Figure 3-4. Embedding Scheme: Approximating the rest (CM region) of the rod beyond the pseudo-atoms by dipoles. The problem of large forces was solved by including the two dipoles in TH-NDDO calculation to incorporate the effects of polarization of the ring due to the dipoles. Figures 3-5(A) and 3-5(B) show the significant improvement in the values of forces and charge densities with the incorporation of dipoles respectively. From Figure 3-5(B), it is seen that the normalized difference of charge density is reduced by an order of magnitude. 0100200300400500600700800pseudo-atomsdipolesdipoles inNDDOinequlibriumforce (in kcal/mol/A) 0.0010.010.110.10.30.50.70.9distance from plane of ringnormalized difference Without dipoles With dipoles (A)(B) Figure 3-5. Effect of inclusion of the dipoles in the TH. (A) Improvement in the value of forces. (B) Reduction in the value of normalized difference of charge densities Open symbols, without dipoles; squares, with dipoles.

PAGE 58

42 So we observe that the pseudo-atoms together with the dipoles in the TH-NDDO can reproduce the charge densities and the forces in the QM domain accurately as in the presence of the bulk CM. This scheme was found to be applicable to strained configurations as well. We studied several non-equilibrium cases including: i) When the ring (the QM domain) was radially expanded by 5%, ii) A distorted ring, in which one Si atom is radially pushed out and one Si pushed in (Figure 3-6), iii) A uniaxially strained rod, and iv) A longer rod with ten rings in equilibrium and cases i) and ii) in this longer rod. We have compared the force on one of the Si nuclei (for the distorted ring case we have chosen the force on the Si atom pushed out) in the QM domain (i) with our scheme, (ii) with LA plus dipoles and (iii) the actual force in the molecule. Figure 3-7 compares the forces for all the cases mentioned above. Similarly the charge densities have been compared near the plane of the ring in the same fashion as described above. Table 3-2 shows that the percentage charge density with respect to bulk calculated with our method in each case lies below 1%. Sipushed outSipushed in(A)(B) Figure 3-6. Cross-sectional view. (A) Distorted ring. (B) Ring in equilibrium.

PAGE 59

43 0100200300400500600700 80 0 rod in equi5% expandedringdistorted ringuniaxiallystrained rodlonger rod -10rings5% exp ring inlong roddist ring in longrod3-memberedringforces (kcal/mol/A ) actual our scheme with links Figure 3-7. Comparison of forces on Si nuclei with our scheme and LA to those in the actual molecule for various cases studied. Table 3-2. Percentage charge difference with respect to bulk calculated with our method for various cases studied Various cases studied % charge density difference with respect to bulk calculated with our method I For rod with 6 rings. a) Equilibrium 0.1 b) 5% radially expanded ring 0.5 c) Distorted ring 0.2 d) Uniaxially strained rod 0.7 II. For the longer rod with 10 rings a) Equilibrium 0.1 b) 5% radially expanded ring 0.7 c) Distorted ring 0.5 III. 3-membered ring 0.8 As a further test, we have used this method to partition a highly strained small molecule, which is a 3-membered silica ring [Figure 3-8]. These types of ringed structures are found abundantly on the surface of amorphous silica [Du and Cheng,

PAGE 60

44 2003]. These strained molecules readily hydrolyze in the presence of water and are used to study hydrolysis for a-SiO 2 [Du et al., 2004]. Again, the results for forces and charge density in the QM domain are found to be quite satisfactory (last row of Table 3-2). OSiHQMCM Figure 3-8. Partitioning of 3-membered silica ring. Other Options for Links In the results discussed above, we found that the pseudo-atoms and the dipoles reproduce the force in the QM domain within 1% error but when we placed the LA at a fixed position of 0.97 from the terminating Si atom and along the Si-O bond, there is a large force on the Si atom (Figure 3-3). One might argue here that the results with the LA could be improved if positions of Las are allowed to vary. In our next analysis, the Si-H bondlength as well as the alignment of the LA was allowed to vary to give minimum force on the Si atom. It was found that the minimum force on the Si atom is obtained when the LA is placed at a distance of 1.45 from the Si atom but is not aligned along the Si-O bond. Instead it is about 5 degrees deviated from the Si-O bond (Figure 3-9). The blue atoms are the LAs; the green atoms are the O atoms that are still included in the picture to show the alignment of the Si-O bond. It is seen that the alignment of the Si-H bond is a little bit deviated from the original Si-O bond.

PAGE 61

45 Though the LA at this position gives forces in the QM domain comparable to our scheme with the pseudo-atoms plus the dipoles, it fails to reproduce the correct charge densities as given by the pseudo-atoms. Figure 3-10 shows the charge density with the LA at these deviated positions. The plane at which the charge density is calculated is same for which charge density was calculated with the pseudo-atoms earlier. Here, we do not see any overlapping region between the six spots on the contour. This indicates the LA cannot account for the bonding as successfully as pseudo-atoms. 1.45AoOSiH Figure 3-9. Deviated position of optimized LA. Figure 3-10. Charge density with LA at optimized positions.

PAGE 62

46 LA on O atoms instead of on Si : Another option for putting the links on O atoms was tried. These link atoms were placed at a distance of 0.97 beyond the six intermediate O instead of on Si in the nanorod (Figure 3-11). This O-H bond essentially behaves as an F atom having a valency of -1. It was seen that the LA at this position reproduced the forces in the QM quite correctly. However the role of the external field in this case then becomes unclear as a reduction in the magnitude of forces on the Si atom was observed when the two dipoles were added. We consider again the example of a 5% expanded ring and compare the forces when the links are placed on the O atoms and then with the dipoles. Figure 3-12 shows that the difference of force between when the LA is placed on O atoms and that of the actual molecule is only ~ 52 kcal/mol/ but when the dipoles are included there is a significant decrease in the magnitude of force on Si and the difference with the actual molecule becomes more than 170 kcal/mol/ HOSi Figure 3-11. Position of LA beyond intermediate oxygen atoms. This is somewhat artificial because the CM environment beyond the next neighbor also plays an important role in determining the correct properties of the QM domain, an effect that is anticipated to become more prominent in bigger and unsymmetrical systems. Thus we note that in none of the cases studied does the LA represent a realistic picture.

PAGE 63

47 0100200300400500600700force (kcal/mo/A) actualour schemelinks on Si+dipoleslinks on Olinks on O +di p oles actual our scheme links on Si+dipoles links on O links on O + dipoles Figure 3-12. Comparison of force on Si nuclei with LA on oxygen to LA on silicon and our embedding scheme. Conclusions A tool for embedding a QM domain in its classical environment was developed in this chapter. This procedure requires incorporation of the information about the CM region into the QM domain in a simple fashion without making the quantum calculations more complicated. Both kinds of interactions, short-range as well as long-range, between the QM and CM were taken into consideration. The short-range exchange interactions were modeled by replacing the atoms at the boundary of QM/CM by pseudo-atoms constructed to describe the chemistry of the nearest neighbor interactions. The pseudo-atoms are F

PAGE 64

48 atoms with electronic properties adjusted to yield the correct charge density and forces in the QM part and these pseudo-atoms are placed at the same distance and angle as original O atoms from the Si-O bond to be cut. The second kind of interactions i.e., the long-range many-body Coulombic interactions were represented to a sufficient accuracy using dipoles. The generality of the proposed scheme extends to systems such as strained rods, longer nanorods and small molecule such as a 3-membered ring as well. The pseudo-atom method is a better way for saturating the dangling bond than the conventional LA method. As can be seen from the charge density plots, the pseudo-atoms take chemical bonding into account more or less as in the actual system unlike the behavior of the LA. Several other possibilities for the alignment of LA were tried, but each of those failed to represent the actual charge density in the QM domain.

PAGE 65

CHAPTER 4 CONSTRUCTION OF A NEW CLASSICAL POTENTIAL FOR THE CM DOMAIN Introduction After developing the tools to satisfy the first criterion of multi-scale modeling in the previous chapter, we focus on the second criterion, that is, to obtain a classical potential for the CM environment which yields the same geometry and non-equilibrium elastic properties as predicted by the QM calculations. Independent treatment of the QM and the CM regions can result in mismatching of many physical quantities such as forces, charge densities, and elastic constants across the QM/CM boundary. To avoid this mismatch, the interatomic forces used in the molecular dynamics (MD) simulation in the CM region must be chosen to preserve selected, crucial properties of the underlying quantum theory used in the QM domain. The process of constructing such a potential is demonstrated here. It is developed using ab initio data on the equilibrium structure and on weakly strained configurations calculated from the quantum description (TH), rather than the more usual approach of fitting to a wide range of empirical data. The latter approach is appropriate when a classical description of the entire material is reasonable (e.g., near equilibrium) so that maximum correlation between the material simulated and the real material is attained. In contrast, for multi-scale modeling of strain to failure, the primary constraints on the potential are those of internal theoretical consistency between the quantum and classical descriptions. The final accuracy of the model is therefore set by the quality of the quantum description, not by the fitting of the potential to experiments. The difficulty in extending the conventional approach to the multi-scale 49

PAGE 66

50 modeling problem has been studied by colleagues [Al-Derzi et al., 2004; Flocke et al., 2005; Zhu et al., 2005] The potential is chosen to have the same functional form as TTAM (equation 4.6). A initial guess for parameters is made using a genetic algorithm with force data obtained from a quantum calculation (TH-NDDO). The algorithm effectively locates a domain of the relevant parametric minimum for the given constraints; the parameters found are then checked for stability, followed by a global scaling (described in Appendix C), resulting in the final choice for the parameters that determine all the pair potentials. It will be shown that the Youngs modulus (Y) obtained from this classical potential matches closely with that obtained from the QM method for strains up to 6%, unlike the standard TTAM for which Y differs by almost 20% from the TH. Furthermore, the bond lengths and bond angles in the rod are an order of magnitude more accurate for the new parameterization in comparison to those from the current TTAM or BKS potential parameters. In the first section we discuss briefly the basic concepts of MD used in the program to obtain the stress-strain curve. Allen and Tildesley [1987] and Haile [1992] detail the MD simulation methods. Next we show the failure of standard pair potentials to give the correct value of Youngs modulus for the nanorod. The method for developing a new potential followed by the structure and elastic properties obtained from the new potential are then described. The possible wider applications of this potential to silica systems other than the nanorod are considered in chapter 7. Molecular Dynamics Corrector-Predictor Algorithms Molecular Dynamics allows us to study phenomena under extreme environmental conditions, which may be difficult to study experimentally. In molecular dynamics,

PAGE 67

51 successive configurations of the system are generated by integrating Newtons laws of motion. iiirVtrmtF )()( (4-1) where (t) is the force acting on the i iF th particle of mass m due to all the other particles at time t. The potential, V describes the interaction among the particles. These equations cannot be solved analytically for a large number of particles, and instead are typically solved using finite difference algorithms [Lamberti, 2000]. The most commonly used finite difference methods are the Verlet method [1967] and Predictor-Corrector algorithms. All algorithms assume that the dynamical properties (velocities and accelerations etc.) can be approximated as a Taylor series expansion in time. There are many algorithms for finding the motion of the particle. We used the Gear predictor-corrector algorithm [Gear, 1971]. The predictor corrector algorithm consists of three basic steps. First, new positions, velocities, accelerations and higher-order terms are predicted according to the fifth-order Taylor expansion equations. In the second stage, the forces are evaluated at the new positions to give the accelerations. These accelerations are then compared with the accelerations predicted from the Taylor expansion. The differences between the predicted and calculated accelerations are then used to correct the positions, velocities etc. in the correction step. The molecular dynamics simulations were done at constant temperature. There are many methods for controlling the temperature for example velocity rescaling [Woodcock, 1971], Andersons method [Anderson, 1980], and the Nose-Hoover thermostat [Nose, 1984; Hoover, 1985]. Velocity rescaling is used in our calculations.

PAGE 68

52 Velocity Rescaling This is the simplest and widely used method for maintaining temperature in a MD simulation. The temperature of a system is related to its average kinetic energy given by TNkEKBNVT23.. (4-2) The obvious way to change the temperature is thus to change the velocity. If the temperature at a time t is and the velocities are multiplied by a factor )(tT then the corresponding change in temperature is NiNiBiiBiiNkvmNkvmT112232213)(221 )()1(2tTT )(/tTTnew (4-3) So the temperature can be controlled by multiplying the velocities at each time step by a factor of currreqTT/ where is the current temperature calculated from the kinetic energy and is the desired temperature. currT reqT This procedure is implemented in the MD code used for obtaining stress-strains curves, the results of which are discussed in the following sections. Failure of Standard Potentials It is found that standard classical pair potentials for silica fail to yield the correct value of elastic constants as predicted by the quantum mechanical TH method for even small strains, let alone predicting fracture in materials. This is not surprising, as these potentials have been parameterized using data for equilibrium states only. For example, Figure 4-1 shows the stress-strain curve for the silica nanorod from MD for the entire rod

PAGE 69

53 using the standard TTAM [Tsuneyuki et al., 1988] and BKS [van Beest et al., 1990] pair potentials (equation 4-6). Also shown for comparison are the reference results obtained from the TH method for the entire rod. It is observed that the stress-strain curves given by the TTAM and the BKS potentials differ significantly from the results of the quantum theory. Examination of the curvature of the stress-strain curves indicates that for strains up to 5%, the standard classical potentials predict a stiffer rod than the TH rod. Above 5% strain the nanorod undergoes more of a plastic deformation (characteristic of ductile fracture) for these than the TH rod, for which the stress-strain relation remains linear even for higher strains. The value of Youngs modulus (given by the initial slopes of these curves) is therefore considerably different for each potential. One can see that for strains up to 4%, the Youngs modulus for the TTAM (Y=1214 GPa) rod differs from that of the TH rod (Y=1026 GPa) by 18% while the difference between the BKS (Y=1516 GPa) and TH rod is even more than 47%. These stress-strain curves were generated by the MD simulation in which the nanorod was subjected to longitudinal strain along the z direction by pulling the end caps with a fixed velocity. The remaining atoms in the rod were allowed to relax with their positions determined by MD and the temperature was kept constant by the velocity rescaling method. The stress was calculated by taking the sum of the forces parallel to the loading direction on the constrained atoms divided by the projected cross-sectional area of the nuclei in the end caps. The cross-sectional area was updated after every step to account

PAGE 70

54 for any motion in x or y direction. The simulation was done at a temperature of 10K with a strain rate of 25 m/s and time steps of 2 fs. Figure 4-1. Comparison of stress-strain curves for the standard pair potentials (BKS, TTAM) with that of Transfer Hamiltonian. Objective and Method for Constructing a New Potential The foregoing results confirm that existing potentials cannot be used for consistent embedding of the QM domain in a CM representation of the nanorod. It is proposed here that each development of a multi-scaling simulation on a specific system should include the construction of its own classical potential. That construction should be governed primarily by two components: the chosen quantum method used in the QM domain, and the particular physical properties requiring the QM treatment. In addition, it is essential that there be a limiting domain (usually, near equilibrium) for which the CM and QM calculations of these properties are indistinguishable.

PAGE 71

55 The proposed method is as follows. A QM method (TH in the present case; DFT in chapter 6) is chosen to represent the essential physical mechanisms not explicitly accounted for in a CM representation. Selected properties of the system are then calculated directly using the QM method. Here, these properties are the forces on nuclei incorporating non-local electronic configuration and exchange effects, both for equilibrium and weakly strained conditions. Next, a functional form for the potential energy is chosen for classical point particles that substitute for the atomic constituents. The parameters of the potential energy function are determined by some optimization method based on fitting the data calculated from the QM method. A key issue is the choice of potential energy function. For practical MD simulation it is useful to restrict attention to two-and three-body potentials, but the functional form of these few-body potentials is still optional. Researchers have constructed potentials having as many as 42 parameters [Watanabe et al., 1999]. Other simpler silica potentials having three-body terms exist, Tersoff [1998], Vashishta et al. [1990], Jin et al. [1993], Biswas and Hanmann, [2001], Stillinger and Weber [1985]. There are a few review articles on modeling of silica using different potentials [Schaible, 1999; Brenner, 2000]. Brenner discussed two aspects for the development of an effective interatomic potential energy functions. The first is the derivation of sound functional forms that are motivated by quantum-mechanical bonding principles. The second is the development of empirical corrections and fitting parameters that are often necessary to make a potential function practical for specific applications. Our method is in fact based on both of these aspects. The functional forms that are considered frequently have their genesis in the theory of intermolecular interactions. However, their application to materials simulation is a

PAGE 72

56 rather speculative endeavor. In fact, one can argue that there is no fundamental origin for such potentials since the underlying electronic origin is both non local and many-body in nature. Instead, the potentials are viewed as constructs of the multi-scale model trained to give selected material properties but with no unique status otherwise. Here, for simplicity, only pair potentials are considered and their functional forms are taken to be the same as the phenomenological forms currently in use for silica but with undetermined parameters. Al-Derzi et al., [2004] showed based on high-level QM calculations on several H x Si y O z clusters and -quartz that the parameterization logic behind TTAM parameters is not capable of providing a consistent embedding. We have parameterized the potential only on this particular nanorod as we interested in its properties. However, in chapter 7 we shall see that this potential does have applicability on other silica material. The objective here is to develop a new classical potential that yields the correct equilibrium structure (bond lengths, bond angles) and small strain elastic properties, as given by the quantum TH theory. Such a potential then can be used to meet the third criterion of our multi-scale method, which is to construct the composite nanorod (i.e., the QM and CM regions joined together) with compatible small-strain elastic properties throughout. We emphasize that an attempt to fit the entire stress-strain curve is not reasonable as the physics of large-strain domains is inherently quantum mechanical, thus beyond the limitations of a simple classical potential. Classical potentials cannot account for charge transfer or distinguish between different paths of dissociations: neutral or ionic mode; both of these phenomena are essential to bond breaking.

PAGE 73

57 We assume a classical pair potential that has the same functional form as the TTAM and BKS potentials, but with new parameters determined as described below. Despite its known limitations, this form is chosen because of its simplicity in implementation and widespread use. The potential energy between iand j ions is given by 6)exp()(ijijijijijijjiijrcrbarqqrV (4-4) The first term is the Coulomb interaction for ions of charge and The remaining two terms are collectively called the Buckingham term. They model the short-range repulsive and dipole dispersion, or van der Waals, interactions, respectively. The Coulomb term is determined by a single parameter, the charge on silicon or oxygen, since there must be charge neutrality for the SiO iq jq 2 molecule (i.e.,). There are all together 10 parameters, consisting of the charge q Sioqq2 Si and three pair parameters a ij b ij and c ij for each pair of interactions (O-O, Si-O and Si-Si). In earlier work, the parameters for the TTAM and BKS potentials have been chosen by fitting ab initio Hartree-Fock potential energy surfaces of SiO4(4-)+4e + and SiO 4 H 4 respectively Furthermore, the Si-Si interactions in the BKS potential are arbitrarily set equal to zero. The resulting parameter values are given in Table 4-2. However it is noted that this fitting was done with a low level Hartree-Fock approximation and the energy surface was explored in equilibrium regions only. Here, we choose the parameters in a quite different way. Since MD is using only the forces between ions, the forces determined directly from the TH quantum mechanics of the entire nanorod are used as the primary reference data. Ideally, these forces should be zero for the correct equilibrium

PAGE 74

58 structure and its stability. Furthermore, to assure good elastic properties, information about these forces for small strained conditions also is used in the search for appropriate parameters. The standard approach for parameterization consists of following steps: (i) choose properties of interest (total forces on atoms), (ii) choose reference data (values of forces for equilibrium and strained configurations), (iii) perform some estimation of the parameters, (iv) test against reference data with some error function and then (v) choose a variational method for the next estimate and convergence of the procedure. In our case, the reference quantities are the TH quantum forces on the atoms of the rod in equilibrium, and when strained from its equilibrium structure. The error function for testing the quality of any chosen set of parameters is the weighted sum of absolute difference between the forces calculated from the trial potential and the reference forces, denoted by and respectively for the forces on the ion of the rod. L )(ig )(if thi kkkikifigaL)()( (4-5) The index i runs over all ions of the rod, while k indexes the sum over data sets for different structural configurations of the rod. The numbers denote the weights assigned to the different data sets with ka kka1 The magnitude of the weight decreases with the increasing strain, with the maximum weight assigned to the equilibrium case. A typical distribution of weights was = (0.5, 0.2,0 .1, 0.1, 0.05, 0.05) for percent strains (0, 1, 2, 3, 4, 5), respectively. It is to be noted that this error function is not limited to the form shown; one could choose other forms such as taking the square of difference or normalizing the difference etc. but it was found that these forms gave similar results. This ka

PAGE 75

59 error function is minimized in a space whose dimension is the number of adjustable parameters. The problem of such a high dimensional domain is the existence of several local minima and maxima. Gradient optimizing routines [Nocedal, 1999] seek extrema in the neighborhood of an initial chosen point without sampling the whole space and hence are only local in scope. This characteristic is illustrated for a simple two-dimensional domain in Figure 4-2, where a search near the second minimum would miss the global minimum. zf(x,y) Figure 4-2. Presence of local and global minima. Further improvement can be sought only through some sort of random or unbiased restart. One can repeatedly restart the search algorithm at other points but as the dimension of the space increases such a process becomes highly inefficient. Therefore a genetic algorithm [Goldberg, 1989] has been used since it does not suffer from such disadvantages. Furthermore, the genetic algorithm (GA) does not require calculation of

PAGE 76

60 derivatives and hence is more efficient. We have used the standard PIKAIA code [Charbonneau, 1995] for the GA. Appendix B details the GA method. However, GAs are not without problems. Once a neighborhood of an extremum is reached the GA converges slowly because it relies primarily on the mutation rate to generate small incremental changes in the population. If the global minimum has not been located the search will begin to stagnate. This is due to the successive removal of comparatively less fit individuals from the population, which, after many generations, results in a population of parameter sets, which are all more or less equally fit hence the evolution of the parameters stalls. This problem is analogous to that of the gradient methods getting trapped in local minima. At this point one might think of switching to Newton methods, that is to construct hybrid schemes of GA and gradient-based methods using GA results as a starting guess. Instead we have used an overall spatial rescaling (Appendix C) of the GA parameters to improve the predicted elastic properties (e.g., Youngs modulus) as the final step in determining the choice of potential parameters. Results Finding Parameters for the Potential The parameters obtained from the GA are checked to see whether they produced a stable rod. This is done using a Broyden-Fletcher-Goldfarb-Shanno [Zhu, 1994] energy minimization simulation at 0 K temperature. In this method, the function to be minimized (energy in our case) is expanded to second order in terms of parameter set p (different geometries of the rod in our case) as )()(21)()()(oTooToppHppppgpfpf (4-6)

PAGE 77

61 where g is the gradient of the function and H is the Hessian. The derivative of the function is calculated numerically. To have the geometry with minimum energy requires the determinant of Hessian to be positive. The optimization proceeds with the generation of a new geometry p given an old geometry using This process continues until goes below some acceptable tolerance. As anticipated, it was found that there are numerous local minima in the parameter space, and several possible sets of parameters were found from the GA that yielded a stable rod. Table 4-1 lists these different parameters. op gHppo1)( f Table 4-1. Different values of potential parameters obtained from GA which gave a stable rod Parameters I II III IV V VI q Si 2.36 2.78 2.68 2.22 2.4 2.4 q O -1.18 -1.39 -1.34 -1.11 -1.2 -1.2 a oo 1.78E+03 2.01E+03 1.87E+03 1.9E+03 4.31E+04 1.45E+05 b oo 2.80 2.48 2.68 2.64 4.27 5.91 c oo 239.05 558.28 326.69 366.13 11.64 0.19 a osi 1.06E+04 1.24E+04 1.10E+04 1.02E+04 2.62E+05 4.11E+03 b osi 4.86 4.13 4.65 5.04 7.19 4.32 c osi 64.55 200.46 87.83 50.26 3.87 22.4 a sisi 8.66E+08 1.02E+09 9.06E+08 8.3E+08 2.34E+04 1.86E+05 b sisi 15.33 13.07 14.66 15.6 14.28 6.6 c sisi 22.08 67.6 30.18 16.4 1.05 8.15 The units of a ij b ij and c ij are eV, () -1 and eV-( ) 6 respectively. Table 4-1 shows that the c ij s vary immensely from one set to another, while the other coefficients are of the same order of magnitude for set I-IV. Also the a ij s are more than an order of magnitude different for set V-VI. All this shows that though these parameters correspond to very different points in the parameter space, they all gave rods

PAGE 78

62 with similar geometry. The length of the rods obtained from these parameters differed from the TH rod by 7-12%. Set VI was found to give least structural error in bond lengths and bond angles. It is noted that for this set, the a ij and c ij differ from the other sets by an order of magnitude or more while the values of the b ij and charges are similar to those of other sets. Set VI above was found to give most accurate rescaled parameters. Table 4-2. Parameters for New Potential in comparison to TTAM and BKS potential parameters Parameters TTAM BKS New potential q Si 2.4 2.4 2.24 q O -1.2 -1.2 -1.12 a oo 1756.98 1388.773 133037.6 b oo 2.846 2.760 6.18 c oo 214.75 175.00 0.133 a osi 10722.23 18003.757 3767.451 b osi 4.796 4.873 4.514 c osi 70.739 133.538 15.754 a sisi 8.73E+08 0 1.698E+05 b sisi 15.22 0 6.90 c sisi 23.265 0 5.70 The units of a ij b ij and c ij are eV, () -1 and eV-( ) 6 respectively Table 4-2 shows the values of these rescaled parameters, defining the new potential developed here, in comparison to the usual TTAM and BKS values. Note that the values of the charges and b ij for the new potential are similar to those of TTAM or BKS, but that the c ij are always much smaller for the new potential. This suggests a more repulsive character of this potential. Figure 4-3 compares the potentials for pairwise Si-O interactions.

PAGE 79

63 TTAMBKSTH Classical potential R (angstroms) V (eV) Figure 4-3. Comparison of Si-O pairwise interactions for the three potentials. (The inset shows the zero crossings of the potential) Table 4-3 shows the comparison of several bond-lengths for the three potentials and those for the reference TH rod. Table 4-3. Comparison of structure of the rod obtained from the different potentials. Bond lengths() Bond angles (degrees) TH New Potential % error TTAM % error BKS % error In silica planes Si-O 1.641 1.642 0.04 1.65 0.76 1.616 1.5
PAGE 80

64 Stress-Strain Curves This new potential was used in MD simulations to obtain the stress-strain curves for a classical model of the nanorod. The simulation was carried out under the same conditions as for the TH (temperature of 10 K and pulling the end-caps at 25 m/s). Figure 4-4 shows the results for this new potential, together with those for the TTAM and BKS potentials, in comparison with the TH results. The stress-strain curve obtained from the potential agrees well with the TH curve up to about 5% strain. Accordingly, the Young's modulus obtained from the new potential agrees very well with that from the reference TH. The MD simulations were repeated at lower temperatures and also at lower strain rates. Figure 4-5 plots the stress-strain curve for several temperatures using the new potential. It is observed that the wiggles gradually decrease with lower temperature and then vanish at 0.2K. This result was confirmed by doing an adiabatic expansion of the rod (0.0K) as well, leading to the same results. Strain Stress (GPa)Value of Y in (GPa) Figure 4-4. Comparison of Youngs modulus at low strains for different classical potentials with that of TH.

PAGE 81

65 Strain Stress (GPa) Figure 4-5. Reduction of the wiggles with temperature. Strain Stress (GPa) Figure 4-6. Stress-strain curves at lower strain rate of 2m/s and a temperature of 0.2K.

PAGE 82

66 Absence of wiggles in the TH curve : Classical potentials give numerous wiggles in the stress-strain curves while that for TH is smooth at any temperature. No satisfactory reasons were found for the cause of wiggles. The first anticipated reason was the possibility of large fluctuations in the radius of the nanorod while pulling, as the x and y coordinates of atoms were allowed to relax. However it was found that the fluctuation in the radius was less than 0.03. The second reason thought was the presence of many lower lying phonon modes for the classical potential compared to TH. Thus, the phonon spectra for all the three classical pair potentials (new, BKS and TTAM) were compared to that of TH method. For classical potentials, the phonon spectrum is obtained by solving the eigenvalue problem for the harmonic oscillator while for the TH the eigenvalues are obtained by diagonalizing the Hessian matrix (finite differences of gradients to get the second derivatives of the Hamiltonian). The following histogram presents the number of modes vs. wavenumber in an interval of 100cm -1 (up to1600 cm -1 ) for the TH (green), TTAM (red), BKS (pink) and new potential (blue). At the lowest wavenumber 100cm -1 we see that number of modes of TH is highest (~28) compared to BKS (~ 12), TTAM (~17) or new potential (~27). This is contradictory to our assumption of lower-lying phonon modes for the classical potential. Hence we see neither fluctuation in the radius of the ring nor the phonon spectra are the causes for wiggles in the stress-strain curves for the classical potential. The exact cause of wiggles yet remains an open problem to study.

PAGE 83

67 TTAMTHNEW CLASSICALBKS Figure 4-7. Histogram of number of modes vs. wave number (cm -1 ). Improving Stress-Strain Curves for Higher Strains Figure 4-6 shows that the classical potential works only up to 5% strain. Beyond that, the stress curve for the classical rod becomes more ductile in nature as is evident from the bending of the curve. As mentioned, the classical potential is not expected to fit the entire TH stress strain curve, since the physics of large-strained domains is inherently quantum mechanical. However, it is desirable for multi-scale modeling that the CM domain should grow less ductile with strain than the QM domain so that the material fails only in the QM domain. There are two different methods to make this classical potential more rigid than the TH. The first technique is to induce a charge-transfer term in the potential, and the second is to go back and re-parameterize the potential by including the information about the forces for higher expansions of the nanorod in the GA. We shall

PAGE 84

68 discuss only the results of the latter technique in this chapter; the charge-transfer potential is described in Appendix D. We repeated the process of parameterization using QM force data up to 10% expansions of the nanorod in the GA. The data for forces at the end caps are excluded as the bonding at end-caps became significantly different from the bulk at higher strains and it ends up in misleading the GA, which gives equal weight to bulk as well as end cap data. It is observed that the information about the bulk force data is self sufficient to generate the end caps. The disadvantage of this approach is that incorporation of more strain data sets in the GA results in a delay in its convergence and although it gives a rod with higher value of Y, the accuracy in the bond lengths is lost (by 3%) so the rod has to be rescaled more than once. The following rescaled parameters (Table 4-4) were found to give a rod with similar geometry (Table 4-5) as that of the TH rod and a stress-strain curve (Figure 4-8) that matches with the TH curve well up to 5% but becomes stiffer afterwards. Table 4-4. New parameters for classical potential-2 (NTH-2 Classical Potential) Parameters New-TH Classical potential # 2 Old-TH classical potential # 1 q Si 2.05 2.24 q O -1.025 -1.12 a oo 39439.87 133037.6 b oo 4.66 6.18 c oo 5.85 0.133 a osi 240101.9 3767.451 b osi 7.88 4.514 c osi 2.06 15.754 a sisi 27530.45 1.698E+05 b sisi 21.42 6.90 c sisi 1.55 5.70

PAGE 85

69 The units of a ij b ij and c ij are eV, () -1 and eV-( ) 6 respectively We term this new potential the NTH-2 classical potential (New-TH-Classical potential # 2). Table 4-5. Structure obtained from NTH-2 Classical Potential in comparison to that of the TH nanorod Bond lengths ( ) and bond angle (degrees) TH NTH-2 Classical Potential % error In Silica Planes Si-O 1.641 1.642 0.04
PAGE 86

70 As can be observed from Table 4-5, the accuracy of the bond-lengths and bond angles (compare Table 4-3) is not lost but the rod becomes much stiffer than the TH after 5% strain. This stiffness will guarantee that in the multi-scale treatment the large strain behavior will be localized in the QM domain where it is treated properly with accurate electronic methods. Conclusions The key to successful multi-scale modeling lies in the consistent embedding of a QM domain (TH in our case) in its classical MD region. This requires the CM region to have the same structure and elastic properties as the QM domain. A method is presented for finding a classical potential that would meet the above criteria for the silica nanorod. A pair potential having the same functional form as the TTAM potential is chosen but unlike the standard method of fitting the parameters of the potential to empirical data, the parameters are fitted to ab initio QM force data for the equilibrium as well as for small strained configurations of the particular system under study. A GA followed by a spatial scaling procedure is used to determine these parameters. The potential obtained this way is found to be more ductile compared to the TH method beyond 5-6% strain. The classical potential can be made much stiffer than the TH by reparameterizing the GA using QM force data until 10% strain. Appendix D describes an alternative approach to stiffen the classical rod. The NTH-2 Classical Potential found using the strain data up to 10% is a good candidate for multi-scale modeling of this rod where a seamless coupling between the QM and CM regions of the nanorod is required. Since this potential is much stiffer than the QM method, the phenomena involving large strain will be localized only in the QM

PAGE 87

71 domain. We have essentially found a classical representation of the TH silica nanorod that includes small strain behavior. Further application of this potential to other silica systems can be found in chapter 7.

PAGE 88

CHAPTER 5 CONSTRUCTION OF THE COMPOSITE NANOROD Introduction We developed a classical pair potential (NTH-2 Classical Potential) for the CM region (chapter 4) that predicted the structure and elastic properties in agreement with those of the QM domain obtained from the TH. The potential thus found is a good candidate for multi-scale modeling of this rod. In this chapter we will use this potential to construct a composite rod that is indistinguishable from the TH or the classical rod in terms of equilibrium geometry and Youngs modulus. This indistinguishability is the final criterion of our proposed multi-scale modeling method. Such a composite rod allows us to study multi-scale phenomena for the events involving large-strains localized to the QM domain. Multi-scale modeling done in this fashion will have much of computational efficiency of classical pair potentials but will retain the accuracy of ab initio methods. The composite rod is built by embedding the QM region in its CM environment as described in chapter 3. The forces on the atoms in the QM domain are then given by the TH method while the forces on the atoms in the CM region are calculated by taking the gradient of the NTH-2 Classical Potential. While calculating the TH forces for the atoms in the QM region, the effects of the CM environment are incorporated with pseudo-atoms and dipoles as described in chapter 3. Results for the elastic properties of the composite rod are reported in this chapter. 72

PAGE 89

73 Results for the Composite Rod Figure 5-1 shows the stress-strain behavior of the nanorods obtained from three different methods: (i) the full quantum mechanics TH method, (ii) full NTH-2 Classical Potential and (iii) the composite rod constructed as described above. The three overlaying curves (measured at 0.01K) indicate that the composite rod is identical to the rod obtained from TH or the NTH-2 nanorod in terms of small strain elastic properties and structure. The stress-strain results shows the success of our multi-scale method indicating that the composite rod is indistinguishable from the underlying quantum mechanics for states near equilibrium.. However the elastic properties of the composite rod obtained from both of these potentials do not agree beyond 10% strain. This is because the pseudo-atoms, trained at regions only close to the equilibrium configuration, fail to give the correct charge densities at such high strains. This diagnosis was checked by comparing the charge density in the QM domain of the 12% strained rod with that of equilibrium configuration and a difference of 6% was found. It is conceivable that retraining of the pseudo-atoms may improve the stress-strain performance of the composite rod beyond a strain of 10% but our next goal is to find an application with the method developed so far in some real silica systems: for example cristobalites and glasses. It is be noted that the 10% strain here corresponds to approximately 10% strain in the bond-lengths, so the failure of pseudo-atoms at high strains does not pose any constraint so far as application of the proposed multi-scale scheme to glass is concerned since glass fails much before 10% strain in bond length is reached. The reason is that real systems like glasses have many inherent defects which act as stress concentrators that cause the material to break at much lower strains than

PAGE 90

74 observed for the nanorod. One should note that the yield stress of this nanorod (Figure 5-1) from the TH method is around 190 GPa, which is notably larger than any observed silica system probably because of its high aspect ratio (one-dimensional system) and defect free nature. NTH-2 Figure 5-1. Comparison of stress-strain of the composite rod with those obtained from NTH-2 Classical Potential and TH. The satisfactory construction of such a composite rod completes our recipe for the multi-scale modeling for all practical purposes. Conclusions In this chapter, we have satisfied the last criterion of our proposed multi-scale modeling scheme, namely to construct a composite nanorod that would be indistinguishable from the TH rod for all relevant properties (structure and elasticity) in states near equilibrium. Such a composite rod can be used for multi-scale studies with substantially reduced computational times.

PAGE 91

CHAPTER 6 INDEPENDENCE OF MULTI-SCALE MODELING WITH RESPECT TO THE CHOICE OF UNDERLYING QM Introduction As mentioned in chapter 3, a proper multi-scale modeling procedure should be independent of the choice of underlying quantum mechanical method. To test the method we have developed, we have essentially repeated the work done in chapter 3 and 4 (embedding of the QM domain and finding a classical potential) with Density Functional Theory (DFT) instead of the Transfer Hamiltonian as the quantum mechanical approximation. In the field of DFT, 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 beyond the scope of this thesis. Hence only a short review of the basic concepts relevant to this chapter, such as Born-Oppenheimer approximation, Hohenberg-Kohn theorem, Kohn-Sham equations, pseudopotential, and dual space formalism will be presented. More detail on DFT can be found in Parr and Yang (1989), Wimmer (1997), Jones and Gunnarsson (1989), and Dreizler and Gross (1990). The primary results of this chapter are a confirmation that the isolation of the QM domain with pseudo-atoms and dipoles, plus the construction of a classical potential based on the DFT forces leads to accuracies of the same quality as those described already using the TH quantum mechanics. Hence, the multi-scale modeling scheme is faithful to the chosen form for the underlying quantum mechanics in both cases. 75

PAGE 92

76 Review of DFT Method Born-Oppenheimer Approximation This approximation follows from the fact that the ions are much heavier than the electrons, so the electrons react instantaneously to the position of the ions. Thus the ions can be treated adiabatically, leading to the separation of the electronic and ionic coordinates in the many-body wave functions. This reduces the many-body problem to the solution of the dynamics of the electrons in a frozen configuration of ions. The total energy of a system consisting of ions and electrons can be reduced to )]([||/||)})({},({221rERRZZRmrREelecJIJIJIIIIItotal (6-1) where, and are the position, mass and charge of the IR Im IZ I th ion. The first two terms correspond to the kinetic and Coulomb interaction energy respectively. )]([rEelec is the ground state energy of the electrons evaluated for a particular configuration of ions which we consider next. }{IR Hohenberg-Kohn Theorems The DFT approach is based on the first Hohenberg-Kohn (HK) theorem [1964], which states that all the ground state properties of a system are functionals of the electron density. The second Hohenberg-Kohn theory enables us to write the total energy as a variational functional of electron density )( elecE (6-2) ][)()(Fdrrvr where ][][][ eeVTF (6-3)

PAGE 93

77 ][ T is the kinetic energy ][ eeV is the electron-electron interaction,. The ground state electron density is the density that minimizes ][ E subject to the condition (6-4) Ndrr)( N being the total number of electrons. ][ F is a functional of the density and is independent of the external potential but no general exact expressions for it are known. Kohn and Sham (1965) invented an ingenious approach to the kinetic energy functional ][ T that turned the Hohenberg-Kohn variational functional into a practical tool for rigorous calculations. The next section shows the derivation of the Kohn-Sham equations from the HK theorems. Kohn-Sham Equations The Functional ][ F can be decomposed into three different components as ][][][][ xcsEJTF (6-5) ][ sT is the kinetic energy of a non-interacting inhomogeneous electron gas in its ground state, ][ J is the classical coulomb energy (or Hartree energy) and the functional ][ xcE is the exchange-correlation energy. ][ xcE in fact includes exchange, self-interaction, Coulomb correlation terms as well as the KE increment ][][ sTT Using equations 6-2 and 6-3 the total energy functional can be written as ][][][)()(][xcselecEJTdrrrvE (6-6) Minimizing this expression subject to equation 6-4 gives 0)(][rEelec being a Lagrange multiplier that enforces 6-4

PAGE 94

78 )(][rTvseff (6-7) where is called the effective KS potential defined as effv )('|'|)'()()()(][)(rvdrrrrrvrErJrvvxcxceff (6-8) From the definitions of the non-interacting system 2|)(|)( iiirfr ( is the electron occupation number) and if iiiisfT2*21 the variation in is equivalent to variation of the i hence one gets (in Hartree) iiieffrv)](21[2 (6-9) These are known as Kohn-Sham equations [1965]. They have to be solved self consistently as depends on effv )(r .Though the KS equations incorporate the exact non-interacting kinetic energy ][ sT the exact form of the exchange correlation functional is still unknown. Many approximations for ][ xcE have been developed. The ones used for our calculations are discussed below. Local Density Approximation The first such approximation proposed by Kohn and Sham is known as the Local Density Approximation (LDA) which uses the basic idea of using the exchange correlation energy of the homogenous electron gas to evaluate ][ xcE Thus (6-10) drrrExcXC)()(

PAGE 95

79 where xc is the exchange correlation energy per particle of a uniform electron gas of density This model assumes the electron density to be locally equivalent to a uniform gas, hence the corresponding exchange-correlation contribution to equation (6-10) becomes )()())(()()(xcLDAxcLDAxcrrrErv (6.11) Considering the inexact nature of the LDA, it is remarkable that calculations using LDA have been so successful. Generally, the total energy differences between related structures are accurate to within a few percent. However there are many failures of LDA also. For example, it can predict incorrect cohesive energies, often over binds molecules and solids, and it fails to explain certain classes of molecular bonds (e.g., hydrogen bonds) for which the charge density in the binding region is very small. Therefore we have used a more refined approximation for known as the Generalized Gradient Approximation. xcE Generalized Gradient Approximation The LDA is based on the expressions for the energy density of a homogeneous electron gas. For slowly varying densities, the energy functional can be expanded as a Taylor series in terms of the gradient of the density, with the first term being the LDA term. Retention of the first two terms is known as the gradient expansion approximation. However the series is shown to be nonconvergent [Kleinman and Sahni, 1990]. The approximation can be made physically sensible, however, and several such approximations have been proposed which now are known as generalized gradient approximations (GGA) (e.g., the PBE named after Perdew, Burke and Ernzerhof [1992]).

PAGE 96

80 In many applications, the GGA shows improved results over LDA in total energies (Perdew et al., 1992), atomization energies (Becke, 1992; Perdew et al., 1992, 1996; Proynov et al., 1995), energy barriers (Hammer et al., 1993; Hammer and Scheffler 1995; Philipsen et al., 1994) etc. Typically, GGA favors the density inhomogeneity more than LDA. However, there are notable exceptions [Trickey, 1997]. Throughout this discussion we have neglected the spin of the electron. The code we used employs Local Spin Density Functional Theory (LSDFT). In this theory, electron and the spin densities are fundamental quantities with net spin density being the difference between the up-spin and down-spin electrons: )()()(rrr (6-11) The total electron density is just the sum of the densities for the two types of electron. The exchange correlation functional is different for the two cases, leading to a set of spin polarized Kohn-Sham equations ,)()(],[)(211121222rrrVdrrrRZiiiXCIII (6-12) So we have two sets of wavefunctions, one for each spin. Troullier-Martin Pseudo-Potential Implementation of all-electron DFT with a plane wave basis set is computationally too intensive. This section introduces the concept of pseudopotential, which reduces this computational cost. The pseudopotential replaces the effect of the chemically inert core states upon valence states. This replacement saves the effort of using a large number of plane waves to expand both the tightly bound core orbitals and the rapid oscillations of

PAGE 97

81 the valence Kohn-Sham orbitals in the core region. The pseudopotential can be obtained from all-electron atomic calculations by solving the radial KS equations self consistently. )()(];[2)1(21222rrRrrRrVrlldrdnlnlnl (6-13) where )(];[];[LDAxcHvrVrZrV ];[rVH being the Hartree potential The pseudopotential has to meet certain conditions: (i) the pseudo-wave function generated from this potential should not contain any nodes (it should be smooth and not have any wiggles), (ii) the normalized radial function of the pseudo-wave function should be identical to the normalized radial real wavefuntion outside a cut-off radius (iii) the charge enclosed within should be equal for the two wave function, and (iv) the valence of all electron and pseudopotential eigenvalues must be equal. A pseudopotential which meets all the above conditions is said to be a norm conserving potential [Troullier and Martin, 1991]. Once the pseudowave function is obtained, the screened pseudopotential can be obtained by inversion of equation ctr ctr )]([)(212)1()(222,rrRdrdrrRrllrVPPlPPllPPlscr (6-14) The screening of the valence electrons is removed and an ionic potential is generated by subtracting the Hartee and the exchange correlation potentials calculated from the valence pseudo-wave functions from the screened potential )(rVPPH )(rVPPXC (6-15) )()()()(,,rVrVrVrVPPXCPPHPPlscrPPlion

PAGE 98

82 Also an accurate pseudopotential should reproduce the scattering of the ionic core. This requires a nonlocal pseudopotential, that is, one which uses a different potential for each angular momentum component of the wave function. The pseudopotential operator can be decomposed into local and nonlocal parts, (6-16) lnlcPPllcPPPPPrVVrV)()(,, where is the local potential and is the non-local potential. is the operator that projects out the l )(,rVlcPP )(,rVnlcPPl lP th angular momentum component from the wave function. In principle, the local potential can be chosen arbitrarily but since the summation in equation 6-16 has to be truncated for some value of l, the local potential should be chosen such that it adequately reproduces the atomic scattering for all the higher angular momentum channels. These forces are then used at each MD step to move the ions on the BO potential energy surface. Dual Space Formalism The computer code we used employs the dual space formalism for calculation of the DFT energy. This method evaluates the kinetic energy in momentum space and potential energy in real space (Martin and Cohen, 1988; Wentzcovitch and Martin, 1991). It is found that this method is computationally more efficient than the pure momentum-space formalism (Ihm et al., 1979) in which all calculations are performed in momentum space. Computational Details All the theory discussed above is implemented into a parallel multi-scale program package, known as Born-Oppenheimer molecular dynamics BOMD [Barnett and Landman, 1993]. The separation of time scales between the electron and the ionic motion

PAGE 99

83 in the BO approximation restricts the dynamical evolution of the ionic system to a single electronic potential surface (PES). For a given ion configuration the PES and forces on the ions are calculated. These then are used to govern the classical equations of motion of ions advancing them by a small time increment to a new configuration. Evaluation of Hamiltonian matrix elements and the operations on the wave functions are performed using the dual-space formalism. A plane wave basis set and cut-off energy of 30.84 Rydbergs is chosen. Results and Discussion We next use the DFT to study the multi-scale modeling of the nanorod. The accuracy our multi-scale scheme is tested against charge densities and forces obtained using DFT on the whole nanorod. In the next section, we discus the training of the DFT-pseudo-atoms that model the nearest neighbor interactions. Training of DFT-Pseudo-atoms A pseudo-atom is constructed based on the Troullier-Martin (TM) pseudopotential discussed above. A cut-off radius of 1.5 is used. Once again the pyrosilicic acid molecule is chosen for parameterization of the F atom and its position is constrained to be at the same place as the O atom (Figure 3-1.). But unlike the TH-NDDO method in which we could change the electron-ion as well as electron-electron interaction parameters, in DFT we can modify only the electron-ion interactions. The three options that can be varied to alter the electron-ion interactions for the F atom using the TM pseudopotential are: (i) the core charge on F, (ii) omission of the non-local part in the potential, and (iii) switching the local and non-local part between s and p orbitals. We tried all three possible choices and their combinations to find the best alternative in the sense of reproducing the

PAGE 100

84 force on the terminated Si atom in the pyrosilicic acid as in the presence of the whole molecule. Table 6-1. Magnitude of force on the terminated Si atom in pyrosilicic acid with changing parameters of the dft-pseudo-atom zcore Only changing zcore Without non-local part Switching between s & p 7.0 0.88 2.9E-04 0.08 6.8 0.97 0.04 1.56 6.2 1.15 1.15 1.87 7.2 0.72 0.81 0.16 7.4 0.62 0.78 0.13 7.6 0.76 0.88 0.74 the units of forces are in Rydberg a.u This method of parameterizing a pseudo-atom by trial and error might seem arbitrary but as it is often said parameterization is not a science but an art. Our case of obtaining a pseudo-atom truly exemplifies this quote. We started with the same core charge (7.0) as in real F, then changed the core charge in steps of 0.2 while neglecting the non-local part of the potential in the second case and interchanging s and p orbitals in the third case until the force on the Si atom is minimized. Table 6-1 shows that the minimum force is obtained when the core charge is 7.0 and the non-local part is omitted. The force for this particular set of combination is lower by two orders of magnitude than for the other cases. This pseudo-atom is then transferred to the nanorod. The nanorod is partitioned in QM and CM regions and the CM portions are approximated by dipoles as before (Figure 3-4). The values of the dipoles are calculated from the charge density obtained from the DFT calculation for the CM portions of the rod. We have computed the force on the Si atom (in the QM domain with

PAGE 101

85 this DFT pseudo-atom plus the dipoles) for the rod in equilibrium and all the strained cases considered in chapter 3. The results obtained from a DFT calculation on the whole rod are taken to be the actual forces. Once again, we have compared our results to conventional link atoms (LA), which are placed at 0.97 from Si. Figure 6-1 shows the magnitude of the force on the Si atom in all the cases studied. 0200400600800100012001400rod in equi5% expandedringdistorted ring3-memberedringforce (kcal/mol/Ao) actual our scheme with links Figure 6-1. Comparison of force on the Si nucleus with different termination schemes and dipoles for the DFT case. Next, we compare the charge densities for different strained cases in the nanorod using DFT-pseudo-atoms and dipoles. Table 6-2 shows the percentage difference of the charge densities of the whole nanorod and that obtained from our scheme, in a plane located at 0.2 from the plane of the ring. Table 6-2. Percentage charge density difference obtained with our scheme versus the charge density in the actual molecule for different cases studied in the DFT method Cases Studied % Charge density difference In equilibrium 0.34 5% expanded ring 0.86 Distorted ring 0.95 3-membered ring 0.82

PAGE 102

86 Like the TH method, here too we observe that the forces and the charge densities in the QM domain when isolated from its bulk can be generated using DFT pseudo-atoms and dipoles within 1% accuracy. Construction of Classical-DFT Potential for the CM Region Having completed the embedding of the QM domain, focus is shifted to finding a classical potential for the CM region. A classical potential having the same form as TTAM is found that will predict the same geometry and elastic properties as the QM domain. A GA with DFT force data up to 4% expansion followed by the scaling procedure is used as before (chapter 4) to find the parameters for the potential. Table 6-3 shows the resulting parameters in comparison to the standard TTAM and BKS parameters. Again, we see that the charge on the ions is lower, in this case almost halved, from that given in the TTAM and BKS potentials. Also, the van der Waals interaction, the c ij s, are much less than in the published potentials. (One might argue that in TTAM potential the information that should be carried by term, is not accurately represented by DFT. The DFT result as atoms get farther apart is not accurate, it falls off exponentially not like However recall that BKS and TTAM use a term but neither uses any information far enough from equilibrium to be considered a van der Waals interaction). 6/rc 6/1r 6/1r In table 6-4, the bond lengths and bond angles of the rod obtained from this potential are compared to the DFT treatment of the rod. Clearly the agreement between the classical and DFT rods is quite good. We have also presented the results for the BKS and TTAM rods. The BKS rod shows a close resemblance in structure with the DFT rod.

PAGE 103

87 Table 6-3. Parameters for New-Classical-DFT Potential in comparison to standard TTAM and BKS potential parameters Parameters New -Classical DFT potential BKS TTAM q Si 1.4 2.4 2.4 q O -0.7 -1.2 -1.2 a oo 1281.153 1388.773 1756.98 b oo 3.655 2.760 2.846 c oo 3.899 175.00 214.75 a osi 7846.292 18003.757 10722.23 b osi 6.1194 4.873 4.796 c osi 1.295 133.538 70.739 a sisi 6.28E+08 0 8.73E+08 b sisi 25.414 0 15.22 c sisi 0.595 0 23.265 Table 6-4. Structure of the nanorod obtained using New-Classical-DFT, BKS potentials versus structure of DFT method Bond lengths () and angles DFT New DFT Classical % error BKS % error TTAM % error In planar silica rings Si-O 1.62 1.62 1.62 1.65 1.85 Si-O-Si 158.02 157.96 0.04 161.1 1.9 162.90 3.09 Between planes and intermediate O O-Si-O 107.15 103.74 3.0 104. 2.9 104.55 2.43 At the end-caps Si-O 1.64 1.63 0.6 1.62 1.2 1.67 1.83 Si-O-Si 92.57 100.5 8.0 100.8 8.9 100.55 8.6 Approximate length 15.79 16.1 1.96 15.86 0.44 16.31 3.29 Approximate diameter 6.44 6.44 6.41 0.47 6.57 2.02 We next want to compare the Youngs Modulus of the rod obtained from this potential with that obtained using DFT. However unlike an MD using TH forces, a stress strain curve from MD using DFT forces cannot be calculated, as the DFT treatment is computationally too intensive to study the dynamics of a system of such size (containing 108 atoms). Taylor [2004] reported that the computer time (CPU seconds) necessary to

PAGE 104

88 perform one energy plus gradient (force) step calculation (IMB Power 2S chip) for a pyrosilicic molecule (15 atoms) is 375 seconds for DFT whereas its only 0.17 seconds for the TH. Since in the DFT case the QM calculations are considerably more time intensive, only equilibrium and selected adiabatic strain configurations were calculated with DFT. The equilibrium structure was determined by sequential DFT calculations and nuclear relaxation to find the minimum energy configuration. The strained configurations were obtained from an affine transformation of the minimum energy configuration by 1, 2, 3, and 4 %, with a single DFT calculation of forces at each of the expanded configurations. Then the average force on the Si atom in each ring of the DFT nanorod is computed for these four cases. So we have compared the stresses for adiabatic configurations for DFT, new-DFT-classical and TTAM potentials using the procedure just described. The values of stress obtained using BKS potential are similar to those of TTAM potential. The stress at each configuration was calculated as before by dividing the average force on end-cap atoms by the projected area of the end cap. Strain in Figure 6-2 corresponds to affine expansions. Figure 6-2. Comparison of adiabatic stress curves for DFT, new-DFT-classical and TTAM potential.

PAGE 105

89 Since the forces computed from this potential are approximately equal to the DFT forces, this potential can be used for multi-scale modeling of the rod using DFT forces for the atoms in the QM domain. Conclusions A DFT pseudo-atom is generated by modifying the electron-ion interaction parameters of F atom based on TM pseudopotential. The core charge was 7.0 and the non-local part of the TM pseudopotential was omitted. The pseudo-atom is placed at the same position as the O atom in the Si-O bond being cut. This is found to be a better termination scheme than the LA atoms The DFT pseudo-atoms together with the dipoles calculated from DFT charges can reproduce the charge densities and dipoles in the QM within 1% accuracy, behavior much like the TH method. A classical potential of the same form as the TTAM is parameterized to give geometry and elastic properties up to 4% stains in agreement with the DFT method. The rod was made further stiffer using a charge transfer term in the potential (Appendix D). The potential found this way is a good candidate for multi-scale modeling of the nanorod using DFT forces in the QM domain. These results show that our multi-scale modeling method is independent of the underlying quantum mechanical method used.

PAGE 106

CHAPTER 7 APPLICATIONS OF THE MULTI-SCALE MODELING PROCEDURE TO OTHER SILICA SYSTEMS Introduction In chapter 5 we illustrated the success of our multi-scale method for the silica nanorod using the TH quantum mechanics and the corresponding NTH-2 Classical Potential. In this chapter, we extend the application of this multi-scale scheme to other silica systems with the long-term objective of studying the effects of hydrolysis on mechanical properties. The systems considered are: (i) a 10-membered silica nanoring, (ii) a notched silica nanorod containing one oxygen less than the nanorod and (iii) a bulk glass containing 2940 atoms. Each of these are susceptible to rapid hydrolyzation (the nanorod without notch is stable, defect free and chemically inert, so it is not appropriate for such investigations of chemo-mechanical properties). Before continuing to study hydrolysis or any other complex phenomena in these systems, it is necessary first to see if the NTH-2 Classical Potential yields the precise structure and elastic properties of the target systems. This is an open question since the potential has not been trained for these systems. Here we show that although the NTH-2 Classical Potential was trained on the nanorod, nevertheless it predicts the structure and elastic properties for these other silica systems within 4% accuracy. Also, the composite systems constructed using the TH quantum mechanics and this potential for the 10-membered nanoring and the notched rod were found to be indistinguishable from their corresponding full quantum systems for up 90

PAGE 107

91 to 9% strain. These results indicate the success of our multi-scale scheme. In the final section of this chapter, we report similar results for the 10-membered nanoring and the glass obtained using the New-DFT-Classical Potential (developed in chapter 6). Only the structure comparisons are done in this case, because the calculations for the elastic properties using DFT were limited by computer and time constraints. Systems Studied Silica Nanorings Ring structures such as that shown in Figure 7-1 have been shown to exist on the surface of amorphous silica as highly strained molecules [Bromley et al., 2003]. Consequently these rings act as preferential sites for water attack and get hydrolyzed very quickly. It was found that such ring structures are energetically more stable than the corresponding linear chain structures and that these rings display frequency modes that fit well with the infrared bands measures on dehydrated silica surfaces. These features indicate their potential as models of strained extended silica systems. QM Figure 7-1. The 10-membered silica nanoring. In the original work, the structure was found from DFT calculations using the Becke-type three parameter hybrid exchange-correlation functional B3LYP [Becke,

PAGE 108

92 1993]. Here, however we chose to use the structure obtained from the TH as our reference data. Table 7-1 compares the structure obtained from the TH, the NTH-2 Classical Potential, and the TTAM potential. As in the original reference, we too found that the internal strain of the molecule leads to some folding of the two rings along their Si-Si axes and a mismatch between the inner and outer Si-O-Si internal two-ring angles and the Si-O bond lengths. The inner Si-O-Si internal two-ring angles are slightly larger than the outer internal two-ring Si-O-Si angles and the inner Si-O bond lengths are slightly smaller than their outer counterparts as can be seen from the table 7-1. It is surprising to note that though the NTH-2 Classical Potential was originally trained with the force data for the nanorod, it can predict the structure of the nanoring successfully with respect to that given by TH. The percentage error in bond angles and bond lengths lies within 5%, considerably better agreement than that obtained from the TTAM potential. Table 7-1. Comparison of structure of the nanoring obtained from TH, NTH-2 Classical Potential and TTAM potential Bond lengths () and angles TH NTH-2 Classical Potential % error TTAM % error Si-O (outer) 1.79 1.73 3.3 1.66 6.7 Si-O (inner) 1.78 1.72 3.3 1.65 7.3 O-O (outer) 2.425 2.324 4.2 2.173 10.4 O-O (inner) 2.421 2.321 4.1 2.171 10.3 Si-Si 2.56 2.49 2.7 2.44 4.7 < Si-O-Si (outer) 90.84 92.83 2.2 95.80 5.5 < Si-O-Si (inner) 91.79 93.3 1.6 96.25 4.9
PAGE 109

93 We next determine the stress-strain curve using a MD simulation as in chapter 4. In this simulation, the strain is imposed by giving all the Si atoms a fixed radial velocity while the rest of the atoms are allowed to relax. The MD was done at a temperature of 0.01K as before. The stress-strain curve obtained using the TH and NTH-2 Classical Potential is shown in Figure 7-2. From the figure, it is seen that the Youngs modulus obtained from NTH-2 Classical Potential matches well with the TH curve up to 9% strain. Hence we see that the NTH-2 potential gives both the correct equilibrium structure and the correct elastic properties. Finally a composite nanoring was constructed, in which one of the silica units is considered to be the QM domain (Figure 7-1). The two O atoms on each of the terminated Si nuclei were replaced by the pseudo-atoms described in chapter 3. The stress-strain behavior of this composite nanoring is shown in Figure 7-2. It matches with the TH curve well up to a strain of 9%. Beyond that it disagrees with the TH result due to failure of pseudo-atoms at this high strain. NTH-2 Strain Stress (GPa) Figure 7-2. Stress-strain curve for the nanoring obtained from TH and NTH-2 Classical Potential.

PAGE 110

94 We have also compared the structure of 2-membered (Si 2 O 2 ), and 3-membered (Si 3 O 3 ) silica rings (Figure 7-3) using TH quantum mechanics, and the NTH-2 and TTAM potentials. In these small clusters also the NTH-2 potential shows a better agreement with TH geometry than does the TTAM (Table 7-3). Table 7-2. Structure comparison of 2-membered and 3-membered silica rings Bond lengths () and angles TH NTH-2 % error TTAM % error Si 2 O 2 Si-O 1.73 1.65 4.0 1.61 6.9
PAGE 111

95 of these molecules, the percentage error in bondlengths and bondangles was found to be less than 5% for the NTH-2 Classical Potential. AB C D Figure 7-4. Silica Nanoclusters. A) Si 12 O 24 B) Si 18 O 36 C) Si 24 O 48 -Cage. D) Si 24 O 48 -Fullerene. Notched Nanorod A defect notch was placed in the 108 atom nanorod by removal of an oxygen atom as shown in Figure 7-5. The MD stress-strain curves for this notched rod were found using TH quantum mechanics and the NTH-2 Classical Potential. The two stress-strain curves again matched up to 9% strain (Figure 7-6). It is interesting to note how the presence of only one small defect can significantly reduce the yield stress of the material and make it prone to fracture. The TH curve for the defect-free rod is plotted in the same

PAGE 112

96 figure to contrast the value of the yield stress. As can be seen there is a reduction of ~60GPa in the yield stress. A composite rod was constructed; however in this case, the QM domain was chosen to consist of 2 silica planes and the intermediate 5 oxygen atoms (Figure 7-5), so that the defect could be located in the QM region. The Y for this composite notched rod matched satisfactorily with that for the TH notched rod for up to 9% strain (Figure 7-6). This composite notched rod is an ideal candidate for studying hydrolysis since the water molecules are known to attack this notched portion which is a stress concentrator and therefore is the weak portion in the rod. Missing Oxygen QM Figure 7-5. Front view of the 107 atom notched nanorod showing one defect (missing oxygen). It is important to note that at higher strains, the stress-strain curve for the composite notched nanorod follows that of the TH nanorod, instead that of the NTH-2 Classical Potential nanorod showing that the large strain behavior is localized to the QM domain. This is exactly what we expected from our multi-scale modeling. The inset in Figure 7-6

PAGE 113

97 shows the zoomed view of the stress-strain curve for the composite rod. We see that the multi-scale modeling studied following our scheme is accurate because the phenomena involving higher strains occur in the QM domain (where it is treated properly with electrons) and computationally efficient as a classical potential is used for the bulk. Figure 7-6. Comparison of stress-strain curves for the notched rod obtained from different methods, showing the significant reduction in the yield stress with the presence of a single defect. Bulk Glass The next application is to a more realistic and bigger system namely a bulk glass containing 2940 atoms with periodic boundary conditions. (This work is done in collaboration with Dr. Krishna Muralidharan). Such a glass is too big to have a quantum calculation on the whole sample; a classical description for the bulk domain is necessary. A glass is formed using the NTH-2 Classical Potential based on the recipe proposed by Huff et al. [1999]. The basic procedure is to melt the crystal structure of a -cristobalite

PAGE 114

98 containing 2940 atoms with periodic boundary conditions, and then to quench it to form a glass. -cristobalite is chosen as the initial structure as its density is very close to experimental glass density (2.2 g/cc). A typical MD simulation for silica glass begins with a NVT simulation at some high temperature, typically 6000 K or higher, for tens of picoseconds. This is done to remove any bias toward of the starting arrangement of the atoms. The temperature of the system is then reduced in steps of 100 2000 K performing NVT simulations at each step for typically tens of picoseconds until ambient temperatures are reached. Following this recipe, a MD simulation was done with the NTH-2 Classical Potential in which the temperature was reduced from 8000K to 300K in steps of 125K/ps. A cubic cell of 34.43 was used for periodic boundary conditions, and the temperature was controlled by a Nose-Hoover thermostat. Figure 7-7 shows the structure of the glass that emerged. Figure 7-7. Structure of glass obtained from the NTH-2 Classical Potential.

PAGE 115

99 To test the accuracy of the NTH-2 classical potential for this system, properties such as the radial distribution function (RDF), bond angle distribution (BAD), density etc., obtained using the NTH-2 Classical Potential were compared to those from experimental results. In fact the properties obtained using the phenomenological BKS potential also fit the experimental results well. Radial distribution functions (RDF) are a useful way to describe the structure of a system. They are defined as the probability of finding two particles separated by a distance between r and r r These can be determined easily from MD simulation results since the exact location of all the particles is known. The RDF is given by )4/()(2)(2drrNVrnrgijij (7-1) where V is the volume of the cell, n ij is the number of pairs of type ij found between r and r+dr and N ij is the possible ij pairs in accumulation. The RDF can be measured experimentally using X-ray or neutron diffraction. The arrangement of atoms gives a characteristic diffraction pattern, which then can be analyzed to calculate an experimental distribution function and compared with that obtained from the simulation. The RDFs obtained from the glass prepared with the NTH-2 classical potential for each pair of Si-O, Si-Si and O-O are shown below in Figure 7-8. The bond angle distribution (BAD) is calculated by finding the first-neighbor of each atom and then computing the bond angle between a central atom and its neighbors defined as )(cos1ikijikijikijikijijkrrzzyyxx (7-2)

PAGE 116

100 where i denotes the apex atom, j and k denote the end atoms. The BAD calculated using this definition for
PAGE 117

101 Table 7-3. Comparison of different properties of the glass obtained from the NTH-2 Classical Potential to those of the BKS glass. (The numbers in the RDF and BAD correspond to the peaks in their corresponding distribution functions) NTH-2 Classical glass Experimental/ BKS glass Si-O RDF 1.64 () 1.62 () O-O RDF 2.64 () 2.63 () Si-Si RDF 3.26 () 3.08 () Si-O-Si BAD 139 o -160 o 140 o -155 o O-Si-O BAD 105 o 109 o Density (g/cc) 2.4 2.29 2-membered ring 4 2 1-coordinated O 6 18 3-coordinated O 12 10 3-coordinated Si 0 2 5-coordinated Si 6 2 Table 7-3 shows the comparison of some of the thermodynamic properties obtained using the NTH-2 Classical Potential with the BKS potential or experimental glass [Wright, 1994]. The density, the RDFs and the BADs for the NTH-2 Classical Potential glass match the experimental results within a precision of 2%. We have also studied the presence of two-membered rings in each case as these are an indicator of the stability of the glass structure. The two-membered rings have very high internal strains, so a stable glass is marked by a lower number of 2-membered rings and as we can see from Table 7-3, the glass obtained using NTH-2 Classical Potential satisfies this condition. It is also interesting to look at the number of non-regular coordinated atoms. Non-regular here refers to Si and O atoms which are not in their normal coordination (the normal coordination are 4 and 2 for Si and O atom respectively). The formation of two-membered rings involves formation of these non-regular coordinated species, so it is an important feature to study the presence of such atoms as well. We observe that the

PAGE 118

102 number of such non-regular coordinated atoms for the NTH-2 Classical Potential is also close to those found in the experimental glass. We see from Table 7-3 that almost all the properties of this NTH-2 Classical Potential glass agree well with those of the experimental glass. A real glass consists of ring structures. Sixand five-membered rings are the dominant species but there are also some smaller ones such as four-, three-, and two-membered rings on the surface. The ring analysis for the NTH-2 Classical Potential is compared to the experimental results and again a good agreement between the two results is observed. Next, a stress-strain study for the NTH-2 Classical Potential Glass was performed, and compared with those of BKS potential. In this stress-strain MD simulation, a uniform uniaxial strain at a rate of 0.1/ps was applied throughout the sample. BKSNTH-2Region IRegion IIRegion IIIRegion IV Figure 7-10. Stress-strain curves for NTH-2 and BKS glass.

PAGE 119

103 There are four characteristic regions in the stress-strain curve (Figure 7-10). Region I corresponds to the elastic deformation of Si-O bonds, region II corresponds to the plastic deformation, region III the failure and surface reconstruction, and region IV separation of the material. Muralidharan [2004] gives the detailed results for BKS glass. Here we see the NTH-2 Classical Potential glass matches with BKS until ~ 10% strain. At larger strains it is no longer clear that a classical pair potential can represent the correct physics of the problem, and the reliability of both the BKS and the NTH-2 potential is uncertain. Still it is interesting to note the qualitative features in regions II IV. The NTH-2 is much weaker than the BKS glass in region II for two reasons: (i) the charges on Si and O atoms in BKS potential are much higher than the NTH-2 potential thereby making the Coulombic interaction stronger and restraining the stretch in the Si-O bond. (ii) the Si-O interaction in the BKS potential becomes even stronger to compensate for the fact that there is no Si-Si interaction. In regions III and IV the system has failed in quite different ways and there is probably no simple relationship between the systems formed with the two different potentials. Unfortunately there are no experimental data or full quantum calculations available for the stress-strain curve for this glass to discern the correct classical potential. Instead, it should be quite instructive to apply our multi-scale modeling to see the form of the stress-strain curve in regions II-IV when regions of highest strain concentration are treated quantum mechanically. The status of this investigation is discussed below. We note from Figure 7-10 that the yield stress (~ 20 GPa for the BKS potential) is an order of magnitude smaller than the artificially high yield stress of the nanorod (~190GPa). Therefore, although our multi-scale scheme fails for the composite nanorod

PAGE 120

104 after 9% strain, nevertheless it is expected to work for bulk glass which should fail at a lower strain (as we predicted in chapter 5). The success of the NTH-2 Classical Potential, originally trained on the nanorod, in predicting the correct results for the glass and the nanoring is somewhat surprising since the local structure of the glass and nanoring are quite different from that of the rod. This suggests a new phenomenology whereby ab initio data for a specific system and specific properties are used to construct a classical potential, followed by an empirical exploration of its range of validity elsewhere. Results Obtained from New-DFT-Classical-Potential In this section we present the structure analysis using the New-DFT-Classical Potential. The 10-Membered Silica Nanoring For comparing the structure of the 10-nanoring obtained from New-DFT-Classical-Potential, the structure obtained using DFT on the whole sample is chosen to be the reference system. Table 7-4 shows the comparison of the structure. Clearly there is relatively large percentage error for the structure of the 10-membered ring than compared to the nanorod (Table 64). This could be because the New-DFT-Classical Potential is trained using expansion data only until 4% strain. This nanoring is an internally very strained molecule, hence training with higher strain data might give better performance (the nanoring is unlike the stable nanorod in which training of the New-DFT-Classical Potential with small strain data was sufficient to give a reasonable structure). It is to be noted that the NTH-2 Classical Potential which was trained on data through 10% strain gives more precise structure (percentage error is always less than 5%) for both the nanorod as well as the nanoring (Table 7-1) compared

PAGE 121

105 to the New-DFT-Classical Potential. Thus we can conclude that the incorporation of higher strain data, in the search algorithm for finding a classical potential, might be a significant factor Table 7-4. Comparison of structure of the nanoring obtained from DFT and New-DFTClassical Potential Bond lengths () and angles DFT New DFT Classical Potential % error Si-O (outer) 1.69 1.64 2.96 Si-O (inner) 1.67 1.63 2.39 O-O (outer) 2.35 2.31 1.7 O-O (inner) 2.36 2.31 1.7 Si-Si 2.31 2.42 4.76 < Si-O-Si (outer) 86.29 95.55 10.7 < Si-O-Si (inner) 89.3 96.5 8.06
PAGE 122

106 The numbers in the RDF and BAD correspond to the peaks in their corresponding distribution functions A glass having the same number of atoms (2940) was simulated in an identical manner as before but using the New-DFT-Classical Potential parameters. Table 7-5 compares the properties of this glass with respect to BKS Glass. We see a good agreement in the properties of the two glasses Conclusions This chapter shows the wider applicability of the NTH-2 potential to two entirely different silica systems namely a 10-membered nanoring and an amorphous glass. The structure predicted for these two systems matches with those predicted from TH within 4% accuracy. A notched rod containing one less oxygen atom was studied to illustrate reduction in the value of yield stress with the presence of defects. We saw that the elastic properties of the composite notched rod agreed with those obtained from TH or pure NTH-2 classical potential for strains up to 9%. This result was found to be true in the case of the 10-membered nanoring as well. The structure analysis for the 10-membered ring and the glass was repeated with the New-DFT-Classical Potential as well. We found that though the properties of the glass created from this potential match quite well with the BKS glass, but the structure in the nanoring shows a relatively higher percentage error. Recalling that the NTH-2 Potential gave better results for the 10-membered nanoring as well, we infer that training the potential with higher strain data is necessary to find a better potential.

PAGE 123

CHAPTER 8 COMPARISON OF THE DIFFERENT QUANTUM MECHANICAL METHODS Introduction The problem posed in this thesis is how to construct a multi-scale model for a solid whose quantum mechanical description was given. In practice, many different choices are made for that description governed by the specific solid considered and practical considerations of computational time. In chapter 6, we showed the independence of the multi-scale modeling scheme for two different choices of underlying quantum mechanical (QM) methods: TH and DFT. Ultimately, of course, the accuracy of prediction from the modeling scheme is linked to the accuracy of the QM description. The assessment of QM descriptions can be made for special cases, such as specific small molecules, but there is no clear guidance for larger systems such as those of interest here. The two different benchmarks for tests of our modeling scheme have provided the results needed to compare the two quantum methods as a side product of this study. In this chapter some comparisons of the results obtained here for the cases of TH and DFT methods are provided, along with some comments. In summary, although the modeling scheme was accurate for both TH and DFT QM descriptions, the quantum solid being studied in each case turned out to be quite different. This raises an important question to be faced in the general approach to multi-scale modeling. An accurate multi-scale modeling scheme is ultimately conditioned by the accuracy of the QM description. 107

PAGE 124

108 The Nanorod In previous sections the 108 atom SiO 2 nanorod has been described using both TH and DFT quantum mechanics. The accuracy of each is uncertain at the outset, but there are reasons to be optimistic about both. Thus a priori expectations were that the equilibrium structure and elastic properties would be quite close. This was found to be the case for structure but not for the elastic properties, as is described quantitatively in the next section. Comparison of Structure The structure of the nanorod obtained by applying TH on the whole rod is compared to that obtained from DFT on the full rod in Table 8-1. Table 8-1. Comparison of structure of the nanorod obtained from DFT and TH quantum methods Bond lengths () and angles DFT TH In planar silica rings Si-O 1.62 1.64 Si-O-Si 158.02 170.06 Between planes and intermediate O O-Si-O 107.15 103.82 At the end-caps Si-O 1.64 1.71 Si-O-Si 92.57 102.03 Approximate length 15.79 16.5 Approximate diameter 6.44 6.55 We see the two rods have noticeable differences in the structure: the TH rod is about 4% longer and 2% wider than the DFT rod; the Si-O-Si bond angle in-plane is about 7% larger for the TH rod. However, these relatively small differences suggest that the two methods are describing the same electronic bonding responsible for the stable

PAGE 125

109 structure obtained in each case. Indirectly, this is provisional support for the accuracy of each. This deduction is proven false by a further comparison of the elastic properties. Elastic Properties While the differences in the structure are not large, there is a striking difference between the elastic properties of the two rods. As mentioned before, the DFT is computationally too intensive to do a MD simulation hence we have compared only equilibrium and selected adiabatic strain configurations for DFT and TH. The stress was calculated in the same manner as explained in chapter 4. Figure 8-1 shows the elastic properties of the TH and DFT nanorods. Comparing the initial slopes of the two curves, we note the TH rod is about three times stronger than the DFT rod. (GPa) Figure 8-1. Comparison of elastic properties for DFT and TH nanorods. The equilibrium structures for the nanorod are local minima in the electronic potential energy surface as a function of the atomic coordinates. Evidently both TH and DFT are locating similar minima, and hence similar structures. However, the small

PAGE 126

110 quadratic deviations on this surface away from the minimum are quite different. This is a small strain effect so it is showing a significant difference in the way the electrons are being described in the two methods. Since the difference is so large, at least one of the two descriptions (and possibly both) is clearly wrong. Comparison of Different Types of DFT In this section, we compare the results obtained using two different types of DFT in the 10-membered nanoring. In the original work on the 10-membered nanoring, Bromley et al., [2003] had used DFT method with a Gaussian basis set and the Becke-type three parameter hybrid exchange-correlation functional B3LYP [Becke, 1993]. In Table 8-2, we compare the bond lengths and bond angles of the 10-membered nanoring obtained using Bromleys method to those obtained from the BOMD method (plane wave basis set, PBE exchange correlation functional). Table 8-2. Comparison of the 10-membered nanoring using two different types of DFT methods Bond lengths () and angles DFT (PBE plane waves) DFT (B3LYP, Gaussian 98) Si-O (outer) 1.69 1.69 Si-O (inner) 1.67 1.67 O-O 2.36 2.79 Si-Si 2.31 2.35 < Si-O-Si (outer) 86.29 87.9 < Si-O-Si (inner) 89.3 89.6
PAGE 127

111 limitation; that is both of these exchange correlation functionals are local. We can only conclude that the differences between the TH and DFT cannot be accounted for based solely on the basis or exchange correlation functional chosen for the DFT calculations. Bulk Glass We also consider the structural and elastic properties of the bulk glass. However since a full quantum calculation cannot be done on bulk glass, we compare instead the structure (Table 8-3) elastic properties (Figure 8-2) obtained using the NTH-2 Classical Potential and New-DFT-Classical Potential. Table 8-3. Comparison of structure of bulk glass using two different classical potentials NTH-2 Classical glass New-DFT-Classical-glass Si-O RDF 1.64 1.6 O-O RDF 2.64 2.62 Si-Si RDF 3.26 3.16 Si-O-Si BAD 139 o -160 o 147 o -160 o O-Si-O BAD 105 o 109 o Density (g/cc) 2.4 2.24 2-membered ring 4 4 1-coordinated O 6 16 3-coordinated O 12 5 3-coordinated Si 0 6 5-coordinated Si 6 5 Table 8-3 shows that the structure of the bulk glass obtaining using NTH-2 and New-DFT Classical Potentials are quite similar. However, from Figure 8-2 one sees that in the elastic region (region I) the glass obtained from the NTH-2 potential is about 3 times stronger than the DFT-classical glass. Of course, this is expected based on the difference between these potentials when applied to the nanorod.

PAGE 128

112 Strain (rate @ 0.1/ps) Figure 8-2. Comparison of elastic properties of the two glasses in the elastic region. Discussion The question now arises as to which of the two methods TH or DFT is more accurate as only one of the values for the Youngs modulus can be correct. The answer to this question is not straightforward because of the lack of experimental or more sophisticated (e.g., coupled cluster) theoretical data on the systems modeled. From a formal point of view too, one cannot easily justify the preference of using one theory over the other as there are advantages and disadvantages associated with both methods. For example, the advantage of the TH method is that it can predict both ionic and radical modes of Si-O bond dissociation (chapter 3) while DFT can predict only the radical dissociation mode. It is known that both of these bond dissociation modes are essential to predict fracture properly. However for the TH method, the parameters optimized (using high quality CC theory) on a small pyrosilicic acid molecule (chapter 3) were assumed to

PAGE 129

113 be transferable (because of its NDDO form) to the bigger nanorod; but there are no data available from the CC theory on this nanorod that can ensure this justification. Moreover the TH method was parameterized only on a single molecule (pyrosilicic acid) and its accuracy has not been tested on other molecules. In the DFT method, no such transferability of parameters was used but then DFT also uses an uncontrolled approximation namely the exchange-correlation functional (PBE in our case). These observations highlight the need for more detailed characterization of the accuracy and domain of applicability of candidates for the QM description in multi-scale modeling. In practice, the size of the quantum domain chosen is determined both by the physical mechanisms deemed important and the practical limitations of the computation. As mentioned several times here, the TH method is computationally much more efficient compared to the DFT method. However, there is no point in performing an efficient calculation whose accuracy is limited. Conversely, a more intensive DFT calculation can not be justified until its accuracy is assured. In our opinion, modeling combined with experiment could possibly reveal the accuracy between the two QM methods for their applications to multi-scale modeling. The proposal is to compare the Youngs modulus for a realistic material like glass. Multi-scale modeling of glass would be required because of its large size. The two QM methods and their corresponding classical potentials would then be used to predict the yield stress of glass. The method with the best agreement with the experimental value of the yield stress should be the preferable QM method for calculating the elastic properties. However, this is easier said than done as there are additional uncertainties associated with assuring a faithful correspondence of the theory to experiment. For example, one has to

PAGE 130

114 make sure that the quenching process in the preparation and measurement of the stress for the simulated glass and the experimental glass are done in an identical manner. Only then would the experiment data be meaningful in discerning the accurate QM method (in terms of predicting the correct elastic properties). At this point, the test of accuracy of these two QM methods remains as an open question but it is not an issue of concern for this dissertation. One of the primary conclusion from the dissertation is that given a particular QM method, our scheme can construct an equivalent composite quantum-classical system that would be faithful (have the same physical properties) to the underlying choice for the QM method for states near equilibrium.

PAGE 131

CHAPTER 9 SUMMARY AND FUTURE RECOMMENDATIONS Conclusions The technique of multi-scale modeling has emerged as a powerful tool in the last decade for studying and predicting complex phenomena. However there remain a few issues yet to be solved. These include the problem of consistent embedding of a QM domain in its CM bulk, and the handshaking region at the boundary of QM/CM domains where researchers have artificially matched the difference in the forces between the QM and CM domain by taking some sort of average of the two forces. Moreover, the question of what variables are most appropriate for control in multi-scale modeling still remains fuzzy. We have addressed these problems in this project. This dissertation provides a prescription for multi-scale modeling and a first example of its implementation and critical assessment that is internally self-consistent and based at all scales on the fundamental underlying quantum mechanics. We have developed the techniques for a consistent embedding of a QM domain inside its CM environment. This required the incorporation of information about the state of the CM environment into the QM domain by approximating its effects in terms of pseudo-atoms and dipoles. For the examples considered, SiO 2 structures, the pseudo-atoms are fluorine atoms with electronic and ionic parameters modified to account for the short-range exchange interactions between the QM and the CM regions. It has been shown that these pseudo-atoms are more accurate to saturate the dangling bonds resulting 115

PAGE 132

116 from the isolation of the QM domain than the more commonly used hydrogen (LA) atoms. Our characterization of the embedding was based on the predicted quality of the exact charge densities and forces (which implies correct geometry as well) in the QM domain as the fundamental variables of interest, rather than the more commonly used tests for fidelity of structure and proton affinities. The results for the charge densities and the forces have been verified by comparing them with the data obtained from direct quantum calculations. The choice of the charge density as a property to assess the quality of such QM/CM partitioning has also been shown to be a more physically reasonable test. We have confronted the problem of the handshaking region by constructing a classical pair potential that predicts the same structure and Youngs modulus for small strains as given by the quantum method. This potential is built using ab initio data on the equilibrium and weakly strained configurations calculated from the quantum description, rather than the more usual approach of fitting a wide range of empirical data. In fact, we have proposed that each development of a multi-scale algorithm should include the construction of its own classical potential. The functional form of the classical potential can be chosen to be the same as any reasonable phenomenological form available but the parameters should be determined using data from the underlying quantum calculations only. In our case we have determined the potential parameters from a genetic algorithm using the quantum force data followed by a spatial scaling. Finally using our embedding tools and the derived new classical potential, we constructed a composite system that is indistinguishable from the one obtained by applying quantum mechanics on the whole system for states near equilibrium. This procedure guarantees the multi-scale phenomena to be: (i) accurate by localizing the large

PAGE 133

117 strain behavior to the QM domain (where it is treated electronically) and (ii) computationally efficient with the use of new classical potential developed for the rest of the system. Although the study has been limited to silica based materials, this scheme for multi-scale modeling should provide the guidance for extension to other materials as well. In addition, this multi-scale modeling method is shown to be independent of the choice of quantum mechanical method used by applying it to the quite different methods of TH and DFT quantum mechanics. Our project is a part of the National Science Foundation ITR program [Bartlett et al., 2004] which focuses on the science, algorithm, and software necessary for the predictive simulations of the material behavior in which chemo-mechanical interactions are critical like crack propagation under the presence of water. (Predictive simulation means accurate results without parameterization using experimental data). Such chemo-mechanical phenomena depend upon the systematic identification of the chemically active region that does not require intervention by the researcher. A wavelet transform analysis is under construction to identify the candidate QM regions [Muralidharan, 2005], An efficient quantum chemistry code is developed [Taylor, 2003] that can handle more than a hundred atoms fast enough to update forces in an MD code. In this way, the simulations are made chemically realistic by obtaining the forces driving the MD from highly accurate electronic structure methods. This dissertation is relevant for the last steps of this ITR project, which after locating the position of the candidate quantum region embeds it in the classical environment with techniques discussed above for studying the multi-scale phenomena.

PAGE 134

118 Recommendations for Future Work We have shown the success of our multi-scale scheme to be independent of the choice of the quantum method used. However the issues of the relative accuracy of the two quantum methods (TH and DFT) used and the reasons for obtaining a nanorod having different structures and elastic properties using the two different methods have not been addressed in this work. It is therefore necessary to calibrate the two methods on the same molecule. The parameterization of the TH method has been done on the pyrosilicic acid for two different reaction pathways of Si-O bond dissociation as explained in chapter 3. A similar calculation has to be done for the same reaction pathways using DFT to calibrate the two methods. For finding the elastic properties of the nanorod using DFT, only adiabatic stresses have been used due to our limited computer resources. A future project could be developing efficient methods using DFT that would allow us to execute the MD simulation in a feasible amount of computational time. While training the electronic properties of the pseudo-atom the Coulombic effect of part of the molecule beyond the right hand side of the F atom (Figure 3-1) was ignored. It would be interesting to approximate this part by a lower order multipole followed by the retraining of the parameters for the pseudo-atom. Moreover the pseudo-atoms have been trained only at regions very close to equilibrium. It is therefore needed to retrain these parameters for higher strains in the Si-O bond to explore their domain of applicability. Finally the techniques for multi-scale modeling developed in this dissertation need to be extended to (i) metals, where the whole approach for partitioning would be entirely different as metals have non-localized electrons [Kroger, 2003], (ii) microstructure and dislocation such as point and line defects, crystallographic texture and grain boundaries,

PAGE 135

119 phase interfaces in solids, (iii) predict chemical kinetics in different reaction pathways and (iv) studying biological systems, integrating quantum to continuum models for sub-cellular process.

PAGE 136

APPENDIX A NEGLECT OF DIATOMIC DIFFERENTIAL OVERLAP (NDDO) MODEL In the restricted Hartree-Fock theory, each electron is assumed to move in an average field of all other electrons. For each electron, the Schrdinger equation can be written as iiiif (A-1) if is the one-electron Fock operator defined as )}1()1({)1()1(1jjcoreNjifKJ (A-2) The core Hamiltonian operator MAAArZ1121core21H (A-3) This operator corresponds to the motion of a single electron moving in the field of the bare nuclei. The Coulomb operator, is jJ )2(1)2(122jjrdJ (A-4) This operator corresponds to the average potential due to an electron in j ( 212drdrd is the volume element). The exchange operator is )1(jK )1()2(1)2(d21(1))1(122jjijrK (A-5) 120

PAGE 137

121 The orbitals s are written as a linear combination of single-electron functions Kiic1 (A-6) The functions are called a basis set and may have the form of atomic orbitals. The Fock operator in terms of these basis functions becomes KKiiiiccf11)1()1()1( (A-7) Premultiplying each side by )1( and integrating we get )1()1()1()1()1(2112dcfdcKiiKii (A-8) We introduce a charge density matrix 2/12NiiiccP The Fock matrix then becomes KKcoreiPHfdF112)]|(21)|[()1()1()1( (A-9) where we have used short-hand notation for the integrals )2()2(1)1()1()|(1221rdd )1()1(1d is the overlap integral between the basis functions and written as S Equation (A-9) is the standard form of the Fock matrix called the Roothaan-Hall Form. In ab initio calculations all elements of Fock matrix are calculated using equation (A-9) irrespective of whether the basis functions , and are on the same atoms,

PAGE 138

122 on atoms that are bonded, or not. The Fock matrix can be partitioned into three elements (the diagonal elements); ( F F and are on the same atom) and ( F and are on the different atoms). In Hartree-Fock SCF calculations most of the time is spent in calculating the integrals. Semi-empirical methods reduce this computational cost by approximating some of the integrals. A feature common to all semi-empirical methods is that the overlap matrix S is set to the identity matrix I. The Neglect of Diatomic Differential Overlap (NDDO) theory further neglects the differential overlap between the atomic orbitals on different atoms. Thus all of the two-electron, two-center integrals of the form )|( where and are on the same atom (say A) and and are also on the same atom (say B, with A B) are retained. In this notation the NDDO Fock matrix elements becomes onAonAABonBonBcorePPHF)|()]|(21)|[( (A-10) where BBcoreVUH, onAABonBonBonAcorePPPHF)|()]|(21)|([ (A-11) and both on A onBonAcorePHF)|(21 on A and on B (A-12) The following terms that appear in the NDDO Fock matrix have been parameterized for our application using CCSD data as described in chapter 3:

PAGE 139

123 (i) One-center, one-electron energies which represent the sum of the kinetic energy of an electron in orbital U at atom A and its potential energy due to the attraction by the core of atom A. (ii) One-center, two-electron repulsion integrals, i.e., Coulomb integrals g )|( and exchange integrals h )|(. (iii) Two-center one electron core resonance integrals (iv) Two-center one electron attractions between an electron in the distribution BV, at atom A and core of atom B. (v) Two-center two-electron repulsion integrals )|(

PAGE 140

APPENDIX B GENETIC ALGORITHM The Genetic algorithms (GAs) are stochastic processes based on Darwin's principle of evolution where the fittest individuals survive and are allowed to reproduce through exchange of chromosomes. In the present context a GA is a computational manifestation of the same paradigm of survival of the fittest' where the fittest parameter sets, i.e. parameters with least error function (e.g., L in equation 4.7) survive and produce new parameter sets. This is done by creating a 'population' of 'individuals', defining rules to measure the fitness of individuals, and defining operators for reproduction' (cross-over) and 'mutation' and for evolving the population until the process converges to the fittest population. The selection procedure is such that the probability of an individual being selected for breeding is proportional to that of the individual's fitness. This is done using the Roulette wheel algorithm [Goldberg, 1989]. The process of reproduction and mutation is simulated in a following way. Suppose that two parents (i.e., two (x,y) pairs) have been selected for breeding. The first step is to construct an encoding of each parent, in the form of a string-like structure or chromosome. This is accomplished by fragmenting each parameter into simple decimal integers, and then splicing the resulting integers strings one behind the other. The number of digits, n d to be retained in the encoding must be specified externally. For example, with n d = 8 one would obtain (x,y) = (0.123456789, 0.987654321) 1234567898765432 The length of the chromosome is then n d times the number of parameters, n, defining the solution. In this example n=2 so the chromosome has length 16. Each digit can be thought of as a gene 124

PAGE 141

125 occupying a chromosomal site for which there exist 10 possible gene values. Decoding is simply the reverse process: 1234567898765432 (0.12345678, 0.98765432) = (x,y). The first step of the breeding process proper is the application of the crossover operation to the pair of parent chromosomes defined as follows. Consider now two parent chromosomes produced by the encoding process described in the preceding section: 1234567898765432 (Parent #1) and 7654321023456789 (Parent #2). One of the 16 sites is selected at random, both chromosomes are cut at that site, and the fragment right of and including the cutting site is interchanged. CUT: 123456789 8765432 765432102 3456789 SWAP: 123456789 3456789 765432102 8765432 SPLICE: 1234567893456789 Offspring #1 7654321028765432 Offspring #2 The results of this process are two new offspring chromosomes, each containing intact chunks of chromosomal material from each parent. In the GA literature this is called a one-point uniform crossover. The second step of the breeding process is the application of the mutation operator to each offspring chromosome. For each gene of each offspring chromosome, a random number between 0 and 1 is generated, and if this number is smaller than a preset mutation probability P m then the gene value is randomly changed to any other legal value. For example the following is a mutant at site 2:

PAGE 142

126 7654321028765432 754321028765432 7154321028765432 In the GA literature this is called one-point uniform mutation; one point because there is only one gene affected at a time, uniform because each gene is equally likely to be subjected to mutation, independently of the site it occupies along the chromosome. The effects of crossover and mutation on the decoded version of a parameter set can be drastic (if one of the leading digits is affected) or quite imperceptible (if one of the least significant digits is affected). The result is that the breeding process can cause large jumps in parameter space, as well as slight displacements; the resulting search algorithm can explore the space efficiently. In this example we would have gone from two parents: (x,y) = (0.12345678,0.98765432) and (x,y) = (0.76543210,0.23456789) to two offsprings: (x,y) = (0.12345678,0.93456789) and (x,y) = (0.71543210,0.28765432). Clearly the offsprings have taken sizable jumps through the parameter space, as compared to their parents. The crossover operation is what distinguishes genetic algorithms from other adaptive stochastic techniques. Consider again one of the 16-gene chromosomes defined previously. Replace now each gene by one of three possible symbols (-, +, or X) according to the following coding: if the gene value is good, in the sense of giving its bearer (on average) above average fitness, replace the gene by the symbol +; if its bearer has below average fitness, replace it by; if the gene value does not seem to matter, replace it by X (this would usually be the case for gene encoding least significant digits early in an evolutionary run). For example, two possible results of this process for a 16-gene chromosome: ++--XXXX---+XXXX Parent #1

PAGE 143

127 ---++XXX+++XXXXX Parent #2 Crossover occurring at site #9 (without any subsequent mutation) would produce ++--XXXX+++XXXXX Offspring #1 ---++XXX---+XXXX Offspring #2 Clearly, Offspring #1 is doing a lot better than any of its parents while Offspring #2 will be mostly below average fitness, and so will have a low probability of being selected for breeding at the next generational iteration. Observe that crossover can combine advantageous chunks of good chromosomal material in a single offspring. With natural selection operating, these chunks will end up (on average) being copied into the next generation rather frequently, as the lucky offspring finds itself above average. This is the essence of genetic algorithm theory. For the problem considered here, the GA was initialized with 100 trial parameters and evolved through 500 generations.

PAGE 144

APPENDIX C SCALING METHOD The parameters of Table 4-1 obtained from the GA represent stable equilibrium structures, but these structures still differ somewhat from those of the underlying quantum TH calculation. To improve the potential further, new potential parameters are determined from these by mapping the GA structures onto the correct TH structure through an affine scaling of the coordinates. Let p ij = (q i, q j a ij b ij c ij ) denote one set of parameters from Table 4-1 and let r ij be the separation for the pair i,j at equilibrium. Also, let be the corresponding equilibrium separation from the TH. The objective is to find a new set of parameters p 'ijr ij for the pair potential such that the new equilibrium configuration is at. 'ijr As a first step require that the pair forces under change of parameters and rescaling of the separation are proportional to the original forces ),()','(ijijijijijijprprFF (C-1) The sum of the forces over i for each j is (nearly) zero as a result of the GA search. This condition assures that the new parameters will yield a new equilibrium at in agreement with the TH calculation. This is strictly true only if the r ),(ijijijprF 'ijr ij and are related by the same scale transformation for every pair, which is observed to be approximately the case. Next, the stability of this new equilibrium condition is assured by the second constraint 'ijr 128

PAGE 145

129 ),()','('ijijijijijijijijprdrdprdrdFF (C-2) For the chosen pair potential 6exp)(ijijijijijijjiijrcrbarqqrV (C-3) the magnitude of the force and its derivative are 726)exp()(ijijijijijijijjiijijijrcrbbarqqrVrF (C-4) 82342)exp(2ijijijijijijijjiijijrcrbbarqqrF (C-5) Equating the force derivatives for the original and rescaled potentials gives jiijijjiqqrrqq3''' )exp()''exp(''22ijijijijijijijijrbbarbba ijijijijcrrc8'' (C-6) while the condition for the forces to be proportional gives in addition '/'ijijijijrrbb

PAGE 146

APPENDIX D CHARGE TRANSFER POTENTIAL The potentials used in chapter 4 assume that the effective charges associated with Si and O atoms are fixed and do not vary as function of their immediate environment. Alavi [1992] proposed that there is some charge transfer between a positive and a negative ion leading to a formation of covalent bond between the two ions thereby making the material stronger. According to his method, if be the coordinates of particles i of type A carrying a positive charge q Air Bjr A and particles j of type B carrying a negative charge q B respectively, then the net charge on each atom consists of the bare charge plus the amount of charge transferred to or from particles of other type given by the formula (D-1a) jijABAiirfqqq)( (D-1b) iijABBjjrfqqq)( The summation runs over the nearest neighbors of opposite particle type. In these equations, is the amount of maximum charge transferred and is an empirical charge transfer function given by q )(ijABrf (D-2) ]}/)tanh[(1{*5.0)(ABijijABRrrf where is the range over which the charge transfer takes place and is the separation at which the charge transferred is ABR 2/q 130

PAGE 147

131 Charge-Transfer Potential for Old Transfer Hamiltonian Classical Potential #1 The parameters used for charge-transfer classical potential, to stiffen the old TH classical potential #1, are chosen to be q Si = 2.84, q o = 1.42, 15.0 q and oABAR0.2 1.0 such that in the fully coordinated state Si and O atoms have a charge of 2.24 and -1.12 respectively and in the zero-coordinated state Si has a charge of 2.84 and O has a charge of -1.42. The amount of charge transferred from Si to O as a function of Si-O distance typically looks like Figure D-1. Distance (Ang) Charge transferred Figure D-1. Charge transferred from Si to O as a function of distance. We used this formulation in conjunction with our previously derived old-TH-classical potential #1 (chapter 4) to obtain the new charge transfer potential. Figure D-2 shows the change in the stress-strain curve with the introduction of the charge transfer term.The nanorod becomes more rigid with the introduction of the charge-transfer term in the potential as we wanted. This is because the magnitude of charge on Si and O increases with increasing Si-O distance, making the Coulombic interaction stronger, thereby resisting the straining of Si-O bonds.

PAGE 148

132 Figure D-2. Increase in the stiffness with introduction of the charge transfer term. Although this charge potential serves our purpose of making the potential stiffer, it makes the functional expression for the potential and hence the MD code somewhat complicated to use. Charge on Silicon 0.052.24 2.84 Figure D-3. Increase in the charge of Si with strain. Figure D-3 shows the increase in the amount of charge on Si as a function of the strain.

PAGE 149

133 Charge-Transfer Classical Potential for New-Density Functional Theory The same charge-transfer method could be applied to stiffen the New-DFT-Classical Potential (chapter 6) as well. The parameters for charge transfer potential in this case are chosen to be, and 15.0q oABAR9.1 8.0 Figure D-4 shows in the increase in the Youngs modulus with the introduction of charge transfer for New DFT Classical Potential. Figure D-4. Increase in the stress with the charge transfer DFT potential.

PAGE 150

LIST OF REFERENCES Alavi, A., Alvarez, L.J., Elliot S.R., and McDonald, I.R., 1992, Philos. Mag. B., 65(3), 489. Al-Derzi, A.R., Cory, M.G., Runge, K., and Trickey, S.B., 2004, J. Phys. Chem. A., 108, 11679 Allen, M.P., and Tildesley, D.J., 1987, Computer Simulation of Liquids, (Clarendon Press). Amara, P., Volbeda, A., Fontecilla-Camps, J.C., and Field, M.J., 1999, J. Am. Chem. Soc., 121, 4468. Anderson, H. C., 1983, J. Chem. Phys., 72, 2389. Antes, I., and Thiel, W., 1999, J. Phys. Chem. A., 103, 9290. Aqvist, J., and Warshel, A., 1993, Chem. Rev., 93, 2523. Assfeld, X., and Rivail, J.L., 1996, Chem. Phys. Lett., 263, 100. Bartlett, R.J., 1995, D.R Yarkony (Ed.), Modern Electronic Structure Theory Part II, World Scientific, 1047. Bartlett, R.J., Trickey, S.B., Cheng, H-P, Harris, F.E and Dufty, J.W., 2004 http://www.qtp.ufl/kdi/ (Last accessed Jan 2005). Barnett, R.N., and Landman, U., 1993, Phys. Rev. B., 48, 2081. Becke, A.D., 1993, J. Chem. Phys., 98, 5648. Biswas, R., and Hamann, D.R., 1985, Phys. Rev. Lett., 55, 2001. Brenner, D.W., 2000, Phys. Stat. Sol. (b), 217, 23. Bromley, S.T., Zwijnenburg, M.A. and Th. Maschmeyer, 2003, Phys. Rev. Lett., 90, 035502. Broughton, J. Q., Abraham, F. F., Bernstein, N., and Kaxiras, E., 1999, Phys. Rev. B., 60, 2391. Chabana, G.M., and Gerber, R.B., 2001, J. Chem. Phys., 115(3), 1340. 134

PAGE 151

135 Charbonneau, P., 1995 Astrophysical J. Supplement series, 101, 309. Clementi, E., 1988, Philos. Trans. R. Soc. London, A 326, 445. Crespo, A., Scherlis, D.A., Mart, M.A., Ordejon, P., Roitberg, A.E., and Estrin, D.A., 2003, J. Phys. Chem. B., 107, 13728. Cui, Q., and Karplus, M., 2002, J. Phys. Chem. B., 106, 1768. Dapprish, S., Kormaromi, I., Byun, K.S., Morokuma, K., and Frish, M.J., 1999, J. Mol. Struct: THEOCHEM, 1, 461. Das, D., Eurenius, K.P., Billings, E.M, Chatfield, D.C., Hodos, M., and Brooks, B.R., 2002, J. Chem. Phys., 117(23), 10534. De Vries, A.H., and. van Duijnen, P.Th, 1996, Internat. J. Quant. Chem., 57, 1067. Dewar, M.J.S., and Thiel, W., 1976, J. Am. Chem. Soc., 4899. DiLabio, G.A., Hurley, M.M., and Christiansen, P.A., 2002, J. Chem. Phys., 116, 9578. Dreizler, R.M., and Gross, E K.U., 1990, Density Functional Theory (Springer-Verlag). Du, M-H., and Cheng, H-P., 2003, J. Chem. Phys., 119, 6418. Du, M-H., Kolchin, A., and Cheng, H-P., 2004, J. Chem. Phys., 120, 1044. Eichinger, M., Tavan, P., Hutter, J., and Parrinello, M., 1999, J. Chem. Phys., 110(21), 10452. Flocke, N., Zhu, W., and Trickey, S.B., 2005, J. Phys. Chem. B., 109, 468. Foulkes, W.M.C., Mitas, L., Needs, R.J., and Rajagopal, G., 2001, Rev. Mod. Phys.,73, 33. Gao, J., 1996, Acc. Chem. Res., 29, 298. Gao, J., Amara, P., Alhambra, C., and Field, M.J., 1998, J. Phys. Chem., 102, 4714. Gao, J., and Freindorf, M., 1997, J. Phys. Chem. A., 101, 3182. Gear, C.W., 1971, Numerical Initial Value Problems in Ordinary Differential Equations, (Prentice Hall). Gogonea, V., 2002, Internal electronic J. of Molecular Design, 1, 173. Goldberg, D.E., 1989, Genetic Algorithms in Search, Optimization, and Machine Learning, (Addison-Wesley).

PAGE 152

136 Guo, J., and Xia, X., 1992, Science, 258, 613. Haile, J.M., 1992, Molecular Dynamics Simulation, (John Wiley and Sons). 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. Hohenberg, H., and Kohn, W., 1964, Phys. Rev., 136B, 864. Houjou, H., Inoue, Y., and Sakurai, M., 2001, J. Phys. Chem. B., 105, 867. Hsiao, Y-W, Runge, K., Cory, M.G., and Bartlett, R.J., 2001, J. Phys. Chem. A., 105, 704. Huff, N.T., Demiralp, E., Cagin, T., and Goddard, W.A., III, 1993, J. Non-Cryst. Solids, 253, 133. Hughes, T.J.R., 1987, The Finite Element Method (Prentice Hall, Englewood Cliffs) Hunbel, S., Sieber, S., and Morokuma, K., 1996, J. Chem. Phys., 16, 1959. Ihm, J., Zunger A., and Cohen, M. L., 1979, J. Phys. C., 12, 4409 Jin, W., Vashishta, P., Kalia, R.K., and Rino J.P., 1993, Phys. Rev. B., 48, 9369. Jones, R.O., and Gunnarsson, O., 1989, Rev. Mod. Phys., 61, 689. Jung, Y., Choi, C.H., and Gordon, M.S., 2001, J. Phys. Chem. B., 105, 4039. Kleinman, L., and Sahni, V., 1990, Adv Quant. Chem. 21, 201. Kohn, W., and Sham, L.J., 1965, Phys. Rev., 140A, 1133. Kroger, M. Stankovic, I., and Hess, S., 2003, Multiscale Model.Simul. 1, 25. Lamberti, V.E., Fosdick, L.D., Jessup, E.R., and Schauble, C., 2002, J. Chem. Ed., 79, 601. Loferer, M.J., Loeffler, H.H., and Liedl, K.R., 2003, J. Comput. Chem., 24,1240. Mallik, A., Runge, K., Cheng, H-P., and Dufty, J.W., 2005, Molecular Simulations, accepted.

PAGE 153

137 Mallik, A., Taylor, C.E., Runge, K., and Dufty, J.W., 2004, Internat. J. Quan. Chem., 100, 1019. Martins, J. L., and Cohen, M. L., 1988, Phys. Rev. B., 37, 6134. Maseras, F., and Morokuma, K., 1995, J. Comp. Chem., 16, 1170. Monard, G., and Merz, Jr. K.M., 1999, Acc. Chem. Res., 32, 904. Mulliken, R.S., 1955, J. Chem. Phys., 23, 1833. Muralidharan. K., 2004, The Molecular Dynamics Simulation of Fracture in Amorphous Silica, Ph.D. Dissertation, University of Arizona, 2004 Murphy, R.B., Philipp, D.M., and Friesner, R.A., 2002, J. Comp. Chem., 21(16), 1442. Nicoll, R.A., Hindle, S.A., MacKenzie, G., Hillier, I.H., and Buton, N.A., 2001, Theor. Chem. Acc., 106, 105. Nocedal, J., and Wright. S.J., 1999, Numerical Optimization, (Springer). Nos, S., 1984, Mol. Phys. 52, 255, W.G. Hoover, 1985, Phys. Rev. A., 31, 1695. Odette, G.R., Wirth, B.D., Bacon, D.J., and Ghoniem, N.M., 2001, MRS Bull., March, 176. Ogata, S., and Belkada, R., 2004, Comp. Mat. Sc., 30, 189. Parr, R.G., and Yang, W., 1989, Density Functional Theory of Atoms and Molecules (Oxford). Perdew, J.P., Burke, K., and Ernzerhof, M., 1996, Phys. Rev. Lett., 77, 3865. 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. Poteau, R., Ortega, I., Alary, F., Solis, Barthelat, and J-C., Daudey, 2001, J. Phys. Chem., 105, 198. Pu, J., Gao, J., and Truhlar, D.G., 2004, J.Phys.Chem A., 108, 632. Proynov, E. I., Ruiz, E., Vela, A., and Salahub, D. R., 1995, Internat. J. Quan. Chem,. 29, 61. Reuter, N., Dejaegere, A., Maigret, B., and Karplus, M., 2000, J. Phys. Chem. A., 104, 1720

PAGE 154

138 Rudd, R.E., and Broughton, J.Q., 2000, Phys. Stat. Sol., 217, 251. Sauer, J., and Sierka, M., 2000, J. Comp. Chem., 21(16), 1470. Schaible, M., 1999, Critical Reviews in Solid State and Material Sciences, 24(4), 265. Singh, U.C., and Kollman, P. A., 1986, J. Comp. Chem., 7, 718. Singh, R.K., and Bajaj, R., 2002, MRS Bulletin, Oct, 743. Stewart, J., 1990 Reviews in Computational Chemistry Volume1, (VCH Publishers). Stillinger, F.H., and. Weber, T.A., 1985, Phys. Rev. B., 31, 5262. Svensson, M., Humbel, S., Froese, R.D.J., Matsubara, T., Sieber, S. and Morokuma, S. J. 1998, Phys. Chem., 100 (50), 19357. Swart, M., 2003, Inernat. J. Quan. Chem., 91, 177. Taylor C.E., 2004, The Transfer Hamiltonian: A Tool for Large Scale, Accurate, Molecular Dynamics Simulations Using Quantum Mechanical Potentials, Ph.D. Dissertation, University of Florida, 2004. Taylor, C.E., Cory, M., and Bartlett, R.J., 2003, Comp. Mat, Sci., 27, 204. Tersoff, J., 1998, Phys. Rev. B., 37, 6991. Thompson, M.A., 1996, J. Phys. Chem., 100 (34), 14492. Titmuss, S.J., Cummins, P.L., Bliznyuk, A.A., Rendell, A.P., and Gready, J.E., 2000, Chem. Phys. Lett., 300,169. Trickey, 1997, Internat. J. Quant. Chem., 61, 641. Troullier, N., and Martins, J.L., 1991a, Phys. Rev. B., 43, 1993. Troullier, N., and Martins, J.L., 1991b, Phys. Rev. B., 43, 8861. Tsuneyuki,S., Tsukada, M.,Aoki, H.,and Matsui, Y., 1998, Phys. Rev. Lett., 61, 869. Van Beest B.W.H., Kramer G.J., and van Santen R.A., 1990, Phys. Rev. Lett., 64, 1955. Vashishta, P., Kalia, R.K., Rino, J.P. and Ebbsj, I., 1990, Phys. Rev. B., 41, 12197. Verlet, L., Phys. Rev., 1967, 159, 98. Warshel, A. and Levitt, M., 1976, J. Mol. Biol., 103, 227.

PAGE 155

139 Watanabe, T., Fujiwara, H., Noguchi, H., Hoshino, T., and Ohdomari, I., 1999, Jpn Journal of Applied Physics, 38, L366. Wentzcovitch, R. M., and Martins, J. L., 1991, Solid State Commun. 78, 831 Wimmer, E., 1997, Electronic Structure Methods. In Catlow C R A and A K Cheetham (Editors) New trend in material chemistry, NATO ASI series C 498. Dordrecht, Kluwer. Woodcock, L.V., 1971, Chem. Phys. Lett. 10, 257. Wright, A.C., 1994, J. Non-Cryst. Solids, 179, 84. Zhang, Y., Lee, T-S, and Yang, W., 1999, J. Chem Phys, 110(1), 46. Zwijnenburg, M.A., Bromley S.T., Flikkema, E., Maschmeyer, T., 2004, Chem. Phys. Lett., 385 389. Zhu, C., Byrd, P.Lu., and Nocedal, J.,1994, Tech Report NAM-II, EFCS Department, (Northwestern University). Zhu, T., Li, J., Yip, S., Bartlett, R.J., Trickey, S.B., and deLeeuw, N., 2003, Molecular Simulations, 29, 671. Zhu, W., Taylor, D.E., Al-Derzi, A.R., Runge, K., Trickey, S.B., Zhu, T., Li, J., and Yip, S., 2005, Comp.Mat.Sc. (submitted).

PAGE 156

BIOGRAPHICAL SKETCH Aditi Mallik was born on the 21 st of January, 1977 (India). She grew up in Kanpur, where she completed high school. After high school, she joined the Presidency College of Calcutta in 1995, with physics as a major and chemistry and mathematics as subsidiary subjects. She graduated with first Class Honors in 1998 and joined the prestigious Physics Department of Indian Institute of Technology (IIT), Kanpur to pursue her M.Sc. level studies. In the course of her masters program she developed a special interest in condensed matter theory. Encouraged by the Professors there, she decided to come to USA to join the Physics Department at the University of Florida (UF) for Ph.D. studies in Fall 2000. Professor James W. Dufty was one of her course instructors. She took two graduate courses on classical and statistical mechanics from him. In summer 2001, she joined his research group where she worked on the problem of embedding scheme in the multi-scale modeling (which is also a part of the Quantum Theory Project). She expects to graduate with a doctorate degree (Ph.D.) in spring 2005. She has plans to stay in academics in the near future, and to pursue post doctoral research in her field of specialization. 140


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

Material Information

Title: Multi-Scale Modeling of Solids as a Composite of Quantum Mechanical (QM) and Classical Mechanical (CM) Domains
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: UFE0009740:00001

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

Material Information

Title: Multi-Scale Modeling of Solids as a Composite of Quantum Mechanical (QM) and Classical Mechanical (CM) Domains
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: UFE0009740:00001


This item has the following downloads:


Full Text











MULTI-SCALE MODELING OF SOLIDS AS A COMPOSITE OF QUANTUM
MECHANICAL (QM) AND CLASSICAL MECHANICAL (CM) DOMAINS















By

ADITI MALLIK


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


2005

































Copyright 2005

by

Aditi Mallik
































This dissertation is dedicated to my parents.















ACKNOWLEDGMENTS

This dissertation owes its existence to my supervisory committee chair Professor

James W. Dufty of University of Florida from whom I have learnt to perceive a research

problem in proper context, translate rudimentary ideas into formal analysis, and check the

logical consistency of the presentation. I cannot adequately thank him for his attention,

emotional support, and inspiration that I have received ever since I joined the University

of Florida, Physics Department in 2000. I would also like to thank my other supervisory

committee members, Professors Samuel B. Trickey, Hai-Ping Cheng, John Mecholsky,

Chris J. Stanton, and Keith Runge. Their valuable time and suggestions helped me to

complete the work within a reasonable time frame.

I am deeply indebted to Dr. Keith Runge for his unstinted interest in my research

progress and especially his remarkable patience in answering all my naive questions on

Quantum Chemistry. Thanks go to Dr. DeCarlos Taylor for helping me with the

computer programming. The implication of this support will only be appreciated by those

who frequently got lost in the programming-maze as I did. I am also grateful to Professor

Hai-Ping Cheng for teaching me the parallel programming code on Density Functional

Theory. Dr. Krishna Muralidharan is gratefully acknowledged for collaboration in his

work on glass which developed into a part of a chapter in this dissertation.

The Physics Department has been a wonderful place for my learning experience. It

was a great experience to be a part of the Quantum Theory Project (QTP) group with

people who have been good friends, sources of discussion and technical support. I wish to









thank all former and current members of the QTP group especially Dr. Ling-Ling Wang,

Dr. Moahua Du, Chao Cao, and Ying-Xao Wang for their help and mentorship. I am also

grateful to my friend Aparna Baskaran who offered pertinent comments on a few

chapters in this dissertation.

I am deeply indebted to those who sustained me at a personal level. My heart-felt

appreciation goes to my parents and all my family members back home in India. Without

their constant emotional support, I simply could not have come this far. I wish to thank

my husband, Samir, for making my stay in the US so enjoyable; not to mention the

painstaking job of proof-reading this dissertation.

This work was support by ITR grant 0325553.
















TABLE OF CONTENTS

page

A C K N O W L E D G M E N T S ......... .................................................................................... iv

LIST OF TABLES ..................................... ....... .... ..... .......... ix

LIST OF FIGURES ......... ......................... ...... ........ ............ xi

GLOSSARY ................................... ........................ xiv

CHAPTER

1 IN TR OD U CTION ............................................... .. ......................... ..

D definition of M ulti-Scale M odeling .................................. .................. ..................
Exam ples of M ulti-Scale M odeling.................................... ........................... ......... 2
M ethods for Different Length Scales ....................... ...... ........... ................... 5
Quantum Mechanical (QM)/Classical Mechanical (CM) Simulations ....................6
Proposed Schem e .................................... ................................ ......... 7
M odel System ................................................................. ...... ......... 8

2 THE MATERIAL AND ITS FORMAL PARTITIONING ..................................... 11

In tro du ctio n ......................................................... ............. ................ 1 1
The Idealized Q uantum Solid ............................................................... ............... 12
The Form al Partition and Com posite Solid ............................................................. 16
The Representative Classical Solid ............................................................16
The Composite Quantum/Classical Solid................................. ............... 18
Description of the Classical Subsystem ............................................. ...............21
Description of the Quantum Subsystem .......................................... ............... 21
Coulomb Effects of Environment.................................................................24
Electron Exchange with Environm ent.............................................................. 25

3 DESCRIPTION OF THE QUANTUM DOMAIN..................................................28

Introdu action ...................................... ................................................. 2 8
Q uantum M methods .................. ....... .......................... .... .. .... ...............29
D description of the Q M D om ain ...................................................................... .. .... 32
D different Term nation Schem es ........................................ ....... ............... 32
P seudo-atom M ethod ................................................................ .....................35









Calculation of Charge D ensities.................................................................. 37
Em bedding Schem e ........................................... .. .. .. .... ................. 40
O their O options for Links.................................... ...... ................. ............... 44
C o n c lu sio n s........................................................................................................... 4 7

4 CONSTRUCTION OF A NEW CLASSICAL POTENTIAL FOR THE CM
D O M A IN .............................................................................................................. 4 9

Intro du action ........................................... .. ....49..........
M olecular D ynam ics...... .... ........ ........................................................ .... 50
Corrector-Predictor A lgorithm s...................................... ......................... 50
V elocity R escaling .......... ..... ........................................................ .. .... .. ... .. 52
Failure of Standard Potentials............................. ................................................ 52
Objective and Method for Constructing a New Potential ........................... ........54
R results .............. ......... ................ ................. ............ ................ 60
Finding Param eters for the Potential ....................................... ............... 60
Stress-Strain C urves .......................... .... ...... .............................. .. .......... 64
Improving Stress-Strain Curves for Higher Strains ........................... .........67
C o n clu sio n s..................................................... ................ 7 0

5 CONSTRUCTION OF THE COMPOSITE NANOROD ......................................72

In tro d u ctio n ......................................................... ............. ................ 7 2
R results for the C om posite R od ............................................................................. ...73
C o n clu sio n s..................................................... ................ 7 4

6 INDEPENDENCE OF MULTI-SCALE MODELING WITH RESPECT TO THE
CHOICE OF UNDERLYING QM...................................... .......................... 75

In tro d u ctio n ................... .................................................... ................ 7 5
R review of D FT M ethod ............................................................................76
Born-Oppenheimer Approximation...................... .......................... 76
H ohenberg-K ohn Theorem s ........................................... ......................... 76
K ohn-Sham Equations.............. ............... ................ ................. ............... 77
Local Density Approximation ............................. ................................. 78
Generalized Gradient Approximation ...................................... ............... 79
Troullier-Martin Pseudo-Potential.................. .............. ..... ............... 80
D u al Sp ace F orm alism ............................................................. .....................82
Com putational D details ........................ ................ ................... ..... .... 82
Results and D discussion ........................... ..... ..... ...... .. ............. 83
Training of DFT-Pseudo-atoms...................... ..... .................. 83
Construction of Classical-DFT Potential for the CM Region .............................86
C o n clu sio n s..................................................... ................ 8 9









7 APPLICATIONS OF THE MULTI-SCALE MODELING PROCEDURE TO
O TH ER SIL IC A SY STE M S ........................................................... .....................90

In tro d u ctio n ........................................................................................................... 9 0
System s Studied.................................................91
Silica N anorings ........................................ .... ....... .... .. .....9 1
Silica N anoclusters ...................... ................ ................... ..... .... 94
N otched N anorod ............ .... ............................................................. ..... ..... .. 95
B ulk G lass .............. ........................... ...................................... 97
Results Obtained from New-DFT-Classical-Potential ..........................................104
The 10-Membered Silica Nanoring ......... ......... .....................104
Bulk-Glass ............... ......... ................. 105
C o n c lu sio n s......................................................................................................... 1 0 6

8 COMPARISON OF THE DIFFERENT QUANTUM MECHANICAL
M E TH O D S ...................................................... ...........................107

Introduction ......... .... ..............................................................................107
The Nanorod ......................................... ........ ....................... 108
C om prison of Structure...................................................................... ... ..... 108
Elastic Properties .................. ............................. .. .. .... ................. 109
Comparison of Different Types of DFT ........................................ ..............110
B u lk G la ss ................................................................................................ ..... 1 1 1
D iscu ssion ................................................................................................ ..... 1 12

9 SUMMARY AND FUTURE RECOMMENDATIONS ......................................115

C o n c lu sio n s............................... ........................................................................ 1 1 5
Recommendations for Future Work ........................... ........... ...............118

APPENDIX

A NEGLECT OF DIATOMIC DIFFERENTIAL OVERLAP (NDDO) MODEL......120

B G EN E TIC A L G O R ITH M ............................ ................................. .....................124

C SCA LIN G M ETH OD .. ................................................................. ............... 128

D CHARGE TRANSFER POTENTIAL ........................................ ............... 130

Charge-Transfer Potential for Old Transfer Hamiltonian Classical Potential #1.....131
Charge-Transfer Classical Potential for New-Density Functional Theory..............133

L IST O F R E FE R E N C E S ....................................................................... .................... 134

BIOGRAPHICAL SKETCH ............................................................. ............... 140
















LIST OF TABLES


Table p

3-1. Comparison of charge densities of the ring with LA and ring with pseudo-atoms
with that of the ring in the bulk calculated at planes parallel to the plane of the
ring ................................................................................. 39

3-2. Percentage charge difference with respect to bulk calculated with our method for
various cases studied ...................... .................. ............................ 43

4-1. Different values of potential parameters obtained from GA which gave a stable
ro d .............. .................... ..................................... ......... ...... 6 1

4-2. Parameters for New Potential in comparison to TTAM and BKS potential
p aram eters ........................................................................... 62

4-3. Comparison of structure of the rod obtained from the different potentials ..............63

4-4. New parameters for classical potential-2 (NTH-2 Classical Potential)...................68

4-5. Structure obtained from NTH-2 Classical Potential in comparison to that of the
T H nanorod ........................................... ........................... 69

6-1. Magnitude of force on the terminated Si atom in pyrosilicic acid with changing
param eters of the dft-pseudo-atom ........................................... .......................... 84

6-2. Percentage charge density difference obtained with our scheme versus the charge
density in the actual molecule for different cases studied in the DFT method ........85

6-3. Parameters for New-Classical-DFT Potential in comparison to standard TTAM
and BK S potential param eters........................................................ ............... 87

6-4. Structure of the nanorod obtained using New-Classical-DFT, BKS potentials
versus structure of D FT m ethod ................................................................... ......87

7-1. Comparison of structure of the nanoring obtained from TH, NTH-2 Classical
Potential and TTA M potential.......................................... ............................ 92

7-2. Structure comparison of 2-membered and 3-membered silica rings.........................94

7-3. Comparison of different properties of the glass obtained from the NTH-2
Classical Potential to those of the BK S glass ............... .............. ..................... 101









7-4. Comparison of structure of the nanoring obtained from DFT and New-DFT-
C classical P potential ................................................................... .. .... .. 105

7-6. Comparison of different properties of the glass obtained from New-DFT
Classical Potential to those of BK S glass......... ................................... .............105

8-1. Comparison of structure of the nanorod obtained from DFT and TH quantum
m eth o d s ......................................................................... 1 0 8

8-2. Comparison of the 10-membered nanoring using two different types of DFT
m methods ...................................................... ........ .. ......... ........... 110

8-3. Comparison of structure of bulk glass using two different classical potentials ......111
















LIST OF FIGURES


Figure page

1-1.Use of different methods at different length scales in crack propagation. ...................2

1-2. Chemo-mechanical polishing of silica wafers ........................................................3

1-3. Structure of the nanorod. A) Front view. B) Top view. The red spheres are the Si
atom s and the blue spheres are atom s................................... ....... ............... 9

1-4. Partitioning of the nanorod into QM and CM domains. Localized (blue) electron
cloud of the CM region shows the appropriateness of such partitioning.................10

3-1. Training of Pseudo-atom on smaller molecules, here Pyrosilicic acid....................36

3-2. Charge densities with different termination schemes. (A) Isolated ring (B) with
links (C) with pseudo-atoms (D) in actual nanorod ..............................................38

3-3. Comparison of magnitude of force on Si nuclei with different termination
sch em e s. .......................................................... ................ 4 0

3-4. Embedding Scheme: Approximating the rest (CM region) of the rod beyond the
pseudo-atom s by dipoles. ............................................... .............................. 41

3-5. Effect of inclusion of the dipoles in the TH. (A) Improvement in the value of
forces. (B) Reduction in the value of normalized difference of charge densities
Open symbols, without dipoles; squares, with dipoles. ........................................41

3-6. Cross-sectional view. (A) Distorted ring. (B) Ring in equilibrium.........................42

3-7. Comparison of forces on Si nuclei with our scheme and LA to those in the actual
molecule for various cases studied ......................... ......... .. .............. 43

3-8. Partitioning of 3-membered silica ring........ ............... ........ ..... ............ 44

3-9. Deviated position of optimized LA. ............................................... ............... 45

3-10. Charge density with LA at optimized positions. .............................................. 45

3-11. Position of LA beyond intermediate oxygen atoms. .................. ....... ............46









3-12. Comparison of force on Si nuclei with LA on oxygen to LA on silicon and our
em bedding scheme e. ................................................ ...............47

4-1. Comparison of stress-strain curves for the standard pair potentials (BKS, TTAM)
w ith that of Transfer H am iltonian ........................................ ........................ 54

4-2. Presence of local and global m inim a ........................................ ...... ............... 59

4-3. Comparison of Si-O pairwise interactions for the three potentials. ..........................63

4-4. Comparison of Young's modulus at low strains for different classical potentials
w ith that of TH ........................................................................64

4-5. Reduction of the wiggles with temperature..........................................................65

4-6. Stress-strain curves at lower strain rate of 2m/s and a temperature of 0.2K ...........65

4-7. Histogram of number of modes vs. wave number (cm). ............................... 67

4-8. Behavior of NTH-2 Classical Potential near fracture..............................................69

5-1. Comparison of stress-strain of the composite rod with those obtained from NTH-
2 Classical Potential and TH ............................................................................. 74

6-1. Comparison of force on the Si nucleus with different termination schemes and
dipoles for the D FT case. ............................................... .............................. 85

6-2. Comparison of adiabatic stress curves for DFT, new-DFT-classical and TTAM
potential ..................................... .................. .............. ........... 88

7-1. The 10-m embered silica nanoring. ........................................ ........................ 91

7-2. Stress-strain curve for the nanoring obtained from TH and NTH-2 Classical
Potential ..................................... .................. .............. ........... 93

7-3. The 2- and 3-m embered silica rings. ................... .......... .......... ............ .....94

7-4. Silica Nanoclusters. A) Si12024. B) Si18036. C) Si24048-Cage. D) Si24048-
F ullerene ................... .................................... ...........................95

7-5. Front view of the 107 atom notched nanorod showing one defect (missing
oxygen) ................... ....... ......................... ....... ................. 96

7-6. Comparison of stress-strain curves for the notched rod obtained from different
methods, showing the significant reduction in the yield stress with the presence
of a single defect. ................................................... ................. 97

7-7. Structure of glass obtained from the NTH-2 Classical Potential.............................98









7- 8. Radial distribution function of the Si-O, 0-0 and Si-Si pairs in the glass.............100

7-9. Bond angle distribution of O-Si-O and Si-O-Si ................................ ............... 100

7-10. Stress-strain curves for NTH-2 and BKS glass. ............................................... 102

8-1. Comparison of elastic properties for DFT and TH nanorods............................. 109

8-2. Comparison of elastic properties of the two glasses in the elastic region .............12

D-1. Charge transferred from Si to O as a function of distance. .................................... 131

D-2. Increase in the stiffness with introduction of the charge transfer term................ 132

D-3. Increase in the charge of Si with strain..........................................................132

D-4. Increase in the stress with the charge transfer DFT potential..............................133















GLOSSARY

BAD Bond angle distribution

CC Coupled Cluster

CM Classical mechanics / mechanical

DFT Density Functional Theory

FE Finite Element

GA Genetic Algorithm

LA Link Atoms

MD Molecular dynamics

NDDO Neglect of Diatomic Differential Overlap

NTH-2 New Transfer Hamiltonian Classical Potential # 2

QM Quantum mechanics / mechanical

RDF Radial distribution function

TB Tight Binding

TH Transfer Hamiltonian

TM Troullier-Martin

Y Young's modulus















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

MULTI-SCALE MODELING OF SOLIDS AS A COMPOSITE OF QUANTUM
MECHANICAL (QM) AND CLASSICAL MECHANICAL (CM) DOMAINS

By

Aditi Mallik

May 2005

Chair: James W. Dufty
Major Department: Physics

The simulation of non-equilibrium phenomena such as fracture in a macroscopic

sample using ab initio quantum methods is not feasible because of the high

computational time required. The alternative approach of a classical mechanical method

does not have this limitation, but an accurate description cannot be made for important

states far from equilibrium. A promising tool for studying such phenomena is the

application of quantum techniques only to the small reactive regions and classical

mechanics to describe the near-equilibrium bulk. Such methods are referred to here as

quantum mechanical/classical mechanical (QM/CM) methods in multi-scale modeling.

Our study addressed the problem of how to bridge these QM and CM domains by

constructing an approximate semi-classical representation of the quantum (actual) solid

as an equivalent composite system comprising a small QM domain plus the CM bulk

region. Considering each domain as an open system, it is necessary to incorporate the

relevant information about the state of the complimentary domain in a simple and









accurate manner. For the QM domain, the CM environment is reduced in terms of trained

pseudo-atoms and dipoles to account for the short-range chemical bonding and long-

range Coulombic interactions, respectively; for the CM region a new classical pair

potential has been constructed using quantum force data, such that this potential yields

the same equilibrium structure and elastic properties as the QM region.

This approach to multi-scale modeling has been applied to a SiO2 nanorod and SiO2

nanoring for critical testing. It is shown that the modeling reproduces the forces, charge

densities, and elastic properties to within a few percent of the actual quantum values.

Applications to hydrolytic weakening and fracture in bulk glass are discussed. Although

our study was limited to silica-based materials, it should provide guidance for extension

to other materials as well.














CHAPTER 1
INTRODUCTION

Definition of Multi-Scale Modeling

The field of multi-scale modeling has opened a new era to computational science

for studying complex phenomena such as fracture, hydrolysis, enzymatic reactions,

solute-solvent studies, hydrogen embrittlement, and many other chemo-mechanical

processes in macroscopic samples. Often these phenomena require a very accurate

description at one scale, while at another scale, one can apply a much less sophisticated

method to get satisfactory results [Broughton et al., 1999]. In fact, it is necessary to have

the less sophisticated description for the bulk sample because the higher accuracy

methods are computationally too intensive to be applied to the bulk. These scales may be

length or time or a combination of both. This scheme of combining different models at

different scales to achieve a balance of accuracy, efficiency and realistic description is

known as multi-scale modeling. It is accomplished by applying the very high accuracy

method only in a small domain where it is needed the most and more approximate

methods for the rest of the bulk where they are appropriate. The next example shows the

need to combine different models.

Crack propagation: For the case of crack propagation (Figure 1-1), one has to apply

detailed quantum mechanics at the tip of the crack where the bonds are breaking between

the atoms, causing a marked deformation of the electron cloud and charge transfer among

ions. Classical mechanical methods cannot be applied in this region because classical

Hamiltonians forgo any electronic treatment. But far from the crack tip, where the atoms









are deformed less, atom interactions can be described well by classical potentials. It

would be computationally infeasible to apply the detailed quantum mechanical methods

to the whole sample; as the application of these methods are limited to only very small

molecules (except for crystals). Thus, one must link different models for the different

length scales; as no single theory is sufficient to describe the process of crack propagation

in the bulk sample. Broughton et al., [1999] described multi-scale simulations in detail.














Figure 1-1.Use of different methods at different length scales in crack propagation
[Broughton et al., 1999].


Examples of Multi-Scale Modeling

(i) Chemo-mechanical polishing (CMP): combines chemical and mechanical

interactions to planarize silica wafer surfaces, using a slurry composed of solvents,

wetting agents, other compounds and submicron-sized particles [Singh and Bajaj, 2002].

The chemical compounds interact with the material to form a chemically modified

surface and the abrasives in the slurry mechanically interact with the chemically modified

surface layers, resulting in material removal. Here interactions occur at both micro and

nano-scale, so these two scales have to be treated generally with two different methods.

At the microscale, the pad carrying the abrasive particles interacts with the surface of the









wafer (Figure 2-2). At the nanoscale, the kinetics of the formation and removal of the thin

surface layer controls CMP output parameters (such as removal rate, surface planarity,

density of surface defects etc.).




Subsma slur


Pad

SuIbstate Chemically pa I aling







Figure 1-2. Chemo-mechanical polishing of silica wafers [Singh and Bajaj, 2002].



(ii) Tidal wave prediction: Clementi [1988] studied the interactions of several water

molecules on an angstrom scale using quantum calculations. From this database, he found

an empirical potential for molecular dynamics simulation, which was then used to

determine the viscosity and density of water. Finally they used the density and viscosity

of water in computational fluid dynamics to model tidal circulation in Buzzard's Bay of

Massachusetts.

(iii) Atmospheric science: In this type of multi-scale modeling one has to combine

different length as well as time scales. High accuracy computational methods are used to

determine reaction barriers for simple chemical reactions which are then used in large

rate equations coupled to grid computation to predict chemical meteorology.









(iv) Embrittlement of pressure-vessels: Another multi-scale process in which the

phenomenon occurs over different time scales is embrittlement of pressure-vessel steels

of nuclear reactors caused by prolonged exposure to nuclear radiation [Odette et al.,

2001].

In the proceeding examples of tidal waves and atmospheric science, the various

scales are weakly coupled i.e., the behavior at one scale does not influence the

phenomena at the others strongly. Rudd and Broughton [2000] described such models as

"serial multi-scale-modeling" which means the computation of parameters at smaller

scale for its use in more phenomenological models at a larger scale. The other type of

multi-scale problem, in which the scales are strongly coupled is known as "concurrent

multi-scale modeling". In this type, all the different models applied on different scales

must be nested with proper boundary conditions. For example in brittle fracture, the

fracture takes place on a macroscopic scale with a propagation direction determined by

local stresses in the material. At the atomic scale, the chemical bonds are taken far from

the equilibrium. So the overall fracture mechanism depends on the instantaneous atomic

bond lengths, at the time of strain application, yet the crack changes character as it

propagates over macroscopic distances. Similarly, the CMP process can be

conceptualized as the altered reactivity of a strained bond. This alteration occurs in

combination with effects of local surface, neighboring atoms, and strain distribution. In

the presence of a chemically reactive species, a stress that normally would produce a

local elastic strain instead can lead to bond breaking.

Our multi-scale modeling is more of the "concurrent" type, though some of the

work can be considered as the "'t a/i, type. Chapter 3 shows that it is necessary to









include the information about other scales while studying the phenomena at one

particular scale as in the "concurrent" method. However the parameterization of the

quantum method (Transfer Hamiltonian) used is initially done on smaller molecules. The

parameterization then is transferred to bigger molecules as in the "serial" method. Such

distinctions are discussed more in chapter 3.

Methods for Different Length Scales

One can distinguish three different length scale levels in the context of multi-scale

modeling of materials

(1) The nanoscale: Here electrons play an important role so the phenomena are

described using quantum mechanical methods like Density Functional theory (DFT) [Parr

and Yang, 1989], post-Hartree-Fock methods [Stewart, 1990] and quantum Monte Carlo

(QMC) [Foulkes et al., 2001].

(2) The atomic or microscale: Here phenomena are described in terms of classical

interactions among atoms via Molecular Dynamics (MD) [Haile, 1992] or Monte Carlo

(MC) simulations.

(3) The macroscale: Here continuum fields such as density, temperature, velocity,

stress-fields etc., are important players, and therefore finite-element (FE) methods

[Hughes, 1987] are used to examine the large-scale properties of the material.

We studied only one particular subclass of multi-scale modeling, known as hybrid

quantum mechanical (QM)/classical mechanical (CM) problems, and their application to

solid insulators, silica in particular. This scheme is sometimes referred to as QM/MM

(molecular mechanics) simulation. In the literature, the classical MD region is further

embedded in the FE continuum region [Broughton et al., 1999]. However this MD/FE

coupling is unbalanced because of the differences in the dispersion relations and the









underlying elastic properties between the MD and FE continuum region. The work on

MD/FE coupling is beyond the scope of our study and is not discussed here.

Quantum Mechanical (QM)/Classical Mechanical (CM) Simulations

Over the past decade a large body of work has been done on this QM/CM topic.

Warshel and Levitt [1976] laid out their first formalism on the QM/CM theory in their

seminal paper. Aqvist and Warshel [1993] and Gao [1996] gave a complete overview of

QM/CM simulation methods. Applications of QM/MM include study of

(i) Biological systems like enzymes, DNAs and proteins [Amara et al., 1999;

Monard and Merz, 1999; Hall et al., 2000; Titmuss et al., 2000; Houjou et al., 2001;

Nicoll et al., 2001; Cui and Karplus, 2002; Dilabio et al., 2002; Gao and Truhlar, 2002;

Murphy et al., 2002; Crespo et al., 2003; Lofere et al., 2003]. For studying enzymes, one

applies the quantum method to the active sites in the enzyme molecules where ligands are

binding, while the rest of the molecule that does not take part in the chemical reaction is

treated classically. Gogonea [2002] reviewed applications of QM/CM methods to

biological systems.

(ii) Vibrational spectroscopy of complexes [Chabana and Gerber, 2001].

(iii) Electronic excitations [Thompson, 1996; Gao et al., 1996; Vries and Duijen,

1996].

(iv) Hydrolysis of silica [Jung et al., 2001; Du et al., 2003].

(v) Solute-solvent problems [Eichinger et al. 1999; Guo and Xia, 1992; Gao and

Freindorf, 1997] in which the solute is described quantum mechanically and solvent that

plays the role of an external electrostatic environment is treated classically.

However if the QM/CM partitioning has to be done in a single molecule (as in our

case) the partitioning becomes more complicated because one or more covalent bonds









have to be cut between the QM and CM regions. Then a proper treatment has to be done

to saturate the dangling bonds because these dangling bonds can give rise to incorrect

charge densities in the QM domain (chapter 3) and to virtual energy bands in the energy

gap [Sauer and Sierka, 2000], electronic states around the Fermi level which are localized

at the peripheral region of the cluster containing dangling bonds that are otherwise not

observed in bulk Si [Ogata and Belkada, 2004].

Proposed Scheme

The challenge of the QM/CM simulation is to move from one length scale to

another, as smoothly as possible. We have proposed a scheme for multi-scale modeling

that would give a "seamless" coupling between the QM and CM regions. This is based on

the following three criteria:

(i) Modeling environmental effects for the quantum domain in an accurate and

simple manner. In QM/CM partitioning methods for multi-scale modeling, one is often

forced to introduce uncontrolled phenomenological effects of the environment (CM) in

the QM domain. Chapter 3 shows that one should consider the precise effects of the CM

region in the QM domain to get an accurate description of the QM domain when

separated from its CM bulk. Here, accurate means that the quantities under study are

indistinguishable whether determined from a length-scale coupled system or from a

system comprising only the finest description scale. Our characterization is based on

predicted quality of charge densities and forces (which implies correct geometry) in the

QM domain [Mallik et al., 2004]. Previous research focused on obtaining only the correct

geometry; so far, no studies have reported on getting the correct charge densities to the

best of author's knowledge. The rationale behind selecting the charge density as a

criterion to assess the quality of QM/CM embedding is explained in chapter 2. The









details of implementations of embedding the QM domain in its classical environment are

given in chapter 3

(ii) Developing a new classical potential for the classical region. Regarding the

coupling to the CM region, most researchers consider some kind of handshaking region

in which they average the force between two regions [Rudd and Broughton, 2000]. This

technique of force matching by choosing the average of two forces is however very

artificial. We avoided this obstacle by designing a new potential for the entire CM region,

such that forces obtained from this potential are the same as those given by the QM

domain [Mallik et al., 2005]. Chapter 4 proposes that each development of a multi-

scaling algorithm should include the construction of its own classical potential. The

potential is constructed using the information about the forces of the QM domain when

strained as well as in equilibrium such that the structure and elastic properties predicted

by this classical potential agree with those of the QM domain. The use of this new

potential will prevent any mismatch of mechanical properties across the QM/CM

boundary.

(iii) Building a composite material by joining the QM and CM regions with suitable

boundary conditions. The composite material is required to be indistinguishable from the

one obtained from full quantum or classical calculation on the entire system for small

strain elastic response.

The accuracy of the proposed method is demonstrated in detail for a benchmark

system: a silica nanorod.

Model System

The complete formalism has been illustrated in detail for a model system: a silica

nanorod [Zhu et al., 2003]. This silica nanorod containing 108 nuclei has been chosen as









our model system since a quantum treatment of the whole sample is possible to assess the

quality and the physical basis for the proposed scheme and yet the system is big enough

to yield bulk properties such as the Young's modulus. The nanorod has the proper

stoichiometric ratio of silicon to oxygen observed in real silica (1:2) and is considered a

viable model for studying silica. Figure 1-3 shows front and top views of this nanorod.

The nanorod consists of six Si606 rings and each of the four planar ring shares,

both above and below, a ring of oxygen atoms with adjacent six member rings. To

alleviate any dangling bonds, the nanorod is terminated with capped rings whereby each

silicon atom of the terminating ring is connected by bridging oxygen or two interstitial

oxygens. The size of this nanorod can be readily adjusted by adding or removing Si606

planar rings with corresponding O atoms. In the present study, one of the rings near the

center of the rod is chosen to be the QM domain and rest of the rod is treated classically

(Figure 1-4A).


















Figure 1-3. Structure of the nanorod. A) Front view. B) Top view. The red spheres are the
Si atoms and the blue spheres are O atoms.









The localized nature of the valence electron charge density (the blue colored cloud

in Figure 1-4) of the CM regions ensures the appropriateness of such a partitioning. The

lighter shades represent higher electron densities and the black region in between shows

that no overlap exists between the two clouds.


CM


CM


Figure 1-4. Partitioning of the nanorod into QM and CM domains. Localized (blue)
electron cloud of the CM region shows the appropriateness of such
partitioning.

The first, second and third criteria of multi-scale modeling are detailed in chapters

3, 4, and 5 respectively.














CHAPTER 2
THE MATERIAL AND ITS FORMAL PARTITIONING

Introduction

In this chapter, we present the theoretical basis for multi-scale modeling. The

objective is to provide a means to assess the approximations made more directly, and also

to provide the basis for their extension to other systems (e.g., metals). Here, the problem

of multi-scale modeling is introduced by stating what variables are most appropriate for

control in the multi-scale modeling. First, the quantum description in terms of ions and

electrons is given and the limitations for practical implementation are noted. Next,

reduced self-consistent descriptions are given for the ions and electrons, each coupled to

the other through their average charge density. In this form, the classical limit for the ions

can be taken, allowing the use of MD simulation methods for their dynamics. The

equation for the electrons is simplified exploiting the different time scales for electron

and nuclear motion (the Born-Oppenheimer limit). The first description constitutes the

idealized quantum solid for which the subsequent multi-scale modeling scheme is

proposed. This idealized quantum solid is then represented as a composite of two

domains, the larger bulk in which a classical representation is used and a smaller

"reactive" domain in which the original quantum domain is retained. The objective of the

modeling is then to construct a potential function that gives an accurate description in

both the small and the classical domain. This has two components: the determination of a

pair potential for the forces on ions in the classical domain and an accurate calculation of

the charge density for the ions in the quantum domain. For the latter calculation,









approximations are made in term of "pseudo-atoms" and dipoles to include the relevant

effects of the classical part onto the quantum domain.

The Idealized Quantum Solid

At the fundamental level a simple atomic or molecular solid can be described in

terms of N, "ions" with charge number Z for the corresponding atoms of species, and

a set of N, electrons, with overall charge neutrality (-N,Z, = N, ). To introduce the


various levels of approximation we start with the density operatorD for the system as a

whole. Properties of interest A are given by the expectation value

(A) = Tr,,DA (2-1)

where the trace is taken over all electron and ion degrees of freedom. The density

operator obeys the Liouville-von Neumann equation


,D+ -[H,D]= 0 (2-2)
h

The Hamiltonian operator H is comprised of the Hamiltonians H, and He for the

isolated systems of ions and electrons, respectively, and their interaction U,e

H =H, +He + U,, (2-3)


U f dr dr-d-' P (f)p (') (2-4)
r r'

Here p, (f) and Pe (r) are the ion and electron charge density operators. This is the most

fundamental level for a quantum description of the system. For a pure state, D is a

projection onto a state Vf which is governed by the corresponding Schrodinger equation.

For a small system, this can be solved using high-quality quantum methods; however this









is not a practical method for application to larger systems. Here only non-crystalline

systems are considered in which one cannot apply periodic boundary conditions.

In many cases the properties of interest are functions only of the ion degrees of

freedom (e.g. structure), A A, or they are purely electronic (e.g., optical), A A,.

Then a description is possible in terms of the reduced density operators for the ions and

for the electrons, D, andDe, resulting from appropriate partial traces over all the

electrons or ions, respectively

(A,) =TrDA,, A,) = TrD,A,, (2-5)

D, TrD, D, Tr D. (2-6)

Their equations of motion follow directly from equation 2-2.

a 1
t,D, +-[H,,D, ]+I(U,D, D,U,) = 0, (2-7)
h h

I 1
8AD +[HD,]+ (UD, -D,U, ) = 0. (2-8)


where the potential energy operators coupling the nuclear and electronic degrees of

freedom are


U, = fdif r- p,()pel'), e(f?)= (Trpe(F')D)(TreD) -, (2-9)


Ue = fcr- (f) p,('), p(f ) (Tr p, (')D)(Tr~D) 1 (2-10)
r- r- |

This is similar to the microscopic ion electron coupling of 2- 4 except that now the

electron charge density operator pe (f) is replaced by its conditional average, Pe (r') in

the equation for D, and the ion charge density p, (r)is replaced by its conditional

average in the equation for De. The equations 2-7 and 2-8 are still exact but formal since









these average charge densities are not determined by these equations alone. The simplest

realistic approximation (mean field) is to neglect the direct correlations in the charge

densities p, and ,e i.e. make the replacement D DD, in equations 2-7 and 2-10 to get

p, -Tr p, (F)D, p, -> TePe (F)D, (2-11)

As a consequence, the potentials U, and Ue become Hermitian and the average charge

densities are now self-consistently determined by 2-7 and 2-8. Self-consistency is

required since De and D, are functionals of the average charge densities p, and pe,

respectively.

The advantage of this reduced description is that the ions and electrons are

described by separate equations, reducing the degree of difficulty of each. The large

differences in electron and nuclear masses imply corresponding differences in time scales

and thermal de Broglie wavelengths. Consequently, the equation for the ions admits a

classical limit for the conditions of interest. The classical limit of equation 2-7 becomes

,D, + {(H, +U,[pe]),DD,} 0 (2-12)

where { now denotes a Poisson bracket operation, andD, is a function of the ion

positions and moment D, D, ({R, }, {P1 }, t). This equation now can be solved

accurately and efficiently by molecular dynamics simulation methods, even for large

systems.

Implementation of equation 2-12 still requires calculation of the potential energy

U, [Pe ], which in turn requires determination of the electronic charge density from De.

However, the general solution to 2-8 is a formidable problem: determination of the

dynamics of electrons self-consistently in the presence of a changing ion charge density.









Two simplifications are made to bring this problem under control. First, is the Born-

Oppenheimer approximation and the second assumption is that only the lowest energy

state contributes to D, at any given time. In principle, 2-12 is solved in time steps At. At

each time step the electronic charge density is computed for the ion configurations at that

time step in order to recompute new forces for the next time step. The electronic charge

density calculation is a ground state eigenvalue problem for the given ion configuration.

While simple in principle, the number of time steps and number of electrons for each

calculation easily can be unreasonably large for practical implementation.

To summarize; the final description of this idealized solid consists of a set of point

ions governed by the classical equation 2-12 for their probability distribution

D, D,{R,),P,),t)

8,D, +fH, +U,[pe], D,= 0 (2-13)





The average electron density pe (r)is determined from the ground state solution to 2-8.
= 1_





[He + U,,De] = 0, (2-15)


U, f= r r- (r, t)p (f'), p,(f,t)= Tr, p,() D, [p(, t),t] (2-16)
Sr- r'|

The analysis proceeds stepwise. The classical equations 2-13 are solved

analytically in discrete time steps for the atomic coordinates and moment. At each step

the electron problem 2-15 is solved for the electronic ground state using the ion

configuration at the previous time step. From this ground state the electron charge density

is determined. This gives the potential energy function U, and consequently the forces









required to change the ion positions and moment at the next time step. The process is

repeated with a new electron charge density calculated at each time step using the new

atomic configurations. All electron correlations are accounted for quantum mechanically

in the eigenvalue problem; all ionic correlations are determined classically through direct

solution of Newton's equations. All structural properties of interest can be calculated

since the phase points for the ions are known at all times. The dynamics of the electrons

is only coarse-grained as they are "slaved" to the time dependence of the ions.

The foregoing treatment defines a "quantum molecular dynamics" representation of

the idealized quantum solid. The resulting description allows an accurate classical

treatment of the atomic structure while retaining relevant quantum chemistry for the

interatomic forces due to electrons. If the scheme could be implemented in practice for

bulk systems of interest over reasonable time intervals there would be little need for

multi-scale modeling. With MD simulation the solution to 2-13 once U, has been

provided is straightforward. So the problem has been reduced to a determination of the

electron charge density. Unfortunately, the solution to 2-15 using realistic ab initio

quantum methods for even a few hundred ions at each time step becomes prohibitively

time intensive.

The Formal Partition and Composite Solid

The Representative Classical Solid

In many cases of interest (e.g., equilibrium structure, thermodynamics) the

computational limitations of ab initio quantum methods can be avoided through a purely

classical representation of the solid. This approach entails an idealization that has many

variants. It consists of the representation of the true potential energy function U({R, })









by a suitably chosen function Uc ({R }) Consequently, its form does not need to be

computed at each time step and the speed and efficiency of classical molecular dynamics

is not compromised.

The problem with this approach lies in the choice for Uc ({R,, }). In principle, an

exact mapping for some fundamental property such as the free energy can be imposed

F[U]= F[U] (2-17)

where F5 is the corresponding classical calculation. Generally, such exact methods can

be inverted only perturbatively and lead to a sequence of effective many ion interactions

involving increasingly more particles. A more practical method is to assume pairwise

additivity for effective point "atoms"


Uc ({R, ) -> Z Z- V (I R, R, ). (2-18)
2 ,j a 8

The exact determination of the pair potentials V, (| R, Rp I) is now much more

restricted as not all properties of interest have such a representation. Nevertheless pair

properties such as the radial distribution functions g, (| r |) might be used to determine

the pair potentials. As the g, (| r |) are unknown and difficult to calculate, the inversion

is again difficult and not practical. Instead, experimental data often is used to fit a

parameterized functional form chosen for the pair potentials. At this stage control over

the approximation is lost and the method becomes phenomenological.

This phenomenological approach has been and remains a valuable tool of materials

sciences. However, in the context of multi-scale modeling it must be reconsidered

carefully. Instead of comparing bulk properties of a solid and its classical representation,









multi-scale modeling constructs classical and quantum models of subsystems and then

requires fidelity at their interface. This is a severe test of the modeling assumptions in

each subsystem. The approach proposed here confronts this issue directly in the

construction of appropriate pair potentials for the problem considered.

The Composite Quantum/Classical Solid

The use of parameterized classical potentials Uc ({R }) to avoid the time

intensive quantum electron calculations cannot be justified for domains far from

equilibrium, where an accurate solution to the quantum description for the electrons is

required. It is presumed that there is some method for identifying such domains within

the solid under given conditions for which quantum chemical effects should be treated in

detail. The quantum solid then is represented as a composite of two domains, the larger

bulk in which a classical representation is used and a smaller "reactive" domain in which

the original quantum description is retained. The objective of the modeling described

here is therefore to construct the potential function U, such that it gives an accurate

description in both the reactive and non-reactive domains. This has two components, the

determination of a pair potential for the forces on ions in the classical domain, and an

accurate calculation of the charge density for the forces on ions in the quantum domain.

The classical and quantum domains are defined by labeling all ions as classical or

quantum, and spatial domains associated with the coordinates of each, denoted

{Rc, } and {Rq, } respectively. The quantum domains are assumed small, to allow

practical calculation of the electronic structure. In principle there could be several

disconnected quantum domains, but to simplify the discussion we consider only one. It is

assumed initially that the two sets are contiguous with a smooth interface and that









diffusion or migration between them is not significant over the times of interest. In

addition to the ions in the quantum domain, there are n, electrons where n, is

determined by the condition that the quantum domain be charge neutral. The total

average electron charge density then is decomposed as

P, (f, t) = Pe (f, t)(zQ (+) + XZ ()) Peq (, t) + Pec (, t), (2-19)

where Z (?) and Zc (f) are characteristic functions for Q and C. The boundaries of the

quantum spatial domain Q are fixed by the choice of ions and

nD = dr Pe (, t), (2-20)


where Q encloses all {R,,,}. It is to be noted that the boundary of Q is not unique. To

make the calculations simple, we chose the boundary of Q to be some smooth surface

such as a plane, beyond which the charge density falls below a particular value, e.g. 105.

The complement of this domain is the classical domain C. This gives a corresponding

decomposition of the total ion potential energy function for the ions V, (ion -ion

Coulomb interactions) plus U, of 2-14 so that 2-13 becomes

atD + {K, + + ,q),D, = 0, (2-21)

where K, is the ion kinetic energy and


c If dFp '1 ((p) ') (r-+ ('t))+ 1(p')

(2-22)

and

2 r --' r -' (, (p,(y)( )+pqt))+I a (P(2-23)
2 QT \ Q` T-T ~f-/









The first terms of the integrands on the right sides of 2-22 and 2-23 represent the

interactions of the ions with a given subsystem in the presence of the average electronic

charge density of that subsystem. The second terms represent the interactions of those

ions with their complementary subsystem. Half of the ion-ion potential energy between

the two subsystems has been associated with each potential in this decomposition so that

the total force acting on the quantum domain by the classical domain is equal and

opposite to that on the classical domain due to the quantum domain.

The potential V, is due to ions in the classical domain. By definition these ions are

in near equilibrium states and therefore this part of the potential should be represented

well by classical pair potentials of the form 2-18, with appropriately chosen parameters.

The potential V, is due to ions in the quantum domain. As this is the domain that can be

far from equilibrium the average electronic charge densities must be calculated in detail

from the quantum description 2-15. There are two parts to this charge density affecting

the ions in the quantum domain, that due to the nD electrons of the quantum domain

peq (i, t) and that due to the surrounding classical domain p, (r) + pec (f, t). The equation

governing the nD electrons of the quantum domain is coupled to this same classical

domain average charge density. The scheme proposed here consists of modeling this

charge density p'e (f, t) as an accurate and practical representation for the environment

of the quantum domain, allowing calculation of peq (f, t) and therefore determining both

V and V,.









Description of the Classical Subsystem

The steps in constructing the pair potentials are the following. First the specific

quantum method to be used in the multi-scale modeling is identified and applied to a

large cluster or sample of the solid to be modeled. The forces on atoms for both

equilibrium and near equilibrium states are then calculated quantum mechanically. Next,

a simple functional form for the pair potentials is chosen and the parameters are adjusted

for a best fit to the quantum force data. Chapter 4 details the construction of such pair

potentials for the silica nanorod.

This method of construction describes the modeling of the first term on the right

side of equation 2-22 where both the ions and the average electron charge density refers

to the classical domain, assumed in a near equilibrium state. The second term involves a

coupling to ions and average charge density in the quantum domain peq (F, t). This

charge density is computed for the forces in the quantum domain (next section) and hence

the second term of 2-22 is known as well.

In summary, the potential energy for the ions of the classical domain is determined

from synthesized pair potentials among the ions and electrons of the classical domain,

plus an interaction with the electrostatic potential of the ions and average electron charge

density of the quantum domain.

Description of the Quantum Subsystem

The quantum domain is a subsystem of nD electrons localized about the designated

ions defining that domain. Consider the reduced density operator for nD electrons

defined by

(nD) =TrNe fDD (2-24)








More specifically, this partial trace is defined in coordinate representation [Coleman and

Yukalov, 2000] by

(r ,...,ir |D |r',,.P..r' )=DD D | D, |r,. .,dr r'' e ID eN
(2-25)

Clearly the full exchange symmetry among all electrons is preserved. However, this

reduced density operator is not specific solely to the quantum domain defined above. For

example, the diagonal elements (r,..., f, D, ,..,?, ) give the probability density to

find nD electrons at the specified positions, and the latter can be chosen anywhere inside

the system. Thus, only when the positions are restricted to the quantum domain Q does

this reduced density operator represents the electrons of that domain. Similarly, if

p"D) () is the charge density operator for nD electrons its average is

p (rD) (F)= Tre pD (F)D, = TrenD p(F)D) D (2-26)
= = Tre Pe (r)De (2-26)

where the trace in the second equality is over nD degrees of freedom. This average

charge density represents the average contribution of electrons at any point r. Both

(r,.., D I ,.., r,,) and pe)(?) change with r since there is an absolute reference

background set by the functional dependence on the ion charge density. Consequently, in

all of the following discussion of this section (f,..," D D, Ii,.., ) and pe"nD)() are

considered only for positions within the chosen quantum domain. Accordingly this

coordinate representation { ~,..., ? )} defines an nD electron Hilbert space of functions

defined over the volume of quantum domain Q. In this context D nD) becomes the









reduced density operator for the quantum domain and pe(D') () its charge density

operator

nD) D,, p"D () ()> p (r). (2-27)

The equation determining the reduced density operator for the quantum subsystem

follows directly from this definition and equation 2-15 for De

I
S[(K, + Ve),Dq = 0. (2-28)


where Ke is the kinetic energy for the nD electrons and


Tv 1 rpe) (P) ( o(Y' 0) + (P, (d')+ P (' 0)
1 1 I 1
2 Q r | er Ir- r
(2-29)

The first term on the left side of 2-29 describes the isolated quantum subsystem,

composed of the Coulomb interactions among the electrons and their coupling to the

average charge density of the ions in the quantum domain. The second term is the

interaction of these electrons with their environment, the total average charge density of

the classical domain.

There are two distinct types of contributions from this charge density of the

classical domain. The first is associated with a subset of ions at the border of the QM/CM

domains which describe chemical bonds in the full quantum solid. These ions locate

regions where there is a highly localized electron charge density shared with the quantum

domain, including both strong correlation and exchange effects, localized and non-

uniform. The second type of contribution is associated with the atoms beyond next

nearest neighbors for which the quantum correlation and exchange effects are much









weaker, and the dominant effect is that of a polarized charge density due to the ions and

electrons.

So the total charge density for the classical region can be written as

PC() -/Pec (F')+P, (F) = Z (| k, I)P (f)+ Apc (f) (2-30)
Rcv

where v denotes the set of border ions in the environment for which a bond has been

broken in identifying the quantum domain. Also, ( r'- R, ) is a characteristic

function specifying a domain centered on R, such that it does not overlap neighboring

ions. Its size is taken large enough to incorporate the bound electrons forming the "atom"

for this nucleus in the quantum solid. The second term Apc (F) is the charge density for

all remaining ions and electrons of the classical environment

Coulomb Effects of Environment

Since by definition, Apc (F') does not include any of the chemical bonding with the

quantum domain required for its valency, the electrostatic potential associated with it can

be expected to have a regular multiple expansion

__1 r .d
fr' Apc () r-> f + (2-31)
rr r

The leading monopole term is zero due to charge neutrality of the classical subsystems

(the integral excludes the border ions). The dipole moment d for the entire environment

is

d = IJi'Ao (F') (2-32)

These results are still quite formal but they provide the basis for the

phenomenology proposed for this part of the modeling: replace all effects of the classical









environment on the quantum system, exclusive of bonding, by an effective dipole

representing the polarization of the medium by the quantum domain. The origin of this

dipole in the above analysis shows that in general it will depend on the geometry and the

state of both the classical and quantum domains (e.g., the dipole will change under

conditions of strain). In the phenomenological application of this prescription the dipole

must be supplied by some simpler means since the electronic contribution to A,p () is

not known.

Electron Exchange with Environment

The contributions from pc (r) in the regions where bonds have been cut require a

more detailed treatment. Clearly, a necessary condition is that valence saturation should

be restored in the quantum domain. However, that is not sufficient to assure the correct

charge density or the correct forces within that domain. Instead, the charge density near

the border ions responsible for bonding should induce a realistic charge density within

the quantum domain. To see how this can be done first write the contribution from one

such ion as


i f (-I f'-R a )P () -p (f) (2-33)


The potential ,pa () represents the actual electrostatic and exchange effects of the

ion of the border on the quantum domain. These effects are comprised largely of the ion

at the site plus its closed shell electrons distorted by exchange and correlation effects of

the site in the quantum domain with which it is bonding. Consequently, an appropriate

pseudo-potential is introduced at each such border ion whose behavior in the direction of









the quantum system is the same as that of the actual ion. These ions plus their pseudo-

potentials will be called "pseudo-atoms".

To accomplish this, an appropriate candidate for the pseudo-atom is chosen based

on valency of the particular pair of border ion and its neighbor in the quantum domain.

Next, a small cluster containing that bonding pair is chosen for training the pseudo- atom.

The bond is then broken and the relevant member of the pair replaced by the pseudo-

atom. The training consists of parameterizing the pseudo potential to give the same forces

as in the original cluster. For example, the pseudo-atom might consist of an ion plus its

closed shell electrons chosen to satisfy the valency for the bond broken. The adjustable

parameters refer to a characterization of the closed shell electron distribution. In this way,

it is assured that the pseudo-atom not only gives the correct saturation of the dangling

bond but also reproduces the forces within the cluster and hence gives a realistic

representation of the charge distribution between the chosen pair. The primary

assumption is that the bonding of interest is a local effect, so that the pseudo atom trained

in the cluster can be used in the bulk solid with a similar accuracy. The training depends

on the particular quantum method used to describe the solid, and is considered in more

detail in chapters 3 and 6.

In summary, the environmental effects on the quantum system are accounted for

approximately by a dipole representing its polarization and a collection of pseudo-atoms

located at the sites of ions where bonds have been cut


Vq ( d Pe(r) d r'(r)-(r- r) + P, ()) +r2 + I p ( (2-34)
2 1 I f rI aiv









This model now allows solution to 2-28 for the electron distribution in the quantum

domain, including its coupling to the classical environment. Then peq (r) is calculated

from 2-26 and 2-27. Finally, the desired potentials V,, and V of 2-22 and 2-23 are fully

determined


S( R )i N, (p (') + p-eq ('',t)), (2-35)
V,, a I I f, r


and


V, 1f d ( ) (p, (F')- 3( -') +eq(,t))+ p(r) (2-36)
2 | r -r' r acv

The solution to the classical ion motion 2-21 can then proceed via MD simulation.

This completes the proposed scheme for multi-scale modeling of the idealized quantum

solid. The verification of this modeling in the silica nanorod and related structures is

considered in chapters 3 through 6.














CHAPTER 3
DESCRIPTION OF THE QUANTUM DOMAIN

Introduction

This chapter illustrates the technique to meet the first criterion for the multi-scale

modeling, proposed in chapter 1, which is to obtain the correct charge densities and

forces in the QM domain when isolated from its CM bulk. This requires incorporating the

information about the state of the CM environment into the QM domain. The CM region

has two kinds of effects on the QM domain: a short-range electronic exchange interaction

and long-range Coulombic interactions. Most previous studies of QM/CM simulations

neglected the long-range interactions. However it will be shown here that both kinds of

interactions must be taken into account for an accurate description of the QM domain.

Also the effects of the CM environment should be modeled in a simple enough manner

that it can be implemented in the quantum calculations without significantly increasing

the computational time. A method for such modeling is shown in this chapter. In doing

so, the chapter describes: (i) the quantum mechanical tool employed to describe the QM

domain, (ii) an overview of various termination schemes used in the literature to treat the

bond cutting region, (iii) the method proposed here for construction of a pseudo-atom to

saturate dangling bonds and to account for short-range exchange interactions, and (iv)

modeling the rest of the CM environment beyond the next neighbor interactions with

lower order multipoles. Finally we compare the charge densities and forces in the QM

domain obtained from our scheme of pseudo-atoms and dipoles with those obtained from

the conventionally used link atoms (bond saturation with hydrogen atoms) and dipoles.









The next section describes the different choices for the quantum mechanical

method used to obtain the reference data against which our method of multi-scale

modeling is to be judged.

Quantum Methods

At the quantum level, each possible choice for the ab initio method entails different

approximations used to optimize the accuracy and speed of the calculations. Each

approximation has its own advantages and limitations, but a general embedding scheme

should be insensitive to these. To test how robust the proposed embedding scheme is, we

chose two different methods for the underlying quantum mechanical description. The first

method is the Transfer Hamiltonian (TH) method [Taylor et al., 2003] based upon the

Neglect of Diatomic Differential Overlap (NDDO) approximation. The second method is

the Born-Oppenheimer density functional theory (DFT) using the generalized gradient

approximation (GGA) [Barnett and Landman, 1993]. Chapters 3 through 5 give the

results obtained from the TH method; chapter 6 shows results of the DFT method. The

next section briefly introduces the TH strategy as the details are published elsewhere

[Taylor et al., 2003]. The DFT method is summarized in chapter 6.

Most of the previous work on multi-scale modeling is based on tight binding (TB)

models [Broughton et al. 1999; Rudd and Broughton, 2000] because of their ability to

treat hundreds of atoms in the QM domain. However, the Hamiltonian in TB models,

once formed, is not solved self-consistently. It is oversimplified, and in particular, it

cannot account for charge transfer, which is essential to bond breaking. The underlying

Hamiltonian in TB models is an extended Huckel Hamiltonian to which a short-ranged

repulsion has been added.









Transfer Hamiltonian (TH): The goal of the TH strategy is to provide forces for

realistic MD simulations of the quality of the coupled cluster singles and doubles (CCSD)

method [Taylor et al., 2003] in computationally accessible times. Bartlett [1995] gave a

comprehensive review of coupled cluster theory and its application to chemistry. The

relevant aspects here are that CC theory is the present state-of-the-art method for high

accuracy predictions on small molecules and that the method is computationally very

costly. The TH is a low-rank, self-consistent, quantum mechanical single-particle

operator whose matrix elements are given in terms of parameterizable functions. The

functional form of the TH can be chosen to be of any semi-empirical functional form but

here the NDDO (Neglect of diatomic differential overlap) form is chosen because of its

better qualitative description over Intermediate Neglect of Differential overlap (INDO)

Hamiltonians [Hsiao et al., 2001]. The NDDO approximation restricts charge distribution

to basis functions on the same center, so this has at most two-atom interactions. Appendix

A details the NDDO method. Because of the NDDO Hamiltonian form it is anticipated

that the parameters are saturated for small cluster size and that these parameters can be

then transferred to extended systems.

An overview of the derivation of TH from coupled cluster theory is as follows. In

coupled cluster (CC) theory, the wavefunction is written as

IF)=e' e0) (3-1)

where is a single Slater determinant and T is an excitation operator that creates single,

double, and so on up to nth excited states. T is given by


T = T, + T2 +.... = K tr a +1 a tab +... (3-2)
,a 2 l,j,a,b









where a, and a, are the standard creation/annihilation operators. In principle, the CC

theory considers all possible electron excitations and correlations like post-Hartree Fock

theories. However, for simplicity we retain only terms with single and double excitations.

This is a standard working approximation in contemporary small molecule computation.

The independent particle TH can be connected to many-body CC theory by inserting

equation 3-1 into the Schrodinger equation, yielding

e THe To)= E o) (3-3)


Ht o)= E o) with H= e HeT (3-4)

The energy and the forces given by the gradient of the energy obtained from equation 3-4

are exact up to the level of correlation included (here single and double excitations). The

ground state eigenfunction of the effective Hamiltonian H is a single Slater determinant,

but it has electron-electron correlations corresponding to the truncation of 3-2 employed.

To apply the TH to silica system for studying complex processes such as fracture,

or hydrolysis where the role of electrons is important, Taylor et al. [2003] parameterized

the TH for the following reactions.

Si2H6 -- 2SiH3

}" "Si(OH)3 + SiO(OH)3 (radical mode of dissociation)
Si207H6 -
S Si(OH)3 + -OSi(OH)3 (ionic mode of dissociation)


Since the TH models two different dissociation pathways for Si207H6, electron

state specificity is thus systematically introduced into TH. Therefore MD driven by TH

forces is expected to predict fracture realistically because in fracture a variety of

dissociations always occur. The TH thus is particularly trained to drive high-quality MD









results while requiring only the computational intensity of semi-empirical methods which

are orders of magnitude less demanding than CC calculations on the same systems.

We will assume that this TH strategy serves as the framework within which we

investigate our scheme for the embedding a QM region inside a CM region. In the current

application, we have presumed that the forces that arise from the TH are the reference

forces by which our proposed method for multi-scale modeling is to be judged. The next

section shows the procedure for accurately representing the QM domain when separated

from its CM environment.

Description of the QM Domain

Description of the QM domain requires incorporating the relevant effects of the

bulk CM environment on the reactive QM region. There are two qualitatively different

effects of the environment. The first is the short-range interaction at the QM/CM

boundary, where the partitioning involves cutting covalent bonds between the two

regions. The second is the long-range Coulomb interaction between the rest of the bulk

and the QM region. As will be seen, both effects must be accounted for to describe the

QM region accurately. A successful QM/CM treatment needs a strategy to treat: (i) the

bond-cutting region, called termination schemes, and (ii) the long-range interactions

through embedding.

Different Termination Schemes

Over the past decade there have been numerous proposals for different types of

termination schemes. The description of all of these is beyond the scope of this work.

Only a few of the more relevant ones are discussed in some detail.

(i) Link Atom method: The predominant scheme used is the link atom (LA)

method, first presented by Singh and Kollman [1986]. In this method, hydrogen atoms are









added to the CM side of broken covalent bond to satisfy the valency of the QM system.

There are many variations within the implementation of the LA method, for example, the

double Link Atom method [Das et al., 2002], the Add-Remove Link Atom method

[Swart, 2003] or the Scaled-Position-Link-Atom Method (SPLAM) [Eichinger et al.,

1999]. Reuter et al. [2000] reported that in all cases the LA should be initially placed at

approximately 0.97 A from the QM atom. One can choose from two different options

depending on whether the link interacts with the CM atoms or not. The type of links that

do not interact with the CM, referred to as QQ Links give large errors in energy

computation. The second type of link called HQ Links, when geometry optimizations are

performed leads to incorrect bond length for the frontier QM/CM bond, which in turn can

affect the rotational barrier around the frontier bond. So, one has to remove the

interactions between the link and the classical group at the frontier. But then, to avoid the

electric field problems mentioned above, it is necessary to prevent the frontier classical

group from interacting with all QM atoms. This procedure is recommended only if the

charge of the classical group at the frontier is small. In spite of all the above-mentioned

difficulties, the LA method is still widely used in many types of QM/CM applications

because of its simplicity in implementation.

For the comparison with the method proposed here, we only studied consequences

of the QQ types of LA in detail where the hydrogen atom is placed at a distance of 0.97 A

from the QM atom at the boundary (which is a Si atom in our case) and the interactions

of the LA with the CM atoms are ignored. We also studied the effects of changing the

position and alignment of LA. These results will be discussed in the last section.









(ii) Frozen Orbital method: This approach also is known as the Local-Self-

Consistent-Field (LCSF) method [Assfeld and Rivail, 1996]. In this method, the electron

density in the cut bond is precalculated and then the electronic density along the frontier

bond is represented by a strictly localized frozen atomic orbital, which has a preset

geometry and electronic population. These frozen orbitals are not allowed to adapt during

the QM/CM calculations. The frozen orbitals are defined by their hybridization

coefficients and their electron population. These two parameters are determined by

quantum chemical calculations on small model systems and are assumed to be

transferable to a bigger system. This method avoids the use of any additional atoms but

Reuter found it to be less robust compared to the LA method. It was seen that energy

calculations are very sensitive to QM size and values of CM charges. He compared the

results for proton affinities and deprotonation enthalpies for a number of molecules using

LA and LSCF methods and found the results to be comparable.

(iii) Generalized Hybrid Oribital method: Gao et al. proposed a refinement of this

LSCF concept known as the Generalized Hybrid Orbital (GHO) Method [Gao et al.,

1998; Pu et al., 2004] in which the electron density in the cut bond is allowed to readjust.

In this method, the hybrid orbitals are divided into auxiliary and active orbitals, the latter

of which are then optimized in the QM calculations. But both of these frozen orbital

methods require substantial computational work and hence are difficult to implement.

(iv) Connection Atom method: In this method [Antes and Thiel, 1999] a connection

atom is developed to saturate a C-C bond, such that the connection atom mimics the

effect of a methyl group. This connection atom interacts with the QM atom quantum

mechanically and the interactions with the CM atoms are handled classically using a









carbon force field. The parameters of the connection atom are determined using semi-

empirical methods [Stewart, 1990] such as AM1, MNDO or PM3 designed to reproduce

theoretical QM data for energies, geometry and net charges. About 30 different methyl

hydrocarbons were used as reference molecules. The parameters adjusted are the orbital

exponent,;, one-center one-electron energy U,,, one-center two-electron integral g,,,

resonance parameter/P and repulsion term a. The mathematical functional forms of

these parameters can be found in any book on semi-empirical theory [Stewart, 1990] or

papers by Dewar and Theil [1976].

There are numerous other termination schemes for the QM/CM boundary like the

'pseudobond' scheme [Zhang et al., 1999], 'IMOMM' [Meseras and Morukuma, 1995;

Humbel et al., 1996] and 'ONIOM' procedures [Svesson et al., 1996; Maseras et al.,

1999; Dapprish et al., 1999], 'effective group potential (EGP)' [Poteau, 2001]. However

these methods are not discussed here.

Pseudo-atom Method

These are developed on the basis of the previously described TH method for

saturating a bond terminating in Si (in silica systems). In keeping with the TH strategy,

we have based the parameterization of the pseudo-atom by matching forces for the QM

portion of the system. We intend that the short-range interactions, particularly the

exchange interaction, will be taken into account using the pseudo-atom. The following

NDDO parameters were modified to account for the near environment, based on a

molecular model, pyrosilicic acid: one-center-one-electron integrals U,,, U,; Coulomb

integrals g,., g,, g,; exchange integral h,; two-center-one-electron resonance

integrals fP, tp [Dewar and Theil, 1976]. As already mentioned, because of the NDDO









form it is anticipated that the parameters are saturated for small cluster size and that these

parameters can be then transferred to extended systems.



.. .... .................

*/ ,... Modified

.. .- ... .... ...






Figure 3-1. Training of Pseudo-atom on smaller molecules, here Pyrosilicic acid.


Figure 3 -1 shows the parameterization scheme for the pseudo-atom on pyrosilicic

acid. The part of the molecule within the dotted lines is replaced by a Fluorine (F) atom

whose NDDO parameters are then adjusted to give the correct QM forces (which implies

correct geometry as well) and charge density in other parts of the molecule (outside the

dotted lines). The F atom is placed at the same position as the neighboring CM atom (O

in this case) in the bond being cut and hence geometry optimization as in the LA method

is not required.

This method is similar to the previously described connection atom developed by

Antes and Thiel [1999]. However the advantage of using pseudo-atoms is that those are

trained to give the correct QM forces and charge densities for both the equilibrium Si-0

bond being cut and for small stretches (up to 3-4% from the equilibrium) in the Si-O

bondlength as well. This fact will allow us to use the pseudo-atoms even while studying

dynamics in the system. The actual charge density from the TH method is calculated

instead of the more commonly used Mulliken charge densities [Mulliken, 1955], which









are known to have several common errors ( e.g., a Mulliken charge distribution always

has equal apportioning of electrons between pairs of atoms, even if their

electronegetavities are very different, which leads to quite unrealistic net atomic charge.

In extreme cases some orbitals may contain a negative number of electrons and others

more than two electrons). Also in all the other termination methods reported, the role of

the more remote CM environment is not quite clear.

Calculation of Charge Densities

In the nanorod (1-4), the charge densities of the (i) ring with LA, and (ii) ring with

pseudo-atoms are compared to the charge densities of the ring in bulk (i.e., the charge

density for QM domain resulting from the benchmark QM calculation for the entire rod).

The density is calculated on a planar grid placed parallel to the plane of ring (which is

taken to be at z = 0). Visualization of the charge densities provides some insights into the

descriptive differences between the pseudo-atom and LA. Figures 3-2(A), 3-2(B), 3-2(C)

and 3-2(D) show the charge density of the isolated ring, the ring with LA, the ring with

pseudo-atoms, and the ring in bulk respectively in the plane of the ring. The six red spots

(high density) on the contour of Figure 3-2(A) correspond to the six oxygen nuclei of the

ring. It is noted that the charge density found from superposition of single atom

calculations of oxygen and silicon atoms is not same as that of the isolated ring. The six

blue spots (low density) in the contour of Figure 3-2(B) correspond to the LAs. In the

contour of Figure 3-2(C), an overlap between the spots indicates the bonding between the

nuclei similar to the ring in the bulk [Figure 3-2(D)]. These plots also indicate that

pseudo-atoms are a more realistic representation of the bulk compared to the LAs

because they take into account the bonding of the ring.












































Figure 3-2. Charge densities with different termination schemes. (A) Isolated ring (B)
with links (C) with pseudo-atoms (D) in actual nanorod.


Table 3-1 shows the normalized square difference of the charge density for these

two terminations with respect to the ring in bulk. Note that as the z coordinate gets larger

we approach the pseudo-atom location where we would not expect the termination

scheme to reproduce the 'exact' result, hence we stop at z = 0.8 A. The right-most

column of the Table 3-1 shows that, while the difference is least for the ring with pseudo-









atoms, still the discrepancy from the 'exact' result is not less than 1% even in the best

case and the relative agreement gets worse the further away from the ring we look.

Hence, even though the pseudo-atom provides reasonable chemical bonding behavior,

there is some significant effect from the rest of the system.



Table 3-1. Comparison of charge densities of the ring with LA and ring with pseudo-
atoms with that of the ring in the bulk calculated at planes parallel to the plane
of the ring.
Density in Normalized difference
Distance ring in
from plane of bulk Ring with
ring (in A) elections / Ring with Pseudo-
Link-atoms
Volume) atoms
Z=0.2 942 0.05 0.02
Z=0.4 507 0.12 0.05
Z=0.6 246 0.32 0.19
Z=0.8 171 0.42 0.42


The forces on the nuclei of the ring with LA and pseudo-atoms were compared to

those of the ring in bulk (the benchmark). Figure 3-3 shows the magnitude of the forces

on one of the Si nuclei of the ring resulting from the LA or pseudo-atom termination

schemes. Since the ring in the rod is in equilibrium, there should not be any forces as is

evident from the graph. It is seen that neither the LA nor the pseudo-atom method alone

is sufficient to represent the bulk, and additional contributions from the rest of the rod are

necessary. These are the long-range Coulomb interactions for which the detailed quantum

mechanics is not expected to be important.













1000
:Z 900
800
700
600
500
400
S300
0 200
100

Ring with links Ring with pseudo- Ring in rod
atoms


Figure 3-3. Comparison of magnitude of force on Si nuclei with different termination
schemes.




Embedding Scheme

To improve our embedding procedure, we describe the remainder of the

environment as two dipoles for the top and bottom portions of the rod (Figure 3-4). The

values of the dipole have been calculated using the TH-NDDO charge density for these

two portions of the rod. These domains are taken to be charge neutral, but are polarized

by the presence of the QM domain. The validity of the approximation of the rod by

dipoles has been checked by finding the difference between the force on Si nuclei due to

all the charges and that due to the dipole and it is found to be very small of (about 6.25

kcal/mol/A and is within our computational error bar which is about 50 kcal/mol/A). It is

found that even with dipoles placed externally above and below the ring, the force on Si

nuclei is 364.7 kcal/mol/A. Hence, it is clear that using pseudo-atoms or external dipoles

alone will not lead to a satisfactory embedding.










+ Rest of
Short range by Pseudo-atoms the rod








Y CM Long range interactions
modeled by dipoles t Rest o

-i- the rod

Figure 3-4. Embedding Scheme: Approximating the rest (CM region) of the rod beyond
the pseudo-atoms by dipoles.


The problem of large forces was solved by including the two dipoles in TH-NDDO

calculation to incorporate the effects of polarization of the ring due to the dipoles. Figures

3-5(A) and 3-5(B) show the significant improvement in the values of forces and charge

densities with the incorporation of dipoles respectively. From Figure 3-5(B), it is seen

that the normalized difference of charge density is reduced by an order of magnitude.


800 -*-Without dipoles
6 700 (B)
S- With dipoles (B)
o 600 (A)
500 -
S400 O

S200 0.01
S100 "
0 0 _I0.001 .
pseudo- dipoles dipoles in in c 0.1 0.3 0.5 0.7 0.9
atoms NDDO equlibrium
distance from plane of ring


Figure 3-5. Effect of inclusion of the dipoles in the TH. (A) Improvement in the value of
forces. (B) Reduction in the value of normalized difference of charge densities
Open symbols, without dipoles; squares, with dipoles.









So we observe that the pseudo-atoms together with the dipoles in the TH-NDDO

can reproduce the charge densities and the forces in the QM domain accurately as in the

presence of the bulk CM. This scheme was found to be applicable to strained

configurations as well. We studied several non-equilibrium cases including:

i) When the ring (the QM domain) was radially expanded by 5%, ii) A distorted

ring, in which one Si atom is radially pushed out and one Si pushed in (Figure 3-6), iii) A

uniaxially strained rod, and iv) A longer rod with ten rings in equilibrium and cases i) and

ii) in this longer rod. We have compared the force on one of the Si nuclei (for the

distorted ring case we have chosen the force on the Si atom pushed out) in the QM

domain (i) with our scheme, (ii) with LA plus dipoles and (iii) the actual force in the

molecule. Figure 3-7 compares the forces for all the cases mentioned above.

Similarly the charge densities have been compared near the plane of the ring in the

same fashion as described above. Table 3-2 shows that the percentage charge density

with respect to bulk calculated with our method in each case lies below 1%.


Figure 3-6. Cross-sectional view. (A) Distorted ring. (B) Ring in equilibrium.












700

600

500
0

U 400
U
300

200

100


rod in equal 5% expanded distorted ring uniaxially longer rod 5% exp ring in dist ring in long 3-membered
ring strained rod 10rings long rod rod ring
Figure 3-7. Comparison of forces on Si nuclei with our scheme and LA to those in the
actual molecule for various cases studied.


Table 3-2. Percentage charge difference with respect to bulk calculated with our method
for various cases studied
% charge density
difference with
Various cases studied respect to bulk
calculated with our
method
I For rod with 6 rings.
a) Equilibrium 0.1
b) 5% radially expanded ring 0.5
c) Distorted ring 0.2
d) Uniaxially strained rod 0.7
II. For the longer rod with 10
rings
a) Equilibrium 0.1
b) 5% radially expanded ring 0.7
c) Distorted ring 0.5
III. 3-membered ring 0.8


As a further test, we have used this method to partition a highly strained small

molecule, which is a 3-membered silica ring [Figure 3-8]. These types of ringed

structures are found abundantly on the surface of amorphous silica [Du and Cheng,









2003]. These strained molecules readily hydrolyze in the presence of water and are used

to study hydrolysis for a-Si02 [Du et al., 2004]. Again, the results for forces and charge

density in the QM domain are found to be quite satisfactory (last row of Table 3-2).



QMA


'- CM
Si





Figure 3-8. Partitioning of 3-membered silica ring.



Other Options for Links

In the results discussed above, we found that the pseudo-atoms and the dipoles

reproduce the force in the QM domain within 1% error but when we placed the LA at a

fixed position of 0.97 A from the terminating Si atom and along the Si-O bond, there is a

large force on the Si atom (Figure 3-3). One might argue here that the results with the LA

could be improved if positions of Las are allowed to vary. In our next analysis, the Si-H

bondlength as well as the alignment of the LA was allowed to vary to give minimum

force on the Si atom. It was found that the minimum force on the Si atom is obtained

when the LA is placed at a distance of 1.45 A from the Si atom but is not aligned along

the Si-O bond. Instead it is about 5 degrees deviated from the Si-O bond (Figure 3-9).

The blue atoms are the LAs; the green atoms are the O atoms that are still included in the

picture to show the alignment of the Si-O bond. It is seen that the alignment of the Si-H

bond is a little bit deviated from the original Si-O bond.









Though the LA at this position gives forces in the QM domain comparable to our

scheme with the pseudo-atoms plus the dipoles, it fails to reproduce the correct charge

densities as given by the pseudo-atoms. Figure 3-10 shows the charge density with the

LA at these deviated positions. The plane at which the charge density is calculated is

same for which charge density was calculated with the pseudo-atoms earlier. Here, we do

not see any overlapping region between the six spots on the contour. This indicates the

LA cannot account for the bonding as successfully as pseudo-atoms.


Figure 3-9. Deviated position of optimized LA.


Figure 3-10. Charge density with LA at optimized positions.


L









LA on O atoms instead of on Si: Another option for putting the links on O atoms

was tried. These link atoms were placed at a distance of 0.97 A beyond the six

intermediate O instead of on Si in the nanorod (Figure 3-11). This O-H bond essentially

behaves as an F atom having a valency of -1. It was seen that the LA at this position

reproduced the forces in the QM quite correctly. However the role of the external field in

this case then becomes unclear as a reduction in the magnitude of forces on the Si atom

was observed when the two dipoles were added. We consider again the example of a 5%

expanded ring and compare the forces when the links are placed on the O atoms and then

with the dipoles. Figure 3-12 shows that the difference of force between when the LA is

placed on O atoms and that of the actual molecule is only 52 kcal/mol/ A, but when the

dipoles are included there is a significant decrease in the magnitude of force on Si and the

difference with the actual molecule becomes more than 170 kcal/mol/ A.












Figure 3-11. Position of LA beyond intermediate oxygen atoms.


This is somewhat artificial because the CM environment beyond the next neighbor

also plays an important role in determining the correct properties of the QM domain, an

effect that is anticipated to become more prominent in bigger and unsymmetrical

systems. Thus we note that in none of the cases studied does the LA represent a realistic

picture.














700


600


500- actual
our scheme
0 links on Si+dipoles
0
E 400- o links on 0
links on + dipoles

2 300


200


100


0
actual our scheme links on Si+dipoles links on 0 links on 0 +
diDoles
Figure 3-12. Comparison of force on Si nuclei with LA on oxygen to LA on silicon and
our embedding scheme.




Conclusions

A tool for embedding a QM domain in its classical environment was developed in

this chapter. This procedure requires incorporation of the information about the CM

region into the QM domain in a simple fashion without making the quantum calculations

more complicated.

Both kinds of interactions, short-range as well as long-range, between the QM and


CM were taken into consideration. The short-range exchange interactions were modeled

by replacing the atoms at the boundary of QM/CM by pseudo-atoms constructed to

describe the chemistry of the nearest neighbor interactions. The pseudo-atoms are F









atoms with electronic properties adjusted to yield the correct charge density and forces in

the QM part and these pseudo-atoms are placed at the same distance and angle as original

O atoms from the Si-O bond to be cut. The second kind of interactions i.e., the long-range

many-body Coulombic interactions were represented to a sufficient accuracy using

dipoles.

The generality of the proposed scheme extends to systems such as strained rods,

longer nanorods and small molecule such as a 3-membered ring as well.

The pseudo-atom method is a better way for saturating the dangling bond than the

conventional LA method. As can be seen from the charge density plots, the pseudo-atoms

take chemical bonding into account more or less as in the actual system unlike the

behavior of the LA. Several other possibilities for the alignment of LA were tried, but

each of those failed to represent the actual charge density in the QM domain.














CHAPTER 4
CONSTRUCTION OF A NEW CLASSICAL POTENTIAL FOR THE CM DOMAIN

Introduction

After developing the tools to satisfy the first criterion of multi-scale modeling in

the previous chapter, we focus on the second criterion, that is, to obtain a classical

potential for the CM environment which yields the same geometry and non-equilibrium

elastic properties as predicted by the QM calculations. Independent treatment of the QM

and the CM regions can result in mismatching of many physical quantities such as forces,

charge densities, and elastic constants across the QM/CM boundary. To avoid this

mismatch, the interatomic forces used in the molecular dynamics (MD) simulation in the

CM region must be chosen to preserve selected, crucial properties of the underlying

quantum theory used in the QM domain. The process of constructing such a potential is

demonstrated here. It is developed using ab initio data on the equilibrium structure and

on weakly strained configurations calculated from the quantum description (TH), rather

than the more usual approach of fitting to a wide range of empirical data. The latter

approach is appropriate when a classical description of the entire material is reasonable

(e.g., near equilibrium) so that maximum correlation between the material simulated and

the real material is attained. In contrast, for multi-scale modeling of strain to failure, the

primary constraints on the potential are those of internal theoretical consistency between

the quantum and classical descriptions. The final accuracy of the model is therefore set

by the quality of the quantum description, not by the fitting of the potential to

experiments. The difficulty in extending the conventional approach to the multi-scale









modeling problem has been studied by colleagues [Al-Derzi et al., 2004; Flocke et al.,

2005; Zhu et al., 2005]

The potential is chosen to have the same functional form as TTAM (equation 4.6).

A initial guess for parameters is made using a genetic algorithm with force data obtained

from a quantum calculation (TH-NDDO). The algorithm effectively locates a domain of

the relevant parametric minimum for the given constraints; the parameters found are then

checked for stability, followed by a global scaling (described in Appendix C), resulting in

the final choice for the parameters that determine all the pair potentials. It will be shown

that the Young's modulus (Y) obtained from this classical potential matches closely with

that obtained from the QM method for strains up to 6%, unlike the standard TTAM for

which Y differs by almost 20% from the TH. Furthermore, the bond lengths and bond

angles in the rod are an order of magnitude more accurate for the new parameterization in

comparison to those from the current TTAM or BKS potential parameters.

In the first section we discuss briefly the basic concepts ofMD used in the program

to obtain the stress-strain curve. Allen and Tildesley [1987] and Haile [1992] detail the

MD simulation methods. Next we show the failure of standard pair potentials to give the

correct value of Young's modulus for the nanorod. The method for developing a new

potential followed by the structure and elastic properties obtained from the new potential

are then described. The possible wider applications of this potential to silica systems

other than the nanorod are considered in chapter 7.

Molecular Dynamics

Corrector-Predictor Algorithms

Molecular Dynamics allows us to study phenomena under extreme environmental

conditions, which may be difficult to study experimentally. In molecular dynamics,









successive configurations of the system are generated by integrating Newton's laws of

motion.

av
IF(t)= mr,(t) (4-1)


where F (t) is the force acting on the ith particle of mass m due to all the other particles at

time t. The potential, V describes the interaction among the particles.

These equations cannot be solved analytically for a large number of particles, and

instead are typically solved using finite difference algorithms [Lamberti, 2000]. The most

commonly used finite difference methods are the Verlet method [1967] and Predictor-

Corrector algorithms. All algorithms assume that the dynamical properties (velocities and

accelerations etc.) can be approximated as a Taylor series expansion in time. There are

many algorithms for finding the motion of the particle. We used the Gear predictor-

corrector algorithm [Gear, 1971]. The predictor corrector algorithm consists of three

basic steps. First, new positions, velocities, accelerations and higher-order terms are

predicted according to the fifth-order Taylor expansion equations. In the second stage, the

forces are evaluated at the new positions to give the accelerations. These accelerations are

then compared with the accelerations predicted from the Taylor expansion. The

differences between the predicted and calculated accelerations are then used to correct the

positions, velocities etc. in the correction step.

The molecular dynamics simulations were done at constant temperature. There are

many methods for controlling the temperature for example velocity rescaling [Woodcock,

1971], Anderson's method [Anderson, 1980], and the Nose-Hoover thermostat [Nose,

1984; Hoover, 1985]. Velocity rescaling is used in our calculations.









Velocity Rescaling

This is the simplest and widely used method for maintaining temperature in a MD

simulation. The temperature of a system is related to its average kinetic energy given by


(K.E.)N = NkBT (4-2)
NW 2

The obvious way to change the temperature is thus to change the velocity. If the

temperature at a time t is T(t) and the velocities are multiplied by a factor A, then the

corresponding change in temperature is

A 1 2m,(2v,)2 1 2m v2
2,1 3Nk, 2, 3 Nk,

AT = (22 1)T(t)

A= lw / T(t) (4-3)

So the temperature can be controlled by multiplying the velocities at each time step by a

factor of A = JT / T where T,,,, is the current temperature calculated from the

kinetic energy and T,, is the desired temperature.

This procedure is implemented in the MD code used for obtaining stress-strains

curves, the results of which are discussed in the following sections.

Failure of Standard Potentials

It is found that standard classical pair potentials for silica fail to yield the correct

value of elastic constants as predicted by the quantum mechanical TH method for even

small strains, let alone predicting fracture in materials. This is not surprising, as these

potentials have been parameterized using data for equilibrium states only. For example,

Figure 4-1 shows the stress-strain curve for the silica nanorod from MD for the entire rod









using the standard TTAM [Tsuneyuki et al., 1988] and BKS [van Beest et al., 1990] pair

potentials (equation 4-6). Also shown for comparison are the reference results obtained

from the TH method for the entire rod.

It is observed that the stress-strain curves given by the TTAM and the BKS

potentials differ significantly from the results of the quantum theory. Examination of the

curvature of the stress-strain curves indicates that for strains up to 5%, the standard

classical potentials predict a stiffer rod than the TH rod. Above 5% strain the nanorod

undergoes more of a plastic deformation (characteristic of ductile fracture) for these than

the TH rod, for which the stress-strain relation remains linear even for higher strains. The

value of Young's modulus (given by the initial slopes of these curves) is therefore

considerably different for each potential. One can see that for strains up to 4%, the

Young's modulus for the TTAM (Y=1214 GPa) rod differs from that of the TH rod

(Y=1026 GPa) by 18% while the difference between the BKS (Y=1516 GPa) and TH rod

is even more than 47%.

These stress-strain curves were generated by the MD simulation in which the

nanorod was subjected to longitudinal strain along the z direction by pulling the end caps

with a fixed velocity. The remaining atoms in the rod were allowed to relax with their

positions determined by MD and the temperature was kept constant by the velocity

rescaling method.

The stress was calculated by taking the sum of the forces parallel to the loading

direction on the constrained atoms divided by the projected cross-sectional area of the

nuclei in the end caps. The cross-sectional area was updated after every step to account









for any motion in x or y direction. The simulation was done at a temperature of 10K with

a strain rate of 25 m/s and time steps of 2 fs.















0 / r m A, nAL
0-






0.00 0.02 0.04 0.06 0.
Stra in

Figure 4-1. Comparison of stress-strain curves for the standard pair potentials (BKS,
TTAM) with that of Transfer Hamiltonian.


Objective and Method for Constructing a New Potential

The foregoing results confirm that existing potentials cannot be used for consistent

embedding of the QM domain in a CM representation of the nanorod. It is proposed here

that each development of a multi-scaling simulation on a specific system should include

the construction of its own classical potential. That construction should be governed

primarily by two components: the chosen quantum method used in the QM domain, and

the particular physical properties requiring the QM treatment. In addition, it is essential

that there be a limiting domain (usually, near equilibrium) for which the CM and QM

calculations of these properties are indistinguishable.









The proposed method is as follows. A QM method (TH in the present case; DFT in

chapter 6) is chosen to represent the essential physical mechanisms not explicitly

accounted for in a CM representation. Selected properties of the system are then

calculated directly using the QM method. Here, these properties are the forces on nuclei

incorporating non-local electronic configuration and exchange effects, both for

equilibrium and weakly strained conditions. Next, a functional form for the potential

energy is chosen for classical point particles that substitute for the atomic constituents.

The parameters of the potential energy function are determined by some optimization

method based on fitting the data calculated from the QM method. A key issue is the

choice of potential energy function. For practical MD simulation it is useful to restrict

attention to two-and three-body potentials, but the functional form of these few-body

potentials is still optional. Researchers have constructed potentials having as many as 42

parameters [Watanabe et al., 1999]. Other simpler silica potentials having three-body

terms exist, Tersoff [1998], Vashishta et al. [1990], Jin et al. [1993], Biswas and

Hanmann, [2001], Stillinger and Weber [1985]. There are a few review articles on

modeling of silica using different potentials [Schaible, 1999; Brenner, 2000]. Brenner

discussed two aspects for the development of an effective interatomic potential energy

functions. The first is the derivation of sound functional forms that are motivated by

quantum-mechanical bonding principles. The second is the development of empirical

corrections and fitting parameters that are often necessary to make a potential function

practical for specific applications. Our method is in fact based on both of these aspects.

The functional forms that are considered frequently have their genesis in the theory

of intermolecular interactions. However, their application to materials simulation is a









rather speculative endeavor. In fact, one can argue that there is no fundamental origin for

such potentials since the underlying electronic origin is both non local and many-body in

nature. Instead, the potentials are viewed as constructs of the multi-scale model "trained"

to give selected material properties but with no unique status otherwise. Here, for

simplicity, only pair potentials are considered and their functional forms are taken to be

the same as the phenomenological forms currently in use for silica but with undetermined

parameters. Al-Derzi et al., [2004] showed based on high-level QM calculations on

several HxSiyOz clusters and a -quartz that the parameterization logic behind TTAM

parameters is not capable of providing a consistent embedding. We have parameterized

the potential only on this particular nanorod as we interested in its properties. However,

in chapter 7 we shall see that this potential does have applicability on other silica

material.

The objective here is to develop a new classical potential that yields the correct

equilibrium structure (bond lengths, bond angles) and small strain elastic properties, as

given by the quantum TH theory. Such a potential then can be used to meet the third

criterion of our multi-scale method, which is to construct the composite nanorod (i.e., the

QM and CM regions joined together) with compatible small-strain elastic properties

throughout. We emphasize that an attempt to fit the entire stress-strain curve is not

reasonable as the physics of large-strain domains is inherently quantum mechanical, thus

beyond the limitations of a simple classical potential. Classical potentials cannot account

for charge transfer or distinguish between different paths of dissociations: neutral or ionic

mode; both of these phenomena are essential to bond breaking.









We assume a classical pair potential that has the same functional form as the

TTAM and BKS potentials, but with new parameters determined as described below.

Despite its known limitations, this form is chosen because of its simplicity in

implementation and widespread use. The potential energy between i and j ions is given

by

q, q c
V (r) =q +a exp(-bj)- 6(4-4)
r r

The first term is the Coulomb interaction for ions of charge q, and q,. The

remaining two terms are collectively called the "Buckingham" term. They model the

short-range repulsive and dipole dispersion, or van der Waals, interactions, respectively.

The Coulomb term is determined by a single parameter, the charge on silicon or oxygen,

since there must be charge neutrality for the SiO2 molecule (i.e., 2q, = -q, ). There are

all together 10 parameters, consisting of the charge qs, and three pair parameters a,, b,

and c, for each pair of interactions (0-0, Si-O and Si-Si).

In earlier work, the parameters for the TTAM and BKS potentials have been chosen

by fitting ab initio Hartree-Fock potential energy surfaces of SiO4(4-)+4e+ and Si04H4

respectively. Furthermore, the Si-Si interactions in the BKS potential are arbitrarily set

equal to zero. The resulting parameter values are given in Table 4-2. However it is noted

that this fitting was done with a low level Hartree-Fock approximation and the energy

surface was explored in equilibrium regions only. Here, we choose the parameters in a

quite different way. Since MD is using only the forces between ions, the forces

determined directly from the TH quantum mechanics of the entire nanorod are used as the

primary reference data. Ideally, these forces should be zero for the correct equilibrium









structure and its stability. Furthermore, to assure good elastic properties, information

about these forces for small strained conditions also is used in the search for appropriate

parameters.

The standard approach for parameterization consists of following steps: (i) choose

properties of interest (total forces on atoms), (ii) choose reference data (values of forces

for equilibrium and strained configurations), (iii) perform some estimation of the

parameters, (iv) test against reference data with some error function and then (v) choose a

variational method for the next estimate and convergence of the procedure.

In our case, the reference quantities are the TH quantum forces on the atoms of the

rod in equilibrium, and when strained from its equilibrium structure. The error function

L, for testing the quality of any chosen set of parameters is the weighted sum of absolute

difference between the forces calculated from the trial potential and the reference forces,

denoted by g(i) and f(i) respectively for the forces on the ith ion of the rod.


L = Z ak (k (i) (4-5)
k i

The index i runs over all ions of the rod, while k indexes the sum over data sets for

different structural configurations of the rod. The numbers ak denote the weights

assigned to the different data sets with ak = 1. The magnitude of the weight decreases
k

with the increasing strain, with the maximum weight assigned to the equilibrium case. A

typical distribution of weights was ak = (0.5, 0.2,0 .1, 0.1, 0.05, 0.05) for percent strains

(0, 1, 2, 3, 4, 5), respectively. It is to be noted that this error function is not limited to the

form shown; one could choose other forms such as taking the square of difference or

normalizing the difference etc. but it was found that these forms gave similar results. This









error function is minimized in a space whose dimension is the number of adjustable

parameters. The problem of such a high dimensional domain is the existence of several

local minima and maxima. Gradient optimizing routines [Nocedal, 1999] seek extrema in

the neighborhood of an initial chosen point without sampling the whole space and hence

are only local in scope. This characteristic is illustrated for a simple two-dimensional

domain in Figure 4-2, where a search near the second minimum would miss the global

minimum.













f(x,y)






Figure 4-2. Presence of local and global minima.



Further improvement can be sought only through some sort of random or unbiased

restart. One can repeatedly restart the search algorithm at other points but as the

dimension of the space increases such a process becomes highly inefficient. Therefore a

genetic algorithm [Goldberg, 1989] has been used since it does not suffer from such

disadvantages. Furthermore, the genetic algorithm (GA) does not require calculation of









derivatives and hence is more efficient. We have used the standard PIKAIA code

[Charbonneau, 1995] for the GA. Appendix B details the GA method.

However, GAs are not without problems. Once a neighborhood of an extremum is

reached the GA converges slowly because it relies primarily on the mutation rate to

generate small incremental changes in the population. If the global minimum has not

been located the search will begin to stagnate. This is due to the successive removal of

comparatively less fit individuals from the population, which, after many generations,

results in a population of parameter sets, which are all more or less equally fit hence the

evolution of the parameters stalls. This problem is analogous to that of the gradient

methods getting trapped in local minima. At this point one might think of switching to

Newton methods, that is to construct hybrid schemes of GA and gradient-based methods

using GA results as a starting guess. Instead we have used an overall spatial rescaling

(Appendix C) of the GA parameters to improve the predicted elastic properties (e.g.,

Young's modulus) as the final step in determining the choice of potential parameters.

Results

Finding Parameters for the Potential

The parameters obtained from the GA are checked to see whether they produced a

stable rod. This is done using a Broyden-Fletcher-Goldfarb-Shanno [Zhu, 1994] energy

minimization simulation at 0 K temperature. In this method, the function to be minimized

(energy in our case) is expanded to second order in terms of parameter set p (different

geometries of the rod in our case) as

1 (4-6
f(p) ;f(po)+g, x(p po)+-(p-po)T x H x (p po) (4-6)
2









where g is the gradient of the function and H is the Hessian. The derivative of the

function is calculated numerically. To have the geometry with minimum energy requires

the determinant of Hessian to be positive. The optimization proceeds with the generation

of a new geometry p given an old geometry p, using (p po) = H -g This process

continues until f goes below some acceptable tolerance. As anticipated, it was found

that there are numerous local minima in the parameter space, and several possible sets of

parameters were found from the GA that yielded a stable rod. Table 4-1 lists these

different parameters.



Table 4-1. Different values of potential parameters obtained from GA which gave a stable
rod
Paramet
I II III IV V VI
ers
qsi 2.36 2.78 2.68 2.22 2.4 2.4
qo -1.18 -1.39 -1.34 -1.11 -1.2 -1.2
aoo 1.78E+03 2.01E+03 1.87E+03 1.9E+03 4.31E+04 1.45E+05
boo 2.80 2.48 2.68 2.64 4.27 5.91
Coo 239.05 558.28 326.69 366.13 11.64 0.19
aosi 1.06E+04 1.24E+04 1.10E+04 1.02E+04 2.62E+05 4.11E+03
bosi 4.86 4.13 4.65 5.04 7.19 4.32
Cosi 64.55 200.46 87.83 50.26 3.87 22.4
asisi 8.66E+08 1.02E+09 9.06E+08 8.3E+08 2.34E+04 1.86E+05
bsisi 15.33 13.07 14.66 15.6 14.28 6.6
Csisi 22.08 67.6 30.18 16.4 1.05 8.15


The units of ai, bi and c, are eV, (A)-1 and eV-( A)6 respectively.

Table 4-1 shows that the c, s vary immensely from one set to another, while the

other coefficients are of the same order of magnitude for set I-IV. Also the a, s are more

than an order of magnitude different for set V-VI. All this shows that though these

parameters correspond to very different points in the parameter space, they all gave rods









with similar geometry. The length of the rods obtained from these parameters differed

from the TH rod by 7-12%. Set VI was found to give least structural error in bond

lengths and bond angles. It is noted that for this set, the a, and c, differ from the other sets

by an order of magnitude or more while the values of the b, and charges are similar to

those of other sets. Set VI above was found to give most accurate rescaled parameters.


Table 4-2. Parameters for New Potential in comparison to TTAM and BKS potential
parameters
Parame New
TTAM BKS
tersM potential
qsi 2.4 2.4 2.24
qo -1.2 -1.2 -1.12
aoo 1756.98 1388.773 133037.6
boo 2.846 2.760 6.18
Coo 214.75 175.00 0.133
aosi 10722.23 18003.757 3767.451
bosi 4.796 4.873 4.514
Cosi 70.739 133.538 15.754
asisi 8.73E+08 0 1.698E+05
bsisi 15.22 0 6.90
Csisi 23.265 0 5.70


The units of aj, by and c, are eV, (A)-1 and eV-( A)6 respectively

Table 4-2 shows the values of these rescaled parameters, defining the new potential

developed here, in comparison to the usual TTAM and BKS values. Note that the values

of the charges and b, for the new potential are similar to those of TTAM or BKS, but that

the c, are always much smaller for the new potential. This suggests a more repulsive

character of this potential. Figure 4-3 compares the potentials for pairwise Si-O

interactions.











40"



3C








C



C _

-1


TTAM
STH Claial oe
TH rldsisiil potential


- I I I I


Figure 4-3. Comparison of Si-O pairwise interactions for the three potentials. (The inset
shows the zero crossings of the potential)


Table 4-3 shows the comparison of several bond-lengths for the three potentials

and those for the reference TH rod.


Table 4-3. Comparison of structure of the rod obtained from the different potentials.
Bond lengths(A) New
Bond angles TH PtTTAM BKS
(degrees)Potential error error error
(degrees)
In silica planes
Si-O 1.641 1.642 0.04 1.65 0.76 1.616 1.5
Between planes and
intermediate O's
In the end caps
Si-O 1.71 1.68 1.75 1.67 2.5 1.62 5.1
Approx. length 16.49 16.69 1.19 16.31 1.1 15.86 3.8
Approx. diameter 6.55 6.57 0.35 6.57 0.35 6.41 2.0


R I.III,'2 11lllls)










Stress-Strain Curves

This new potential was used in MD simulations to obtain the stress-strain curves

for a classical model of the nanorod. The simulation was carried out under the same

conditions as for the TH (temperature of 10 K and pulling the end-caps at 25 m/s). Figure

4-4 shows the results for this new potential, together with those for the TTAM and BKS

potentials, in comparison with the TH results. The stress-strain curve obtained from the

potential agrees well with the TH curve up to about 5% strain.

Accordingly, the Young's modulus obtained from the new potential agrees very

well with that from the reference TH. The MD simulations were repeated at lower

temperatures and also at lower strain rates. Figure 4-5 plots the stress-strain curve for

several temperatures using the new potential. It is observed that the wiggles gradually

decrease with lower temperature and then vanish at 0.2K. This result was confirmed by

doing an adiabatic expansion of the rod (0.0K) as well, leading to the same results.

60-


50


40-
a3_

S30- Value of Y in (GPa)
V ^ TH, Y=1026.1
20 0 New Potential, Y = 1021.8
A TTAM, Y = 1214
x BKS, Y= 1516



0-
0.00 0.01 0,02 0.03 0,04 0.05
Strain

Figure 4-4. Comparison of Young's modulus at low strains for different classical
potentials with that of TH.




















W",




Pr: + 50K
10K
5K
4+ -o 2K

0.2K


I
0.00


0.01


' I
0.02
Strain


.03
0.03


' I
0.04


Figure 4-5. Reduction of the wiggles with temperature.


0.00
0.00


0.02 0.04 0.06


Strain
Figure 4-6. Stress-strain curves at lower strain rate of 2m/s and a temperature of 0.2K.


0.05
0.05


60






-- 40

20
0
(0

65

20









Absence of wiggles in the TH curve: Classical potentials give numerous wiggles in

the stress-strain curves while that for TH is smooth at any temperature. No satisfactory

reasons were found for the cause of wiggles. The first anticipated reason was the

possibility of large fluctuations in the radius of the nanorod while pulling, as the x and y

coordinates of atoms were allowed to relax. However it was found that the fluctuation in

the radius was less than 0.03A.

The second reason thought was the presence of many lower lying phonon modes

for the classical potential compared to TH. Thus, the phonon spectra for all the three

classical pair potentials (new, BKS and TTAM) were compared to that of TH method.

For classical potentials, the phonon spectrum is obtained by solving the eigenvalue

problem for the harmonic oscillator while for the TH the eigenvalues are obtained by

diagonalizing the Hessian matrix (finite differences of gradients to get the second

derivatives of the Hamiltonian).

The following histogram presents the number of modes vs. wavenumber in an

interval of 100cm-1 (up tol600 cm-1) for the TH (green), TTAM (red), BKS (pink) and

new potential (blue). At the lowest wavenumber 100cm-1 we see that number of modes of

TH is highest (-28) compared to BKS (- 12), TTAM (-17) or new potential (-27). This

is contradictory to our assumption of lower-lying phonon modes for the classical

potential. Hence we see neither fluctuation in the radius of the ring nor the phonon

spectra are the causes for wiggles in the stress-strain curves for the classical potential.

The exact cause of wiggles yet remains an open problem to study.











50
50 -------------------------------------
TTAM
TH
45 NEW CLASSICAL
BKS
40

35






20 -

15 :: ,' ,' I
I s --- i i .......... ...i

[ ]I I ii i I I
I :




Ii I I I I I

0 200 400 600 800 1000 1200 1400 1800 1800


Figure 4-7. Histogram of number of modes vs. wave number (cm-1).


Improving Stress-Strain Curves for Higher Strains

Figure 4-6 shows that the classical potential works only up to 5% strain. Beyond

that, the stress curve for the classical rod becomes more ductile in nature as is evident

from the bending of the curve. As mentioned, the classical potential is not expected to fit

the entire TH stress strain curve, since the physics of large-strained domains is inherently

quantum mechanical. However, it is desirable for multi-scale modeling that the CM

domain should grow less ductile with strain than the QM domain so that the material fails

only in the QM domain. There are two different methods to make this classical potential

more rigid than the TH. The first technique is to induce a charge-transfer term in the

potential, and the second is to go back and re-parameterize the potential by including the

information about the forces for higher expansions of the nanorod in the GA. We shall









discuss only the results of the latter technique in this chapter; the charge-transfer potential

is described in Appendix D.

We repeated the process of parameterization using QM force data up to 10%

expansions of the nanorod in the GA. The data for forces at the end caps are excluded as

the bonding at end-caps became significantly different from the bulk at higher strains and

it ends up in misleading the GA, which gives equal weight to bulk as well as end cap

data. It is observed that the information about the bulk force data is self sufficient to

generate the end caps. The disadvantage of this approach is that incorporation of more

strain data sets in the GA results in a delay in its convergence and although it gives a rod

with higher value of Y, the accuracy in the bond lengths is lost (by 3%) so the rod has to

be rescaled more than once. The following rescaled parameters (Table 4-4) were found to

give a rod with similar geometry (Table 4-5) as that of the TH rod and a stress-strain

curve (Figure 4-8) that matches with the TH curve well up to 5% but becomes stiffer

afterwards.



Table 4-4. New parameters for classical potential-2 (NTH-2 Classical Potential)

New-TH Old-TH
Parameters Classical classical
potential # 2 potential # 1
qsi 2.05 2.24
qo -1.025 -1.12
aoo 39439.87 133037.6
boo 4.66 6.18
Coo 5.85 0.133
aosi 240101.9 3767.451
bosi 7.88 4.514
Cosi 2.06 15.754
asisi 27530.45 1.698E+05
bsisi 21.42 6.90
Csisi 1.55 5.70







69


The units of aj, bi and c, are eV, (A)- and eV-( A)6 respectively

We term this new potential the NTH-2 classical potential (New-TH-Classical

potential # 2).


Table 4-5. Structure obtained from NTH-2 Classical Potential in
TH nanorod


25


20


1C


NTH-2
Bond lengths (A) and H Class %
TH Classical
bond angle (degrees) Potential error
Potential
In Silica Planes
Si-O 1.641 1.642 0.04
Between Planes
End Caps
Si-O 1.71 1.67 2.5
Length 16.49 16.34 0.9
Diameter 6.55 6.54 0.15



0-


0 -
/ m





M fall OM(TH)
0 -X- I ITH-2





0
,"


comparison to that of the


0.00 0.05 0.10


0.15 0.20
Strain


0.25 0.30 0.35


Figure 4-8. Behavior of NTH-2 Classical Potential near fracture.


15









As can be observed from Table 4-5, the accuracy of the bond-lengths and bond

angles (compare Table 4-3) is not lost but the rod becomes much stiffer than the TH after

5% strain. This stiffness will guarantee that in the multi-scale treatment the large strain

behavior will be localized in the QM domain where it is treated properly with accurate

electronic methods.

Conclusions

The key to successful multi-scale modeling lies in the consistent embedding of a

QM domain (TH in our case) in its classical MD region. This requires the CM region to

have the same structure and elastic properties as the QM domain. A method is presented

for finding a classical potential that would meet the above criteria for the silica nanorod.

A pair potential having the same functional form as the TTAM potential is chosen

but unlike the standard method of fitting the parameters of the potential to empirical data,

the parameters are fitted to ab initio QM force data for the equilibrium as well as for

small strained configurations of the particular system under study. A GA followed by a

spatial scaling procedure is used to determine these parameters.

The potential obtained this way is found to be more ductile compared to the TH

method beyond 5-6% strain. The classical potential can be made much stiffer than the TH

by reparameterizing the GA using QM force data until 10% strain. Appendix D describes

an alternative approach to stiffen the classical rod.

The NTH-2 Classical Potential found using the strain data up to 10% is a good

candidate for multi-scale modeling of this rod where a seamless coupling between the

QM and CM regions of the nanorod is required. Since this potential is much stiffer than

the QM method, the phenomena involving large strain will be localized only in the QM






71


domain. We have essentially found a classical representation of the TH silica nanorod

that includes small strain behavior.

Further application of this potential to other silica systems can be found in chapter

7.














CHAPTER 5
CONSTRUCTION OF THE COMPOSITE NANOROD

Introduction

We developed a classical pair potential (NTH-2 Classical Potential) for the CM

region (chapter 4) that predicted the structure and elastic properties in agreement with

those of the QM domain obtained from the TH. The potential thus found is a good

candidate for multi-scale modeling of this rod. In this chapter we will use this potential to

construct a composite rod that is indistinguishable from the TH or the classical rod in

terms of equilibrium geometry and Young's modulus. This indistinguishability is the

final criterion of our proposed multi-scale modeling method. Such a composite rod

allows us to study multi-scale phenomena for the events involving large-strains localized

to the QM domain. Multi-scale modeling done in this fashion will have much of

computational efficiency of classical pair potentials but will retain the accuracy of ab

initio methods.

The composite rod is built by embedding the QM region in its CM environment as

described in chapter 3. The forces on the atoms in the QM domain are then given by the

TH method while the forces on the atoms in the CM region are calculated by taking the

gradient of the NTH-2 Classical Potential. While calculating the TH forces for the atoms

in the QM region, the effects of the CM environment are incorporated with pseudo-atoms

and dipoles as described in chapter 3. Results for the elastic properties of the composite

rod are reported in this chapter.









Results for the Composite Rod

Figure 5-1 shows the stress-strain behavior of the nanorods obtained from three

different methods: (i) the full quantum mechanics TH method, (ii) full NTH-2 Classical

Potential and (iii) the composite rod constructed as described above. The three overlaying

curves (measured at 0.01K) indicate that the composite rod is identical to the rod

obtained from TH or the NTH-2 nanorod in terms of small strain elastic properties and

structure. The stress-strain results shows the success of our multi-scale method indicating

that the composite rod is indistinguishable from the underlying quantum mechanics for

states near equilibrium..

However the elastic properties of the composite rod obtained from both of these

potentials do not agree beyond 10% strain. This is because the pseudo-atoms, trained at

regions only close to the equilibrium configuration, fail to give the correct charge

densities at such high strains. This diagnosis was checked by comparing the charge

density in the QM domain of the 12% strained rod with that of equilibrium configuration

and a difference of 6% was found.

It is conceivable that retraining of the pseudo-atoms may improve the stress-strain

performance of the composite rod beyond a strain of 10% but our next goal is to find an

application with the method developed so far in some real silica systems: for example

cristobalites and glasses. It is be noted that the 10% strain here corresponds to

approximately 10% strain in the bond-lengths, so the failure of pseudo-atoms at high

strains does not pose any constraint so far as application of the proposed multi-scale

scheme to glass is concerned since glass fails much before 10% strain in bond length is

reached. The reason is that real systems like glasses have many inherent defects which

act as stress concentrators that cause the material to break at much lower strains than










observed for the nanorod. One should note that the yield stress of this nanorod (Figure 5-

1) from the TH method is around 190 GPa, which is notably larger than any observed

silica system probably because of its high aspect ratio (one-dimensional system) and

defect free nature.

250





150-


100 ll .A1
10o I -x- NTH-2
,4' --A-- composite
50 -, -"-N



0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
strain

Figure 5-1. Comparison of stress-strain of the composite rod with those obtained from
NTH-2 Classical Potential and TH.

The satisfactory construction of such a composite rod completes our recipe for the

multi-scale modeling for all practical purposes.

Conclusions

In this chapter, we have satisfied the last criterion of our proposed multi-scale

modeling scheme, namely to construct a composite nanorod that would be

indistinguishable from the TH rod for all relevant properties (structure and elasticity) in

states near equilibrium. Such a composite rod can be used for multi-scale studies with

substantially reduced computational times.














CHAPTER 6
INDEPENDENCE OF MULTI-SCALE MODELING WITH RESPECT TO THE
CHOICE OF UNDERLYING QM

Introduction

As mentioned in chapter 3, a proper multi-scale modeling procedure should be

independent of the choice of underlying quantum mechanical method. To test the method

we have developed, we have essentially repeated the work done in chapter 3 and 4

(embedding of the QM domain and finding a classical potential) with Density Functional

Theory (DFT) instead of the Transfer Hamiltonian as the quantum mechanical

approximation. In the field of DFT, 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

beyond the scope of this thesis. Hence only a short review of the basic concepts relevant

to this chapter, such as Born-Oppenheimer approximation, Hohenberg-Kohn theorem,

Kohn-Sham equations, pseudopotential, and dual space formalism will be presented.

More detail on DFT can be found in Parr and Yang (1989), Wimmer (1997), Jones and

Gunnarsson (1989), and Dreizler and Gross (1990).

The primary results of this chapter are a confirmation that the isolation of the QM

domain with pseudo-atoms and dipoles, plus the construction of a classical potential

based on the DFT forces leads to accuracies of the same quality as those described

already using the TH quantum mechanics. Hence, the multi-scale modeling scheme is

faithful to the chosen form for the underlying quantum mechanics in both cases.









Review of DFT Method

Born-Oppenheimer Approximation

This approximation follows from the fact that the ions are much heavier than the

electrons, so the electrons react instantaneously to the position of the ions. Thus the ions

can be treated adiabatically, leading to the separation of the electronic and ionic

coordinates in the many-body wave functions. This reduces the many-body problem to

the solution of the dynamics of the electrons in a frozen configuration of ions. The total

energy of a system consisting of ions and electrons can be reduced to


Etotal({R },{p(r)}) = I -m, I R I2 + Z Zz, Z/ R,- R, I+Eeec[p(r)] (6-1)
i I>J

where R, m and Z1 are the position, mass and charge of the Ith ion. The first two

terms correspond to the kinetic and Coulomb interaction energy respectively. Eeze[p(r)]

is the ground state energy of the electrons evaluated for a particular configuration {R, } of

ions which we consider next.

Hohenberg-Kohn Theorems

The DFT approach is based on the first Hohenberg-Kohn (HK) theorem [1964],

which states that all the ground state properties of a system are functionals of the electron

density. The second Hohenberg-Kohn theory enables us to write the total energy as a

variational functional of electron density p.

Eeec (P)= p(r)v(r)dr + F[p] (6-2)

where


F[p]= T[p]+V,,[p]


(6-3)









T[p] is the kinetic energy V,,[p] is the electron-electron interaction,. The ground state

electron density is the density that minimizes E[p] subject to the condition

Sp(r)dr= N (6-4)

N being the total number of electrons.

F[p] is a functional of the density and is independent of the external potential but

no general exact expressions for it are known. Kohn and Sham (1965) invented an

ingenious approach to the kinetic energy functional T[p] that turned the Hohenberg-

Kohn variational functional into a practical tool for rigorous calculations. The next

section shows the derivation of the Kohn-Sham equations from the HK theorems.

Kohn-Sham Equations

The Functional F[p] can be decomposed into three different components as


F[p]= T, [p]+ J[p] + Ex [p] (6-5)

T, [p] is the kinetic energy of a non-interacting inhomogeneous electron gas in its ground

state, J[p] is the classical coulomb energy (or Hartree energy) and the functional Ex [p]

is the exchange-correlation energy. Ex [p] in fact includes exchange, self-interaction,

Coulomb correlation terms as well as the KE increment T[p] T [p]. Using equations 6-

2 and 6-3 the total energy functional can be written as

Eee[P] = v(r)p(r)dr + T [p] + J[p] + Ex [p] (6-6)

Minimizing this expression subject to equation 6-4 gives

ee- u = 0 p being a Lagrange multiplier that enforces 6-4
5p(r)









S= v T + (6-7)
e p(r)

where veff is called the effective KS potential defined as

dJ[p] +E,,
Vef = v(r) + +
dp(r) dp(r) (6-8)
(6-8)
= v(r)+ p(r') dr'+vx (r)
SIr-r'|

From the definitions of the non-interacting system p(r) = 1f |, (r) 12 (f is the


electron occupation number) and T, = f, *V20 the variation in p is equivalent


to variation of the hence one gets (in Hartree)

1
[- V2 +efr)]V g (6-9)
2

These are known as Kohn-Sham equations [1965]. They have to be solved self

consistently as vef depends on p(r) .Though the KS equations incorporate the exact non-

interacting kinetic energy T [p], the exact form of the exchange correlation functional is

still unknown. Many approximations for E, [p] have been developed. The ones used for

our calculations are discussed below.

Local Density Approximation

The first such approximation proposed by Kohn and Sham is known as the Local

Density Approximation (LDA) which uses the basic idea of using the exchange

correlation energy of the homogenous electron gas to evaluate Ex[p]. Thus


(6-10)


Exc = Exc (r)p(r)dr









where Ex, is the exchange correlation energy per particle of a uniform electron gas of

density p. This model assumes the electron density to be locally equivalent to a uniform

gas, hence the corresponding exchange-correlation contribution to equation (6-10)

becomes

LDA 8ELDA a (p)
v r (r) = -(p(r)) + p(r) (6.11)
8p(r) Op

Considering the inexact nature of the LDA, it is remarkable that calculations using

LDA have been so successful. Generally, the total energy differences between related

structures are accurate to within a few percent. However there are many failures of LDA

also. For example, it can predict incorrect cohesive energies, often over binds molecules

and solids, and it fails to explain certain classes of molecular bonds (e.g., hydrogen

bonds) for which the charge density in the binding region is very small. Therefore we

have used a more refined approximation for Ex known as the Generalized Gradient

Approximation.

Generalized Gradient Approximation

The LDA is based on the expressions for the energy density of a homogeneous

electron gas. For slowly varying densities, the energy functional can be expanded as a

Taylor series in terms of the gradient of the density, with the first term being the LDA

term. Retention of the first two terms is known as the gradient expansion approximation.

However the series is shown to be nonconvergent [Kleinman and Sahni, 1990]. The

approximation can be made physically sensible, however, and several such

approximations have been proposed which now are known as generalized gradient

approximations (GGA) (e.g., the PBE named after Perdew, Burke and Ernzerhof [1992]).









In many applications, the GGA shows improved results over LDA in total energies

(Perdew etal., 1992), atomization energies (Becke, 1992; Perdew etal., 1992, 1996;

Proynov et al., 1995), energy barriers (Hammer et al., 1993; Hammer and Scheffler 1995;

Philipsen et al., 1994) etc. Typically, GGA favors the density inhomogeneity more than

LDA. However, there are notable exceptions [Trickey, 1997].

Throughout this discussion we have neglected the spin of the electron. The code we

used employs Local Spin Density Functional Theory (LSDFT). In this theory, electron

and the spin densities are fundamental quantities with net spin density being the

difference between the up-spin and down-spin electrons:

(r)= p(r) -p(r) (6-11)

The total electron density is just the sum of the densities for the two types of

electron. The exchange correlation functional is different for the two cases, leading to a

set of spin polarized Kohn-Sham equations


v- dr, + [r, :,+(r )= (r) = ac,/ (6-12)
S2 \iRr)2dr

So we have two sets of wavefunctions, one for each spin.

Troullier-Martin Pseudo-Potential

Implementation of all-electron DFT with a plane wave basis set is computationally

too intensive. This section introduces the concept of pseudopotential, which reduces this

computational cost. The pseudopotential replaces the effect of the chemically inert core

states upon valence states. This replacement saves the effort of using a large number of

plane waves to expand both the tightly bound core orbitals and the rapid oscillations of









the valence Kohn-Sham orbitals in the core region. The pseudopotential can be obtained

from all-electron atomic calculations by solving the radial KS equations self consistently.


Sd + ) +V[p; r] rR,,(r)= srR,2(r) (6-13)
2 dr2 2r2

where

-Z
V[p;r] = +V [p;r]+ vj (p), VH [p;r] being the Hartree potential
r

The pseudopotential has to meet certain conditions: (i) the pseudo-wave function

generated from this potential should not contain any nodes (it should be smooth and not

have any wiggles), (ii) the normalized radial function of the pseudo-wave function should

be identical to the normalized radial real wavefuntion outside a cut-off radius r,, (iii) the

charge enclosed within r, should be equal for the two wave function, and (iv) the

valence of all electron and pseudopotential eigenvalues must be equal. A pseudopotential

which meets all the above conditions is said to be a "norm conserving" potential

[Troullier and Martin, 1991]. Once the pseudowave function is obtained, the screened

pseudopotential can be obtained by inversion of equation

V, () 1(1 + 1) 1 d2
V (r) = R(r1 + )Pd [rRP (r)] (6-14)
2r 2rRP (r) dr

The screening of the valence electrons is removed and an ionic potential is generated by

subtracting the Hartee v PP(r) and the exchange correlation VPP (r) potentials calculated

from the valence pseudo-wave functions from the screened potential

V,~ .= V, () PP (r) ( () (6-15)









Also an accurate pseudopotential should reproduce the scattering of the ionic core.

This requires a nonlocal pseudopotential, that is, one which uses a different potential for

each angular momentum component of the wave function. The pseudopotential operator

can be decomposed into local and nonlocal parts,

VPP (r) = VPP, +c + iVPP, nic (r)P (6-16)


where VPP'c () is the local potential and VIPPnlc (r) is the non-local potential. P, is the

operator that projects out the th angular momentum component from the wave function.

In principle, the local potential can be chosen arbitrarily but since the summation in

equation 6-16 has to be truncated for some value of the local potential should be chosen

such that it adequately reproduces the atomic scattering for all the higher angular

momentum channels. These forces are then used at each MD step to move the ions on the

BO potential energy surface.

Dual Space Formalism

The computer code we used employs the dual space formalism for calculation of

the DFT energy. This method evaluates the kinetic energy in momentum space and

potential energy in real space (Martin and Cohen, 1988; Wentzcovitch and Martin, 1991).

It is found that this method is computationally more efficient than the pure momentum-

space formalism (Ihm et al., 1979) in which all calculations are performed in momentum

space.

Computational Details

All the theory discussed above is implemented into a parallel multi-scale program

package, known as Born-Oppenheimer molecular dynamics "BOMD" [Barnett and

Landman, 1993]. The separation of time scales between the electron and the ionic motion









in the BO approximation restricts the dynamical evolution of the ionic system to a single

electronic potential surface (PES). For a given ion configuration the PES and forces on

the ions are calculated. These then are used to govern the classical equations of motion of

ions advancing them by a small time increment to a new configuration. Evaluation of

Hamiltonian matrix elements and the operations on the wave functions are performed

using the dual-space formalism. A plane wave basis set and cut-off energy of 30.84

Rydbergs is chosen.

Results and Discussion

We next use the DFT to study the multi-scale modeling of the nanorod. The

accuracy our multi-scale scheme is tested against charge densities and forces obtained

using DFT on the whole nanorod. In the next section, we discus the training of the DFT-

pseudo-atoms that model the nearest neighbor interactions.

Training of DFT-Pseudo-atoms

A pseudo-atom is constructed based on the Troullier-Martin (TM) pseudopotential

discussed above. A cut-off radius of 1.5 A is used. Once again the pyrosilicic acid

molecule is chosen for parameterization of the F atom and its position is constrained to be

at the same place as the O atom (Figure 3-1.). But unlike the TH-NDDO method in which

we could change the electron-ion as well as electron-electron interaction parameters, in

DFT we can modify only the electron-ion interactions. The three options that can be

varied to alter the electron-ion interactions for the F atom using the TM pseudopotential

are: (i) the core charge on F, (ii) omission of the non-local part in the potential, and (iii)

switching the local and non-local part between s and p orbitals. We tried all three possible

choices and their combinations to find the best alternative in the sense of reproducing the









force on the terminated Si atom in the pyrosilicic acid as in the presence of the whole

molecule.



Table 6-1. Magnitude of force on the terminated Si atom in pyrosilicic acid with
changing parameters of the dft-pseudo-atom
Only
Snin Without non- Switching
zcore changing
z e local part between s & p
zcore
7.0 0.88 2.9E-04 0.08
6.8 0.97 0.04 1.56
6.2 1.15 1.15 1.87
7.2 0.72 0.81 0.16
7.4 0.62 0.78 0.13
7.6 0.76 0.88 0.74
the units of forces are in Rydberg a.u



This method of parameterizing a pseudo-atom by trial and error might seem

arbitrary but as it is often said "parameterization is not a science but an art". Our case of

obtaining a pseudo-atom truly exemplifies this quote. We started with the same core

charge (7.0) as in real F, then changed the core charge in steps of 0.2 while neglecting the

non-local part of the potential in the second case and interchanging s and p orbitals in the

third case until the force on the Si atom is minimized.

Table 6-1 shows that the minimum force is obtained when the core charge is 7.0

and the non-local part is omitted. The force for this particular set of combination is lower

by two orders of magnitude than for the other cases. This pseudo-atom is then transferred

to the nanorod. The nanorod is partitioned in QM and CM regions and the CM portions

are approximated by dipoles as before (Figure 3-4). The values of the dipoles are

calculated from the charge density obtained from the DFT calculation for the CM

portions of the rod. We have computed the force on the Si atom (in the QM domain with