UFDC Home  myUFDC Home  Help 



Full Text  
MULTISCALE 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, HaiPing 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 programmingmaze as I did. I am also grateful to Professor HaiPing 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. LingLing Wang, Dr. Moahua Du, Chao Cao, and YingXao 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 heartfelt 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 proofreading 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 ultiScale M odeling .................................. .................. .................. Exam ples of M ultiScale 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 seudoatom 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 CorrectorPredictor 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 StressStrain C urves .......................... .... ...... .............................. .. .......... 64 Improving StressStrain 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 MULTISCALE 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 BornOppenheimer Approximation...................... .......................... 76 H ohenbergK ohn Theorem s ........................................... ......................... 76 K ohnSham Equations.............. ............... ................ ................. ............... 77 Local Density Approximation ............................. ................................. 78 Generalized Gradient Approximation ...................................... ............... 79 TroullierMartin PseudoPotential.................. .............. ..... ............... 80 D u al Sp ace F orm alism ............................................................. .....................82 Com putational D details ........................ ................ ................... ..... .... 82 Results and D discussion ........................... ..... ..... ...... .. ............. 83 Training of DFTPseudoatoms...................... ..... .................. 83 Construction of ClassicalDFT Potential for the CM Region .............................86 C o n clu sio n s..................................................... ................ 8 9 7 APPLICATIONS OF THE MULTISCALE 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 NewDFTClassicalPotential ..........................................104 The 10Membered Silica Nanoring ......... ......... .....................104 BulkGlass ............... ......... ................. 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 ChargeTransfer Potential for Old Transfer Hamiltonian Classical Potential #1.....131 ChargeTransfer Classical Potential for NewDensity 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 31. Comparison of charge densities of the ring with LA and ring with pseudoatoms with that of the ring in the bulk calculated at planes parallel to the plane of the ring ................................................................................. 39 32. Percentage charge difference with respect to bulk calculated with our method for various cases studied ...................... .................. ............................ 43 41. Different values of potential parameters obtained from GA which gave a stable ro d .............. .................... ..................................... ......... ...... 6 1 42. Parameters for New Potential in comparison to TTAM and BKS potential p aram eters ........................................................................... 62 43. Comparison of structure of the rod obtained from the different potentials ..............63 44. New parameters for classical potential2 (NTH2 Classical Potential)...................68 45. Structure obtained from NTH2 Classical Potential in comparison to that of the T H nanorod ........................................... ........................... 69 61. Magnitude of force on the terminated Si atom in pyrosilicic acid with changing param eters of the dftpseudoatom ........................................... .......................... 84 62. 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 63. Parameters for NewClassicalDFT Potential in comparison to standard TTAM and BK S potential param eters........................................................ ............... 87 64. Structure of the nanorod obtained using NewClassicalDFT, BKS potentials versus structure of D FT m ethod ................................................................... ......87 71. Comparison of structure of the nanoring obtained from TH, NTH2 Classical Potential and TTA M potential.......................................... ............................ 92 72. Structure comparison of 2membered and 3membered silica rings.........................94 73. Comparison of different properties of the glass obtained from the NTH2 Classical Potential to those of the BK S glass ............... .............. ..................... 101 74. Comparison of structure of the nanoring obtained from DFT and NewDFT C classical P potential ................................................................... .. .... .. 105 76. Comparison of different properties of the glass obtained from NewDFT Classical Potential to those of BK S glass......... ................................... .............105 81. Comparison of structure of the nanorod obtained from DFT and TH quantum m eth o d s ......................................................................... 1 0 8 82. Comparison of the 10membered nanoring using two different types of DFT m methods ...................................................... ........ .. ......... ........... 110 83. Comparison of structure of bulk glass using two different classical potentials ......111 LIST OF FIGURES Figure page 11.Use of different methods at different length scales in crack propagation. ...................2 12. Chemomechanical polishing of silica wafers ........................................................3 13. 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 14. Partitioning of the nanorod into QM and CM domains. Localized (blue) electron cloud of the CM region shows the appropriateness of such partitioning.................10 31. Training of Pseudoatom on smaller molecules, here Pyrosilicic acid....................36 32. Charge densities with different termination schemes. (A) Isolated ring (B) with links (C) with pseudoatoms (D) in actual nanorod ..............................................38 33. Comparison of magnitude of force on Si nuclei with different termination sch em e s. .......................................................... ................ 4 0 34. Embedding Scheme: Approximating the rest (CM region) of the rod beyond the pseudoatom s by dipoles. ............................................... .............................. 41 35. 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 36. Crosssectional view. (A) Distorted ring. (B) Ring in equilibrium.........................42 37. Comparison of forces on Si nuclei with our scheme and LA to those in the actual molecule for various cases studied ......................... ......... .. .............. 43 38. Partitioning of 3membered silica ring........ ............... ........ ..... ............ 44 39. Deviated position of optimized LA. ............................................... ............... 45 310. Charge density with LA at optimized positions. .............................................. 45 311. Position of LA beyond intermediate oxygen atoms. .................. ....... ............46 312. Comparison of force on Si nuclei with LA on oxygen to LA on silicon and our em bedding scheme e. ................................................ ...............47 41. Comparison of stressstrain curves for the standard pair potentials (BKS, TTAM) w ith that of Transfer H am iltonian ........................................ ........................ 54 42. Presence of local and global m inim a ........................................ ...... ............... 59 43. Comparison of SiO pairwise interactions for the three potentials. ..........................63 44. Comparison of Young's modulus at low strains for different classical potentials w ith that of TH ........................................................................64 45. Reduction of the wiggles with temperature..........................................................65 46. Stressstrain curves at lower strain rate of 2m/s and a temperature of 0.2K ...........65 47. Histogram of number of modes vs. wave number (cm). ............................... 67 48. Behavior of NTH2 Classical Potential near fracture..............................................69 51. Comparison of stressstrain of the composite rod with those obtained from NTH 2 Classical Potential and TH ............................................................................. 74 61. Comparison of force on the Si nucleus with different termination schemes and dipoles for the D FT case. ............................................... .............................. 85 62. Comparison of adiabatic stress curves for DFT, newDFTclassical and TTAM potential ..................................... .................. .............. ........... 88 71. The 10m embered silica nanoring. ........................................ ........................ 91 72. Stressstrain curve for the nanoring obtained from TH and NTH2 Classical Potential ..................................... .................. .............. ........... 93 73. The 2 and 3m embered silica rings. ................... .......... .......... ............ .....94 74. Silica Nanoclusters. A) Si12024. B) Si18036. C) Si24048Cage. D) Si24048 F ullerene ................... .................................... ...........................95 75. Front view of the 107 atom notched nanorod showing one defect (missing oxygen) ................... ....... ......................... ....... ................. 96 76. Comparison of stressstrain 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 77. Structure of glass obtained from the NTH2 Classical Potential.............................98 7 8. Radial distribution function of the SiO, 00 and SiSi pairs in the glass.............100 79. Bond angle distribution of OSiO and SiOSi ................................ ............... 100 710. Stressstrain curves for NTH2 and BKS glass. ............................................... 102 81. Comparison of elastic properties for DFT and TH nanorods............................. 109 82. Comparison of elastic properties of the two glasses in the elastic region .............12 D1. Charge transferred from Si to O as a function of distance. .................................... 131 D2. Increase in the stiffness with introduction of the charge transfer term................ 132 D3. Increase in the charge of Si with strain..........................................................132 D4. 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 NTH2 New Transfer Hamiltonian Classical Potential # 2 QM Quantum mechanics / mechanical RDF Radial distribution function TB Tight Binding TH Transfer Hamiltonian TM TroullierMartin 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 MULTISCALE 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 nonequilibrium 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 nearequilibrium bulk. Such methods are referred to here as quantum mechanical/classical mechanical (QM/CM) methods in multiscale modeling. Our study addressed the problem of how to bridge these QM and CM domains by constructing an approximate semiclassical 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 pseudoatoms and dipoles to account for the shortrange 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 multiscale 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 silicabased materials, it should provide guidance for extension to other materials as well. CHAPTER 1 INTRODUCTION Definition of MultiScale Modeling The field of multiscale modeling has opened a new era to computational science for studying complex phenomena such as fracture, hydrolysis, enzymatic reactions, solutesolvent studies, hydrogen embrittlement, and many other chemomechanical 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 multiscale 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 11), 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 multiscale simulations in detail. Figure 11.Use of different methods at different length scales in crack propagation [Broughton et al., 1999]. Examples of MultiScale Modeling (i) Chemomechanical polishing (CMP): combines chemical and mechanical interactions to planarize silica wafer surfaces, using a slurry composed of solvents, wetting agents, other compounds and submicronsized 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 nanoscale, 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 22). 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 12. Chemomechanical 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 multiscale 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 pressurevessels: Another multiscale process in which the phenomenon occurs over different time scales is embrittlement of pressurevessel 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 multiscalemodeling" which means the computation of parameters at smaller scale for its use in more phenomenological models at a larger scale. The other type of multiscale problem, in which the scales are strongly coupled is known as "concurrent multiscale 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 multiscale 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 multiscale 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], postHartreeFock 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, stressfields etc., are important players, and therefore finiteelement (FE) methods [Hughes, 1987] are used to examine the largescale properties of the material. We studied only one particular subclass of multiscale 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) Solutesolvent 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 multiscale 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 multiscale 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 lengthscale 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 13 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 14A). Figure 13. 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 14) 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 14. 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 multiscale 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 multiscale 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 multiscale modeling is introduced by stating what variables are most appropriate for control in the multiscale modeling. First, the quantum description in terms of ions and electrons is given and the limitations for practical implementation are noted. Next, reduced selfconsistent 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 BornOppenheimer limit). The first description constitutes the idealized quantum solid for which the subsequent multiscale 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 "pseudoatoms" 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 (21) where the trace is taken over all electron and ion degrees of freedom. The density operator obeys the Liouvillevon Neumann equation ,D+ [H,D]= 0 (22) 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,, (23) U f dr drd' P (f)p (') (24) 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 highquality quantum methods; however this is not a practical method for application to larger systems. Here only noncrystalline 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,, (25) D, TrD, D, Tr D. (26) Their equations of motion follow directly from equation 22. a 1 t,D, +[H,,D, ]+I(U,D, D,U,) = 0, (27) h h I 1 8AD +[HD,]+ (UD, D,U, ) = 0. (28) 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) , (29) Ue = fcr (f) p,('), p(f ) (Tr p, (')D)(Tr~D) 1 (210) 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 27 and 28 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 27 and 210 to get p, Tr p, (F)D, p, > TePe (F)D, (211) As a consequence, the potentials U, and Ue become Hermitian and the average charge densities are now selfconsistently determined by 27 and 28. Selfconsistency 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 27 becomes ,D, + {(H, +U,[pe]),DD,} 0 (212) 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 212 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 28 is a formidable problem: determination of the dynamics of electrons selfconsistently 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, 212 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 212 for their probability distribution D, D,{R,),P,),t) 8,D, +fH, +U,[pe], D,= 0 (213) The average electron density pe (r)is determined from the ground state solution to 28. = 1_ [He + U,,De] = 0, (215) U, f= r r (r, t)p (f'), p,(f,t)= Tr, p,() D, [p(, t),t] (216) Sr r' The analysis proceeds stepwise. The classical equations 213 are solved analytically in discrete time steps for the atomic coordinates and moment. At each step the electron problem 215 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 coarsegrained 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 multiscale modeling. With MD simulation the solution to 213 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 215 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] (217) 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, ). (218) 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 multiscale modeling it must be reconsidered carefully. Instead of comparing bulk properties of a solid and its classical representation, multiscale 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 nonreactive 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), (219) 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), (220) 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 214 so that 213 becomes atD + {K, + + ,q),D, = 0, (221) where K, is the ion kinetic energy and c If dFp '1 ((p) ') (r+ ('t))+ 1(p') (222) and 2 r ' r ' (, (p,(y)( )+pqt))+I a (P(223) 2 QT \ Q` TT ~f/ The first terms of the integrands on the right sides of 222 and 223 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 ionion 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 218, 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 215. 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 multiscale 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 222 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 222 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 (224) 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 (225) 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 (226) = = Tre Pe (r)De (226) 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). (227) The equation determining the reduced density operator for the quantum subsystem follows directly from this definition and equation 215 for De I S[(K, + Ve),Dq = 0. (228) 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 (229) The first term on the left side of 229 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) (230) 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 + (231) 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') (232) 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) (233) 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 pseudopotential 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 "pseudoatoms". To accomplish this, an appropriate candidate for the pseudoatom 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 pseudoatom 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 pseudoatom 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 pseudoatoms located at the sites of ions where bonds have been cut Vq ( d Pe(r) d r'(r)(r r) + P, ()) +r2 + I p ( (234) 2 1 I f rI aiv This model now allows solution to 228 for the electron distribution in the quantum domain, including its coupling to the classical environment. Then peq (r) is calculated from 226 and 227. Finally, the desired potentials V,, and V of 222 and 223 are fully determined S( R )i N, (p (') + peq ('',t)), (235) V,, a I I f, r and V, 1f d ( ) (p, (F') 3( ') +eq(,t))+ p(r) (236) 2  r r' r acv The solution to the classical ion motion 221 can then proceed via MD simulation. This completes the proposed scheme for multiscale 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 multiscale 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 shortrange electronic exchange interaction and longrange Coulombic interactions. Most previous studies of QM/CM simulations neglected the longrange 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 pseudoatom to saturate dangling bonds and to account for shortrange 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 pseudoatoms 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 multiscale 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 BornOppenheimer 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 multiscale 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 selfconsistently. 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 shortranged 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 stateoftheart method for high accuracy predictions on small molecules and that the method is computationally very costly. The TH is a lowrank, selfconsistent, quantum mechanical singleparticle operator whose matrix elements are given in terms of parameterizable functions. The functional form of the TH can be chosen to be of any semiempirical 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 twoatom 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) (31) 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 +... (32) ,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 postHartree 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 manybody CC theory by inserting equation 31 into the Schrodinger equation, yielding e THe To)= E o) (33) Ht o)= E o) with H= e HeT (34) The energy and the forces given by the gradient of the energy obtained from equation 34 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 electronelectron correlations corresponding to the truncation of 32 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 highquality MD results while requiring only the computational intensity of semiempirical 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 multiscale 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 shortrange interaction at the QM/CM boundary, where the partitioning involves cutting covalent bonds between the two regions. The second is the longrange 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 bondcutting region, called termination schemes, and (ii) the longrange 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 AddRemove Link Atom method [Swart, 2003] or the ScaledPositionLinkAtom 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 abovementioned 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 LocalSelf ConsistentField (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 CC 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,;, onecenter oneelectron energy U,,, onecenter twoelectron integral g,,, resonance parameter/P and repulsion term a. The mathematical functional forms of these parameters can be found in any book on semiempirical 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. Pseudoatom 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 pseudoatom by matching forces for the QM portion of the system. We intend that the shortrange interactions, particularly the exchange interaction, will be taken into account using the pseudoatom. The following NDDO parameters were modified to account for the near environment, based on a molecular model, pyrosilicic acid: onecenteroneelectron integrals U,,, U,; Coulomb integrals g,., g,, g,; exchange integral h,; twocenteroneelectron 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 31. Training of Pseudoatom on smaller molecules, here Pyrosilicic acid. Figure 3 1 shows the parameterization scheme for the pseudoatom 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 pseudoatoms is that those are trained to give the correct QM forces and charge densities for both the equilibrium Si0 bond being cut and for small stretches (up to 34% from the equilibrium) in the SiO bondlength as well. This fact will allow us to use the pseudoatoms 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 (14), the charge densities of the (i) ring with LA, and (ii) ring with pseudoatoms 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 pseudoatom and LA. Figures 32(A), 32(B), 32(C) and 32(D) show the charge density of the isolated ring, the ring with LA, the ring with pseudoatoms, and the ring in bulk respectively in the plane of the ring. The six red spots (high density) on the contour of Figure 32(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 32(B) correspond to the LAs. In the contour of Figure 32(C), an overlap between the spots indicates the bonding between the nuclei similar to the ring in the bulk [Figure 32(D)]. These plots also indicate that pseudoatoms are a more realistic representation of the bulk compared to the LAs because they take into account the bonding of the ring. Figure 32. Charge densities with different termination schemes. (A) Isolated ring (B) with links (C) with pseudoatoms (D) in actual nanorod. Table 31 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 pseudoatom location where we would not expect the termination scheme to reproduce the 'exact' result, hence we stop at z = 0.8 A. The rightmost column of the Table 31 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 pseudoatom provides reasonable chemical bonding behavior, there is some significant effect from the rest of the system. Table 31. 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 Linkatoms 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 pseudoatoms were compared to those of the ring in bulk (the benchmark). Figure 33 shows the magnitude of the forces on one of the Si nuclei of the ring resulting from the LA or pseudoatom 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 pseudoatom method alone is sufficient to represent the bulk, and additional contributions from the rest of the rod are necessary. These are the longrange 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 33. 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 34). The values of the dipole have been calculated using the THNDDO 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 pseudoatoms or external dipoles alone will not lead to a satisfactory embedding. + Rest of Short range by Pseudoatoms the rod Y CM Long range interactions modeled by dipoles t Rest o i the rod Figure 34. Embedding Scheme: Approximating the rest (CM region) of the rod beyond the pseudoatoms by dipoles. The problem of large forces was solved by including the two dipoles in THNDDO calculation to incorporate the effects of polarization of the ring due to the dipoles. Figures 35(A) and 35(B) show the significant improvement in the values of forces and charge densities with the incorporation of dipoles respectively. From Figure 35(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 35. 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 pseudoatoms together with the dipoles in the THNDDO 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 nonequilibrium 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 36), 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 37 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 32 shows that the percentage charge density with respect to bulk calculated with our method in each case lies below 1%. Figure 36. Crosssectional 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 3membered ring strained rod 10rings long rod rod ring Figure 37. Comparison of forces on Si nuclei with our scheme and LA to those in the actual molecule for various cases studied. Table 32. 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. 3membered ring 0.8 As a further test, we have used this method to partition a highly strained small molecule, which is a 3membered silica ring [Figure 38]. 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 aSi02 [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 32). QMA ' CM Si Figure 38. Partitioning of 3membered silica ring. Other Options for Links In the results discussed above, we found that the pseudoatoms 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 SiO bond, there is a large force on the Si atom (Figure 33). 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 SiH 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 SiO bond. Instead it is about 5 degrees deviated from the SiO bond (Figure 39). 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 SiO bond. It is seen that the alignment of the SiH bond is a little bit deviated from the original SiO bond. Though the LA at this position gives forces in the QM domain comparable to our scheme with the pseudoatoms plus the dipoles, it fails to reproduce the correct charge densities as given by the pseudoatoms. Figure 310 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 pseudoatoms 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 pseudoatoms. Figure 39. Deviated position of optimized LA. Figure 310. 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 311). This OH 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 312 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 311. 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 312. 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, shortrange as well as longrange, between the QM and CM were taken into consideration. The shortrange exchange interactions were modeled by replacing the atoms at the boundary of QM/CM by pseudoatoms constructed to describe the chemistry of the nearest neighbor interactions. The pseudoatoms are F atoms with electronic properties adjusted to yield the correct charge density and forces in the QM part and these pseudoatoms are placed at the same distance and angle as original O atoms from the SiO bond to be cut. The second kind of interactions i.e., the longrange manybody 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 3membered ring as well. The pseudoatom 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 pseudoatoms 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 multiscale 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 nonequilibrium 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 multiscale 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 multiscale modeling problem has been studied by colleagues [AlDerzi 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 (THNDDO). 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 stressstrain 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 CorrectorPredictor 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) (41) 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 higherorder terms are predicted according to the fifthorder 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 NoseHoover 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 (42) 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) (43) 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 stressstrains 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 41 shows the stressstrain 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 46). Also shown for comparison are the reference results obtained from the TH method for the entire rod. It is observed that the stressstrain curves given by the TTAM and the BKS potentials differ significantly from the results of the quantum theory. Examination of the curvature of the stressstrain 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 stressstrain 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 stressstrain 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 crosssectional area of the nuclei in the end caps. The crosssectional 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 41. Comparison of stressstrain 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 multiscaling 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 nonlocal 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 twoand threebody potentials, but the functional form of these fewbody potentials is still optional. Researchers have constructed potentials having as many as 42 parameters [Watanabe et al., 1999]. Other simpler silica potentials having threebody 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 quantummechanical 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 manybody in nature. Instead, the potentials are viewed as constructs of the multiscale 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. AlDerzi et al., [2004] showed based on highlevel 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 multiscale method, which is to construct the composite nanorod (i.e., the QM and CM regions joined together) with compatible smallstrain elastic properties throughout. We emphasize that an attempt to fit the entire stressstrain curve is not reasonable as the physics of largestrain 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(44) 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 shortrange 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 (00, SiO and SiSi). In earlier work, the parameters for the TTAM and BKS potentials have been chosen by fitting ab initio HartreeFock potential energy surfaces of SiO4(4)+4e+ and Si04H4 respectively. Furthermore, the SiSi interactions in the BKS potential are arbitrarily set equal to zero. The resulting parameter values are given in Table 42. However it is noted that this fitting was done with a low level HartreeFock 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) (45) 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 twodimensional domain in Figure 42, where a search near the second minimum would miss the global minimum. f(x,y) Figure 42. 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 gradientbased 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 BroydenFletcherGoldfarbShanno [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 (46 f(p) ;f(po)+g, x(p po)+(ppo)T x H x (p po) (46) 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 41 lists these different parameters. Table 41. 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 41 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 IIV. Also the a, s are more than an order of magnitude different for set VVI. 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 712%. 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 42. 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 42 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 43 compares the potentials for pairwise SiO interactions. 40" 3C C C _ 1 TTAM STH Claial oe TH rldsisiil potential  I I I I Figure 43. Comparison of SiO pairwise interactions for the three potentials. (The inset shows the zero crossings of the potential) Table 43 shows the comparison of several bondlengths for the three potentials and those for the reference TH rod. Table 43. 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 SiO 1.641 1.642 0.04 1.65 0.76 1.616 1.5 intermediate O's SiO 1.71 1.68 1.75 1.67 2.5 1.62 5.1 Approx. diameter 6.55 6.57 0.35 6.57 0.35 6.41 2.0 R I.III,'2 11lllls) StressStrain Curves This new potential was used in MD simulations to obtain the stressstrain 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 endcaps at 25 m/s). Figure 44 shows the results for this new potential, together with those for the TTAM and BKS potentials, in comparison with the TH results. The stressstrain 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 45 plots the stressstrain 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 44. 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 45. Reduction of the wiggles with temperature. 0.00 0.00 0.02 0.04 0.06 Strain Figure 46. Stressstrain 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 stressstrain 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 100cm1 (up tol600 cm1) for the TH (green), TTAM (red), BKS (pink) and new potential (blue). At the lowest wavenumber 100cm1 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 lowerlying 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 stressstrain 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 47. Histogram of number of modes vs. wave number (cm1). Improving StressStrain Curves for Higher Strains Figure 46 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 largestrained domains is inherently quantum mechanical. However, it is desirable for multiscale 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 chargetransfer term in the potential, and the second is to go back and reparameterize 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 chargetransfer 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 endcaps 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 44) were found to give a rod with similar geometry (Table 45) as that of the TH rod and a stressstrain curve (Figure 48) that matches with the TH curve well up to 5% but becomes stiffer afterwards. Table 44. New parameters for classical potential2 (NTH2 Classical Potential) NewTH OldTH 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 NTH2 classical potential (NewTHClassical potential # 2). Table 45. Structure obtained from NTH2 Classical Potential in TH nanorod 25 20 1C NTH2 Bond lengths (A) and H Class % TH Classical bond angle (degrees) Potential error Potential In Silica Planes SiO 1.641 1.642 0.04 SiO 1.71 1.67 2.5 Diameter 6.55 6.54 0.15 0 0  / m M fall OM(TH) 0 X I ITH2 0 ," comparison to that of the 0.00 0.05 0.10 0.15 0.20 Strain 0.25 0.30 0.35 Figure 48. Behavior of NTH2 Classical Potential near fracture. 15 As can be observed from Table 45, the accuracy of the bondlengths and bond angles (compare Table 43) is not lost but the rod becomes much stiffer than the TH after 5% strain. This stiffness will guarantee that in the multiscale 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 multiscale 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 56% 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 NTH2 Classical Potential found using the strain data up to 10% is a good candidate for multiscale 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 (NTH2 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 multiscale 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 multiscale modeling method. Such a composite rod allows us to study multiscale phenomena for the events involving largestrains localized to the QM domain. Multiscale 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 NTH2 Classical Potential. While calculating the TH forces for the atoms in the QM region, the effects of the CM environment are incorporated with pseudoatoms 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 51 shows the stressstrain behavior of the nanorods obtained from three different methods: (i) the full quantum mechanics TH method, (ii) full NTH2 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 NTH2 nanorod in terms of small strain elastic properties and structure. The stressstrain results shows the success of our multiscale 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 pseudoatoms, 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 pseudoatoms may improve the stressstrain 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 bondlengths, so the failure of pseudoatoms at high strains does not pose any constraint so far as application of the proposed multiscale 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 (onedimensional system) and defect free nature. 250 150 100 ll .A1 10o I x NTH2 ,4' A composite 50 , "N 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 strain Figure 51. Comparison of stressstrain of the composite rod with those obtained from NTH2 Classical Potential and TH. The satisfactory construction of such a composite rod completes our recipe for the multiscale modeling for all practical purposes. Conclusions In this chapter, we have satisfied the last criterion of our proposed multiscale 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 multiscale studies with substantially reduced computational times. CHAPTER 6 INDEPENDENCE OF MULTISCALE MODELING WITH RESPECT TO THE CHOICE OF UNDERLYING QM Introduction As mentioned in chapter 3, a proper multiscale 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 BornOppenheimer approximation, HohenbergKohn theorem, KohnSham 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 pseudoatoms 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 multiscale modeling scheme is faithful to the chosen form for the underlying quantum mechanics in both cases. Review of DFT Method BornOppenheimer 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 manybody wave functions. This reduces the manybody 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)] (61) 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. HohenbergKohn Theorems The DFT approach is based on the first HohenbergKohn (HK) theorem [1964], which states that all the ground state properties of a system are functionals of the electron density. The second HohenbergKohn 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] (62) where F[p]= T[p]+V,,[p] (63) T[p] is the kinetic energy V,,[p] is the electronelectron interaction,. The ground state electron density is the density that minimizes E[p] subject to the condition Sp(r)dr= N (64) 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 KohnSham equations from the HK theorems. KohnSham Equations The Functional F[p] can be decomposed into three different components as F[p]= T, [p]+ J[p] + Ex [p] (65) T, [p] is the kinetic energy of a noninteracting inhomogeneous electron gas in its ground state, J[p] is the classical coulomb energy (or Hartree energy) and the functional Ex [p] is the exchangecorrelation energy. Ex [p] in fact includes exchange, selfinteraction, Coulomb correlation terms as well as the KE increment T[p] T [p]. Using equations 6 2 and 63 the total energy functional can be written as Eee[P] = v(r)p(r)dr + T [p] + J[p] + Ex [p] (66) Minimizing this expression subject to equation 64 gives ee u = 0 p being a Lagrange multiplier that enforces 64 5p(r) S= v T + (67) e p(r) where veff is called the effective KS potential defined as dJ[p] +E,, Vef = v(r) + + dp(r) dp(r) (68) (68) = v(r)+ p(r') dr'+vx (r) SIrr' From the definitions of the noninteracting 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 (69) 2 These are known as KohnSham 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 (610) 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 exchangecorrelation contribution to equation (610) 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 upspin and downspin electrons: (r)= p(r) p(r) (611) 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 KohnSham equations v dr, + [r, :,+(r )= (r) = ac,/ (612) S2 \iRr)2dr So we have two sets of wavefunctions, one for each spin. TroullierMartin PseudoPotential Implementation of allelectron 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 KohnSham orbitals in the core region. The pseudopotential can be obtained from allelectron atomic calculations by solving the radial KS equations self consistently. Sd + ) +V[p; r] rR,,(r)= srR,2(r) (613) 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 pseudowave 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 pseudowave function should be identical to the normalized radial real wavefuntion outside a cutoff 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)] (614) 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 pseudowave functions from the screened potential V,~ .= V, () PP (r) ( () (615) 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 (616) where VPP'c () is the local potential and VIPPnlc (r) is the nonlocal 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 616 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 multiscale program package, known as BornOppenheimer 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 dualspace formalism. A plane wave basis set and cutoff energy of 30.84 Rydbergs is chosen. Results and Discussion We next use the DFT to study the multiscale modeling of the nanorod. The accuracy our multiscale 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 pseudoatoms that model the nearest neighbor interactions. Training of DFTPseudoatoms A pseudoatom is constructed based on the TroullierMartin (TM) pseudopotential discussed above. A cutoff 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 31.). But unlike the THNDDO method in which we could change the electronion as well as electronelectron interaction parameters, in DFT we can modify only the electronion interactions. The three options that can be varied to alter the electronion interactions for the F atom using the TM pseudopotential are: (i) the core charge on F, (ii) omission of the nonlocal part in the potential, and (iii) switching the local and nonlocal 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 61. Magnitude of force on the terminated Si atom in pyrosilicic acid with changing parameters of the dftpseudoatom Only Snin Without non Switching zcore changing z e local part between s & p zcore 7.0 0.88 2.9E04 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 pseudoatom 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 pseudoatom 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 nonlocal 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 61 shows that the minimum force is obtained when the core charge is 7.0 and the nonlocal part is omitted. The force for this particular set of combination is lower by two orders of magnitude than for the other cases. This pseudoatom 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 34). 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 