UFDC Home  myUFDC Home  Help 



Full Text  
ADAPTED AB INITIO THEORY: A SIMPLIFIED KOHNSHAM DENSITY FUNCTIONAL THEORY By JOSHUA J. McCLELLAN 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 2008 02008 Joshua J. McClellan To my parents. ACKNOWLEDGMENTS I thank the chair and members of my supervisory committee for their mentoring. I would also like to give special thanks to Mark Ponton, Victor Lotrich, and Dr. Deumens for their extensive help in the last few months. I owe a debt of gratitude to my family, whose support has instilled a deep yearning to pursue my interests. Several people made many other contributions to my life, both personally and academically. Larry urged me to appreciate the finer things life has to offer and kept me sane. I thank the two real chemists, Travis and Lani, for their unending quest to answer the basic question man asks of life: What is fire? The people of the Quantum Theory Project offered a constant source of inspiration, and often a source of morning headaches. I offer many thanks to Andrew, Tom, Dan, ('C!ii i ,, Georgios, Seonah, Martin, Joey, Lena, Julio, and of course Merve. I also acknowledge Prof. Merz for the use of the DIVCON program and Dr. Martin Peters for helping to implement and test some initial ideas for cubic spline interpolation (C'!i pter 4). Special thanks go to Andrew Taube and Tom Hughes for the numerous extensive discussions that were invaluable in my research. Also, I acknowledge Tom Hughes' work in implementing the transfer Hamiltonian parameters for nitromethane clusters (C'!i plter 2). I would like to recognize the help from Dr. Ajith Perera for help throughout my graduate studies. TABLE OF CONTENTS page ACKNOW LEDGMENTS ................................. 4 LIST OF TABLES ....................... ............. 7 LIST OF FIGURES .................................... 8 A BSTR A CT . . . . . . . . .. . 9 CHAPTER 1 INTRODUCTION AND BACKGROUND .................... 10 1.1 Who Ordered a New Semiempirical Theory? ........ ......... 10 1.2 Quantum C('! iiii ry ................... ....... 15 1.3 HartreeFock ............................... 17 1.4 CoupledCluster Theory ............................ 19 1.5 KohnSham Density Functional Theory .................. 22 1.6 Hiickel Method .. ... .............. .......... 23 1.7 Neglect of Diatomic Differential Overlap .................. 25 1.8 TightBinding ............... . . .. ..26 1.9 DensityFunctional TightBinding . . .... ... 27 1.10 SelfConsistentC( rge DensityFunctional TightBinding . ... 29 2 THE TRANSFER HAMILTONIAN .................. ...... .. 35 2.1 Theoretical Method .................. ............ .. 35 2.2 Computational Details .................. .......... .. 37 2.3 NitrogenContaining Energetic Materials .................. .. 37 2.3.1 Nitromethane Monomer. .................. .... .. 40 2.3.2 Nitromethane Dimer .................. ..... .. 41 2.3.3 Nitromethane Trimer ............... ..... .. 42 2.4 Silica and Silicon ............... .......... .. 44 2.4.1 Silicon Clusters ............... ......... .. 45 2.4.2 Silica Polymorphs ............... ........ .. 46 3 IMPETUS FOR A NEW SEMIEMPIRICAL THEORY . . 58 3.1 Neglect of Diatomic Differential Overlap Based Methods . ... 58 3.1.1 Does NDDO Approximate HF? ............ .. .. 58 3.1.2 Artificial Repulsive Bump ..... . . .... 60 3.1.3 Separation of Electronic and Nuclear Energy Contributions ..... .61 3.1.4 Total Energy Expression ............... .. .. 61 3.2 SelfConsistent C' irige DensityFunctional TightBinding . ... 62 3.2.1 SelfConsistent C!i ,i'ge DensityFunctional TightBinding in an AO Framework ...... .............. ........... 62 3.2.2 Mulliken Approximation Connection to SCCDFTB . ... 65 3.2.3 SCCDFTB Repulsive Energy . . . .. ... 66 3.3 Systematic Comparison of HF, NDDO, DFTB, SCCDFTB, KSDFT in a Single Framework .................. ............. .. 67 3.3.1 General Energy Expression . . ........ .. .. 70 3.3.2 Density Independent FockLike Matrix Contribution . . 71 3.3.3 Density Dependent FockLike Matrix Contribution . ... 72 3.3.4 Residual Electronic Energy ................ ..... 73 3.3.5 General Gradient Expression . . ....... ..... 74 3.3.6 Requirements for an Accurate Simplified Method . . ... 75 4 ADAPTED AB INITIO THEORY .................. ..... .. 80 4.1 TwoElectron Integrals ............ . . ... 80 4.2 KohnSham Exchange and Correlation ........... ... .. 89 4.3 Adapted ab initio Theory Model Zero ............. ... .. 92 4.4 Implementation ............... ........... .. 94 4.5 Verification ............... ............ .. 97 4.6 Conclusion ................ ............. .. 101 APPENDIX PARALLEL IMPLEMENTATION IN THE ACES III ENVIRONMENT 117 REFERENCES ..................... ...... .......... .122 BIOGRAPHICAL SKETCH ................... . ... 127 LIST OF TABLES Table page 11 The NDDO onecenter twoelectron integrals ................ 32 12 The NDDO twocenter twoelectron integrals .................. ..33 13 Multipole distributions .................. ............. .. 34 21 Equilibrium geometry for NMT, in A and degrees . . . 49 41 Adapted integrals ............... ............ .. .. 103 42 Average percentage of Fock matrix contribution by adapted approximation 104 LIST OF FIGURES Figure page 21 The NMT force curve for CN rupture .................. .. 50 22 The NMT energy curve for CN rupture .................. ..... 51 23 The NMT dimer energies relative minimum in the intermolecular hydrogen bonding distance . . . . . . . . .. .. . 52 24 Snapshots of the geometry optimizations for NMT dimer and trimer performed with the THCCSD and AM.................. . . ..... 25 Average deviation of SiO stretch...... . . ..... 26 Average deviation of SiSi stretch.... . . ..... 27 Average deviation of OSiO angle.............. . . ..... 28 Density deviation from PW91 DFT... . . ..... 31 Multicenter contribution to the C28N28 matrix element of NMT . . 32 The AM1 artificial repulsive bump for NMT direct bond fission. . . 33 The SCCDFTB total energy breakdown........... . . ..... 34 Energies from Riidenberg and Mulliken approximations .. . ..... 41 Scatter plots of agreement between approximate and exact twoelectron integrals 42 Dissociation curve of CN with 4center Riidenberg approximation . . 43 Orbital products ...................... . . ..... 44 Dissociation curve of CN with GZERO(B3LYP) approximation . . . 45 Dissociation curve of CN with AATMO approximation .. . ..... 46 Timing of ab initio versus Riidenberg fourcenter integrals . . .. 47 Percentage of multicenter terms versus system size .. ............. 48 Timing of Fock build using traditional NDDO and NDDO with cubic splines as a function of system size....... . . ...... 49 Error introduced per atom from cubic splines as a function of system size . 410 Pseudoreactionpath splitting C20 .. .................... 411 Force RMSD from MP2........... . . ..... 412 Force RMSD from B3LYP........ . . ..... 53 54 55 56 57 76 77 78 79 S105 106 107 108 109 110 112 113 114 115 116 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 ADAPTED AB INITIO THEORY: A SIMPLIFIED KOHNSHAM DENSITY FUNCTIONAL THEORY By Joshua J. McClellan May 2008 C'!I ir: Rodney J. Bartlett Major: Chemistry We present work toward a oneelectron Hamiltonian whose solution provides electronic energies, forces, and properties for more than 1000 atoms fast enough to drive large scale molecular dynamics. Ideally, such a method would be as predictive as accurate ab initio quantum chemistry for such systems. To design the Hamiltonian requires that we investigate rigorous oneparticle theories including density functional theory and the analogous wavefunction theory construction. These two complementary approaches help identify the essential features required by an exact oneparticle theory of electronic structure. The intent is to incorporate these into a simple approximation that can provide the accuracy required but at a speed orders of magnitude faster than tod 'v's DFT. We call the framework developed adapted ab initio theory. It retains many of the computationally attractive features of the widelyused neglect of diatomic differential overlap and selfconsistentcharge densityfunctional tightbinding semiempirical methods, but is instead, a .,,1l.:7 1 method as it allows for precise connections to highlevel ab initio methods. Working within this novel formal structure we explore computational aspects that exploit modern computer architectures while maintaining a desired level of accuracy. CHAPTER 1 INTRODUCTION AND BACKGROUND The underlying physical laws for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. P.A.M. Dirac 1.1 Who Ordered a New Semiempirical Theory? Available quantum chemical methods are incapable of achieving socalled chemical accuracy (1 to 2 kcal/mol for bond energies) for a system composed of thousands of atoms within a computationally efficient framework (approximately one d4i on a typical desktop machine). A computational method with this capability would revolutionize the way that matter is studied in biology, chemistry and physics. In the field of quantum chemistry, there are many unexploited connections among wavefunction, density functional and semiempirical theories. Use of such connections will aid in the development of quantum chemical methods with the desired accuracy and computational efficiency. The ultimate goal then is to provide scientists with a computational tool that makes materials simulations predictive, chemically specific, and applicable to a wide variety of mechanical and optical properties of realistic, complex systems. Though highlevel quantum methods approach the desired accuracy, the computational scaling of such methods precludes their use for systems larger than a tens of atoms (e.g., coupledcluster with singles and doubles (CCSD) scales O(N6), where N is the number of basis functions). In contrast, the computational scaling of semiempirical (SE) methods is acceptable (typically scales 0(N3)), though they can be made to scale linearly with N. For instance, if the size of a system is increased by an order of magnitude then a CCSD calculation would incur a factor of 106 more computational time. For the same increase in system size, the computational time of a normal SE method would only increase by a factor of a thousand, but with linear scaling only by ten. In the regime of very large systems (thousands of atoms) it is apparent that a calculation that could be performed in one di, using SE methods might take several years for CCSD. Unfortunately, SE methods or, better, .:,,1 fi.:; methods (SM) are not yet sufficiently reliable, although they do provide an excellent opportunity for improvement. There are two general approaches to improve SE methods: first, existing reference data can be expanded to include additional systems and properties or improved fitting procedures can be used to help locate new minima in the model Hamiltonian's parameter space; second, more rigorous model Hamiltonians can be derived to satisfy more demanding desiderata, such as uniform accuracy of the energy and forces across the entire potential energy surfaces (PES) and improved electronic densities, ionization potentials, and associated firstorder properties. The first route is a mainstay of SE method development, though we will show that the second SM strategy has the potential to be more valuable in the long run. To this end, a new SM is discussed that incorporates rigorous formal connections to highlevel ab initio methods while retaining computational features of commonly used SE methods. An important computational aspect of existing SE methods is that their speed is dictated by the solution of a set of oneparticle eigenvalue equations (typically matrix diagonalization) and not by the construction of those oneparticle equations or any ensuing step. This computational restriction, in addition to a modicum of accuracy, guides the development of the method under discussion, adapted ab initio theory (AAT). AAT not only provides a general theoretical SM framework into which existing SE methods can fit, but also has the advantage that there are welldefined connections between it and highlevel ab initio theories. This novel framework not only allows for improvements to existing SE methods that enhance transferability and extendibility, but also has the flexibility needed for more systematic improvement in future versions of SE methodology. Before we begin the detailed description of our quantum mechanical method the need for a fast and accurate method to describe the properties of large molecules will be demonstrated. Phenomena that involve bondbreaking or bondforming often necessitate a detailed knowledge of the electronic structure and the PES. Such quantum chemical descriptions are important in many complex chemical processes; for example: 1. Gas phase studies. Plasmon absorption of organometallic nanoparticles. Degradation pathiv of chemical warfare agents. Environmentally safe green explosives. 2. Condensed matter studies. Defects in semiconductors. Cracktip propagation and hydrolytic weakening of silica. Bulk tensile properties of silica polymorphs. 3. Biological studies. Bioinformatics: fragmentation patterns of peptides for proteomics. Biocatalysis: enzymatic reactions on oligosaccharides for biofuel generation. Photoreceptors such as rhodopsin. Drug design. Quantum chemical methods are used to some extent in the analysis of these types of problems. The ideal quantum chemical method could be applied to any size system and achieve any desired accuracy, while using minimal computational resources, but, of course, this cannot be achieved. Hence, there are often practical limitations dictated by the desired accuracy and the available computational resources. Our goal is to create a method with chemical accuracy that can be applied to systems with thousands of atoms. A method with these features would be able to readily provide the quantum mechanical contribution to the solution of the systems mentioned. The two main approaches in quantum chemistry are wavefunction theory (WFT) and density functional theory (DFT). Likewise, there are two commonly used SE techniques: neglect of diatomic differential overlap (NDDO) methods, which originated in WFT, and densityfunctional tightbinding (DFTB), including its selfconsistentcharge extension (SCCDFTB), which originated in DFT. Both NDDO and SCCDFTB can be considered semiempirical since they both approximate the formal structure of the original methods and are parameterized to reproduce either experimental or theoretical reference data. Just as the formal connections between WFT and DFT have recently been uncovered and exploited, via ab initio DFT [1], there are similar formal connections between their semiempirical counterparts: NDDO and SCCDFTB. These connections form the basis of AAT, which is a SM rather than a SE method. Existing SE methods are widely used and, when applied judiciously, perform quite well. Regrettably, they are often used in situations for which the results are not reliable, such as locating transition states. To expand the applicability of SE methods, we combine their primary features with those from rigorous WFT or DFT, and define the new AAT framework. Hence, AAT encompasses both NDDO and SCCDFTB. Despite their formal deficiencies, both NDDO and SCCDFTB have been quite successful in a variety of applications [2]. Nonetheless, there remain several aspects of these SE methods that fall short of expectation that could be greatly improved. In addition to obvious formal limitations, existing SE methods were designed to reproduce the energetic of systems at equilibrium and often yield unphysical results near transition states and along dissociation pathi,. Uniform accuracy of the PES is critical for molecular dynamics (\!1)) simulations, for determining the lowest energy (global minimum) structure, and for calculating rate constants. Associated with these failings are unrealistic electronic densities that, in turn, make combining these methods into multiscale modeling schemes difficult. Furthermore, errors in firstorder properties such as the dipole and higherorder moments would be remedied by a method that delivers adequate electronic densities. Attempts to improve such properties have met with modest success, typically at the expense of other aspects of the model. By far, the most common route to incorporate such additional features is to reparameterize the model. Unfortunately, building in accuracy for nonequilibrium structures, excited states, and other properties, by simple reparameterization of existing NDDObased SE methods [3] is not readily achievable, at least not in a way that maintains transferability (similar accuracy regardless of bonding characteristics or system size). This difficulty highlights an important distinction between the current AAT route of designing a new simplified model Hamiltonian based on rigorous ab initio principles versus attempting to extend existing model Hamiltonians to include features that the model was not originally intended to capture. To avoid such complications the AAT SM is primarily designed to yield accurate electronic structure (electronic density) and, secondarily, the total energy and forces. Excited and ionized states are of interest. Clearly, the electronic structure and the energy of a molecule are closely connected. Consequently, the distinction between the two is subtle, nevertheless, their separation is crucial in order to ensure that the AAT achieves uniform accuracy for energies, forces and firstorder properties over an entire PES. The general computational approach is as follows: * Approximate fourcenter terms with a sum of twocenter terms using the Riidenberg approximation [4, 5]. * Determine exchangecorrelation potentials to guarantee the exact result at the separatedatom limit, starting from a sum of atombased densities. * Use cubic splines to interpolate these twocenter potential terms (a lookup table of splines need only be generated once and is very fast). * Solve a set of oneparticle equations with explicit inclusion of orbital overlap. * Include explicit exact (HFlike/ nonlocal) exchange. * At convergence of selfconsistency, make a final energy correction by inserting the density into a DFT functional (B3LYP) and performing a single numeric integration for the exchangecorrelation energy. An important computational consideration is the inclusion of the overlap between orbitals, which is a feature of AAT that distinguishes it from standard NDDO methods. The zerodifferential overlap approximation, which is the fundamental approximation of NDDO, implies that there is no overlap between orbitals. Though this approximation enhances the speed of the method, it neglects an important physical feature of electrons and the basis sets used to describe them. There have been attempts to reintroduce such overlap effects in NDDO methods, most notably Thiel's orthogonalization correction methods (OM1 and OM2) [6]. Still, Elstner's development of SCCDFTB [7] demonstrates that one can include overlap explicitly and still retain a computationally efficient framework. Thus, AAT explicitly includes the overlap between atomic orbitals. The predominant computational aspects of the method have been highlighted; now, consider the formal development of the AAT approach. There is an exact reference to the exchangecorrelation potential with which we can directly compare. The electronic density can be corrected iteratively by an inverse procedure, such as that of Zhao, Morrison, and Parr [8] or that of van Leeuwen and Baerends [9]. These procedures take a highquality reference density and generate the numerical exchangecorrelation potential that best reproduces it. In the limit that the reference density is exact, such procedures produce an exact local exchangecorrelation potential. This exact limit defines a set of oneparticle equations, which in turn, define the firstorder electronic energy. When added to the nuclearnuclear repulsive energy, the electronic energy gives the total energy correct through firstorder in the trial wavefunction. In existing SE methods, the remaining higherorder energy terms are combined with the nuclearnuclear repulsion to give a parameterized atompair based total energy correction that would be predominantly nuclear repulsion. Instead, we choose to retain a separate nuclearnuclear repulsion contribution to ensure that every aspect of our model has a direct onetoone correspondence to ab initio theory. The connections to the exact limits for both density and total energy components of the overall parameterization provide the necessary structure so that every approximation of the model can be rigorously verified. 1.2 Quantum Chemistry Electronic structure theory describes the distribution of electrons around nuclei in atoms, molecules, and solids. Many observable features of a molecule are given by its electronic structure, though computational efficiency requires approximations that can limit the accuracy with which we can describe some chemical phenomena. There is a series of standard approximations that permeate quantum chemistry that restrict its universal applicability, but that are valid for many cases of interest: * Timeindependent Schridinger equation. * BornOppenheimer approximation. * Nonrelativistic Hamiltonian. * Linear combination of atomic orbitals. The timeindependent Schridinger equation provides the solutions needed for the timedependent Schr6dinger equation, and though some problems require a timedependent solution, the vast ii, ii iiy of problems in quantum chemistry can be treated adequately with the timeindependent Schridinger equation. The BornOppenheimer approximation is based on the simple observation that electrons move much faster than nuclei. This difference in speed effectively implies nuclear motion cannot affect the electronic structure significantly and we can solve the Schr6dinger equation for a set of fixed nuclear coordinates. For example in the hydrogen atom the electron is traveling at about two million meters per second, which is orders of magnitude faster than the motion of the nucleus. Also, this speed is less than 1 of the speed of light, which also implies that we need not concern ourselves with relativstic effects that would only become important for atoms with heavier nuclei (typically those with atomic number larger than 25). For such heavy atoms the speed of an electron near the nucleus will approach the speed of light and would require a relativistic Hamiltonian for an accurate description. In such cases effective core potentials that approximate the passive (Darwin and massvelocity) relativistic effects can be used, which still allows one to keep the nonrelativistic Hamiltonian framework. Finally, the linear combination of atomic orbitals (LCAO) approximation introduces basis sets into the problem and gives a general framework for systematically expanding the reference space of the orbitals that the electrons may occupy. The timeindependent Schr6dinger equation within the BornOppenheimer (clamped nucleus) approximation is Hf = ET (11) where E is the energy of the system, T is the wavefunction of the electrons, and the nonrelativistic Hamiltonian in Hartree atomic units is H V Y2 ZA ZAZB rAi rij RAB i A,i i>j A>B (1 2) Sh(i) + + VNN i i>j with i and j referring to electrons and A and B to nuclei. The terms in Equation 12 are the operators that correspond to the kinetic energy of the electrons, electronnuclear attraction, electronelectron repulsion, and nuclearnuclear repulsion. The wavefunction T describes the state of our quantum system. We choose to approximate T within a basis set of atomic orbitals (AOs). Enforcing antisymmetry and satisfying the Pauli exclusion principle can be achieved with an antisymmetric product of orthonormalized AOs within a Slater determinant, the LCAO approximation )=\ IX1(1r)X2(2) ... Xn(rn)) (13) where n is the number of electrons and the spin orbital x(ri) = ((ri)a, and a, denote spin and O(r) is a spatial orbital. 1.3 HartreeFock Focusing on the ground electronic state, the most widely known approximate solution of the electronic Schr6dinger equation is made by a meanfield approximation: HF theory. That is, we insist that a zerothorder approximation to the nelectron wavefunction (0o) be written as an antisymmetrized product of n oneelectron molecular spinorbitals, )o = A[ i(1)Q2(2)... ~(n)]. Here A is the antisymmetrization operator. Then using the RaleighRitz variational principle, EHF > EeaOct, we derive the oneparticle HF equations which yield the energetically best possible single determinant wavefunction. Each of the orbitals in its canonical form is an eigenfunction of the effective oneparticle operator, f. It is fitting to start with the most wellknown approximation in electronic structure theory. HF is an excellent touchstone since correlation energy is defined by it (Ecorreltion Exact EHF) and because many quantum chemists have developed an intuitive understanding of it. Also, HF provides a convenient way to introduce notation that will be used throughout the current work. The oneelectron operator of the full Hamiltonian is h(1)X) ( (_ ZZA)X) (14) 2 A rA and the twoelectron operators are Coulomb J and exchange K, J(1)jXi) XJ2 )2dr2 Xi) (15) iX \ X'(rX 2 Xi (ri2) K(1)j i = dr2 \Xj) (16) j7i tr12 and f(1)Xi) [h(ri) + (J(1) K(1))] X) QXi) (17) denoting the molecular spin orbital Xi as a linear combination of atomic spin orbitals Xi Xi(ri) Cix(ri) (1 8) Using Lagrange's method of undetermined multipliers and enforcing orthonormality of the orbitals gives a set of oneparticle matrix equations, the HartreeFockRoothaan equations: FC SCc (19) where Sp,= J2x(rl) ,(rl)dr, (110) and F, / (ri)f(ri)x (ri)dri (1 11) The construction of a Focklike (or effective oneparticle) matrix is a central theme of this work. A Focklike matrix is an orbital representation of an effective oneparticle operator, such as found in HF or DFT. It is also possible to define others, as in the correlated orbital potential (COP) method [10, 11]. HF illustrates many of the necessary features of a oneparticle theory. One distinct feature, which sets HF apart from most other oneparticle theories, is that HF enforces the Pauli exclusion principle exactly via the exchange operator. Such exact exchange terms are missing in KohnSham density functional theory [12] which leads to the socalled self interaction error. Another benefit of HF is that there is Koopmans' theorem that gives meaning to the eigenvalues of the Pock matrix by relating them to ionization potentials and electron affinities. However, the exact DFT structure can provide a formally correct oneparticle theory that includes electron correlation, which HF cannot do because by definition it does not include electron correlation. COP is another exact alternative that can be defined from WFT. 1.4 CoupledCluster Theory The introduction of electron correlation in wavefunction theory can be achieved by generalizing the meanfield wavefunction to allow electrons to be excited into the orbitals that are unoccupied (virtual, i.e., do not appear in the HartreeFock determinant) orbitals, a,b,c, .... Those additional orbitals allow the electrons to avoid each other in a more detailed way (hence lower the energy) than in the meanfield wavefunction. That is, for single excitations (replacements), double excitations, triple excitations, etc., we define excitation operators Tj, j = 1,2, 3,... T o1 = t + (112) T2+=0 > tt (13) i i with the excited determinants, cb, meaning the antisymmetrized product, j = A[( (1) 2](2)..a(i)..bj). .. (n)], (115) where unoccupied orbitals a and b replace the occupied ones, i and j. The coefficient (amplitude) in front of the excitation introduces its proper weight into the wavefunction. The ultimate version of such a sequence would include all possible excitations up to nfold for n electrons and is called the full configuration interaction or full CI. Full CI is: * The best possible solution in a given basis set. * Invariant to all unitary transformations (rotations) of the orbitals. * Sizeextensive [13] (meaning that it scales properly with the number of electrons.) * An upper bound to the exact energy for a molecular electronic state. The full CI can be written conveniently in the coupledcluster (CC) form, = exp(i + 2+ T3+ ... + T)o. (116) [1316]. Approximations to CC theory are made by restricting the Tp operators to singles and doubles, as in CCSD, or to other approximations that decouple higherexcitations from lower ones, e.g., CCSD(T). Here the triple excitations are approximated from Ti and T2, which can be obtained from CCSD, e.g., without allowing the new T3 estimate then to contribute back to T1 and T2. The details of these methods, which can be found elsewhere [17], are not important to the present discussion. What is important is that all truncated CC methods are sizeextensive, but do not necessarily give upper bounds to the exact energy. The full CC equations provide the exact electronic energy, ECC = (oH o) EHF + (iab)( + ta) t +  f a)t (117) i where the twoelectron and oneelectron integrals are, respectively, (pqlrs) = ()(2)(1 P12)r 1 ( (1)d1d2 (118) (ilf a) J (1)f(1)(1)d1 (119) and P12 is the operator that permutes electrons 1 and 2. The amplitudes, {ti,, f, } come from projection of exp(T)Hexp(T) = H onto single, double, etc. excitations, (NUHo) = 0 (120) (<tH1o) 0 (121) The CC equations are a set of complicated, coupled nonlinear equations. Unlike the truncated CI equations, even when the CC equations are restricted to some subset of excitations, like singles and doubles (CCSD), they still guarantee a sizeextensive result. They converge much more rapidly to the full CI reference solution than does the sequence of partial sums of the CI expansion itself. For these reasons coupledcluster has become the method of choice in highaccuracy quantum chemical applications for small molecules [1719]. CC theory provides a route to exact reference data. It is wellknown that the hierarchy of CC methods are systematically improvable, and lead to the full CI solution for a particular basis set. Hence, for small molecules we have ready access to reference energies, forces, densities, and electronic properties to which we can compare directly. The wealth of information that is afforded by a highlevel ab initio calculation is superior to the typical amount of information available from experiment and provides unique insight into how well an approximation is working. For instance, besides getting energetic information about a molecule we simultaneously have access to the electronic density. This increased level of control, when compared to experimental reference data, allows for systematically verifiable approximations of our effective oneparticle Hamiltonian. The ability to verify every approximation systematically puts the current work into sharp relief from current semiempirical methods that depend intrinsically upon uncontrolled approximations. 1.5 KohnSham Density Functional Theory The KohnSham [12] equations may be written as fKS Xi) = iXi) (122) or, [1V2 + t(r) + vj(r) + Vc(r)]Xi) Ei Xi) (123) 2 where ve r) ZA (124) A RA 7I( Vj(ri) P(2 dr2, (125) J 712 p (r) I i(r) 2, (126) and vxc(r) is the exchangecorrelation potential. Constrainedsearch minimization of the noninteracting kinetic energy alone determines orbitals that yield the ground state density [20]. The other terms do not affect the minimization because they are explicit functionals of the density. The HohenbergKohnSham theorems guarantee that the density obtained from the exact ground state wavefunction will be returned by Equation 126 when using the exact exchangecorrelation potential. In a pragmatic sense, these arguments give a formal proof of existence, though they do not guide the construction of the exact oneparticle operator, or, more specifically, the exchangecorrelation potential. However, in combination with socalled inverse procedures, such as that of Zhao, Morrison and Parr (ZMP) [2123], which have been explored by Tozer et al. [24], or the procedure by Peirs et al. [25], the potential may be refined iteratively such that the solutions to Equation 123 will generate a reference density, presumably a highquality ab initio reference density. This is only part of the overall problem. In addition to the purely electronic properties, we also want total energies and forces. The total energy in KSDFT can be written, E K E i P 1 drd/ v.c(r)p(r)dr + E.[p] (127) where the exchangecorrelation potential is the functional derivative of the exchange correlation energy functional vc (r) = (128) Unfortunately, currently there is no procedure that allows us to take a numerically derived exchangecorrelation potential, one that is determined such that it reproduces a highquality reference density, and generate the exchangecorrelation functional, which is an explicit functional of the density. This inability to determine the exchangecorrelation energy functional from a reference density alone is the first indication of a fundamental problem in semiempirical methods, which will be discussed in further detail later: how does one design a model Hamiltonian that will simultaneously achieve the correct result for both electronic properties and energetic? However, if we are willing to assume the form of a particular exchangecorrelation energy functional then the potential is well defined, as is the effective oneparticle Hamiltonian. 1.6 Hiickel Method The Hiickel relectron method is perhaps the simplest quantum chemical approximation. It is intended to highlight the basic underlying effects of symmetry involved with w electrons for alternant hydrocarbons, though it can be made more general. In relation to HF, the same type of oneparticle matrix equations are solved, though with the overlap matrix set to the identity matrix and the Focklike matrix is composed of only a oneelectron (density independent) contribution, F = H HC = Cc (129) The total energy is E = nici (130) where ni is the occupation number for the ith MO. Note there is no explicit consideration of the nuclear repulsion in Equation 130. The Hiickel approximations are a ifi j Hi = 3 if i are nearest neighbors j (131) 0 otherwise. The extended Hiickel method [26], as its name implies, is merely an extension of the standard Hiickel method with the same total energy formula. The full matrix equation HC = SCc is solved, with overlap explicitly included. And the Hamiltonian is subject to a slightly more sophisticated approximation: a if i = j Hij = KSij(Hii + Hjj) if i and j are nearest neighbors (1 32) 0 otherwise where Hii is based on the ith AO of an isolated atom and K is a dimensionless adjustable parameter for an atompair. Hoffmann ,ii:. 1. 1 a value for K based on its effect on the energy of ethane (C2H6). He argued that, since the energy dependence on K is nearly linear beyond K = 1.5, the MO energies and, in turn, the total energy, will not be too sensitive to K for values greater than that. He then considered the energy difference between the 1 i:. red and eclipsed conformations of ethane and optimized the value of K to minimize the error between calculated and experimental values. The final value he arrived at was K = 1.75, a value is still in common use. The important aspect to note about the Hiickel method, as it pertains to the development of improved semiempirical methods, is that a very simple model can reproduce a few features qualitatively well. The spatial symmetry is correct due to the inclusion of the overlap matrix. However, despite the ability of such basic approximations to capture a large fraction of the total energy, more sophisticated descriptions are needed to guarantee chemical accuracy. 1.7 Neglect of Diatomic Differential Overlap Neglect of diatomic differential overlap (NDDO) is the model that most people associate with semiempirical theory in quantum chemistry todiw. The NDDO approach is meant to approximate the HF solution, not the exact solution, though its parameters can be varied to minimize error to reference data [2732]. It is a minimal basis set, valenceonly approach. NDDO is based on the zerodifferential overlap (ZDO) approximation: products of orbitals with respect the same electron coordinate that are on different centers are defined to be zero, (X B() 0. (133) The orbitals are considered to be L6wdin symmetrically orthogonalized AOs, so that the overlap matrix is the identity matrix, and the set of oneparticle matrix equations solved is FNDDOC = Cc. The Focklike matrix in NDDO contains oneelectron (H) and twoelectron (G) components. The oneelectron part is approximated as S / U, EB#A ZB(ppI BB) if p = v) _1, (/3 +/3,) if/P:/ V where p is on atom A, v is on atom B, and U,, and 3., are adjustable parameters, (ppIBB) models the repulsion between a product of orbitals (p) and a point charge at center B, ((p(1)1 11(1))). Note that the twocenter oneelectron terms include overlap and are therefore not consistent with the ZDO approximation. The onecenter twoelectron integrals are given in Table 11. The parameters Gs, Gpp, Hp,, Gpp, and Gp2 are adjustable. Within a local framework there are twentytwo unique twocenter NDDOtype twoelectron integrals for each atom pair (Table 12). The NDDO twocenter twoelectron integrals are approximated by a multipolemultipole expansion [33]. Each product of orbitals is assigned an approximate charge distribution (Table 13). Adjustable parameters arise from the distances between the point charges used in the multiple expansion and by enforcing the correct limit as R  0 where the equivalent onecenter term should be reproduced. As mentioned above, originally the NDDO Focklike matrix was designed to model the HF oneparticle Hamiltonian. The HF solution neglects all electron correlation effects and on the scale of chemical accuracy one should not a priori expect the NDDO approximation to reproduce experimental results. However, because some parameters are adjustable and fit to experimental data the typical argument is that electron correlation is ,'inl. :./l included. In fact, it is clear that electron correlation is no more implicit in the NDDO model than are relativistic effects, timedependence, nonBornOppenheimer effects, basis set incompleteness, and experimental errors in the reference data. Despite such formal gaps, the NDDO method has enjol, da tremendous amount of success and reproduces heats of formation and other properties with surprising accuracy. The difficulty arises when better accuracy is required in situations for which the model was not originally intended to describe, such as transition states or other nonequilibrium structures. In such cases one quickly approaches a point of diminishing returns for parameterizing to new reference data. The precise origin of these limitations is explored later. 1.8 TightBinding The basic equation for the total energy in tightbinding (TB) theory is assumed to be of the form, 2 ETB = Ci + U(R RB). (135) i= 1 AB where the ci's are the eigenvalues of hff i) = ( V2 + v,,t) i) = i), (1 36) e,,t is the nuclearelectron Coulomb attraction term and U(IRA RBI) is considered to represent all the remaining terms. In total, the latter constitute a repulsive function of the set {RAB} between atoms A and B. For the oneparticle Hamiltonian, the usual LCAO approximation is made and the Qi's are expanded in a set of X,'s giving the equation, HeffC = SC (137) where (Hff)p = (p heff I v) (138) and S,, = (p (139) In practice, these twocenter matrix elements are evaluated using the SlaterKoster parameterization [34]. The TB approximation has been applied successfully to many solid state problems and works especially well when the density of the chemical system is well approximated as a sum of spherically symmetric atomic densities. In the case of molecules, a more advanced approach that I ii the delicate balance of ( i [7] is needed. If we compare TB to HF, we see that U(IRA RBI) has to account for J K and VNN. In DFT it would account for J, VNN, and would combine exchange with correlation by adding vc,. 1.9 DensityFunctional TightBinding To vest the TB energy expression with more rigor, Foulkes and Haydock [35] derived it from firstprinciples DFT considerations. The energy expression in KSDFT, Equation 127, may be written as EDFT et+ iv7 + (r ) ] [p() + R R RE (140) i A,B Equation 140 is exact if p = f where p is the exact density. However, it should be noted that Qi are eigensolutions of h+ (t + vt + f + vc[po]) and not (v + Vext + f' P('). Following Foulkes and H 1i, .d 1 the density p is replaced by a reference density po and a small deviation 6p, such that p = po + 6p. Then we can rewrite the energy expression to be occ V2 P ' EDFT (i + +v',t + Po0 + v~[Po1i) 2Ir r  P(p V[po](po Vp) (141) 2 r r 2 RA B 1 '6p'(po p) + E o + p + A,B to enable us to introduce a form dependent upon the oneparticle Hamiltonian, hs, plus its corrections. The first new term accounts for the doublecounting of the Coulomb potential, the next accounts for the addition of vxc in the first term, and the third new term is the remaining part of the Coulomb potential. Next, the exchangecorrelation contribution Exc is expanded about the reference density, po, with a Volterra (Taylor) expansion for functionals, Exc[po + 6p] Ex[po] + vxc[po]p + f pf p 6P6P' + / / 6" 6p6p'6p" + "" (142) 6 6p6p'5p" PO Substitution of this expression for Exc into the energy equation yields an expression for the energy that is formally exact, yet only includes terms that are secondorder and higher in 6p. E + VEzX[ + Ir rI + v [p Ii) [po] t 6pp + t 62EX 7 1 Poo' + Epo + ZAZB f f r 6p6p' (143) oP " 2 1 r r' 2J 6p6p' P 1 r "3E + i E 6p6pl p" +  6 6p6p'6p" PO All of the terms in this expression that are independent of 6p comprise the standard TB equation. Hence, U(RA RB) = .f f Pfu[po] (Po) + E.,[po + t Z (144) U(IRA RBI) 2 I2 + [ 0+ IRA R:I 2 Ir r'd 2A,BRA RB which obviously is correct to second order in p. Note the inclusion of nuclearnuclear repulsion in this term. One of the great advantages of DFT approaches is that they offer a reasonable zerothorder approximation to the energy simply from knowing some approximation to the reference density, po. However, since E,.[p] is not known, getting the exchangecorrelation potential, vc,(r) = 6EW,/6p(r), and its kernel, fx,(r, r') = 62E,1/6p(r)6p(r'), from first principles, which are required to set up a rigorous DFT, is difficult [1, 36]. However, were we to accept one of the plethora of approximations to Ec (LDA, LYP, PBE, B3LYP, B88, etc.) the equations can be solved for a few hundred atoms, but not on a timescale that can be tied to MD for thousands of atoms. Furthermore, to provide a correct description of bond breaking the electronic charge has to be free to relocate. A common example is LiF, in which at R = Rq we have virtually an exact Li+F form, but at separation in a vacuum, we must have LiF0. Such selfconsistency usually is introduced by solving the AObased DFT equations via diagonalization, but this approach imposes an iterative step that scales as N3 to the procedure. Since the model is meant to be approximate anyway, the alternative charge selfconsistency approach, familiar from iterative extended Hiickel theory [37, 38], can be used more efficiently. This extension leads to the selfconsistentcharge densityfunctional tightbinding (SCCDFTB) approach of Elstner et al. [7]. 1.10 SelfConsistentCharge DensityFunctional TightBinding To introduce charge selfconsistency, the terms that are secondorder in 6p are approximated as a shortrange Coulombic interaction based upon Mulliken atomic charge definitions, 1 'p'p 12E N E2nd = ) + AqAAqBAB. (145) 2 r r' 6p6p' 24 PO A,B The AqA are the differences between some reference charge, qg, and the charge on atom A, occ qA cycisi (146) i pEA v 7AB is a pairpotential, whose common functional forms are taken from traditional semiempirical quantum chemistry schemes such as those due to Klopman [39], Ohno [40], and MatagaNishimoto [41]. For example, the latter form is 7AB = 1/[2(7A + 7B)1 + RAB] (147) where the single index 7A is meant to be a purely atomic twoelectron repulsion term. These are typically chosen as fixed atomic parameters in the calculation. For an assessment of how the SCCDFTB model performs compared to standard semiempirical quantum chemical methods, see Morokuma, et al. [42]. For example, unlike SCCDFTB, standard AM1 does not give the most and least stable isomers of C28 when compared to DFT calculations with the B3LYP Exc and a 631G(d) basis set. In contrast, the relative energies of the isomers of C28 with B3LYP/631G(d)//SCCDFTB are in very good agreement with B3LYP/631G(d), having a linear regression R2 coefficient of 0.9925. Though the geometries generated by SCCDFTB are excellent, the SCCDFTB relative energies are poor, with a linear regression R2 of 0.7571. SCCDFTB is a fairly rigorous semiempirical method. The framework provided by Foulkes and Haydock exposed how higherorder terms could be included. The practical implementation of a wellbehaved approximation of higherorder was subsequently provided by Elstner. Later, we will develop the SCCDFTB formalism from an MO viewpoint to facilitate the comparison between it and other methods. We have discussed the basic features of HF, CC, KSDFT, Hiickel, NDDO, TB, DFTB, and SCCDFTB, and have alluded to the role each pl in the development of an improved semiempirical method. We take two approaches to develop this target method. First, we demonstrate the transfer Hamiltonian, which provides an exact oneparticle framework in WFT. Second, we develop adapted ab initio theory which provides rigorously verifiable approximations to reproduce ab initio reference data. Table 11. The Parameter Gss Gsp Hp Gpp Gp2 NDDO onecenter twoelectron integrals Twoelectron integral (asss) (sp sp) (pp pp) (pp p'p') Table 12. The NDDO twocenter twoelectron integrals Term (ss ss) (ss pp,) (P7P Iss) (PaPa ss) (p,pl,p,) (PwPw Pua,) (ppa lpp,) (Pa\Pa aPa) (spa ss) Term (spa pP,) (spa PaPa) (ss spa) (pPw Spa) (PaPa SPa) (spw spw) (spa Spa) (sP lpWP,) (pPa sp,) (PPa pPa,) (P'Pw' [PPw') Table 13. Multipole distributions Term C'! I'ge distribution Point charges (ss Monopole 1 (sp Dipole 2 (pp Linear Quadrapole 4 (pp' Square Quadrapole 4 CHAPTER 2 THE TRANSFER HAMILTONIAN 2.1 Theoretical Method Historically SE method development has been guided by experiment, the goal being to reproduce a set of experimental quantities. While any QM method should in some limit reproduce experimentally observed phenomena, the exact WFT imposes more demanding and detailed requirements on a model than experiment alone. For instance, the wavefunction is not an observable yet contains all the information about the system. In this sense a model based in WFT aspires to reproduce the exact wavefunction and the effective Hamiltonian that generates it. The expectation value of this exact wavefunction for a particular Hermitian operator will yield a particular experimental observable. Highquality ab initio reference data is readily available from CC theory and provides an unambiguous point of reference for our model. Furthermore, by reproducing the features of an exact WFT we are guaranteed to reproduce all experimental phenomena. However, practical limitations restrict the functional form of the wavefunction and the Hamiltonian. The challenge is to incorporate the i ,inrv body effects of exact WFT into a simple (lowrank) model Hamiltonian. One aspect of the approach taken here is to develop a formally exact oneparticle theory, the transfer Hamiltonian (Th), which includes all manyparticle effects. Another aspect is to connect a simple approximate model Hamiltonian to the exact Th. Our initial work used the wellknown NDDOtype model Hamiltonian. Because these models have limited applicability, it has become evident that more complete models are needed that accurately incorporate the manyparticle effects of the Th. We develop the AAT simplified approach to include the features of the Th in a systematic way, which is discussed in detail later. Here we introduce the Th [43], a oneparticle operator with a correlation contribution, as an extension of the similarity transformed Hamiltonian from CC theory [44]. This approach provides a rigorous formal framework for an exact oneparticle WFT. Also, unlike other methods, such as the specific reaction coordinate approach of Truhlar [45], which attempt to fit a known energy surface to a SE form via reparameterization, the Th fits to the Hamiltonian itself, to provide a form that is simple enough that it can be solved on a time scale consistent with large scale MD. We start from the exponential ansatz in which the exact ground state wavefunction is given by I ,) = eT o) (21) where 4o is a reference function, which can be taken as a single determinant and T is the excitation operator of CC theory. When defined this way, we may write the Schr6dinger equation in terms of the similarity transformed Hamiltonian, H. HIo) = eTHeT ) = E o). (22) The solution of this equation has the same eigenvalues as the bare Hamiltonian. The total energy, E, is exact, thus the forces, i, are exact as well. In secondquantized form H may be written as, 1 H = ,Pq + 2 (p llrs)ptqtsr p,q p,q,r,s + (pqr lstu)ptqtrtuts + (23) p,q,r,s,t,u where we show the inclusion of three and higherbody interaction terms. To include the higherorder terms we renormalize with respect to the one and twobody terms and obtain a generalized, correlated Focklike operator, called the transfer Hamiltonian [10, 43] Th Y[h,, + Y Pxpq, rs(lth (24) PV A,7 where PA, is the single particle density matrix and both h., and (pqllrs) can be approximated with a suitable functional of atombased parameters. We have used the NDDO form in the current work, but the formalism developed thus far is not limited by this choice of Focklike operator. 2.2 Computational Details Parameters were optimized by a qu( iN. .tonRaphson optimization algorithm in tandem with a genetic algorithm to fit Austin Model 1 (AM1) NDDO parameters [27] to reproduce the forces determined with CC theory. The parameterization was performed by a combination of the PIKAIA [46] implementation of a genetic algorithm and the linearized qu, iN. .tonRaphson minimization routine of LBFGS [47]. The reference forces used in this parameterization were obtained with the ACES II [48] program system. We have used CCSD in a triple zeta basis, with a unrestricted reference and HartreeFock stability following (to ensure that we arrive at the proper unrestricted spin solution.) All semiempirical AM1 and OM1 calculations were performed using a modified version of MNDO97 Version 5.0 [3032, 49]. DFT calculations for the nitromethane (NMT) monomer were performed at the B3LYP/631G(d) level with the HONDOPLUS program [50]. The initial NMT dimer and trimer geometries were taken from a DFT study using B3LYP/631++G** [51]. Due to the large number of correlated reference calculations needed in modeling our clusters, we used the Turbomole Version 5 electronic structure package [52, 53] to do allelectron calculations at the resolutionoftheidentity second order M. .11, rPlesset (RIMP2) level in a TZVP basis. We chose an allelectron model in a large auxiliary basis to ensure an accurate treatment of correlation energies and geometries that were demonstrated on a series of compounds representative of atoms in various oxidation states [54]. 2.3 NitrogenContaining Energetic Materials The motivation of this study is to develop a Th that can accurately predict the quantum mechanical behavior of nitrogencontaining energetic materials. We present a new parameterization for a semiempirical NDDO Hamiltonian which improves upon the standard AM1 parameter set. This new parameterization reproduces the decomposition of NMT to the CH3 and NO2 radicals as predicted with CCSD/TZP. We demonstrate transferability to small NMT clusters. For NMT dimers and trimers this model Th predicts rearrangements similar to previously calculated unimolecular mechanisms. The NMT trimer undergoes a concerted rearrangement to the methylnitrite trimer when one of the NMT molecules becomes geometrically strained. NMT has been the subject of many theoretical and experimental studies due to the significant role that nitrogencontaining compounds pl i, in the chemistry of explosives, propellants, and atmospheric pollution [55]. The description of the decomposition of NMT and its products is crucial to the understanding of the kinetics and dynamics of larger energetic materials for which NMT is a model. NMT is small enough to allow detailed investigation of its complex decomposition [56, 57]. Specifically, even for a system of deceiving simplicity, its decomposition to radicals via simple CN bond rupture versus the competing NMTmethylnitrite (\ NT) rearrangement has been the subject of much debate [5659]. The use of classical potentials and standard semiempirical NDDO methods often yield incorrect behavior for energies, forces, and geometries, especially in nonequilibrium cases. Therefore, the use of such techniques to do MD will often lead to incorrect results when bond breaking and forming are studied. Alternatively, ab initio quantum chemistry is known to be predictive, in the sense that it can reproduce most quantities adequately compared to the experimental values in those regions without knowing the result prior to calculation. However, these methods are too expensive and, therefore, currently impractical for MD because of the time and disk space requirements. Previous NDDO parameterizations, such as AM1, can reproduce nearequilibrium structures for a large number of molecules, yet they often fail for nonequilibrium geometries. Thus, for the purposes of doing MD simulations, we are willing to sacrifice the ability to accurately treat a large number of molecules near equilibrium, for the ability to treat a small class of molecules accurately at all geometries. When modeling combustion or detonation it is particularly important to have a method that works equally well for equilibrium and nonequilibrium geometries if the chemistry of bond breaking and forming is to be described. However, this sacrifice means that the new Th parameter set must be applied judiciously. If used for a class of molecules for which it was not trained, then the accuracy of calculated properties could be diminished. If the appropriate set of reference data and a sufficiently complete model are chosen, then the Th will provide an equally appropriate and complete description of the system. Among the first to study NMT theoretically were Dewar et al. [56] who used MINDO/3 and found that the NMTMNT rearrangement occurred with an energy barrier of 47.0 kcal/mol and that MNT dissociated to H2CO+HNO with an energy of 32.4 kcal/mol. Hence they proposed that NMT decomposition occurs via rearrangement. McKee [57] reported the NMTMNT rearrangement barrier to be 73.5 kcal/mol calculated at the MP2 level of theory with a 631G(d) basis and that MNT dissociated to H2CO+HNO with an energy of 44.1 kcal/mol. The high barrier height of the rearrangement I.  that the decomposition pathway is simple bond rupture to NO2 and the CH3 radical. Hu et al. [59] performed calculations with the G2MP2 approximation and found that the NMT dissociation occurs via direct CN bond rupture with an energy of 61.9 kcal/mol. The NMT dissociation is lower than both the NMTMNT rearrangement and the NMTacinitromethane rearrangement by 2.7 and 2.1 kcal/mol, respectively. Nguyen et al. [58] reported that the NMTMNT rearrangement barrier is 69 kcal/mol and NMT direct dissociation to CH3 and NO2 radicals is 63 kcal/mol calculated at the CCSD(T)/ccpVTZ level using CCSD(T)/ccpVDZ geometries. Most recently Taube and Bartlett [60] report the NMTMNT rearrangement at 70 kcal/mol using ACCSD(T). Hu et al. [59] also performed a detailed study of unimolecular NMT decomposition and isomerization channels using DFT with G2MP2//B3LYP level of theory in a 6311++G(2d,2p) basis. Single points were basis set extrapolated along the minimum energy paths to correct for incompleteness using the G2MP2 method. Their article shows that from a particular intermediate state (IS2b) one possible reaction path is the rupture of the NO bond to give the methoxy radical and nitrous oxide. Hu et al. indicate that on this reaction path no transition state could be located. Results from photodissociation and tl!. i i .i.vi experiments .. 1 that one of the primary steps of NMT decomposition is the elimination of oxygen [58, 61, 62]. However, the reliability of both the B3LYP and, to a lesser extent, CCSD(T) at nonequilibrium geometries [59] is limited and the same accuracy should not be expected in regions other than equilibrium. The work of Li et al. [51] ,., i that the threebody interaction energy in bulk NMT is a critical component in accurately describing the potential energy surface. We investigated the use of our Th in predicting how an accurate treatment of threebody effects dictates chemical reactivity. 2.3.1 Nitromethane Monomer We have developed a set of NDDO parameters to reproduce the force curve along the intrinsic reaction coordinate for CN bond rupture and to reproduce the CCSD/TZP equilibrium geometry. As seen in Table 21 the agreement between CCSD and THCCSD is to within a degree for the angles and less than one hundredth of an A for bond lengths. THCCSD is shown in Figure 21 to reproduce the original force curve for CN bond rupture to within an accuracy of 0.005 Hartree/Bohr. The AM1 parameters were designed to reproduce equilibrium structures, and thus they have the correct force at the equilibrium bond length. However, at 0.1 A away from equilibrium the forces are in considerable error. A recurring feature for AM1 force curves, which is problematic for MD, is the unphysical repulsive behavior at large internuclear distances. The transfer Hamiltonian removes this artificial repulsion and gives the qualitatively and quantitatively correct dissociation. One way of determining energy differences in a single reference theory, other than HartreeFock, is by integration of the force curve. For AM1 the error in the forces is reflected in the energy as shown in Figure 22. Also shown in Figure 22 the DFT solution yields a potential energy surface and force curve which are incorrect at nonequilibrium geometries. 2.3.2 Nitromethane Dimer In analogy with how Bukowski et al. [63] have studied the interaction energy of dimers as a function of rigid monomer separation, we calculate the relative energies of the NMT dimer as a function of hydrogen bonding distance. The minimum chosen in our NMT dimer calculations from [51] is stable due to the formation of two hydrogen bonds of equal length and antiparallel arrangement of the molecular dipoles. In Figure 23, we plot the energies relative to the respective methods minima in ROH. Our data indicate that in the case of frozen monomer geometries, AM1 overestimates the repulsive wall at small intermolecular distances. It appears as if the H, C, N, and O parameters in the corecore repulsion function for AM1 were unable to cancel the effect from which MNDO has traditionally suffered, namely overestimation of energies in crowded molecular systems [31, 64]. The THCCSD reproduces the interaction forces quite well for intermolecular distances up to ~3 A, while AM1 is shown to slightly overbind. Thereafter the THCCSD underbinds slightly in disagreement with RIMP2 but in agreement with the OM1 Hamiltonian. The OM1 method has credence as a method which models the dipole and thus hydrogen bonding interaction properly by correcting the overlap matrix [32, 65, 66]. Adiabatic CN bond dissociation for small NMT clusters was investigated starting from published local minima [51] using the THCCSD and AM1. In studying the dimer interaction, for which one of the monomers is undergoing unrestricted CN bond rupture, the THCCSD forces predict that a bimolecular process takes place when the stretched bond is at a distance of 2.42 A (1.6Rgq). The reaction products are nitrosomethane, methoxy radical, and nitrogen dioxide radical. Figure 24 shows snapshots of this particular dimer decomposition channel. The THCCSD geometry optimization shows the monomermonomer distance decreasing from 3.50 A to 2.90 A as the methyl group on monomer one rotates to face an oxygen on monomer two. As this occurs, the CN bond on monomer two increases from 1.5 A to 1.6 A and the length of the NO bond participating in the reaction increases from 1.3 A to 1.7 A. The AM1 Hamiltonian does not predict any chemistry as a result of the CN bond breaking because the forces governing proper dissociation are incorrect, as shown in the case of the NMT monomer (Figure 21). The geometry optimization using AM1 reveals that the monomers actually separate by 0.45 A and that no methyl rotation takes place. The bimolecular dimer reaction predicted by the THCCSD is supported by the unimolecular rearrangement of intermediate IS2b [59] for which no transition state was located. The similarity is that the reaction products both involve the formation of methoxy radical via the cleavage of an NO bond. Unlike the study in [59], our process is bimolecular, where methoxy radical formed is by NMT undergoing oxygen elimination or NO rupture. Contrary to [59], where the methyl group picks up an oxygen from its own monomer, we observe that the oxygen is picked up by a methyl group on a different monomer. Our results for oxygen elimination are in agreement with previous studies [61, 62], which ii:. 1 this is the primary process in detonation. 2.3.3 Nitromethane Trimer The NMT trimer local minimum taken from [51] has a ring structure involving three hydrogen bonds. Performing the same study done for the dimer on the reference structure trimer from [51] shows that at a CN bond distance of 1.94 A or 1.3Req a concerted rearrangement of NMT to MNT occurs. The decomposition of NMT in trimer thus occurs at a CN distance that is 0.5 A less than that in the dimer. This indicates a cooperative reactivity in clusters of NMT initiated by the vibrational excitation of a CN bond. The monomer rearrangement to MNT is widely supported [58, 59] and is believed to be a key step in the detonation of NMT. The initial intermolecular NN distances are 4.5 A, 4.6 A, and 4.8 A and angles are 620, 580, and 600 (Figure 24). During the geometry optimization, because one of the monomers has a stretched CN bond, the system is unable to reach its ringlike hydrogenbonded local minimum, leading to ring strain. The THCCSD predicts that, as a result, the one monomer cleaves at a CN bond distance of 1.8 A and, subsequently, a second monomer dissociates at the same CN bond distance. At this point, all the intermolecular NN bond distances have decreased to 4.3 A, 4.5 A, and 4.6 A; the angles, however, remain unchanged. Geometry optimization with AMI shows that the stretched monomer shifts slightly from the remaining two monomers allowing one of them to rotate to a lower energy structure where it forms a bifurcated hydrogen bond each of which is 2.44 A. We have demonstrated that the THCCSD parameter set trained for NMT monomer transfers to the dimer and trimer calculations. The criteria for transferability used here regards the extent to which these forces yield products that are in accord with theoretically supported decomposition mechanisms of NMT [58, 59]. Note that the training of our Hamiltonian was not tailored to favor any particular rearrangement mechanism, evident by the fact that the dimer decomposition differs markedly from the trimer mechanism. Along the adiabatic surface AM1 does not predict any bond breaking to occur. Therefore, the accuracy of modeling bond breaking and forming processes in bulk simulations using AM1 is questionable, at best. We have shown a parameterization of the NDDO Hamiltonian that reproduces proper dissociation of NMT into radicals, both for the energy and forces. This parameter set reproduces the reference CCSD/TZP geometries and forces for unrestricted CN bond breaking. This NDDO parameterization would appear to potentially offer improved accuracy in MD simulations for clusters of highenergy materials. We also show that concerted reactions for both the dimer (bimolecular rearrangement) and trimer (trimolecular rearrangement) are not generated with an AM1 parameterization, yet are required for the proper description of NMT clusters. This new transfer Hamiltonian will allow for rapid and accurate MD studies due to the inherent speed of NDDO and the fitting of the parameters to highlevel CC reference forces. In this sense we are only limited by the reference data and the quality of the model Hamiltonian. We have demonstrated transferability of monomer parameters to nitromethane dimer and trimer. These forces predict chemical reactivity in clusters that is not observed in AM1 simulations and, in the case of the trimer, we see a threefold rearrangement of NMT to MNT giving further insight into the detonation process. 2.4 Silica and Silicon The TH depends upon the choice of the form of the Hamiltonian. So far we have used the NDDO form [67, 68] because of its wide use and prior successes. However, such methods using standard parameters like AM1 [27], MNDO [3032] and PM3 [28, 29] are also known to have significant failings compared to ab initio methods. Our TH, on the other hand, benefits from the fact that we only require specific parameters [45] for the systems of interest, for example, SiO2 clusters and H20 in our hydrolytic weakening studies [69]. These parameters permit all appropriate clustering and bond breaking for all possible paths for all the components, to allow MD to go where the physics leads. Furthermore, unlike traditional semiempirical methods, we do not use a minimal basis set, preferring to use polarization functions on all atoms, and in some cases, diffuse functions. The form of the Hamiltonian emphasizes twocenter interactions (but is not limited to nearest neighbors as in TB), which allows our parameters to saturate for relatively small clusters. We check this by doing ab initio calculations on larger clusters involving the same units to document that the TH is effectively unchanged with respect to system size. At that point, we have a TH that we can expect to describe extended systems composed of SiO2, and its solutions properly, including all multicenter effects on the energies and forces. Only the Hamiltonian is shortrange. We also check our TH by applying it to related systems. We will show that our TH determined from forces for bond rupture in disiloxane, pyrosilicic acid, and disilane also provide correct structures for several small cluster models for silica and silicon. Moreover, multiple bonds and high coordination Si, though not typical in silica polymorphs, are described well. Finally, we will show that our TH performs well using periodic boundary conditions for silica polymorphs. 2.4.1 Silicon Clusters A method based on rigorous quantum mechanics should be 'ii,, ferable. We define two types of transferability: first, an appropriate method should apply equally well to systems of similar size, but with different bonding characteristics, e.g., single and doublebonds; and second, a method should perform equally well for small and large systems, and even infinite systems as modeled by periodic boundary conditions. Now we will show that the THCCSD parameters for an NDDOtype Hamiltonian trained on a small set of reference molecules are transferable in the first sense. We chose our training set to be reasonable representatives of the local behavior of extended systems: disiloxane (H3SiOSiH3), pyrosilicic acid ((OH)3SiOSi(OH)3) and disilane (H3SiSiH3). The initial parameters used were from MNDO/d [70] because of availability, overall performance on the training set, and inclusion of dorbitals on Si. The fit was based on forces from the CCSD approximation using DZP basis sets. Primarily the THCCSD fit was to the forces along the intrinsic reaction coordinate for SiSi bond breaking in disilane and SiO bond breaking in both disiloxane and pi ilicic acid. Further refinement was achieved by fitting to the forces for the SiOSi bends of disiloxane and pyrosilicic acid. The most reliable error function to use in the type of fitting seems to be the square of the difference between the reference CCSD and THCCSD forces. The fitting procedure itself, though nonlinear, involved iteratively performing a linear minimization of each parameter, starting with those that most reduced the error. To verify the transferability of the parameters one should evaluate properties for molecules outside the reference set. Verification for bond breaking tests were performed on Si(OH)4, Si(SiH3)(OH)3, Si(SiH3)2(OH)2, Si(SiH3)3(OH), Si(SiH3)4, H2Si SiH2, (SiH3)SiSi(SiH3), (SiH3)2Si Si(SiH3)2, and (SiH3)3SiSi(SiH3)3. Verification testing for the OSiO bend involved calculations on Si(OH)4, Si(SiH3)(OH)3, and Si(SiH3)2(OH)2. Figures 25 and 26 show the results relative to DFT (with the B3LYP approximate exchangecorrelation functional) values. It is apparent that the THCCSD parameters outperform many of the commonly used NDDO methods for a variety of different bonding environments. Furthermore, although the parameters were designed to reproduce SiOSi, SiO, and SiSi forces, they transfer well to equilibrium OSiO (Figure 27). 2.4.2 Silica Polymorphs Now we demonstrate transferability of the THCCSD parameters to periodic systems. We choose to do so for < i i ,11ii,. forms of SiO2 first. Note that we need a model that describes amorphous and < ii i 11! phases equally well; if a method is comparatively more accurate (or less accurate) for crystalline or amorphous materials then nonphysical restrictions would be imposed. For example, crack tip propagation in silica is attracted to regions of low potential. Crystalline SiO2 is a good first step because < i ~1 Ii11;!i. phases have been observed at the interface in the otherwise amorphous silica gate oxide 1lv,r in metaloxide semiconductors (\ OS) [71]. Understanding silica and silicon brings us closer to understanding defects in MOS. In them, phenomena such as current leakage, electrical breakdown over time, and the thinning trends of the dielectric 1v.r in gate oxides are strongly influenced by defects [72, 73]. It also should be noted that we focus on training models that perform well primarily for localized units (i.e., clusters) and secondarily for bulk systems with periodic boundary conditions. This is because it is much harder to build local phenomena into a model used primarily for extended systems than it is to construct a model based on local chemistry that is applicable for small molecules, clusters, and bulk. For example, periodic calculations with plane waves require exceedingly large numbers of plane waves to describe oxygen defects in silica [74]. Furthermore, there is evidence that any model that is only dependent upon pairpotentials, (i.e., nearestneighbor distances) such as standard orthogonal TB schemes, lacks the ability to describe some fundamental properties of silica. It is difficult to prove such a claim, as we are able to imagine an expansion of the energy about some parameter, such as nearestneighbor distances. However, Stixrude et al. [75] have shown experimentally, based on the twofold compression of liquid SiO2, that there is essentially no change in nearestneighbor distances though there are significant changes in mediumr r',,' structure. This mediumrange structure is best described by ring and clusterstatistics and cluster geometry [76], which would be not be incorporated in any twobody potential. To avoid the limitations of empirical twobody potentials, one may consider NDDO or nonorthogonal TB methods discussed above, which include some manybody effects due to inclusion of the overlap matrix [77]. Especially promising is the SCCDFTB approach [7] which includes the largest contributions to longrange Coulomb interactions while handling charge in a selfconsistent way. These methods offer a practical balance of chemical accuracy and computational efficiency, as evidenced from their recent increase in popularity. Figure 28 shows the deviations of the THCCSD parameterization, AM1, and MNDO/d NDDO calculated unit cell densities from reference DFT densities. The reference DFT calculations for all silica polymorphs are from the VASP code using the gradientcorrected PW91 GGA exchangecorrelation functional with 3D periodicity, an ultrasoft pseudopotential, and a plane wave basis set (395 eV cutoff). A 2x2x2 kpoint grid is employ, 1 which converges total energies to within 0.05 eV/cell, and the convergence of selfconsistent steps is precise to 104 eV/cell. The optimization of cell parameters and atomic positions is precise to < 103 eV A1. It is clear that the THCCSD parameters trained on small silica clusters are transferable to several silica polymorphs. It is also interesting to note that not only is quantitative improvement achieved but the error from the PW91 DFT calculations is systematic and parallels similar improvement found when this new parameterization is applied to molecular systems. Further refinement of the THCCSD parameters must be made to accurately describe the siliconsilica interface, since the current parameterization fails to converge to the proper topology. The application of the transfer Hamiltonian methodology has been presented for nitrogen containing energetic materials and silicon containing compounds. In both cases we have demonstrated significant improvement over the underlying semiempirical method by reparameterizing with a training set of highquality CC forces for representative molecules. It is apparent that for focused MD studies that use semiempirical model Hamiltonians there is an advantage to creating a special set of parameters that is tailored to the chemical system being studied. However, the transferability of the model is affected by which molecules are selected for the training set and the PES reference points used. These models will often break down when applied to systems to which they have not been parameterized. One route to achieve this level of transferability is to design a simplified method that does not require parameterization, only controlled approximations, which leads to the AAT approach to defining a superior Th. Table 21. Equilibrium geometry for NMT, in A and degrees Method RCN RCH RNO AHCN AONC CCSD/TZP 1.498 1.087 1.219 107 117 THCCSD 1.504 1.084 1.225 108 118 AMI 1.500 1.119 1.201 107 119 0.05 0 0.05 0.1 CCSD/TZP UHF AM1 UHF TH(CCSD) UHF ................ B3LYP/631 G* UHF ........... 1 1.5 2 2.5 3 3.5 CN reaction coordinate (Angstrom) Figure 21. The NMT force curve for CN rupture 80 0 AM1 UHF o TH(CCSD) UHF ................ S60 B3LYP/631G* ...  40 20 0 1 1.5 2 2.5 3 3.5 CN reaction coordinate (Angstrom) Figure 22. The NMT energy curve for CN rupture 15 RIMP2/TZVP RHF A3 AMI RHF o OM 1 RHF ................ 1o TH(CCSD) RHF ...........  O p 5 1 1.5 2 2.5 3 3.5 4 4.5 OH reaction coordinate (Angstrom) Figure 23. The NMT dimer energies relative minimum in the intermolecular hydrogen bonding distance 4 *;:Y A^ V WAV W*. Ork 44*0 ^~ x^ A' 0^0%. Figure 24. Snapshots of the geometry optimizations for NMT dimer and trimer performed with the THCCSD and AM1. A) NMT dimer treated with THCCSD. B) NMT dimer treated with AM1. C) NMT trimer treated with THCCSD. D) NMT trimer treated with AMI D btS. W% 4 <tiO c Y ri I ~tr~T 0.08 0 0.06 0.04 0.02 AM1 PM3 MNDO/d TH(CCSD) Method Figure 25. Average deviation of SiO stretch 0.14 0.12 0.1 0 0.08 AM1 PM3 MNDO/d TH(CCSD) Method Figure 26. Average deviation of SiSi stretch Figure 26. Average deviation of SiSi stretch 6 S5 4 <3 0 AM1 PM3 MNDO/d TH(CCSD) Method Figure 27. Average deviation of OSiO angle 0.9 TH(CCSD)  0.8 AMI MNDO/d  0.7 S 0.6 0.5 0.4  0.3  0.2  0.1 0 Figure 28. Density deviation from PW91 DFT .0j ot^ CHAPTER 3 IMPETUS FOR A NEW SEMIEMPIRICAL THEORY Semiempirical method development has a long history dating back to the 1930's starting with Hiickel theory for describing evenalternant hydrocarbons. In the current study we have been guided by the developments over the last seventy years in regard to both what should and should not be approximated in our model. We focus here on the NDDO and SCCDFTB models. Both models have been widelyused with reasonable success. The relative numerical performance of the most recent versions of each model has been reported recently [2], and extensive comparisons of the formal similarities of each model have been made [78, 79]. It is important to understand the context in which these models were designed to fully appreciate the difference in the design philosophy of the current AAT model. NDDObased methods were originally constructed to reproduce the experimental heats of formation for small organic molecules. Since the heat of formation is defined only for molecules in their equilibrium structures, it is not surprising that NDDO models are parameterized specifically to reproduce equilibrium properties. For nonequilibrium structures the accuracy of NDDO methods varies greatly, as we will show. Similarly, SCCDFTB has been parameterized to reproduce equilibrium structures. The goal of the AAT approach is to achieve uniform accuracy of electronic properties, structures, and energetic over the entire PES to the extent possible for a minimal ba sis set effective oneparticle theory. Uniform accuracy is required to achieve meaningful results when performing MD simulations, in which transition states and other nonequilibrium structures p1l v a critical role. 3.1 Neglect of Diatomic Differential Overlap Based Methods 3.1.1 Does NDDO Approximate HF? A common misconception about NDDObased methods is that the NDDO model Hamiltonian approximates the HF Hamiltonian. This assumption is not the case and it can be shown that the density dependent contribution to the Focklike matrix in NDDO has little to do with the density dependent Fock matrix of HF theory, though they are deceptively similar: HF 1 P/ (pvAU)JABJCD 2 AlV)JADJBC (32) where p, v, A, and a are on centers A, B, C, and D respectively and the terms (puv Aa) are subject to the multipolemultipole expansion approximation. It is true that a typical NDDOtype twoelectron integral encountered is often larger in magnitude than a corresponding nonNDDOtype twocenter, threecenter, or fourcenter twoelectron integral. For molecules with more that three heavy atoms in a valenceonly minimal basis set representation, there are significantly more terms involving three and fourcenter twoelectron integrals. The combined effect of this large number of multicenter terms has a significant effect on the density dependent Fock matrix contributions. We consider a simple example. Shown in Figure 31 are the individual contributions to the NMT density dependent Fock matrix element of the carbon 2stype orbital and nitrogen 2stype orbital in a minimal basis set for the combined three and fourcenter, twocenter, ab initio NDDOtype, and full multicenter contributions. Though only a single matrix element is plotted over this reaction coordinate, the trend is the same for other matrix elements and other molecules. The three and fourcenter contributions are nearly equal and oppo site to the twocenter contributions. This cancellation is a universal trend and holds at all energy scales. And, while it could be said that the NDDOtype class of twoelectron integrals approximates the full set of twocenter twoelectron integrals, it is clear that the NDDOtype integrals alone bear little resemblance to the full twoelectron integral contribution. Instead, it is the delicate balance among all multicenter twoelectron integrals that is critical to obtaining the proper Fock matrix contribution. It has been known since studies were first performed on the NDDO approximation that a parameterized set of NDDO integrals perform better that the ab initio NDDO integrals. However, what has been less clear is the lack of correspondence between the NDDO approximation and the quantity it is supposedly approximating. Instead it seems that the adjustable parameters involved in the multipolemultipole expansion are being fit to the full density dependent part of the Fock matrix, and not the integrals themselves. Furthermore, this fit is nonlinear because of the twoelectron exchange integrals, which means that improving the NDDO approximation by modifying the integrals or trying to add three or fourcenter contributions becomes increasingly convoluted. For the AAT model, as we will show, all integral approximations can be systematically improved and do not involve any adjustable parameters. By avoiding the need to parameterize integral contributions we can retain uniform accuracy of the PES because no structural bias has been introduced. 3.1.2 Artificial Repulsive Bump There are many chemically interesting cases in which existing semiempirical methods seem to work quite well. For example Jorgensen et al. [80] show that their PDDG/PM3 reparameterization of PM3 yields heats of formation with a mean absolute error of 3.2 kcal/mol for a test set of 622 molecules. However, there is little evidence to sIl:. 1 that the same performance should be expected for chemical structures that are significantly distorted from equilibrium structures. Even for simple activation barriers, the error is typically 1015 kcal/mol [80]. For the determination of rate constants and in molecular dynamics studies the PES should be described within 12 kcal/mol. For example, since at 300K an error of 1 kcal/mol in the activation barrier will result in an order of magnitude change in the rate constant, clearly a 1015 kcal/mol error will not yield meaningful rate constants. The situation is especially bad in cases in which direct bond fission is treated with NDDObased model Hamiltonians. In such cases there is an artificial repulsive region near 1.2Req to 1.6Req, as is shown in Figure 32 for NMT UHF direct bond fission. Generally, in cases where it is known that there is barrierless dissociation, one would recognize such features as an error in the method and subsequent results would be questioned. For more complicated chemical processes, for example transition states or concerted reaction mechanisms, such dramatic errors are much more difficult to recognize, though it is precisely those regions that are often of primary interest. 3.1.3 Separation of Electronic and Nuclear Energy Contributions In ab initio theory (without consideration of nonBO effects) inclusion of the nuclearnuclear repulsion term (VNN) is a trivial matter when evaluating total energies and associated forces. While VNN has no part in the determination of purely electronic properties (such as the electronic density, electronic spectra, ionization potentials, or electron affinities) it is critical in the determination of molecular structure, the identification of transition states, and their associated vibration spectra. In simplified theories the role of VNN is obfuscated because the electronic and nuclear repulsion effects are not separated. The corecore repulsion term in NDDO Hamiltonians involves parameters that are designed for energetic at equilibrium. For such a specific set of properties the difference of two large numbers, the attractive (negative) electronic energy and the repulsive (positive) nuclearnuclear repulsion energy, is easier to model than each separately. In our opinion, combining such unrelated terms is the origin of the longstanding contention that electronic properties and total energy properties cannot be described within the same set of parameters. Also, once the corecore terms are parameterized, the correspondence between rigorous theory and semiempirical theory is lost. To solve this problem, we need to search for a better form for our model Hamiltonian that does not require that VNN be combined with terms that are purely electronic. 3.1.4 Total Energy Expression Here we illustrate the origin of the combined electronic and nuclear energy contributions of the total energy expression, the corecore term in NDDO methods, so that we may avoid this pitfall in the development of AAT. Consider the case in which the exact KSDFT exchangecorrelation energy functional and corresponding potential are known, and further are represented in a complete basis set. Then our effective oneparticle Hamiltonian will yield a set of orbitals which yield the exact density. The exact energy is given by the KSDFT total energy expression, Equation 127. Now, if we expect that the NDDO model Hamiltonian will generate the correct electronic density, which is a requirement if we want any of the firstorder electronic properties of our system, then we can relate the NDDO and KSDFT energy expressions. Making the temporary assumption that fNDDO j fKS, and consequently PNDDO PKS, then the total energy difference is NDDO KSNDDO \ Dr [ ZAZB ENDDO EKS= NDDO AB RAB) Exc[p] + 1 Vc(r)p(r)dr (33) A,B>A A,B>A A Now it is clear that if the NDDO oneparticle Hamiltonian were to yield the correct electronic density then the twocenter corecore term (ENDDO) must account for not only the simple nuclearnuclear repulsion term ( ), but also some pairbased contribution, dependent only upon RAB, of the the complicated multicenter exchangecorrelation contribution (Exc[p] f' vxc(r)p(r)). The motivation for approximating this last term is obvious: it drastically simplifies and avoids a complicated and expensive step in the determination of the total energy. However, there is no reason to expect a priori that this approximation would be accurate. For our improved model Hamiltonian we make a special effort to ensure that all approximations are controllable and systematically verifiable. 3.2 SelfConsistent Charge DensityFunctional TightBinding 3.2.1 SelfConsistent Charge DensityFunctional TightBinding in an AO Framework Equations 34 through 38 are the working equations for SCCDFTB which come directly from Elstner's original paper [7]: AqA qA qA (3 4) occ atoms qA i > > cCiS + cj cusS,/) (35) i pEA v M ci(Hs S) 0 (36) 1 atoms Hsc = (x HoIX,) X + S,, (7A + A.)Aq 2 t H0 + / 1 (37) ocC atoms E2TB o I 0I) + S 7ABAqAAqB + Erep(RAB) (3 8) i A,B A#B and the corresponding gradients are oCC OH o H S N a ,A B Erep(RAB) ni ci ( I a' 1 Aq, Aq a  (39) we have used p E A and v C B. To facilitate comparison of the SCCDFTB model Hamiltonian to other methods we force it into a structure that is comprised of a density independent part H and a density dependent part G. Using Mulliken population ,in ,i,~; to determine the net charge associated with atom A yields AqA qA + P ,, (310) pEA,v where the occupation numbers ni = 1 and the MO coefficients are real. FSCCDFTB H0 1+ H S atoms atoms = H,, S,, (7A + .)o + 2 5 (7A + .)PAuSA _i ( SC V^ pL /.^S(5\<7A / (Hs"c),, + PA(Gs ) (311) Aa we have introduced the notation that Satoms (H c),, H Sp, (A + .) (312) ~LV 2 ~L (Gscc)^ 2 SPS, A( + ) (313) Note that unlike HF or NDDO, in SCCDFTB it appears that it is not generally true that GG = Gt, as A E ( and a is unrestricted. If SCCDFTB were to have this property, that would be useful both for the calculation of forces and for comparison among all effective oneparticle methods we consider. To enforce this similarity we make the following observation: 2 1 SPAS AAC 7BC)] C SS(7AC + BC + 7AD + 7BD) (314) A<7 A<7 where p, v, A, and a are on A, B, C, and D respectively, and our working equation for the density dependent contribution in SCCDFTB is (Gsc) Sp, SA (AC + 7BC + 7AD + 7BD) (315) Equation 315 bears a similarity to the Mulliken approximation to Coulomb repulsion, we will address this in further detail later. Now we revisit the total energy expression. Using our new definition of the density dependent and independent terms in the SCCDFTB Focklike matrix. ESCC t [2Hc + P ,(GSC) + ABqAqB + Es c(RAB)) (316) Jv Ac A#B This energy expression closely resembles other oneparticle Hamiltonians, yet is equal to Equation 38. Similarly, after significant manipulation the gradient expression can be shown to be t1 sCC f[2 rSCCah 2 1 a ( 9ca o S cc + + Y ( t 7ABqAq + Escc(RAB)) (317) pv i A#B which again is equal to the original expression, Equation 39, but allows a more systematic comparison with other methods. 3.2.2 Mulliken Approximation Connection to SCCDFTB To facilitate understanding of the densitydependent term in SCCDFTB in relation to concepts in wavefunction theory we must recognize that the term AqA is the charge fluctuation on atom A from some reference charge for that atomtype (qA), with the charge on atom A being defined by the Mulliken population analysis. In the expanded representation in terms of the AO density matrix (P) the additional term to the density dependent part of the Focklike matrix is given by Equation 313, which has been shown to be equivalent to Equation 315, though the latter has a form that some may recognize as the Mulliken expansion for twoelectron integrals, 1 (pV Au7) SSoa(7YAC + 7BC + AD + 7BD) (318) The origin of the Mulliken expansion and the further extension of the basic approximation to the Riidenberg approximation will be explored later. For now it is sufficient that we discuss the numerical properties of this approximation and what advantages and limitations it entails. It is perhaps easiest to consider the nature of the Mulliken approximation by way of a simple example. For an NDDOtype twocenter twoelectron integral with normalized orbitals it follows directly from Equation 318 7AB (PA IABAB) (319) There have been many studies on the various forms for 7, most notably by Klopman [39], Ohno [40], and NishimotoMataga [41]. The particular form within SCCDFTB model is the Coulomb repulsion between two Is Slater type orbitals with an added condition that as the separation between the two centers approaches zero YAB  UA where UA is the Hubbard parameter, which is related to chemical hardness, which in turn is related to the difference between the first IP and EA of that atom. The particular form of the 7AB is only as significant as the underlying approximation in which it is used. In this case the Mulliken approximation ought to be accurate. The Mulliken approximation in SCCDFTB essentially converts three and fourcenter twoelectron integral contributions to the Coulomb repulsion energy term into overlapweighted twocenter repulsive terms. If we apply the Mulliken approximation in standard HF and use the ab initio (ppAA) twocenter integrals, the behavior is quite unphysical. Shown in Figure 34 are errors from the full HF when the Mulliken and Riidenberg approximations are applied to the two classes of threecenter terms, (AAIBC) and (ABIAC), and the fourcenter terms in NMT over the carbonnitrogen reaction coordinate. The error introduced by the Mulliken approximation is drastic, particularly for the (AAIBC) threecenter term, for which the error is over 1000 kcal/mol. Given that the only explicit density dependence of the Focklike matrix in SCCDFTB is introduced by the Mulliken term, and yet the relative energetic of SCCDFTB in practice do not introduce such errors, then some density independent part of the model must be accounting for the difference. We will avoid such uncontrolled approximations in the AAT approach. 3.2.3 SCCDFTB Repulsive Energy The SCCDFTB total energy expression relies on the original TB form for the total energy, in which a sum of repulsive atompair potentials is used. An error is introduced because of this reliance on atompair potentials to describe quantities that are inherently density dependent. Under an ideal parameterization, using the same argument as used earlier for NDDO, where we assume the effective oneparticle model SCCDFTB Hamiltonian returns the correct density, the error in the SCCDFTB energy to the exact form is pscc s KS SCC(R) E., [p] + v1t(r)p(r)dr C ZAZB (2Q E c EE c(RAB) v,([p]+ r)p(r)dr z (320) A,B>A A,B>A The numerical behavior of this error is shown in Figure 33. Clearly, this repulsive potential is carrying the burden of a portion of the electronic energy. This is not by accident; the value of a function (in this case the total energy) that results from the cancellation of two unrelated large functions (the electronic and nuclear repulsion energies) is easier to parameterize against than the values of the larger component functions. Alternatively, if the parameterization is against the larger functions there is a risk that any error introduced in the parameterization will be magnified when they are subsequently combined. If only the total energies or heats of formation (as in NDDO), and hence only some combined function, are of interest then such an average approach would be suitable. However, for our purposes we are interested in two separate and distinct sets of information, the electronic properties and the total energy. Because of this added demand on our semiempirical method we must generate the individual components of the total energy expression with equal accuracy. We recognize three distinct components: firstorder electronic energy, residual electronic energy, and nuclearnuclear repulsion. The firstorder electronic component has been the focus of most semiempirical method development, and yields, in addition to the firstorder electronic energy, the firstorder electronic properties: IPs, EAs, dipole, etc. The residual energy term presents a unique challenge, though we will demonstrate how it can be included in a computationally efficient way. The nuclearnuclear repulsion is trivial to calculate and the only complication arises when screened charges must be used if we use a valenceonly Hamiltonian. 3.3 Systematic Comparison of HF, NDDO, DFTB, SCCDFTB, KSDFT in a Single Framework One of the goals of this research is to bring many seemingly disparate approximate methods under a single common framework. For now, we will not deal explicitly with socalled postHF methods that include correlation based upon a reference single determinant (because of the required computational effort). Instead, we are trying to push the limits of what can be done within a oneparticle theory. The problem we are solving is the set of HartreeFockRoothaan matrix equations, FC = SCc, for a general oneparticle operator f. We use the notation that F = H + G, where H is the density independent contribution and G is the density dependent contribution to F. The spinfree oneparticle density matrix P is defined by a set of coefficients C for the real spinunrestricted case, N P^ z(c a + AC,:c (3 21) This dependence allows us to use a selfconsistent field (SCF) procedure: determine the Focklike matrix from a guess density and determine a new set of coefficients and the corresponding density, repeat until the density output is within some threshold of the input density. Alternatively, FPS PSF and (a lhfj) 0 would be good measures of i Ca convergence. Upon convergence of the SCF equations we have a set of AO coefficients that define the set of MOs = ZC^xv, (3 22) a reference single determinant wavefunction is then generated as an antisymmetrized product of the MOs, =o A{p1(1)02(2).../. (n)} (323) It is useful at this point to define the firstorder electronic energy given a zerothorder wavefunction, IoQ). The exact electronic Hamiltonian is ^ h(i+ (324) i i,j>i +3 where the core Hamiltonian is h(1) ZA (325) 2 A 71A and in its matrix representation the core Hamiltonian is H,`"7 I XJ(r1)h(r1)x,(rdr1 (326) The electronic energy correct through first order is &= (o01H 0o) Tr[PHcore] + (ijij) (327) i,j Note that the firstorder electronic energy is exactly the same as the HF energy expression and furthermore has no explicit dependence on the particular form of the oneparticle operator f. However, the goal is not the firstorder electronic energy but the exact electronic energy. We can consider the CC route, where we know that for a single determinant the exact electronic energy is (with any choice of orbitals) Ec = + Y jab)(t + t t tt ) + (if la)t? (328) i The CC energy expression provides great insight, but the explicit calculation of ta and ta amplitudes requires much more computational effort than we are able to expend for a fast semiempirical model. A more natural environment for a semiempirical methods is given by KSDFT. Assuming that we have the exact exchangecorrelation energy functional, the exact KSDFT energy is (for KSDFT orbitals) Eel =R + (ijji) +Exc[p] i,j = Tr[PHcre] + (ijlij) + E[p] (329) i,j There is another form of the KSDFT energy expression that is more useful for the current development, though it requires that we use a slightly more specific form for our oneparticle operator. If the target is a oneparticle Hamiltonian that returns the correct electronic density, and hence the correct firstorder electronic properties, then the ideal oneparticle operator is indeed the KSDFT oneparticle operator. And, if only the total operator and not its individual components, i.e., V2, Vext(r), VJ(r), and vxc(r), is pertinent then the firstorder energy expression given by Equation 329 would be suitable. By separating terms that are dependent on the density (G) we can write 1 1 EKS Tr[P(2HKS + GKS)] Tr[PG S] + Exc[p] (330) 2 2 where (GKS) c(r) v) (331) and F = H+G. Equation 330 is a central theme of our work, as it gives the unambiguous definition of the exact electronic energy. Tied to it are the oneparticle KS equations that define the orbitals, the density and associated firstorder properties. In this way, the problem for the total energy and gradients is completely specified by * F that generates the ground state oneparticle density matrix P. * F that can be decomposed into density independent and (H) and dependent (G) components. * The exchangecorrelation energy functional. 3.3.1 General Energy Expression We now choose a general electronic energy expression. Including the nuclearnuclear repulsion we have the total energy expression, 1 N ZAZ_ E = Tr [P(2H + G)] + Eridal[P] + ZAZ (332) 2 A,B>A AB We will bring HF, NDDO, DFTB, SCCDFTB and KSDFT into this framework, to make it clear what the three distinct components are: * Density independent H. * Density dependent G. * Residual electronic energy Eresidual[P]. If a method is to give the correct electronic density, and we can know the exact KSDFT HKS and GKS, then the model components for H and G Jl. l. must correspond to them. There is no flexibility with this interpretation. In addition, if we have the F that yields the correct density then to get the exact total energy, Eresidual [P] must be Ersidual [P]. Clearly, we are dealing with approximate models, though instead of defining the target or reference as being a single number, such as an experimental heat of formation, we instead have three separate and distinct functional forms, for which the exact forms are, in principle, known. 3.3.2 Density Independent FockLike Matrix Contribution This component generally accounts for nuclearelectron attraction and the electron kinetic contribution. The following definitions are based on p E A and v E B. * HF: HHF2 A1 1 ) (333) A * NDDO: H S u, 1 YBA ZB(ppBB) if (334) S/I (0/, + 0) if p ' * DFTB: HDFTB H= H (335) Pi/ 2 "P * SCCDFTB: N HSCCDFTB Ho S, (7A (336) Pi/'2 2 * KSDFT: Hv2s 2ZA  HPV (v IV) (337) 2 A r3A The only density independent contribution to the ideal Focklike matrix (the Focklike matrix which returns the correct oneparticle ground state density) is the exact core Hamiltonian, H0re". The 'li. 1 bottleneck for computing contributions for H" r" is the threecenter oneelectron integral, which scale as norbitalsNatoms. Still, they represent a small computational effort for large systems when compared to the orbitals terms for the threecenter twoelectron integrals. Though there is some cancellation, the possibility of complete cancellation among the threecenter oneelectron attraction and twoelectron repulsion integrals is limited. All threecenter oneelectron terms appear only in offdiagonal blocks of the Focklike matrix and are small compared to the dominant threecenter twoelectron terms which appear primarily in the diagonal blocks. 3.3.3 Density Dependent FockLike Matrix Contribution This component generally accounts for electronelectron contributions. * HF: Gcf P,((piv A ) A(IA v)) (3 38) Ac * NDDO: 1 GC DDO Px.((uA )6AB6D 2 (pA AD6Bc) (339) A7 * DFTB: G TB = 0 (340) * SCCDFTB: GSCCDFTB PA[ SpS, (AC + BC + .AD + 7BD) (341) Ac * KSDFT: G'fS = ( j[P] (r) + v.c[P](r) v) (342) The density dependent term of the ideal Focklike matrix is GKS. This choice follows directly from the chosen form for the density independent component and the requirement that we reproduce the exact electronic density. The density dependence of the Focklike matrix accounts for the complex interactions among electrons. This term is dominated by Coulombic repulsion. The second most important term is generally called exchange. The third component is correlation, and can be considered to be the highly complex way in which electrons interact with one another, i.e., their motion is correlated. For large systems, the density dependent part of a Focklike matrix is by far the most computationally costly. Unlike the density independent part, this term must be recalculated for every iteration of the SCF procedure. In addition, the terms involved here are more computationally intensive, e.g., the fourcenter integrals and the numeric integration required for the full treatment of vc. Generally, reduction in the computational cost of the density dependent term translates into an equivalent reduction in cost for the entire calculation. Therefore, the bulk of our attention is on systematic approximations that can be applied that simultaneously reduce the computational cost but maintain a desired accuracy. 3.3.4 Residual Electronic Energy This term places no restriction on the Focklike matrix and hence should have no direct effect on the electronic density. However, its accurate determination is critical for highquality electronic energies and consequently total energies and structures. * HF: E ridual 0 (343) * NDDO: UNDDO \ corecore Z Z\ B S (E DOcore (RAB) (344) residual RAB A,B>A * DFTB: EDFTB ZAZB residual (EA (RAB) AZ) (345) A,B>A AB * SCCDFTB: ESCCDFTB 0 SCC ZAZB ES(DFTB ABqOq +EAs (RAB)_ ZAZB) (346) residual + A) (3 46) A,B>AA * KSDFT: dual = Tr[PGs] + Exc[p] (347) 2 Because we insist that an approximate Focklike matrix reproduce the KS Focklike matrix, and thereby define the firstorder electronic energy, the residual electronic energy must be Ersidual. The residual electronic energy and its combination with the unrelated nuclearnuclear repulsion term in SE methods is notorious. We have attempted to include this term in a simple and computationally efficient way in our method. If our Focklike operator included exact (nonlocal) exchange then this term would introduce only the residual correlation energy contribution. Although the correlation energy is a small component (< 1 .) of the total energy, its accurate inclusion is required to obtain chemical accuracy. 3.3.5 General Gradient Expression One of the overarching objectives of our work is to create a model Hamiltonian to drive largescale MD. Achieving this goal requires an efficient method for obtaining gradients. Given the general energy expression provided by Equation 330, the gradient is [81] dE 2S _H ,,H, 1 8 8Eresid9al da aa aa 2 6a aa i Ac which is based on C C, (F/,, cS,,) = 0 and orthonormality. Note that there is no explicit dependence on !. Equation 348 only holds if terms in G are linear in the density G,, = PcAGZ (349) Ac The derivatives of Eresidual, which are not present in HF, are simple for those methods that only depend on a sum of atompair terms, such as NDDO, DFTB, and SCCDFTB. In KSDFT additional numerical integration is required for the perturbed density. Potential approaches for reducing the computational cost associated with the derivatives in KSDFT will be discussed in further detail later. 3.3.6 Requirements for an Accurate Simplified Method We have discussed the features of a general simplified method aimed at large systems that can be used to reliably obtain electronic properties as well as energetic. It is clear that the framework is essentially that of KSDFT, though we have structured the basic equations in a form amenable to further approximation. The two essential elements for the determination of electronic structure are the terms in the Focklike matrix that are either ,.ii ,,*///. l/ independent (H) or dependent (G) on the electronic density. The final element, which is needed to acquire accurate energetic and structures, is the residual electronic energy term. Furthermore, we have welldefined exact limits for each of these terms. When we discuss the AAT approach we have unambiguous connections to DFT and, from ab initio DFT, an analogous connection to WFT and the transfer Hamiltonian. 6 G3,4C Full 4 G2C 2 \ GNSS DO 6 t 2 12 ' ^ 6 1 1.5 2 2.5 3 CN reaction coordinate (Angstrom) Figure 31. Multicenter contribution to the C28N28 matrix element of NMT. 50 AM1 UHF o 40 30  20 lo 10 0 10 1 1.5 2 2.5 3 CN reaction coordinate (Angstrom) Figure 32. The AM1 artificial repulsive bump for NMT direct bond fission. 0.6 Etot 0.4 Erep E 3 _el ec ..... S0.2 0.2 0.4 0.6 / 0.8 1 1.5 2 2.5 3 CC reaction coordinate (Angstrom) Figure 33. The SCCDFTB total energy breakdown. 1200 1000 A 0 ~ Rudenberg Mulliken 4UU 200 0 1 1.5 2 2.5 3 3.5 CN reaction coordinate (Angstrom) 00 Rudenberg 90 Mulliken 80 70 60 50 40 30 20 10 1 1.5 2 2.5 3 3.5 CN reaction coordinate (Angstrom) Rudenberg Mulliken 1 1.5 2 2.5 3 3.5 CN reaction coordinate (Angstrom) Figure 34. Energies from Riidenberg and Mulliken approximations. A) Threecenter (AABC). B) Threecenter (ABAC). C) Fourcenter (ABCD). CHAPTER 4 ADAPTED AB INITIO THEORY We have specified the target framework of our simplified quantum mechanical method. The task now is to demonstrate simple wellcontrolled approximations that significantly reduce computational effort while maintaining a desired level of accuracy. To this end, a primary focus is on the density dependent component of the Focklike matrix, as it is the most computationally intensive part of the SCF procedure. Furthermore, the residual electronic energy contribution is examined because we have shown that its proper treatment is essential for a SM to achieve uniform accuracy of electronic structure and energetic. Our strategy is to start with a reference theoretical method that is reasonably wellbehaved, in this case KSDFT. Though there is much debate about the best DFT energy functional, we choose a hybrid DFT functional, specifically the widelyused B3LYP, because it includes a portion of exact exchange which adds a substantial degree of flexibility to our simplified method. To simplify the computational procedure we then explore various approximations and verify their accuracy. It may be too much to demand a simplified method that is systematically improvable (in the CC wavefunction theory sense), however, every approximation should be systematically verifiable. We will show that by making a couple of wellcontrolled approximations we are able to reproduce the topographical features of a hybrid DFT PES with a substantial reduction in computational effort. 4.1 TwoElectron Integrals Explicit calculation of twoelectron integrals is a computational bottleneck in DFT because the number of integrals scales as rbitals. For the Coulomb integral this can be reduced to nrbitals by fitting the density prior to the evaluation of the integrals [82], but this cannot be done for the exact exchange, leaving norbitals dependence. Though most quantum chemistry programs employ screening techniques to reduce the number of twoelectron integrals calculated, they remain a significant cost, especially when exact exchange is included. There are two v to lower the computational cost associated with twoelectron integrals: first, the number of integrals can be lowered by applying screening methods; second, the cost associated with the evaluation of each individual integral can be reduced. Several techniques for reducing the number of integrals have been successfully applied. They generally involve a procedure that recognizes a threshold for significant contributions to the Fock matrix. These include fast matrix multiple (FMM) [8385], the linear exact exchange [86], and quantum chemical tree code (QCTC) [87]. These methods still rely on the explicit calculation of integrals using standard techniques for Cartesian Gaussiantype functions. For large systems the most numerous class of integrals are the fourcenter integrals. Furthermore, fourcenter integrals are not only the largest class of integrals, but also the most computational involved integrals. According to Gill et al. [88] a fourcenter integral is typically an order of magnitude more computationally demanding than a twocenter integral. Our focus here is on systematic approximations that can be made to these ab initio multicenter integrals that will reduce the computational effort associated with calculating these terms while retaining acceptable accuracy. The cornerstone of our approach is based on ideas originally presented by Riidenberg in 1951 [5]. The underlying approximation is to expand an orbital on center A in terms of a set of auxiliary orbitals located on center B: 00 A(1) YZSkB, PA M() (4 1) k=1 where SkB, PA kB AA) (42) If the auxiliary orbitals on center B are complete and orthogonal then the expansion is exact. By repeated application of Equation 41 the number of centers that are involved in three and fourcenter twoelectron integrals can be reduced. The possible combinations of this type of expansion has been explored and enumerated by Koch et al. [89], though without comparing the numerical accuracy of such approximations. We will show the numerical precision of the Riidenberg approximation in an effort to avoid the explicit calculation of the numerous fourcenter integrals. By reducing the cost of each fourcenter integral we improve the most computationally expensive step of building the Fock matrix for large systems. First, consider the effect of the orbital expansion on a twocenter product of two orbitals with respect to the same electron coordinate, 10 0 (t(1)0, (1) 2 [SPAk ",. (1)vB(1) + SkAB ;. ( A() (43) k=1 Since a general multicenter twoelectron integral is given as (PAAc 'VBD) (IAB AcD) J B 1 (2)D dld2 (44) J J r12 substituting Equation 43 into Equation 44 for both the AB and CD orbital products yields the following expression, which is the Ridenberg twoelectron integral form, (iAVB ACD) SA[SikB SACD(kBVB ijDU7D) + SPAkB SJCD (kBVB I Acic) k=1 j=1 (45) +SkAV SAjD (/AkA jDc9D) + SkAvB SJCD (1AkA Acjc)] In the limit that the sums over auxiliary orbitals k and j are complete this expansion is exact. In practice, the expansion is approximated by limiting the orbitals over which we are expanding to those included in the AO basis functions, though using a separate auxiliary set of orbitals may also be suitable. The Ridenberg twoelectron integral expression (Equation 45), is more general than the wellknown Mulliken twoelectron integral approximation. The Mulliken approximation is equivalent to restricting the sum in Equation 41 to be only over orbitals of the same type, PA 1VB1 2SPAB (, 1) + B(1)> (1)] (46) which captures a large portion of the Riidenberg expansion and when applied to twoelectron integrals yields, 1 (PAVB AcIOD) SPAB3SACo [(PAPA Ac Ac) + (PAPA D7D) (47) +(VBV IAc AcA) + (VBVB U DU D)  Note the similarity of this expression to that which appears in SCCDFTB, Equation 318, in which 7AB ~ (PAPA VBVB) (48) though in SCCDFTB the orbitals p and v are restricted to be stype functions. We consider here all eight different approximations that result from the application of Equation 41 to three and fourcenter twoelectron integrals. The two classes of threecenter terms are of the form (PAVA IABcc) and (PAVB AAUc), which will be referred to as types A and B respectively. To facilitate our discussion we introduce the numbering convention in Table 41. The Riidenberg and Mulliken formulas can be applied to either class of threecenter or the fourcenter integrals (approximations IVI) using Equations 45 and 47, respectively. For the partial approximations to fourcenter terms (VII and VIII), partial means that the corresponding approximation has been applied only once. The expansion for these partial approximations is then in terms of threecenter twoelectron integrals (instead of only twocenter terms). For the Mulliken partial (VII) we have 1 (PAVB AcI7D) = { SAB [(PAA AAciD) + (VBVBA CUD)1 4 (49) +SACoD[(PAVB AcAc) + (PAVIB UDD)1}. and the Riidenberg partial (VIII) is 1 oo (PAVB AcoD) = 4 [SPAk (kBB IAcUD) + SkAB (pAkA AcUD) k= (410) + SAckD (PA VB I kDUD) + SkD (A VB I Ackc)] With all eight of the twoelectron integral approximations now defined, the accuracy of each will be considered. Since one requirement of our approach is that it reproduce the PES that would arise from a full ab initio treatment, we require that our approximations be uniformly accurate over the entire PES. This is in contrast with many commonly used semiempirical methods, in which the focus is typically on accuracy only about an equilibrium structure. However, uniform accuracy is critical to getting accurate forces to drive largescale MD and to aid in locating transition states for rate constant studies. Unfortunately there is no simple measure for uniform accuracy, due to the complexity of the 3N 6 dimensional PES, so our analysis will be statistical. Note that the total electronic energy is the sum of all the one, two, three, and fourcenter components. And, each of these components individually bear little similarity to the total electronic energy, again consider each contribution to any particular matrix element shown in Figure 31, instead the total energy is the result of the delicate cancellation among these terms. In addition, since only the nonparallelity errors are of importance in reconstructing the PES, the energetic contribution arising from each individual class of terms may be shifted by a constant value independently of one another without affecting the features of a PES. These subtleties, compounded by the 3N 6 degrees of freedom of the PES, make ascribing significance to any single class of terms difficult. Because we want to assess the accuracy of a particular approximation as it applies to a single class of terms we will instead focus on the Fock matrix elements, which contain all the information needed for the zerothorder electronic energy, to overcome some of these difficulties. Given any set of points on the PES (for example, a dissociation curve) and the Fock matrix for each point, each class of multicenter twoelectron integrals will contribute a particular percentage to the full density dependent Fock matrix. The best approximation for a particular class of integrals is then given by how well it reproduces the average percentage of the Fock matrix for the corresponding ab initio set of integrals. Consider the scatter plots of the full densitydependent contribution of the Fock matrix versus the individual contributions arising from the multicenter terms. Performing a leastsquares linear regression yields a line whose slope is a measure for the overall contribution from each term. Summing the one, two, three, and fourcenter slopes yields unity. An individual slope is, in a sense, the average percentage of the densitydependent part of the Fock matrix elements (entries). However, it should not be interpreted as the percentage of the energetic contribution. Because we are interested in not only the energy but also the density (which is given by the Fock matrix) then this average percentage serves as an acceptable measure of the importance of each term. Furthermore, in assessing approximations to individual classes of terms then we have an average way of measuring its significance. This is useful because the contribution from each multicenter component can vary greatly from molecule to molecule and also depending on internuclear separation. Take as an example NMT, and the dissociation of the CN bond. Plotted in Figure 41 are the scatter plots for each multicenter contribution to the Coulomb and exchange contributions of the Fock matrix from HF theory in a minimal basis. The CN reaction coordinate ranges from 1.0 A to 5.0 A in steps of 0.2 A The slopes of the corresponding bestfit lines of the one, two, three and fourcenter contributions are 0.0063, 0.5199, 0.4522, and 0.0216 respectively, note that these slopes represent the fraction of the full contribution, and that their sum is one. Table 42 shows the statistical averages of the densitydependent contributions to the Fock matrix for several molecules along the specified reaction coordinate (Rxn) at the HF level of theory in a minimal basis set. The dissociation in all cases is from 1.0 A to 5.0 A in steps of 0.2 A On average, the twocenter twoelectron integral component accounts for roughly half of the full matrix elements and the threecenter type A integrals recover another third. Again, the corresponding energetic contributions will not exhibit the same ratios. The purpose here is to validate individual approximations made to various multicenter twoelectron integrals. To this end we consider the average contribution: the better the approximation to a particular multicenter class of integrals the closer their average contributions will be. The a priori best approximation to the two types of threecenter contributions is that of Riidenberg. Shown in Figure 34 are the errors introduced by making approximations IVI, for NMT direct bond fission along the carbonnitrogen reaction coordinate. For type A and B threecenter integrals (3C A and 3C B) the average matrix contributions are 35.1'7' and 4'i'. respectively. The Ridenberg approximations II and IV differ by about half a percent at 35.1"' and 5.50'. Though the percent error for both approximations II and IV are roughly the same, the corresponding error introduced energetically is two orders of magnitude larger for the (AAIBC) threecenter terms than for the (ABIAC) threecenter terms. In contrast, the Ridenberg approximation applied to fourcenter terms (VI) is 2.47'. which differs from the fourcenter ab initio (4C) value of 2.3' I' by only 0 (i.' The error introduced by the full Riidenberg fourcenter approximation (VI) is surprisingly small. One would expect that the partial Riidenberg approximation would have the best agreement with the exact value, instead, at least empirically, the full Ridenberg approximation performs best, and has the smallest error relative to the average contribution to the Focklike matrix. Indeed, the effect this 0.0' difference has on the energetic of a dissociation curve is within an acceptable range and reproduces the energetic well for a variety of different dissociation curves. We will focus on the molecules in Table 42 that involve the dissociation of a CN bond. Plotted in Figure 42 are the dissociation curves for NMT, nitroethane (NET), COHNO2, and CH3NH2 at the B3LYP level of theory in a minimal basis set. There is a clear trend for the fourcenter Riidenberg approximation to overestimate the average contribution to the Focklike matrix. The energetic, on the other hand, do not show this clearcut trend since the energy is lowered for NET, raised for CH3NH2, and essentially unchanged for NMT. However, even though the PES are not precisely reproduced they still exhibit all the basic features of the ab initio surface, surprisingly without the explicit calculation of fourcenter twoelectron integrals. We have demonstrated that the Mulliken and Riidenberg approximations do not approximate the threecenter one or twoelectron integrals well. These classes of integrals pose a challenge as they are the only remaining terms that require explicit integration since all onecenter terms can be tabulated and twocenter terms can be interpolated by cubic splines, as will be demonstrated later. However, based on the Gaussian product rule these terms can, without approximation, be reduced to twocenter terms. In principle then cubic spline interpolation could be used to store the intermediates and completely avoid the need to do any integrals explicitly. The Gaussian product rule is a:(rRA)2 a(rRB)2 j (RARB)2 (a+a (IRA + RB)2 ./+.v ./+.v (4 i) We can generate an auxiliary set of orbitals for each Gaussiantype orbitalpair that have the orbital exponent a,,+a,. Also, because products of I= 1 (ptype) orbitals are involved this auxiliary set must also now span I = 2 (dtype). The prefactor is not included in the stored orbitals, instead it is computed onthefly to retain the twocenter character of the terms in which it will be used. There is a significant difficulty with using the Gaussian product in a simple way. Because we work in contracted Gaussiantype orbitals the number of orbital pairs increases as the number of contracted functions squared. The centers of these resulting product Gaussians are not necessarily the same, a A+a .B To use cubic spline interpolation an auxiliary orbital must have a constant exponent, e.g., a. + ac in the case of the simple product of uncontracted Gaussians. We envision an approximation to the exact contracted Gaussian product, which is a sum of several Gaussians with different centers and different exponents, that is a single contracted Gaussian located at a single center. The approximate auxiliary Gaussian could be determined by fitting directly to the exact product. However, this fit will generally yield different exponents and coefficients for two contracted Gaussian functions as there separation changes. If we consider the basis set STOnG, a sum of contracted Gaussians approximate a Slatertype orbital. There is no Slater product rule that reduces a product of two Slaters to a single Slater. Instead we have ajlrRARl9 arRB B a rRA~a &rRB (412) Consider a simple example of the product of two Gaussian functions versus two Slater functions separated by one Bohr, which is represented by the blue curve in Figure 43. The challenge is now clearer, the underlying basis functions we use are fit to Slater functions, and the product of two Slater functions does not resemble a simple Gaussian, or even a sum of contracted Gaussians. The Slater product is defined exactly by the intersection of two functions: ea lrRAleajrRB e ajrRcl e a6lrRnI (4 13) where, for a. < a, and RA < RB 7 a= a + a (414) 6 = a. a (415) Rc RA + aRB (416) ap + a, RD a /RA aVRB (4 17) And for this particular case RA < RC < RB < RD. For the Slater orbital product given in Figure 43 the exponents used were a1 = a, so ac = , as 0, R 0, and RD  o. The purple curve is ealrRcl which is truncated by a horizontal line that results from eas'rRDI. For a, / a, the situation is somewhat more complicated. However we now have a more direct path to twocenter intermediate terms that have a constant exponent for each atompair, and, are therefore, amenable to cubic spline interpolation with universally transferable parameters. The difficulty still remains on how best to incorporate the intersection of these two functions in a practical way. We have developed the strategy and framework for eventually including the threecenter terms. For now we are satisfied with reducing the computational cost through only the fourcenter Riidenberg approximation because the fourcenter terms account for most of the computational effort associated with construction of the Fock matrix for large systems. 4.2 KohnSham Exchange and Correlation If we envision an approximation to KSDFT then a primary focus ought to be avoiding the numerical integration that must be performed for each SCF iteration. The formal scaling of the numerical integration step is norbitls, as it is performed in real space. In practice, the scaling is lower because the exchangecorrelation potential is local and the density decays rapidly. Still, we want to avoid any extra computational effort that contributes little to the overall energetic. To approximate the exchangecorrelation potential contribution consider the expansion of the potential in terms of atomcentered densities, p = LA PA. Then we have the following telescoping series Vxc[ PAI](r) = vYxc[pAI(r) + {vc, [pA+ PBI](r) V,[PAI](r) V.cC[PBl (r)} A A A>B + S {vXc[PA PB+pc](r) A>B>C Vc[PA + PBI ) Vc [PA + c](r) Vc [PB + Pc] (r) + vxc[PA (r) + Vx[PB (r) + xc[pc] (r)} + .. 1center 2 center 3center _18) Vc 1+ Vc + Vzc + (4 8) Equation 418 is exact. It only requires an arbitrary partition of the density into atomic parts, though convergence of the series is affected by the choice of partition. Our target is to include contributions from the exchangecorrelation potential to the KS Focklike matrix in a way that only requires the storage of a limited set of twocenter terms. Such twocenter terms can be interpolated using cubic splines once, then the cubic spline coefficients can be stored in RAM to be efficiently retrieved and computed onthefly. With such a target in mind only terms in Equation 418 that can be retained are vcenter and v center. Furthermore, such contributions ought to be independent of chemical environment and thus completely transferable. Therefore, we insist on densities that are transferable. The natural choice is spherically symmetric atomcentered densities of neutral atoms, PA PA), recognizing we can have a and 3 spin components. Based on these considerations we have an approximate exchangecorrelation potential, xc[ PA](r) vxc[p(r) + [v [p +POB](r) ~Uc[po (r) vc[P](r)] (419) A A A>B We do not want to map and store the realspace exchangecorrelation potential, because of the obvious high degree of complexity. Instead, all that is ever needed are the contributions from the potential projected in the basis set. The Focklike matrix contribution from this approximate exchangecorrelation potential is (Gxc)PA [p3 (PA xc[P]I B) (420) Restricting this to contributions that only depend on two centers yields, for the diagonal block (A = B) (Gxc)PAVA (A vc[ZP1 VA A = (PA I c[AI A) (421) + (PA VLxc[P + PB xcI] IVA) B#A and for the offdiagonal block (A / B) (Gxc) PAVB lA vxc[ P B) A (422) (PA V xc[POA + POBI VB) The offdiagonal block terms can be further improved by applying the Riidenberg approximation. Note that the diagonal block terms do not benefit from this approximation because the orbitals only depend upon one center. The class of neglected terms (pA 1x2 [Pc] IB) becomes (LA \c PcB] = Y [SPAkB (kB ,[pU CVB) + SkA vB(PA c[p1 kA)I (423) C#A,B k The leading order terms are (oPA (PA I c[pOA]I A) if A B (Gc)/ A (4 24) (PA Ic[PA +PII VB) if A/ B and a slight improvement could be achieved with the addition of X)1 B#A(A A xc[POA + POB vxc[[POA VA) if A B S ZCA,B [SPAkB (kBv, [POCI]iB) + SkAVB (PA xc [PO] kA)] if A / B (425) Shown in Figure 44 are results for the GZERO(B3LYP) approximation: G6 is used with the B3LYP energy functional interpolated over all the atom pairs CNOH. The errors introduced are relatively small and GZERO(B3LYP) appears to be generally transferable. In terms of computational feasibility for GZERO we have introduced a single onecenter matrix per atomtype and a single twocenter matrix per atompair. For C, N, O and H atom types the number of onecenter matrices is trivial and does not present a challenge, as they are the dimension of the number of orbitals per atom squared. The twocenter matrices are stored in the local coordinate framework of the atom pair and are roughly the dimension of the number of orbitals on atom A times the number of orbitals on atom B. The actual dimension can be reduced because most of the matrix elements are zero due to orthogonality. To interpolate the values of these twocenter matrix elements cubic splines are used. Cubic splines require five parameters per number interpolated. We used adaptive grids (nonuniformly spaced) that are defined for every pair term from 0.5 A to infinity (effectively). To a good level of numerical precision this fit can be performed with roughly 600 reference points, but more can be added should higher levels of precision be required. For CNOH in a minimal basis set (there are ten unique atom pairs) splines require roughly half a million double precision numbers. Going to a double zeta type basis increases this storage requirement by about a factor of nine. Even for the double zeta case that cost does not present a challenge. 4.3 Adapted ab initio Theory Model Zero To briefly summarize, we have applied two novel approximations to speed up the construction of the Focklike matrix in a KSDFT framework. We have shown that the Riidenberg approximation applied to fourcenter twoelectron integrals retains the features of the ab initio integrals, thus avoiding the most costly class of Cartesian Gaussian integrals. Furthermore, this approximation is readily applied to hybrid KSDFT methods, which include some fraction of exact exchange. In addition to approximations for integrals we also demonstrated approximations that can be applied to the exchangecorrelation potential contribution of the Focklike matrix, and thereby avoid the need to perform numerical integration for every SCF iteration. Extensions to include threecenter contributions from the telescoping series approximation to the potential were also presented. The last component we will consider involves the residual total energy. Though it is worthwhile to avoid numerical integration to generate the potential and subsequent matrix contribution, a final onetime evaluation of Ec[p] upon SCF convergence should not be a bottleneck. Furthermore, as demonstrated by Bartlett and Oliphant [90] the exchangecorrelation energy functional is relatively insensitive to the input density. Indeed, as we have seen already in Figures 42 and 44 the B3LYP energy evaluation was used and introduced no significant error. The current implementation of AAT (AAT Model Zero) is given by HAAT ( Vi + ZZA) (426) A 1A all oneelectron terms are calculated exactly as they represent a small portion of the total computational cost for large systems. Given that we are using a hybrid KSDFT functional we have a fraction of exact exchange in the density dependent Focklike matrix, for B3LYP, a = 0.2. (Ck r I PA,(( A1) 1A )) (427) Af (GIea = ,, center ( + (( t) + ( GAv)) (428) AA where A (,, A V[pI1VA) if A (43 (PA Vxc[A + P VB) if A B AAT (GJaK ( JaK \ i + /(GO \ (43t) PV 1,2,3center)Pv + \(f4center) \Xzcv (4 1 and the residual electronic energy is EAd, = Tr[PG^] + Exc[P] (432) Finally, shown in Figure 45 are UHF dissociation curves for AAT Model Zero (AATMO), where both the Riidenberg and GZERO(B3LYP) approximations have been applied simultaneously. This model reproduces the underlying B3LYP dissociation curves well especially when compared to AM1. The corresponding gradient expression for AATMO is given by Equation 434. The only difference between standard KStype derivatives arise from the modified fourcenter term dependence on overlap. The evaluation of terms is straightforward and only dependent upon derivatives of the overlap, as, and twocenter NDDOtype integrals, 9(PAVA IBOGB) (lAu=) z y[)SkB {SACj (kBVB IjDD) + SjCD (kBVB ACjc)} k=l j=1 Oa (433) + S {S~A cD (AkA DD) + SJgD AkAI AC) +as4JD{SPkB(kBVB JDUD) +SkAVB(AkAiDUD)} (433) a(kBVBujDUD) a(kBVBAcjc) a(P{AklA jD.D) (SBAkA CjC) +SkAlB SACJD aa + SkAVB SJCGD a dE S {HAAT + (Gtc),} = 2 C, iici + P + da aa 8 a + 1 (434) 1 P {(GC)l1,2,3center + (GX)4center]} + Eresidua Y P .( ^ X^ + residual 2 0 a" 70 a Xo 4.4 Implementation The implementation of the approximations outlined above must be efficient because the construction of the Focklike matrix can be a rate limiting step for large systems. The Riidenberg approximation can be implemented in a stored integral or direct SCF procedure. For the Riidenberg approximation to be competitive requires that the contraction of the overlap with the appropriate twocenter twoelectron integrals be significantly faster than the fourcenter twoelectron integrals. Shown in Figure 46 are the relative timings for the evaluation of the ab initio fourcenter integrals versus the Riidenberg approximation. These timings are preliminary and more work could be done to optimize both methods for evaluating this fourcenter contributions. However, it is clear that there is a significant advantage to using the Riidenberg approximation and the possibility of improving the speed over traditional DFT methods by twoorders of magnitude is well within range. It is sufficient to focus only on timings for the fourcenter integrals as they quickly become the dominate contribution for large systems, as shown in Figure 47. Screening of small twoelectron integrals can lead to linear scaling for large systems. Such screening methods can be applied to the Riidenberg approximation very efficiently because of the explicit dependence on overlap. In an effort to minimize computational effort with building the Focklike matrix we have made extensive use of storing transferable onecenter and twocenter terms. In our current implementation all integrals up through threecenters are calculated exactly. There are two reasons for explicitly evaluating these integrals: first, they are a small contribution when compared to the explicit calculation of the fourcenter terms; second, the application of the Riidenberg or Mulliken approximations to threecenter terms (decomposing them into twocenter terms) introduces significant errors. If an accurate, simple and transferable approximation for the threecenter terms were available then all terms in the model could be twocenter. It is a straight forward matter to interpolate functions of two centers using cubic splines. We have implemented such techniques for the evaluation of the matrix elements of the exchangecorrelation: (G0),,, (RAB). Cubic splines are applicable for numerical 2center functions, where the exact functional form is not known, and for complicated analytic 2center functions. In the first case, the benefit of splines is to provide an interpolation scheme for a set of empirical reference points. In the second case splines provide a computationally efficient framework for the evaluation of otherwise prohibitively expensive analytic functions. We are primarily interested in the latter. The ability to store a large number of parameters in computer memory (RAM) is a requirement for the efficient implementation of cubic splines. The memory required is directly determined by the number of cubic spline parameters. A pertinent example is the evaluation of twocenter twoelectron integrals terms in MNDO, specifically, the multipolemultipole expansion of the NDDOtype twoelectron integrals. There are twentytwo unique NDDOtype integrals for every atompair involving two heavy atoms in standard NDDO methods such as AM1, PM3, etc. This means that for CNOH atom types, for which there are ten unique atompairs (CC, CN, CO, CH, etc.), there are fewer than 220 analytic functions. The evaluation of a cubic spline requires fewer operations than even the simplest of the twentytwo unique NDDOtype analytic expressions. As a simple demonstration of the speed of cubic splines we have implemented them to replace the explicit multipolemultipole evaluation for each of the 22 unique NDDOtype integrals. The multipolemultipole approximation is already significantly faster than the ab initio integral evaluation. As shown in Figure 48 a modest speedup is achieved over the already fast NDDO Focklike matrix construction. These timings are for 16 protein systems up to 20000 orbitals. Using cubic splines is up to 15'. faster than the multipolemultipole expansion. There is an error introduced, though it is small and controllable. For example the average error per atom is given in Figure 49. Furthermore, while going to higher multipoles in the expansion of integrals is increasingly expensive, cubic splines cost the same regardless of the underlying function being interpolated. However, the number of parameters involved in a typical NDDO multipolemultipole expansion is relatively small versus the number of parameters needed for cubic splines. The premium placed on RAM (or core) memory when NDDO methods were first developed thirty years ago would have made cubic splines impractical. Now, the amount of memory typically available on most desktop machines is more than sufficient to make cubic splines the method of choice when evaluating nearly any twocenter integral in quantum chemistry. A general cubic spline has the form f(x) = A + Bk(x Xk) + Ck(x xk)2 + Dk(x k xk)3 (4 35) where Xk is a reference point (node), and f(xk) = Ak. The determination of the parameters Bk, Ck, and Dk is straightforward and unambiguous. We have used the freely available Duris routine for discrete cubic spline interpolation and smoothing [91]. The only degrees of flexibility when using cubic spline interpolation are the number of nodes used and if they are evenly spaced (discrete) or based on an adaptive grid. In our case, for both integrals and matrix elements, the functions in all cases go to zero in the separated limit and at smaller length scales the function requires a dense grid to completely capture the subtle changes that occur in that region. A dense node spacing over the entire length scale is unnecessary because the functions change very little at long length scales. To determine a suitable node mapping function we tried several different functional forms and concluded that a simple exponential function provides a good balance between having a dense grid in the more critical bonding region (from around 0.5 A to 5 A) and a sparse grid at long bond lengths (hundreds of A). The actual node mapping function used is node(RAB) Int(1000 e6/(RAB+)) 1. (4 36) All of the functions we have currently implemented, those involving twocenter exchange correlation approximations, are zero beyond 10 A. Using the node mapping function then requires approximately 600 reference points. Since there are 5 variables per reference point we need 3000 parameters that would form a lookup table for a single function. If double precision numbers are used to store the spline parameters the storage requirement is about 0.023 MB per function. Further details of the current implementation of the AAT procedure in the parallel ACES III program system can be found in Appendix A. 4.5 Verification We have shown that AATMO reproduces well the basic features of a PES for a small test set of molecules involving the direct fission of a CN bond. Now we turn our attention to a more thorough comparison to test the behavior of AATMO for a variety of different properties. We use the diverse set of molecules in the socalled G2 test of molecules [92]. The G2 set is made up of 148 neutral molecules with singlet, doublet, and triplet spin multiplicities, though we focus on a subset of these that contain only C, N, 0, and H because the current implementation of AATMO includes only the spline parameters generated for those atom atompairs. There are 71 molecules in this truncated set ranging from 2 to 14 atoms. The separation of electronic and nuclear components in the energy expression of our AAT approach sets it apart from commonly used SE methods. So it follows that we expect more consistent results for purely electronic properties. To this end we have considered the vertical ionization potentials of the relevant molecules in the G2 set. The vertical ionization potential is purely electronic as it does not depend on the nuclearnuclear repulsion. We have used the simple definition Elp = E(N) E(N 1). This requires the total energy of the molecule and, since all molecules in the G2 set are neutral, the corresponding cation with the properly modified multiplicity. We have compared AATMO against the underlying B3LYP minimal basis set results and find a average RMSD from B3LYP for our entire test set of 0.057 Hartree (1.6 eV). The analogous comparison to AM1 yields an RMSD from B3LYP of 0.331 Hartree (9.0 eV). This clearly demonstrates the ability of AATMO to reproduce the electronic structure of the underlying B3LYP. The error of AM1 from B3LYP is substantial, though a more valid comparison for AM1 would be against experimental IPs, the important issue here is the ability our AATMO to reproduce the physics of the method it is approximating. We now consider the dipole moment. The dipole RMSD deviation of AATMO and AMI from minimal basis set B3LYP is 0.219 a.u. (0.56 Debye) and 0.595 a.u. (1.51 Debye) respectively. There were no sign errors observed for the orientation of the dipole moment for either AATMO or AM1. Though the error introduced by AATMO is significantly smaller than that of AMI, it is larger than would be expected given the performance of other properties generated by the AATMO procedure. Another test of a purely electronic property is simply the electronic density. In this case a direct comparison to existing SE methods is not possible in a straight forward way because NDDObased and DFTBbased methods drop the core orbitals since they are using a valenceonly approach. Instead, we have used B3LYP as the reference and compared AATMO, minimal basis HF with the Riidenberg approximation applied to fourcenter twoelectron integrals, and full ab initio minimal basis HF. Considering the oneparticle density matrix (AO density matrix P) from B3LYP the average RMSDs are 0.061, 0.136 and 0.029 for AATMO, HF (with Riidenberg), and full HF, respectively. These results demonstrate some surprising behavior. On average the full HF has the best agreement to B3LYP, at least in the case of a minimal basis set. The AATMO agreement is somewhat worse and HF with the Riidenberg approximation introduces a significant error. Since AATMO includes the Riidenberg approximation, we would expect its deviation from B3LYP to be similar to the deviation of HF with the Riidenberg approximation from the full ab initio HF results, especially since several systems in the test have only two or three atoms and the fourcenter approximation does not affect those cases. We see two probable origins of this unexpected discrepancy; first, there is a cooperative cancellation of errors between the Riidenberg approximation and the exchangecorrelation approximation (Go,); second, because the hybrid B3LYP functional includes only 211'. nonlocal exchange the corresponding error introduced by the Riidenberg approximation is not as significant when compared to the HF treatment in which full nonlocal exchange is used. In addition to the average RMSD of the density of AATMO and full HF we consider the standard deviation of the RMSD of the density for each molecule in the test set. The standard deviations are 0.042 and 0.057 for AATMO and full HF respectively. This  iI. r that while the overall average RMSD is somewhat higher from AATMO the error is slightly more systematic. Turning our attention to energetic properties we consider the atomization energies of the relevant portion of the G2 set. We take as the atomic reference energies the neutral atoms for each method. Since the AATMO is designed to reproduce the exact B3LYP result in the separated atom limit, the atomic energies are identical. So, the error between AATMO and B3LYP is equivalent to the total energy differences. Since we also want to compare B3LYP to AM1 we must also have atomization energy for B3LYP. Similarly, for AMI we perform the neutral atom calculations and subtract the sum of atomic energies from the total energy. The AATMO error from B3LYP is 0.057 Hartree. For AM1 the corresponding error is 0.331 Hartree. This indicates that the AATMO total energies are in better agreement with the underlying B3LYP total energies. Since one goal of our approach is to reproduce a complex PES of large systems we have applied the AATMO approximation to describe a somewhat unphysical reaction mechanism for C20 in order to subject the approximation to an extreme case of bond breaking. The pseudoreactionpath corresponds to splitting the C20 molecule in half, simultaneously breaking ten CC bonds. Shown in Figure 410 are the total energies of B3LYP, AATMO, AM1, and SCCDFTB. The reaction coordinate is the distance between the endcaps and the equilibrium value is 3.223 A, the range is from 2.8 A through 6 A. As expected we see that AM1 exhibits unphysical behavior, and some convergence problems beyond 4.5 A. The SCCDFTB actually performs quite well for C20, this result may not be too surprising since the repulsive energy term in SCCDFTB is parameterized from B3LYP total energies for CC bond breaking and the energy scale is very similar to the sum of CC bond energies. The AATMO overbinds this system significantly and may in fact be converging to a different electronic state, as the eigenvalue spectrum deviates significantly from the B3LYP result. Further work is needed to determine the precise origin of this discrepancy and to correct it. The final comparison we make is to the RMSD of the gradients for the G2 set. We will consider two sets of data, one for the MP2 equilibrium structures and the other for the B3LYP minimal basis set equilibrium structures. Shown in Figure 411 are the average RMSD of the Cartesian gradients from the MP2 reference structures calculated with AATMO, B3LYP, AM1, and SCCDFTB for the G2 set. The standard deviation of the gradients is also shown. AM1 appears to have the best agreement to the MP2 result with an RMSD of about 0.02 Hartree/Bohr, the corresponding standard deviation is also about 0.02 Hartree/Bohr. SCCDFTB does not perform well for the G2 set, however if we eliminate those molecules that contain two or three atoms from the set then the SCCDFTB performance improves significantly and the RMSD is 0.007 Hartree/Bohr, which lower than the AM1 RMSD of 0.017 Hartree/Bohr for this truncated set. In both cases the AATMO and B3LYP RMSD forces are nearly the same at 0.05 Hartree/Bohr, with a standard deviation of 0.01 Hartree/Bohr. A better comparison for how well the AATMO approximates the underlying B3LYP structure is shown in Figure 412. At the B3LYP geometries the RMSD is 0.05 Hartree/Bohr and 0.09 Hartree/Bohr for AM1 and SCCDFTB, respectively. We see the AATMO has an average RMSD of 0.006 Hartree/Bohr with a standard deviation of only 0.003 Hartree/Bohr. This clearly shows that AATMO is reproducing the B3LYP gradients with excellent accuracy. 4.6 Conclusion Current SE methods are not uniformly accurate over an entire PES, owing to severe approximations that are only valid near equilibrium. AAT has been shown to be able to describe bond breaking in its spinpolarized form at least as well as B3LYP, and much better than traditional SE methods, like AM1. The AAT SM framework enables the systematic verification of approximations via comparisons to CC and B3LYP references for prototypical molecules. Furthermore, all approximations can be further refined. Unlike existing SE model Hamiltonians, there are no adjustable parameters. Fast and accurate approximations to fourcenter twoelectron integrals via Riidenberg (for Coulomb or exact exchange) enable substantial savings for very large systems. Since the Riidenberg approximation merely reduces the cost associated with calculating each fourcenter twoelectron integral, methods for screening twoelectron integrals to achieve linear scaling are equally applicable here. The nonlocal exchange that is not considered in SCCDFTB is allowed in AAT, so that B3LYP and other hybrid DFT methods can be a target of the approach. We have also developed a route toward a fully twocenter model based on considerations of the threecenter twoelectron repulsion integrals. The Fock build requires no numerical integration, but the AAT approach does explicitly include exchange and correlation. In the current approach we have used the minimal basis set B3LYP exchangecorrelation potentials as a reference, further improvement is possible by evaluating large basis set ab initio DFT exchangecorrelation potentials and subsequently projecting them in a minimal basis set representation. Extensions of the AATMO that include a more complete expansion of the telescoping series for the exchangecorrelation potential have also been described and their implementation would be a further improvement on the results presented. Table 41. Adapted integrals (PAVA ABec) I Mulliken 3center Type A II Riidenberg (PAVBIAAUC) III Mulliken 3center Type B IV Riidenberg (IAVBIAc D) V Mulliken 4center VI Riidenberg VII Mulliken Partial VIII Riidenberg Partial OC 0 L00 ow a u nOi^o oc gOA~i Co CA O r C NOidLoi 1> o 1 ~~~~ Ln~ L y odLn~ 1 ~00V> 00~0~C> C>  TO ^ CM 1 lO O O C IA  3 ~~\ 7 ~"\ T ~ ^^ ~~\ T ~ cr3~~00 =1~" ti ~~oidoid~Lnoi  oo t t" cd > ,CIA tt 10"t CIA 0 L CIA cid d c 01 00000 d A ' d to B ' d Ci '3 0 D 'C D cd y = 0.0063x R2 = 0.0065 I 10 5 0 5 Full 2e integral (a.u.) 4 3 y, 2 R 0 1 2 3 4 10 4 3 y, 2 R I 0 1 2 3 4 10 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 10 =0.5199x = 0.8988 , r+ 5 0 5 Full 2e integral (a.u.) 0.45: = 0.8 22x 699 7 5 0 5 Full 2e integral (a.u.) y = 0.0216x R2 0.6841 / 5 0 5 Full 2e integral (a.u.) Figure 41. Scatter plots of agreement between approximate and exact twoelectron integrals. A) Onecenter. B) Twocenter. C) Threecenter. D) Fourcenter 60 Rudenberg (4C) S40' B3LYP 20 0 A 1) 20 40 60 1 1.5 2 2.5 3 CN reaction coordinate (Angstrom) 60 Rudenberg (4C)  S40 B3LYP 6 20 0 B 20 40 S60 80 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 60 0 I Rudenberg (4C)  S50 B3LYP 40 30 20 .i 20 C > io 10 S 0 S10 20 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 50 Rudenberg (4C) SB3LYP 0 D50 D I 100 1 1.5 2 2.5 3 3.5 4 CN reaction coordinate (Angstrom) Figure 42. Dissociation curve of CN with 4center Riidenberg approximation. A) NMT. B) NET. C) COHNO2. D) CH3NH2 3 2 1 0 1 2 3 Separation (Bohr) 4 2 0 2 4 Separation (Bohr) Figure 43. Orbital products. A) Gaussian orbital product. B) Slater orbital product 60 GZERO(B3LYP) 40 B3LYP S20 0 A A0 20 40 60 1 1.5 2 2.5 3 CN reaction coordinate (Angstrom) 60 GZERO(B3LYP)  40 B3LYP E 20 o 0 B 20 I 40 S60 80 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 60 GZERO(B3LYP) 50 B3LYP 40 30 20 , 20 10 S0 S10 20 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 50 GZERO(B3LYP) SB3LYP 0 S50 D 1 100 1 1.5 2 2.5 3 3.5 4 CN reaction coordinate (Angstrom) Figure 44. Dissociation curve of CN with AAT approximation. A) NMT. B) NET. C) COHNO2. D) CH3NH2 40 40  B3LYfP 0 A M 1 ................ 20 A _20 .\,"' 40 60 1 1.5 2 2.5 3 CN reaction coordinate (Angstrom) 60 AATMO 40 B3LYP 0 A M 1 ................ C 20 0 B 20 40 T 60 80 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 60 AATMO S50 B3LYP 40 0 AM1 ....... 40 S30 20 20 C 1o 0 S10 20 1 1.2 1.4 1.6 1.8 2 2.2 2.4 CN reaction coordinate (Angstrom) 50 i AATMO C B3LYP o AM1 ..... D 50 D S100 1 1.5 2 2.5 3 3.5 4 CN reaction coordinate (Angstrom) Figure 45. Dissociation curve of CN with AATMO approximation. A) NMT. B) NET. C) COHNO2. D) CH3NH2 40000 ab initio 35000 Rudenberg 2 30000 25000 20000 S 15000 10000 5000 0 400 600 800 1000 1200 1400 1600 Orbitals Figure 46. Timing of ab initio versus Riidenberg fourcenter integrals 1 0 0 ............................. .............................. 100 2center 90 2center 7 0 3center ............ S 760 i4center ................................ 60 50 30 q 20 10 ,, 10 1 0 1 1 1 1 1 1 1 1 0 100 200 300 400 500 600 700 800 9001000 Number of atoms (5 orbitals per atom) Figure 47. Percentage of multicenter terms versus system size 0 5000 10000 15000 20000 25000 Orbitals Figure 48. Timing of Fock build using traditional NDDO and NDDO with cubic splines as a function of system size 1.72e05 1.74e05 1.76e05 1.78e05 1.8e05 1.82e05 1.84e05 1.86e05 1.88e05 1.9e05 1.92e05 Error 0 5000 10000 15000 20000 25000 Orbitals Figure 49. Error introduced per atom from cubic splines as a function of system size 0.5 0 0.5 1.5 t 2 2.5 5 B3LYP 3 AATMO 3.5 A M 1 **................ SCCDFTB ................ 4  2.5 3 3.5 4 4.5 5 5.5 6 R (Angstrom) Figure 410. Pseudoreactionpath splitting C20 RMSD Standard Deviation N f S"cc Figure 411. Force RMSD from MP2 4 RMSD 0 2 Standard Deviation 015 01 S Figure 412. Force RMSD from B3LYP APPENDIX Parallel Implementation in the ACES III Environment The ACES III program system environment has been designed for the efficient parallel implementation of wavefunction theory methods in computational chemistry. The ACES III program is suited to treating large systems, typically 5001000 basis functions and up to 300 electrons with postHF ab initio methods. The super instruction processor (SIP) manages the communication and processing of blocks of data, such blocks are intended to be somewhat large (on the order of 500,000 floating point numbers). To provide a more direct route for implementing new methods the SIP reads the code written in the high level symbolic language, the super instruction ', ,,,11,; ; ,la,;ij,h: (SIAL) (pronounced  i!"). The SIP hides many of the more complicated aspects of parallel programing allowing the SIAL programmer to focus on algorithm and method development. Within SIAL each index of an array is divided into segments, the segments map elements of the array into blocks. An algorithm then involves operations for which the segment is the fundamental unit of indexing. Operations that involve the contraction of two arrays can then be broken into several smaller contractions over two blocks. Each blockpair contraction can then be distributed over many processors, allowing for the parallel implementation of the contraction of two large arrays. A general issue in parallel codes is ensuring that latency does not become a rate limiting step. This can be done by making sure that the amount of computation done is on a par with the latency of any operation. Such balancing is machine dependent, so it requires some amount of testing to determine an optimal solution for a particular machine. For the implementation of the AAT methodology this balancing is critical and requires a different treatment than is usually used for postHF methods. The structure of the AAT equations are quite different than postHF equations. In the former the atomic indices are used as the primary variable in the loops used to construct the Focklike matrix. For postHF methods the orbital indices are the important variable. Moreover, the number of orbitals per atom in AAT is small relative to the number of atoms, the reverse of many postHF calculations, which are generally used to treat systems with fewer atoms but with much larger basis sets. The speed of the AAT approach is achieved by avoiding the explicit calculation of the numerous and costly fourcenter integrals. Instead, these terms are introduced by a sum of contracted twocenter terms. Initially, we used small segment block sizes. The block size ultimately determines the computational efficiency by affecting how the IO, or fetching of data, is done. In these initial studies each block spanned only one atom. In this case, when the twoelectron integrals were evaluated it was straight forward to simply neglect the integrals that involved four atoms (fourcenters.) This procedure worked for systems with few atoms (<20), but, as the system size increased, the number of blocks rapidly increased, and latency became an problem. An improved approach that was implemented, which is more consistent with the original ACES III design philosophy, involved removing the atomspanning block size restriction. This was achieved by rewriting the original routines that returned all twoelectron integrals spanning fourblock units to explicitly exclude integrals involving fourcenters. In addition, the Riidenberg approximation, which relies on contractions between the overlap matrix and NDDOtype twocenter twoelectron integrals ((AAIBB)), required modification to avoid overcounting some terms. Similar to the routine that excludes fourcenter integrals, the modified routine to generate the integrals used by the Riidenberg approximation excluded all integrals that were not of the type (AAIBB). These two modified routines involve some overhead that could be reduced with further refinement. These modified routines are specific to the AAT approach and are not included in the standard release version of ACES III. Since twoelectron integrals possess an eightfold spatial symmetry the overall number of integrals that need to be calculated can be reduced by roughly a factor of eight. To do so requires recognizing the symmetry of every integral type. In ACES III the problem is slightly different. Instead of having this eightfold symmetry manifest itself by AO orbital index, it instead deals with block indices. In the pseudo code for the integral direct implementation given below, p, v, A, and a are block indices and each index generally spans several atoms. The following code for generating the twoelectron integrals evaluates every block of integrals only once, implying that the eightfold degeneracy is optimally incorporated. Do p Do v Ifv I End If (v = ) If v / p (ppI) ^ (ppV), (Iv1PP), (VPPP) Do A If A / p and A > v (pp vA) (pp\Av), (vA\pp), (AvIpp) (pvlpA) (pv\Ap), (vp\pA), (vp1Ap), (pAj\v), (Ap lv), (pAjvp), (Ap1VP) End If (A / p and A > v) End Do (A) End If (v / p) If v > p (pi vv) (vViPP) [Apply Riidenberg approximation to the terms (ppilvv) and Sj,k] (PVPv) (vPjPv), (vP VP), (Pv vP) Do A If A > p and A / v Do a If a > A and a / p and a/ v (pvA,\u) ^ (ivluA), (vpAu), (vpuA), (Au jiv), (u\Ajiv), (Auvp), (7AIVP) End If (a > A and a / p and au v) End Do (a) End If (A > p and A / v) End Do (A) End If (v > p) End Do (v) End Do (p) In addition to the Riidenberg approximation, we have also implemented an interface to the stored cubic spline parameters needed to construct the approximation to the exchangecorrelation contribution to the Focklike matrix. This component is much faster than the twoelectron integral component, therefore it is not a rate limiting step. Still, there are many aspects that need to be implemented efficiently to have a viable, stable code. We have previously discussed the use of cubic splines are an efficient way of balancing the numerical accuracy of the twocenter functions and the number of reference points stored. Another aspect that needs to be considered is the rotation of the twocenter matrix elements in a local atompair coordinate system to the gl..1, l' coordinate system of the molecule. This rotation is unambiguous once the global coordinate system is specified, and although orbitals of the molecule can be rotated simultaneously without affecting the energetic, the corresponding matrix elements will change on rotation. This degree of flexibility must be incorporated to ensure that the Focklike matrix contribution that arises from every atompair corresponds to the orbital orientation of the global coordinate system of the molecule. To achieve this we use a simple procedure that generates the matrix that rotates a pair of atomic coordinates so that they are aligned to the axis of the the stored cubic spline, and then perform the reverse rotation on that subblock of the Focklike matrix. Finally, the third component implemented was the numerical integration that is performed once the SCF procedure has converged. The ACES III program system generates a job archive file (JOBARC) that can be read by ACES II executables, assuming it has been generated on the same computer architecture. To test the efficacy of the AAT approximations, the quickest route was to modify the ACES II xintgrt executable. The xintgrt module performs numerical integration on a density to determine the exchangecorrelation energy that is needed for the residual electronic energy contribution of AAT. Effectively, xintgrt pulls the AO eigenvector matrix from the JOBARC, as well as other information about the system that would typically be generated by the I; . UI executable in ACES II when running a KSDFT calculation (e.g. occupation numbers). Since ACES III does not create all of the necessary records in the JOBARC, some small adjustments to the xintgrt routine led to a new modified executable that performed the needed numerical integration for the density generated by the AAT SCF procedure. The largest modification in this regard was the proper inclusion of the exchange terms. In the original xintgrt routine the Coulomb and exchange energy contributions are recalculated from the integral file (IIII). This integral file is not created by ACES III (because it works in an integral direct environment), removing the dependence on the IIII file and eliminating the revaluation of the Coulomb and exchange terms. Instead those intermediates are stored in the JOBARC by ACES III when they are evaluated in the normal course of evaluating the Focklike matrix. With these modifications all energy contributions needed for the AAT approximations becomes available. REFERENCES [1] R. J. Bartlett, V. F. Lotrich, and I. V. Schweigert, J. C'!, in Phys. 123, 062205 (2005). [2] K. W. Sattelmeyer, J. TiradoRives, and W. L. Jorgensen, J. Phys. ('C!. i, 110, 13551 (2006). [3] J. McClellan, T. Hughes, and R. J. Bartlett, Int. J. Quantum ('!,. i, 105, 914 (2005). [4] R. S. Mulliken, J. Chim. Phys. 46, 497 (1949). [5] K. Ruidenberg, J. C'!, in Phys. 34, 1433 (1951). [6] W. Weber and W. Thiel, Theor. ('I!. i,, Ace. 103, 495 (2000). [7] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998). [8] Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 50, 2138 (1994). [9] R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994). [10] A. Beste and R. J. Bartlett, J. ('C!. i, Phys. 120, 8395 (2004). [11] A. Beste and R. J. Bartlett, J. ('C!. i, Phys. 123, 154103 (2005). [12] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). [13] R. J. Bartlett and G. D. Purvis III, Int. J. Quantum ('!,. i, XIV, 561 (1978). [14] J. Cfiek and J. Paldus, J. C.'!, i Phys. 47, 3976 (1960). [15] J. Paldus and J. Cfiek, J. ('C!. _, Phys. 54, 2293 (1971). [16] R. J. Bartlett, Annu. Rev. Phys. ('C!. i, 32, 359 (1981). [17] R. J. Bartlett and M. Musial, Rev. Mod. Phys. 79, 291 (2007). [18] T. Helgaker, P. Jorgensen, and J. Olsen, Molecular ElectronicStructure Theory (Wiley, New York, 2000). [19] R. J. Bartlett, Ti,. *,y and Applications of Computational C', i,. .:ry (Elsevier, 2005), p. 1191. [20] M. Levy, Proc. of the Nat. Acad. of Sci. 76, 6062 (1979). [21] Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 50, 2138 (1994). [22] Q. Zhao and R. G. Parr, Phys. Rev. A 46, 2337 (1992). [23] Q. Zhao and R. Parr, J. Chem. Phys. 98, 543 (1993). [24] D. J. Tozer, V. E. Ingamells, and N. C. Handy, J. C'!i. in Phys. 105, 9200 (1996). [25] K. Peirs, D. Van Neck, and M. Waroquier, Phys. Rev. A 67, 012505 (2003). [26] R. Hoffmann, J. C'!I. ii Phys. 39, 1397 (1963). [27] M. J. S. Dewar, J. Friedheim, G. Grady, E. F. Healy, and J. Stewart, Organometallics 5, 375 (1986). [28] J. J. P. Stewart, J. Computat. C'!i. i, 10, 209 (1989). [29] J. J. P. Stewart, J. Computat. C'!i. i, 10, 221 (1989). [30] M. J. S. Dewar and W. Thiel, J. Am. C'!., in Soc. 99, 4899 (1977). [31] M. J. S. Dewar, E. G. Zoebisch, E. F. Healy, and J. J. P. Stewart, J. Am. C'I, in Soc. 107, 3902 (1985). [32] M. Kolb and W. Thiel, J. Computat. CI'. in 14, 775 (1993). [33] M. J. S. Dewar and W. Thiel, Theor. C'!., in Ace. 46, 89 (1977). [34] J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954). [35] W. M. C. Foulkes and R. Haydock, Phys. Rev. B 39, 12520 (1989). [36] R. J. Bartlett, I. Grabowski, S. Hirata, and S. Ivanov, J. ('!i, in Phys. 122, 034104 (2005). [37] M. Wolfsberg and L. Helmholz, J. CI. in Phys. 20, 837 (1952). [38] C. J. Ballhausen and H. B. Gray, Inorg. C'!i. i, 1, 111 (1962). [39] G. Klopman, J. Am. C'!., in Soc. 86, 4550 (1964). [40] K. Ohno, Theor. C'!I, i. Ace. 2, 219 (1964). [41] K. Nishimoto and N. Mataga, Z. Physik. ('!I, in (Frankfurt) 12, 335 (1957). [42] G. Z1i. w S. Irle, M. Elstner, and K. Morokuma, J. Phys. C'!, ii. A 108, 3182 (2004). [43] C. E. Taylor, M. G. Cory, R. J. Bartlett, and W. Thiel, Comput. Mat. Sci. 27, 204 (2003). [44] R. J. Bartlett, J. Phys. C'!. i,, 93, 1697 (1989). [45] I. Rossi and D. Truhlar, C'!., i, Phys. Lett. 233, 231 (1995). [46] P. C'!i rbonneau, Astrophys. J. (Suppl.) 101, 309 (1995). [47] R. H. Byrd, LBFGSB version 2.1, University of Colorado. (1997). [48] ACES II is a program product of the Quantum Theory Project, University of Florida. J.F. Stanton, J. Gauss, J.D. Watts, M. Nooijen, N. Oliphant, S. A. Perera, P. G. '. ,1 ,i, W.J. Lauderdale, S.R. Gwaltney, S. Beck, A. Balkov4 D. E. Bernholdt, K.K Baeck, P. Rozyczko, H. Sekino, C. Hober, and R.J. Bartlett. Integral packages included are VMOL (J. Alml6f and P.R. Taylor); VPROPS (P. Taylor) ABACUS; (T. Helgaker, H.J. Aa. Jensen, P. Jorgensen, J. Olsen, and P.R. Taylor) (2006). [49] W. Thiel, MNDO97, Organischchemisches Institut, Universitat Zirich, Winterhurerstrasse 190, CH8057 Ziirich, Switzerland (1997). [50] M. Dupuis, A. Marquez, and E. Davidson, HONDO 99.6, based on HONDO 95.3, Quantum C'!i ii,,I ry Program Exchange (QCPE), Indiana University, Bloomington, IN 47405 (1999). [51] J. Li, F. Zhao, and F. Jing, J. Comput. ('!., i, 24, 345 (2003). [52] R.Ahlrichs, M. Bar, M. Haser, H. Horn, M. Hser, and C. Klemel, ('C!. i. Phys. Lett. 162, 165 (1989). [53] M. Haser and R. Ahlrichs, J. Comput. ('!, ii 10, 104 (1989). [54] F. Weigend, M. Haser, H. Patzelt, and R. Ahlrichs, ('!, in, Phys. Lett. 294, 143 (1998). [55] P. Politzer and J. S. M. Eds., Theoretical and Computational C'i in,.Iry Energetic Materials Part 1. Decomposition, Cr;.ll and Molecular Properties, vol. 12 (Elsevier, Amsterdam, The Netherlands, 2003). [56] M. J. S. Dewar, J. P. Ritchie, and J. Alster, J. Org. ('C!. i. 50, 1031 (1985). [57] M. L. J. McKee, J. Am. ('C!. i,, Soc. 108, 5784 (1986). [58] M. T. Nguyen, H. T. Le, B. H i': tt6, T. Veszprimi, and M. C. Lin, J. Phys. ('!i. i_ A 107, 4286 (2003). [59] W. F. Hu, T. J. He, D. M. C'!i. i, and F. C. Liu, J. Phys. ('C!. i, A 106, 7294 (2002). [60] A. G. Taube and R. J. Bartlett, J. ('C!. i, Phys. accepted (2007). [61] W. D. Taylor, T. D. Allston, M. J. Moscato, G. B. Fazekas, R.Kozlowski, and G. A. Takacs, Int. J. C'!h in Kin. 12, 231 (1980). [62] H. A. Taylor and V. V. Vesselovsky, J. Phys. ('!. ii, 39, 1095 (1935). [63] R. Bukowski, K. Szalewicz, and C. F. ('! ,l. i,!wski, J. Phys. ('!., 11 A 103, 7322 (1999). [64] M. J. S. Dewar and W. Thiel, J. Am. ('C!. i, Soc. 99, 4907 (1977). [65] W. Thiel, J. Mol. Struc. (Theochem) 16, 398 (1997). [66] W. Weber and W. Thiel, Theor. ('!1. ini Acc. 103, 495 (2000). [67] J. Stewart, Reviews in Computational C', i,,.Iry (VCH Publishers, 1990), vol. 1, p. 45. [68] W. Thiel, Adv. Chem. Phys. 93, 703 (1996). [69] D. E. Taylor, K. Runge, and R. J. Bartlett, Mol. Phys. 103, 2019 (2005). [70] W. Thiel and A. Voityuk, Theor. Chem. Acc. 81, 391 (1992). [71] A. Ourmazd, D. W. Taylor, J. A. Rentschler, and J. Bevk, Phys. Rev. Lett. 59, 213 (1987). [72] D. Muller, T. Sorsch, S. Moccio, F. Baumann, K. EvansLutterdodt, and G. Timp, Nature 399, 758 (1999). [73] X. Y. Liu, D. Jovanovic, and R. Stumpf, Applied Physics Letters 86, 82104 (2005). [74] A. A. Demkov, J. Orteg, O. F. Sankey, and M. P. Grumback, Phys. Rev. B 52, 1618 (1995). [75] L. Stixrude and M. S. T. Bukowinski, Science 250, 541 (1990). [76] L. Stixrude and M. S. T. Bukowinski, Phys. Rev. B 44, 2523 (1991). [77] M. Menon and K. R. Subbaswamy, Phys. Rev. B 55, 9231 (1997). [78] M. Elstner, J. Phys. ('!,. i, A 111, 5614 (2007). [79] N. Otte, M. Scholten, and W. Thiel, J. Phys. ('!,. i, A 111, 5751 (2007). [80] M. P. Repasky, J. C'!I iiioh isekhar, and W. L. Jorgensen, J. Comput. ('!I, ii 23, 1601 (2002). [81] P. Pulay, Mol. Phys. 17, 197 (1969). [82] B. Dunlap, J. Connolly, and J. Sabin, J. C'!. in Phys. 71, 3396 (1979). [83] C. A. White, B. G. Johnson, P. M. W. Gill, and M. HeadGordon, C'!. in Phys. Lett. 230, 8 (1994). [84] C. A. White, B. G. Johnson, P. M. W. Gill, and M. HeadGordon, C'!. in Phys. Lett. 253, 268 (1996). [85] M. C. Strain, G. E. Scuseria, and M. J. Frisch, Science 271, 51 (1996). [86] J. C. Burant, G. E. Scuseria, and M. J. Frisch, J. ('!,. ii Phys. 105, 8969 (1996). [87] M. C'!I .!! ...,l>e and E. Schwegler, J. ('!I. ii Phys. 106, 5526 (1997). [88] P. M. W. Gill, A. T. B. Gilbert, and T. R. Adams, J. Computat. C'!. :ii 21, 1505 (2000). [89] W. Koch, B. Frey, J. F. S. Ruiz, and T. Scior, Z. N ,i w fi ch 58a, 756 (2003). [90] N. Oliphant and R. J. Bartlett, J. Chem. Phys. 100, 6550 (1994). [91] C. S. Duris, ACl\ TOMS 6, 92 (1980). [92] L. A. Curtiss, K. Raghavachari, P. Redfern, and J. A. Pople, J. Chem. Phys. 106, 1063 (1997). BIOGRAPHICAL SKETCH I was born in Michigan but have done time in Ohio, Washington, Oregon, and now Florida. My tumultuous career in science began in sixth grade when my science teacher called me a dreamer (I know I am not the only one) for proposing the construction of an elevator to the moon to avoid the costs associated with space flight. In hindsight, I should not have specified my original design as a rigid structure. That same year I was awarded a certificate as the classroom's str ,u. . kid. After making the transition to high school, I doubledup on math and science and finished all available classes in those areas by my Junior year, resulting in a misspent Senior year. Somehow I managed to get into the only college to which I applied. I arrived a Freshman filled with youthful optimism and persuaded my academic advisor to enroll me in a chemistry course 2 years beyond my level. I continued my Sophomore year without the youthful optimism. I participated in two Research Experience for Undergraduates programs sponsored by the National Science Foundation at Cornell and the University of Chicago. In Chicago I worked for Prof. Karl Freed, essentially spending the summer building overly complicated internal coordinate input files. I must have enjoy, 1 it because that same summer I tattooed myself with the mark of a quantum chemist, literally. Eventually, I took an interdisciplinary BA in chemistry and physics from Reed College in beautiful Portland, Oregon, in 2001. PAGE 1 1 PAGE 2 2 PAGE 3 3 PAGE 4 Ithankthechairandmembersofmysupervisorycommitteefortheirmentoring.IwouldalsoliketogivespecialthankstoMarkPonton,VictorLotrich,andDr.Deumensfortheirextensivehelpinthelastfewmonths.Ioweadebtofgratitudetomyfamily,whosesupporthasinstilledadeepyearningtopursuemyinterests.Severalpeoplemademanyothercontributionstomylife,bothpersonallyandacademically.Larryurgedmetoappreciatethenerthingslifehastooerandkeptmesane.Ithankthetworealchemists,TravisandLani,fortheirunendingquesttoanswerthebasicquestionmanasksoflife:Whatisre?ThepeopleoftheQuantumTheoryProjectoeredaconstantsourceofinspiration,andoftenasourceofmorningheadaches.IoermanythankstoAndrew,Tom,Dan,Christina,Georgios,Seonah,Martin,Joey,Lena,Julio,andofcourseMerve.IalsoacknowledgeProf.MerzfortheuseoftheDIVCONprogramandDr.MartinPetersforhelpingtoimplementandtestsomeinitialideasforcubicsplineinterpolation(Chapter4).SpecialthanksgotoAndrewTaubeandTomHughesforthenumerousextensivediscussionsthatwereinvaluableinmyresearch.Also,IacknowledgeTomHughes'workinimplementingthetransferHamiltonianparametersfornitromethaneclusters(Chapter2).IwouldliketorecognizethehelpfromDr.AjithPereraforhelpthroughoutmygraduatestudies. 4 PAGE 5 page ACKNOWLEDGMENTS ................................. 4 LISTOFTABLES ..................................... 7 LISTOFFIGURES .................................... 8 ABSTRACT ........................................ 9 CHAPTER 1INTRODUCTIONANDBACKGROUND ..................... 10 1.1WhoOrderedaNewSemiempiricalTheory? ................. 10 1.2QuantumChemistry .............................. 15 1.3HartreeFock .................................. 17 1.4CoupledClusterTheory ............................ 19 1.5KohnShamDensityFunctionalTheory .................... 22 1.6HuckelMethod ................................. 23 1.7NeglectofDiatomicDierentialOverlap ................... 25 1.8TightBinding .................................. 26 1.9DensityFunctionalTightBinding ....................... 27 1.10SelfConsistentChargeDensityFunctionalTightBinding .......... 29 2THETRANSFERHAMILTONIAN ......................... 35 2.1TheoreticalMethod ............................... 35 2.2ComputationalDetails ............................. 37 2.3NitrogenContainingEnergeticMaterials ................... 37 2.3.1NitromethaneMonomer ......................... 40 2.3.2NitromethaneDimer .......................... 41 2.3.3NitromethaneTrimer .......................... 42 2.4SilicaandSilicon ................................ 44 2.4.1SiliconClusters ............................. 45 2.4.2SilicaPolymorphs ............................ 46 3IMPETUSFORANEWSEMIEMPIRICALTHEORY .............. 58 3.1NeglectofDiatomicDierentialOverlapBasedMethods .......... 58 3.1.1DoesNDDOApproximateHF? .................... 58 3.1.2ArticialRepulsiveBump ....................... 60 3.1.3SeparationofElectronicandNuclearEnergyContributions ..... 61 3.1.4TotalEnergyExpression ........................ 61 5 PAGE 6 .......... 62 3.2.1SelfConsistentChargeDensityFunctionalTightBindinginanAOFramework ................................ 62 3.2.2MullikenApproximationConnectiontoSCCDFTB ......... 65 3.2.3SCCDFTBRepulsiveEnergy ..................... 66 3.3SystematicComparisonofHF,NDDO,DFTB,SCCDFTB,KSDFTinaSingleFramework ................................ 67 3.3.1GeneralEnergyExpression ....................... 70 3.3.2DensityIndependentFockLikeMatrixContribution ......... 71 3.3.3DensityDependentFockLikeMatrixContribution .......... 72 3.3.4ResidualElectronicEnergy ....................... 73 3.3.5GeneralGradientExpression ...................... 74 3.3.6RequirementsforanAccurateSimpliedMethod ........... 75 4ADAPTEDABINITIOTHEORY ......................... 80 4.1TwoElectronIntegrals ............................. 80 4.2KohnShamExchangeandCorrelation .................... 89 4.3AdaptedabinitioTheoryModelZero ..................... 92 4.4Implementation ................................. 94 4.5Verication ................................... 97 4.6Conclusion .................................... 101 APPENDIXPARALLELIMPLEMENTATIONINTHEACESIIIENVIRONMENT 117 REFERENCES ....................................... 122 BIOGRAPHICALSKETCH ................................ 127 6 PAGE 7 Table page 11TheNDDOonecentertwoelectronintegrals .................... 32 12TheNDDOtwocentertwoelectronintegrals .................... 33 13Multipoledistributions ................................ 34 21EquilibriumgeometryforNMT,inAanddegrees ................. 49 41Adaptedintegrals ................................... 103 42AveragepercentageofFockmatrixcontributionbyadaptedapproximation ... 104 7 PAGE 8 Figure page 21TheNMTforcecurveforCNrupture ....................... 50 22TheNMTenergycurveforCNrupture ...................... 51 23TheNMTdimerenergiesrelativeminimumintheintermolecularhydrogenbondingdistance ........................................ 52 24SnapshotsofthegeometryoptimizationsforNMTdimerandtrimerperformedwiththeTHCCSDandAM1. ............................ 53 25AveragedeviationofSiOstretch .......................... 54 26AveragedeviationofSiSistretch .......................... 55 27AveragedeviationofOSiOangle .......................... 56 28DensitydeviationfromPW91DFT ......................... 57 31MulticentercontributiontotheC2sN2smatrixelementofNMT. ........ 76 32TheAM1articialrepulsivebumpforNMTdirectbondssion. ........ 77 33TheSCCDFTBtotalenergybreakdown. ..................... 78 34EnergiesfromRudenbergandMullikenapproximations. ............. 79 41Scatterplotsofagreementbetweenapproximateandexacttwoelectronintegrals. 105 42DissociationcurveofCNwith4centerRudenbergapproximation. ....... 106 43Orbitalproducts. ................................... 107 44DissociationcurveofCNwithGZERO(B3LYP)approximation. ......... 108 45DissociationcurveofCNwithAATM0approximation. .............. 109 46TimingofabinitioversusRudenbergfourcenterintegrals ............ 110 47Percentageofmulticentertermsversussystemsize ................ 111 48TimingofFockbuildusingtraditionalNDDOandNDDOwithcubicsplinesasafunctionofsystemsize ............................... 112 49Errorintroducedperatomfromcubicsplinesasafunctionofsystemsize .... 113 410PseudoreactionpathsplittingC20 114 411ForceRMSDfromMP2 ............................... 115 412ForceRMSDfromB3LYP .............................. 116 8 PAGE 9 WepresentworktowardaoneelectronHamiltonianwhosesolutionprovideselectronicenergies,forces,andpropertiesformorethan1000atomsfastenoughtodrivelargescalemoleculardynamics.Ideally,suchamethodwouldbeaspredictiveasaccurateabinitioquantumchemistryforsuchsystems.TodesigntheHamiltonianrequiresthatweinvestigaterigorousoneparticletheoriesincludingdensityfunctionaltheoryandtheanalogouswavefunctiontheoryconstruction.Thesetwocomplementaryapproacheshelpidentifytheessentialfeaturesrequiredbyanexactoneparticletheoryofelectronicstructure.Theintentistoincorporatetheseintoasimpleapproximationthatcanprovidetheaccuracyrequiredbutataspeedordersofmagnitudefasterthantoday'sDFT. Wecalltheframeworkdevelopedadaptedabinitiotheory.Itretainsmanyofthecomputationallyattractivefeaturesofthewidelyusedneglectofdiatomicdierentialoverlapandselfconsistentchargedensityfunctionaltightbindingsemiempiricalmethods,butisinstead,asimpliedmethodasitallowsforpreciseconnectionstohighlevelabinitiomethods.Workingwithinthisnovelformalstructureweexplorecomputationalaspectsthatexploitmoderncomputerarchitectureswhilemaintainingadesiredlevelofaccuracy. 9 PAGE 10 Theunderlyingphysicallawsforthemathematicaltheoryofalargepartofphysicsandthewholeofchemistryarethuscompletelyknown,andthedicultyisonlythattheexactapplicationoftheselawsleadstoequationsmuchtoocomplicatedtobesoluble.P.A.M.Dirac Thoughhighlevelquantummethodsapproachthedesiredaccuracy,thecomputationalscalingofsuchmethodsprecludestheiruseforsystemslargerthanatensofatoms(e.g.,coupledclusterwithsinglesanddoubles(CCSD)scalesO(N6),whereNisthenumberofbasisfunctions).Incontrast,thecomputationalscalingofsemiempirical(SE)methodsisacceptable(typicallyscalesO(N3)),thoughtheycanbemadetoscalelinearlywithN.Forinstance,ifthesizeofasystemisincreasedbyanorderofmagnitudethenaCCSDcalculationwouldincurafactorof106morecomputationaltime.Forthesameincreaseinsystemsize,thecomputationaltimeofanormalSEmethodwouldonlyincreasebyafactorofathousand,butwithlinearscalingonlybyten.Intheregimeofverylargesystems(thousandsofatoms)itisapparentthatacalculationthatcouldbeperformedinonedayusingSEmethodsmighttakeseveralyearsforCCSD.Unfortunately,SEmethods 10 PAGE 11 Beforewebeginthedetaileddescriptionofourquantummechanicalmethodtheneedforafastandaccuratemethodtodescribethepropertiesoflargemoleculeswillbedemonstrated.PhenomenathatinvolvebondbreakingorbondformingoftennecessitateadetailedknowledgeoftheelectronicstructureandthePES.Suchquantumchemicaldescriptionsareimportantinmanycomplexchemicalprocesses;forexample: 11 PAGE 12 2.Condensedmatterstudies. 3.Biologicalstudies. Quantumchemicalmethodsareusedtosomeextentintheanalysisofthesetypesofproblems.Theidealquantumchemicalmethodcouldbeappliedtoanysizesystemandachieveanydesiredaccuracy,whileusingminimalcomputationalresources,but,ofcourse,thiscannotbeachieved.Hence,thereareoftenpracticallimitationsdictatedbythedesiredaccuracyandtheavailablecomputationalresources.Ourgoalistocreateamethodwithchemicalaccuracythatcanbeappliedtosystemswiththousandsofatoms.Amethodwiththesefeatureswouldbeabletoreadilyprovidethequantummechanicalcontributiontothesolutionofthesystemsmentioned. Thetwomainapproachesinquantumchemistryarewavefunctiontheory(WFT)anddensityfunctionaltheory(DFT).Likewise,therearetwocommonlyusedSEtechniques:neglectofdiatomicdierentialoverlap(NDDO)methods,whichoriginatedinWFT,anddensityfunctionaltightbinding(DFTB),includingitsselfconsistentchargeextension(SCCDFTB),whichoriginatedinDFT.BothNDDOandSCCDFTBcanbeconsideredsemiempiricalsincetheybothapproximatetheformalstructureoftheoriginalmethodsandareparameterizedtoreproduceeitherexperimentalortheoreticalreferencedata.JustastheformalconnectionsbetweenWFTandDFThaverecentlybeenuncoveredandexploited,viaabinitioDFT[ 1 ],therearesimilarformalconnectionsbetweentheir 12 PAGE 13 ExistingSEmethodsarewidelyusedand,whenappliedjudiciously,performquitewell.Regrettably,theyareoftenusedinsituationsforwhichtheresultsarenotreliable,suchaslocatingtransitionstates.ToexpandtheapplicabilityofSEmethods,wecombinetheirprimaryfeatureswiththosefromrigorousWFTorDFT,anddenethenewAATframework.Hence,AATencompassesbothNDDOandSCCDFTB. Despitetheirformaldeciencies,bothNDDOandSCCDFTBhavebeenquitesuccessfulinavarietyofapplications[ 2 ].Nonetheless,thereremainseveralaspectsoftheseSEmethodsthatfallshortofexpectationthatcouldbegreatlyimproved.Inadditiontoobviousformallimitations,existingSEmethodsweredesignedtoreproducetheenergeticsofsystemsatequilibriumandoftenyieldunphysicalresultsneartransitionstatesandalongdissociationpathways.UniformaccuracyofthePESiscriticalformoleculardynamics(MD)simulations,fordeterminingthelowestenergy(globalminimum)structure,andforcalculatingrateconstants.Associatedwiththesefailingsareunrealisticelectronicdensitiesthat,inturn,makecombiningthesemethodsintomultiscalemodelingschemesdicult.Furthermore,errorsinrstorderpropertiessuchasthedipoleandhigherordermomentswouldberemediedbyamethodthatdeliversadequateelectronicdensities.Attemptstoimprovesuchpropertieshavemetwithmodestsuccess,typicallyattheexpenseofotheraspectsofthemodel. Byfar,themostcommonroutetoincorporatesuchadditionalfeaturesistoreparameterizethemodel.Unfortunately,buildinginaccuracyfornonequilibriumstructures,excitedstates,andotherproperties,bysimplereparameterizationofexistingNDDObasedSEmethods[ 3 ]isnotreadilyachievable,atleastnotinawaythatmaintainstransferability(similaraccuracyregardlessofbondingcharacteristicsorsystemsize).ThisdicultyhighlightsanimportantdistinctionbetweenthecurrentAATrouteofdesigninganewsimpliedmodelHamiltonianbasedonrigorousabinitioprinciplesversus 13 PAGE 14 4 5 ]. Animportantcomputationalconsiderationistheinclusionoftheoverlapbetweenorbitals,whichisafeatureofAATthatdistinguishesitfromstandardNDDOmethods.Thezerodierentialoverlapapproximation,whichisthefundamentalapproximationofNDDO,impliesthatthereisnooverlapbetweenorbitals.Thoughthisapproximationenhancesthespeedofthemethod,itneglectsanimportantphysicalfeatureofelectronsandthebasissetsusedtodescribethem.TherehavebeenattemptstoreintroducesuchoverlapeectsinNDDOmethods,mostnotablyThiel'sorthogonalizationcorrectionmethods(OM1andOM2)[ 6 ].Still,Elstner'sdevelopmentofSCCDFTB[ 7 ]demonstrates 14 PAGE 15 Thepredominantcomputationalaspectsofthemethodhavebeenhighlighted;now,considertheformaldevelopmentoftheAATapproach.Thereisanexactreferencetotheexchangecorrelationpotentialwithwhichwecandirectlycompare.Theelectronicdensitycanbecorrectediterativelybyaninverseprocedure,suchasthatofZhao,Morrison,andParr[ 8 ]orthatofvanLeeuwenandBaerends[ 9 ].Theseprocedurestakeahighqualityreferencedensityandgeneratethenumericalexchangecorrelationpotentialthatbestreproducesit.Inthelimitthatthereferencedensityisexact,suchproceduresproduceanexactlocalexchangecorrelationpotential.Thisexactlimitdenesasetofoneparticleequations,whichinturn,denetherstorderelectronicenergy.Whenaddedtothenuclearnuclearrepulsiveenergy,theelectronicenergygivesthetotalenergycorrectthroughrstorderinthetrialwavefunction.InexistingSEmethods,theremaininghigherorderenergytermsarecombinedwiththenuclearnuclearrepulsiontogiveaparameterizedatompairbasedtotalenergycorrectionthatwouldbepredominantlynuclearrepulsion.Instead,wechoosetoretainaseparatenuclearnuclearrepulsioncontributiontoensurethateveryaspectofourmodelhasadirectonetoonecorrespondencetoabinitiotheory.Theconnectionstotheexactlimitsforbothdensityandtotalenergycomponentsoftheoverallparameterizationprovidethenecessarystructuresothateveryapproximationofthemodelcanberigorouslyveried. 15 PAGE 16 ThetimeindependentSchrodingerequationprovidesthesolutionsneededforthetimedependentSchrodingerequation,andthoughsomeproblemsrequireatimedependentsolution,thevastmajorityofproblemsinquantumchemistrycanbetreatedadequatelywiththetimeindependentSchrodingerequation.TheBornOppenheimerapproximationisbasedonthesimpleobservationthatelectronsmovemuchfasterthannuclei.ThisdierenceinspeedeectivelyimpliesnuclearmotioncannotaecttheelectronicstructuresignicantlyandwecansolvetheSchrodingerequationforasetofxednuclearcoordinates.Forexampleinthehydrogenatomtheelectronistravelingatabouttwomillionmeterspersecond,whichisordersofmagnitudefasterthanthemotionofthenucleus.Also,thisspeedislessthan1%ofthespeedoflight,whichalsoimpliesthatweneednotconcernourselveswithrelativsticeectsthatwouldonlybecomeimportantforatomswithheaviernuclei(typicallythosewithatomicnumberlargerthan25).ForsuchheavyatomsthespeedofanelectronnearthenucleuswillapproachthespeedoflightandwouldrequirearelativisticHamiltonianforanaccuratedescription.Insuchcaseseectivecorepotentialsthatapproximatethepassive(Darwinandmassvelocity)relativisticeectscanbeused,whichstillallowsonetokeepthenonrelativisticHamiltonianframework.Finally,thelinearcombinationofatomicorbitals(LCAO)approximationintroducesbasissetsintotheproblemandgivesageneralframeworkforsystematicallyexpandingthereferencespaceoftheorbitalsthattheelectronsmayoccupy. ThetimeindependentSchrodingerequationwithintheBornOppenheimer(clampednucleus)approximationis 16 PAGE 17 2Xir2iXA;iZA withiandjreferringtoelectronsandAandBtonuclei.ThetermsinEquation 1{2 aretheoperatorsthatcorrespondtothekineticenergyoftheelectrons,electronnuclearattraction,electronelectronrepulsion,andnuclearnuclearrepulsion.Thewavefunctiondescribesthestateofourquantumsystem.Wechoosetoapproximatewithinabasissetofatomicorbitals(AOs).EnforcingantisymmetryandsatisfyingthePauliexclusionprinciplecanbeachievedwithanantisymmetricproductoforthonormalizedAOswithinaSlaterdeterminant,theLCAOapproximation wherenisthenumberofelectronsandthespinorbital(r1)=(r1);and;denotespinand(r)isaspatialorbital. Itisttingtostartwiththemostwellknownapproximationinelectronicstructuretheory.HFisanexcellenttouchstonesincecorrelationenergyisdenedbyit(Ecorrelation= 17 PAGE 18 ^h(1)jii=(1 2r21XAZA andthetwoelectronoperatorsareCoulomb^Jandexchange^K, ^J(1)jjii=Zjj(r2)j2 ^K(1)jjii=Xj6=iZj(r2)i(r2) and ^f(1)jii=[^h(r1)+(^J(1)^K(1))]jii=ijii(1{7) denotingthemolecularspinorbitaliasalinearcombinationofatomicspinorbitals~i UsingLagrange'smethodofundeterminedmultipliersandenforcingorthonormalityoftheorbitalsgivesasetofoneparticlematrixequations,theHartreeFockRoothaanequations: where and TheconstructionofaFocklike(oreectiveoneparticle)matrixisacentralthemeofthiswork.AFocklikematrixisanorbitalrepresentationofaneectiveoneparticleoperator, 18 PAGE 19 10 11 ]. HFillustratesmanyofthenecessaryfeaturesofaoneparticletheory.Onedistinctfeature,whichsetsHFapartfrommostotheroneparticletheories,isthatHFenforcesthePauliexclusionprincipleexactlyviatheexchangeoperator.SuchexactexchangetermsaremissinginKohnShamdensityfunctionaltheory[ 12 ]whichleadstothesocalledselfinteractionerror.AnotherbenetofHFisthatthereisKoopmans'theoremthatgivesmeaningtotheeigenvaluesoftheFockmatrixbyrelatingthemtoionizationpotentialsandelectronanities.However,theexactDFTstructurecanprovideaformallycorrectoneparticletheorythatincludeselectroncorrelation,whichHFcannotdobecausebydenitionitdoesnotincludeelectroncorrelation.COPisanotherexactalternativethatcanbedenedfromWFT. PAGE 20 abij=A[1(1)2(2)::a(i)::b(j):::n(n)];(1{15) whereunoccupiedorbitalsaandbreplacetheoccupiedones,iandj.Thecoecient(amplitude)infrontoftheexcitationintroducesitsproperweightintothewavefunction.TheultimateversionofsuchasequencewouldincludeallpossibleexcitationsuptonfoldfornelectronsandiscalledthefullcongurationinteractionorfullCI.FullCIis: 13 ](meaningthatitscalesproperlywiththenumberofelectrons.) ThefullCIcanbewrittenconvenientlyinthecoupledcluster(CC)form, =exp(bT1+bT2+bT3+:::+bTn)0:(1{16) [ 13 { 16 ].ApproximationstoCCtheoryaremadebyrestrictingthebTpoperatorstosinglesanddoubles,asinCCSD,ortootherapproximationsthatdecouplehigherexcitationsfromlowerones,e.g.,CCSD(T).HerethetripleexcitationsareapproximatedfrombT1andbT2;whichcanbeobtainedfromCCSD,e.g.,withoutallowingthenewbT3estimatethentocontributebacktobT1andbT2:Thedetailsofthesemethods,whichcanbefoundelsewhere[ 17 ],arenotimportanttothepresentdiscussion.WhatisimportantisthatalltruncatedCCmethodsaresizeextensive,butdonotnecessarilygiveupperboundstotheexactenergy. ThefullCCequationsprovidetheexactelectronicenergy, 20 PAGE 21 (1{18) (1{19) andP12istheoperatorthatpermuteselectrons1and2.Theamplitudes,ftai;tabij;gcomefromprojectionofexp(bT)^Hexp(bT)= (1{20) (1{21) 17 { 19 ]. CCtheoryprovidesaroutetoexactreferencedata.ItiswellknownthatthehierarchyofCCmethodsaresystematicallyimprovable,andleadtothefullCIsolutionforaparticularbasisset.Hence,forsmallmoleculeswehavereadyaccesstoreferenceenergies,forces,densities,andelectronicpropertiestowhichwecancomparedirectly.Thewealthofinformationthatisaordedbyahighlevelabinitiocalculationissuperiortothetypicalamountofinformationavailablefromexperimentandprovidesuniqueinsightintohowwellanapproximationisworking.Forinstance,besidesgettingenergeticinformationaboutamoleculewesimultaneouslyhaveaccesstotheelectronicdensity.Thisincreasedlevelofcontrol,whencomparedtoexperimentalreferencedata,allowsforsystematicallyveriableapproximationsofoureectiveoneparticleHamiltonian.The 21 PAGE 22 12 ]equationsmaybewrittenas ^fKSjii="ijii(1{22) or, [1 2r2+ext(r)+J(r)+xc(r)]jii="ijii(1{23) where andxc(r)istheexchangecorrelationpotential.Constrainedsearchminimizationofthenoninteractingkineticenergyalonedeterminesorbitalsthatyieldthegroundstatedensity[ 20 ].Theothertermsdonotaecttheminimizationbecausetheyareexplicitfunctionalsofthedensity. TheHohenbergKohnShamtheoremsguaranteethatthedensityobtainedfromtheexactgroundstatewavefunctionwillbereturnedbyEquation 1{26 whenusingtheexactexchangecorrelationpotential.Inapragmaticsense,theseargumentsgiveaformalproofofexistence,thoughtheydonotguidetheconstructionoftheexactoneparticleoperator,or,morespecically,theexchangecorrelationpotential.However,incombinationwithsocalledinverseprocedures,suchasthatofZhao,MorrisonandParr(ZMP)[ 21 { 23 ],whichhavebeenexploredbyTozeretal.[ 24 ],ortheprocedurebyPeirsetal.[ 25 ],thepotentialmayberenediterativelysuchthatthesolutionstoEquation 1{23 willgenerateareferencedensity,presumablyahighqualityabinitioreferencedensity. 22 PAGE 23 2ZZ(r)(r0) wheretheexchangecorrelationpotentialisthefunctionalderivativeoftheexchangecorrelationenergyfunctional Unfortunately,currentlythereisnoprocedurethatallowsustotakeanumericallyderivedexchangecorrelationpotential,onethatisdeterminedsuchthatitreproducesahighqualityreferencedensity,andgeneratetheexchangecorrelationfunctional,whichisanexplicitfunctionalofthedensity.Thisinabilitytodeterminetheexchangecorrelationenergyfunctionalfromareferencedensityaloneistherstindicationofafundamentalprobleminsemiempiricalmethods,whichwillbediscussedinfurtherdetaillater:howdoesonedesignamodelHamiltonianthatwillsimultaneouslyachievethecorrectresultforbothelectronicpropertiesandenergetics?However,ifwearewillingtoassumetheformofaparticularexchangecorrelationenergyfunctionalthenthepotentialiswelldened,asistheeectiveoneparticleHamiltonian. HC=C(1{29) 23 PAGE 24 whereniistheoccupationnumberfortheithMO.NotethereisnoexplicitconsiderationofthenuclearrepulsioninEquation 1{30 .TheHuckelapproximationsare TheextendedHuckelmethod[ 26 ],asitsnameimplies,ismerelyanextensionofthestandardHuckelmethodwiththesametotalenergyformula.ThefullmatrixequationHC=SCissolved,withoverlapexplicitlyincluded.AndtheHamiltonianissubjecttoaslightlymoresophisticatedapproximation: whereHiiisbasedontheithAOofanisolatedatomandKisadimensionlessadjustableparameterforanatompair.HomannsuggestedavalueforKbasedonitseectontheenergyofethane(C2H6).Hearguedthat,sincetheenergydependenceonKisnearlylinearbeyondK=1:5,theMOenergiesand,inturn,thetotalenergy,willnotbetoosensitivetoKforvaluesgreaterthanthat.HethenconsideredtheenergydierencebetweenthestaggeredandeclipsedconformationsofethaneandoptimizedthevalueofKtominimizetheerrorbetweencalculatedandexperimentalvalues.ThenalvaluehearrivedatwasK=1:75,avalueisstillincommonuse. TheimportantaspecttonoteabouttheHuckelmethod,asitpertainstothedevelopmentofimprovedsemiempiricalmethods,isthataverysimplemodelcanreproduceafewfeaturesqualitativelywell.Thespatialsymmetryiscorrectduetotheinclusionoftheoverlapmatrix.However,despitetheabilityofsuchbasicapproximations 24 PAGE 25 27 { 32 ].Itisaminimalbasisset,valenceonlyapproach.NDDOisbasedonthezerodierentialoverlap(ZDO)approximation:productsoforbitalswithrespectthesameelectroncoordinatethatareondierentcentersaredenedtobezero, TheorbitalsareconsideredtobeLowdinsymmetricallyorthogonalizedAOs,sothattheoverlapmatrixistheidentitymatrix,andthesetofoneparticlematrixequationssolvedisFNDDOC=C.TheFocklikematrixinNDDOcontainsoneelectron(H)andtwoelectron(G)components.Theoneelectronpartisapproximatedas 2S(+)if6=(1{34) whereisonatomA,isonatomB,andUandareadjustableparameters,(jBB)modelstherepulsionbetweenaproductoforbitals()andapointchargeatcenterB,(h(1)j1 11 .TheparametersGss,Gpp,Hsp,Gpp,andGp2areadjustable. WithinalocalframeworktherearetwentytwouniquetwocenterNDDOtypetwoelectronintegralsforeachatompair(Table 12 ).TheNDDOtwocentertwoelectronintegralsareapproximatedbyamultipolemultipoleexpansion[ 33 ].Eachproductof 25 PAGE 26 13 ).AdjustableparametersarisefromthedistancesbetweenthepointchargesusedinthemultipoleexpansionandbyenforcingthecorrectlimitasR!0wheretheequivalentonecentertermshouldbereproduced. Asmentionedabove,originallytheNDDOFocklikematrixwasdesignedtomodeltheHFoneparticleHamiltonian.TheHFsolutionneglectsallelectroncorrelationeectsandonthescaleofchemicalaccuracyoneshouldnotaprioriexpecttheNDDOapproximationtoreproduceexperimentalresults.However,becausesomeparametersareadjustableandttoexperimentaldatathetypicalargumentisthatelectroncorrelationisimplicitlyincluded.Infact,itisclearthatelectroncorrelationisnomoreimplicitintheNDDOmodelthanarerelativisticeects,timedependence,nonBornOppenheimereects,basissetincompleteness,andexperimentalerrorsinthereferencedata.Despitesuchformalgaps,theNDDOmethodhasenjoyedatremendousamountofsuccessandreproducesheatsofformationandotherpropertieswithsurprisingaccuracy.Thedicultyariseswhenbetteraccuracyisrequiredinsituationsforwhichthemodelwasnotoriginallyintendedtodescribe,suchastransitionstatesorothernonequilibriumstructures.Insuchcasesonequicklyapproachesapointofdiminishingreturnsforparameterizingtonewreferencedata.Thepreciseoriginoftheselimitationsisexploredlater. 2XA6=BU(jRARBj):(1{35) wherethei'saretheeigenvaluesof 2r2+vext)jii=ijii;(1{36)vextisthenuclearelectronCoulombattractiontermandU(jRARBj)isconsideredtorepresentalltheremainingterms.Intotal,thelatterconstitutearepulsivefunctionofthe 26 PAGE 27 where (Heff)=hjbheffji(1{38) and Inpractice,thesetwocentermatrixelementsareevaluatedusingtheSlaterKosterparameterization[ 34 ].TheTBapproximationhasbeenappliedsuccessfullytomanysolidstateproblemsandworksespeciallywellwhenthedensityofthechemicalsystemiswellapproximatedasasumofsphericallysymmetricatomicdensities.Inthecaseofmolecules,amoreadvancedapproachthat\treatsthedelicatebalanceofcharge"[ 7 ]isneeded.IfwecompareTBtoHF,weseethatU(jRARBj)hastoaccountfor^J^KandVNN.InDFTitwouldaccountforJ,VNN,andwouldcombineexchangewithcorrelationbyaddingvxc. 35 ]deriveditfromrstprinciplesDFTconsiderations.TheenergyexpressioninKSDFT,Equation 1{27 ,maybewrittenas 2Z0(r0) 2NXA;BZAZB Equation 1{40 isexactif=Poccij2ij,whereistheexactdensity.However,itshouldbenotedthatiareeigensolutionsof^hs=(r2 2R0(r0) 27 PAGE 28 2ZZ000(0+) (1{41) +1 2ZZ00(0+) 2NXA;BZAZB 2ZZ02Exc 6ZZ0Z003Exc SubstitutionofthisexpressionforExcintotheenergyequationyieldsanexpressionfortheenergythatisformallyexact,yetonlyincludestermsthataresecondorderandhigherin. 2ZZ0000 2NXA;BZAZB 2ZZ00 2ZZ02Exc +1 6ZZ0Z003Exc PAGE 29 2ZZ0000 2NXA;BZAZB whichobviouslyiscorrecttosecondorderin.Notetheinclusionofnuclearnuclearrepulsioninthisterm. OneofthegreatadvantagesofDFTapproachesisthattheyoerareasonablezerothorderapproximationtotheenergysimplyfromknowingsomeapproximationtothereferencedensity,0:However,sinceExc[]isnotknown,gettingtheexchangecorrelationpotential,vxc(r)=Exc=(r);anditskernel,fxc(r;r0)=2Exc=(r)(r0),fromrstprinciples,whicharerequiredtosetuparigorousDFT,isdicult[ 1 36 ].However,werewetoacceptoneoftheplethoraofapproximationstoExc(LDA,LYP,PBE,B3LYP,B88,etc.)theequationscanbesolvedforafewhundredatoms,butnotonatimescalethatcanbetiedtoMDforthousandsofatoms. Furthermore,toprovideacorrectdescriptionofbondbreakingtheelectronicchargehastobefreetorelocate.AcommonexampleisLiF,inwhichatR=ReqwehavevirtuallyanexactLi+Fform,butatseparationinavacuum,wemusthaveLi0F0:SuchselfconsistencyusuallyisintroducedbysolvingtheAObasedDFTequationsviadiagonalization,butthisapproachimposesaniterativestepthatscalesasN3totheprocedure.Sincethemodelismeanttobeapproximateanyway,thealternativechargeselfconsistencyapproach,familiarfromiterativeextendedHuckeltheory[ 37 38 ],canbeusedmoreeciently.Thisextensionleadstotheselfconsistentchargedensityfunctionaltightbinding(SCCDFTB)approachofElstneretal.[ 7 ]. 29 PAGE 30 2ZZ0(0 2NXA;BqAqBAB:(1{45) TheqAarethedierencesbetweensomereferencecharge,q0A,andthechargeonatomA, 39 ],Ohno[ 40 ],andMatagaNishimoto[ 41 ].Forexample,thelatterformis wherethesingleindexAismeanttobeapurelyatomictwoelectronrepulsionterm.Thesearetypicallychosenasxedatomicparametersinthecalculation.ForanassessmentofhowtheSCCDFTBmodelperformscomparedtostandardsemiempiricalquantumchemicalmethods,seeMorokuma,etal.[ 42 ].Forexample,unlikeSCCDFTB,standardAM1doesnotgivethemostandleaststableisomersofC28whencomparedtoDFTcalculationswiththeB3LYPExcanda631G(d)basisset.Incontrast,therelativeenergiesoftheisomersofC28withB3LYP/631G(d)//SCCDFTBareinverygoodagreementwithB3LYP/631G(d),havingalinearregressionR2coecientof0.9925.ThoughthegeometriesgeneratedbySCCDFTBareexcellent,theSCCDFTBrelativeenergiesarepoor,withalinearregressionR2of0.7571. SCCDFTBisafairlyrigoroussemiempiricalmethod.TheframeworkprovidedbyFoulkesandHaydockexposedhowhigherordertermscouldbeincluded.ThepracticalimplementationofawellbehavedapproximationofhigherorderwassubsequentlyprovidedbyElstner.Later,wewilldeveloptheSCCDFTBformalismfromanMOviewpointtofacilitatethecomparisonbetweenitandothermethods. 30 PAGE 31 31 PAGE 32 TheNDDOonecentertwoelectronintegrals ParameterTwoelectronintegral 32 PAGE 33 TheNDDOtwocentertwoelectronintegrals TermTerm 1(ssjss)12(spjpp)2(ssjpp)13(spjpp)3(ssjpp)14(ssjsp)4(ppjss)15(ppjsp)5(ppjss)16(ppjsp)6(ppjpp)17(spjsp)7(ppjp0p0)18(spjsp)8(ppjpp)19(spjpp)9(ppjpp)20(ppjsp)10(ppjpp)21(ppjpp)11(spjss)22(pp0jpp0) 33 PAGE 34 Multipoledistributions TermChargedistributionPointcharges (ssjMonopole1(spjDipole2(ppjLinearQuadrapole4(pp0jSquareQuadrapole4 34 PAGE 35 However,practicallimitationsrestrictthefunctionalformofthewavefunctionandtheHamiltonian.ThechallengeistoincorporatethemanybodyeectsofexactWFTintoasimple(lowrank)modelHamiltonian.Oneaspectoftheapproachtakenhereistodevelopaformallyexactoneparticletheory,thetransferHamiltonian(Th),whichincludesallmanyparticleeects.AnotheraspectistoconnectasimpleapproximatemodelHamiltoniantotheexactTh.OurinitialworkusedthewellknownNDDOtypemodelHamiltonian.Becausethesemodelshavelimitedapplicability,ithasbecomeevidentthatmorecompletemodelsareneededthataccuratelyincorporatethemanyparticleeectsoftheTh.WedeveloptheAATsimpliedapproachtoincludethefeaturesoftheThinasystematicway,whichisdiscussedindetaillater. HereweintroducetheTh[ 43 ],aoneparticleoperatorwithacorrelationcontribution,asanextensionofthesimilaritytransformedHamiltonianfromCCtheory[ 44 ].ThisapproachprovidesarigorousformalframeworkforanexactoneparticleWFT.Also,unlikeothermethods,suchasthespecicreactioncoordinateapproachofTruhlar[ 45 ], 35 PAGE 36 Westartfromtheexponentialansatzinwhichtheexactgroundstatewavefunctionisgivenby where0isareferencefunction,whichcanbetakenasasingledeterminantandTistheexcitationoperatorofCCtheory.Whendenedthisway,wemaywritetheSchrodingerequationintermsofthesimilaritytransformedHamiltonian, ThesolutionofthisequationhasthesameeigenvaluesasthebareHamiltonian.Thetotalenergy,E,isexact,thustheforces,@E @RA,areexactaswell.Insecondquantizedform 2!Xp;q;r;s 3!Xp;q;r;s;t;u whereweshowtheinclusionofthreeandhigherbodyinteractionterms.Toincludethehigherordertermswerenormalizewithrespecttotheoneandtwobodytermsandobtainageneralized,correlatedFocklikeoperator,calledthetransferHamiltonian[ 10 43 ] wherePisthesingleparticledensitymatrixandbothehand^hpqjjrsicanbeapproximatedwithasuitablefunctionalofatombasedparameters.WehaveusedtheNDDOforminthecurrentwork,buttheformalismdevelopedthusfarisnotlimitedbythischoiceofFocklikeoperator. 36 PAGE 37 27 ]toreproducetheforcesdeterminedwithCCtheory.TheparameterizationwasperformedbyacombinationofthePIKAIA[ 46 ]implementationofageneticalgorithmandthelinearizedquasiNewtonRaphsonminimizationroutineofLBFGS[ 47 ]. ThereferenceforcesusedinthisparameterizationwereobtainedwiththeACESII[ 48 ]programsystem.WehaveusedCCSDinatriplezetabasis,withaunrestrictedreferenceandHartreeFockstabilityfollowing(toensurethatwearriveattheproperunrestrictedspinsolution.)AllsemiempiricalAM1andOM1calculationswereperformedusingamodiedversionofMNDO97Version5.0[ 30 { 32 49 ]. DFTcalculationsforthenitromethane(NMT)monomerwereperformedattheB3LYP/631G(d)levelwiththeHONDOPLUSprogram[ 50 ].TheinitialNMTdimerandtrimergeometriesweretakenfromaDFTstudyusingB3LYP/631++G**[ 51 ].Duetothelargenumberofcorrelatedreferencecalculationsneededinmodelingourclusters,weusedtheTurbomoleVersion5electronicstructurepackage[ 52 53 ]todoallelectroncalculationsattheresolutionoftheidentitysecondorderMllerPlesset(RIMP2)levelinaTZVPbasis.Wechoseanallelectronmodelinalargeauxiliarybasistoensureanaccuratetreatmentofcorrelationenergiesandgeometriesthatweredemonstratedonaseriesofcompoundsrepresentativeofatomsinvariousoxidationstates[ 54 ]. 37 PAGE 38 NMThasbeenthesubjectofmanytheoreticalandexperimentalstudiesduetothesignicantrolethatnitrogencontainingcompoundsplayinthechemistryofexplosives,propellants,andatmosphericpollution[ 55 ].ThedescriptionofthedecompositionofNMTanditsproductsiscrucialtotheunderstandingofthekineticsanddynamicsoflargerenergeticmaterialsforwhichNMTisamodel.NMTissmallenoughtoallowdetailedinvestigationofitscomplexdecomposition[ 56 57 ].Specically,evenforasystemofdeceivingsimplicity,itsdecompositiontoradicalsviasimpleCNbondruptureversusthecompetingNMTmethylnitrite(MNT)rearrangementhasbeenthesubjectofmuchdebate[ 56 { 59 ]. TheuseofclassicalpotentialsandstandardsemiempiricalNDDOmethodsoftenyieldincorrectbehaviorforenergies,forces,andgeometries,especiallyinnonequilibriumcases.Therefore,theuseofsuchtechniquestodoMDwilloftenleadtoincorrectresultswhenbondbreakingandformingarestudied.Alternatively,abinitioquantumchemistryisknowntobepredictive,inthesensethatitcanreproducemostquantitiesadequatelycomparedtotheexperimentalvaluesinthoseregionswithoutknowingtheresultpriortocalculation.However,thesemethodsaretooexpensiveand,therefore,currentlyimpracticalforMDbecauseofthetimeanddiskspacerequirements. PreviousNDDOparameterizations,suchasAM1,canreproducenearequilibriumstructuresforalargenumberofmolecules,yettheyoftenfailfornonequilibriumgeometries.Thus,forthepurposesofdoingMDsimulations,wearewillingtosacricetheabilitytoaccuratelytreatalargenumberofmoleculesnearequilibrium,fortheabilitytotreatasmallclassofmoleculesaccuratelyatallgeometries.Whenmodelingcombustionordetonationitisparticularlyimportanttohaveamethodthatworksequallywellforequilibriumandnonequilibriumgeometriesifthechemistryofbondbreakingandforming 38 PAGE 39 AmongthersttostudyNMTtheoreticallywereDewaretal.[ 56 ]whousedMINDO/3andfoundthattheNMTMNTrearrangementoccurredwithanenergybarrierof47.0kcal/molandthatMNTdissociatedtoH2CO+HNOwithanenergyof32.4kcal/mol.HencetheyproposedthatNMTdecompositionoccursviarearrangement.McKee[ 57 ]reportedtheNMTMNTrearrangementbarriertobe73.5kcal/molcalculatedattheMP2leveloftheorywitha631G(d)basisandthatMNTdissociatedtoH2CO+HNOwithanenergyof44.1kcal/mol.ThehighbarrierheightoftherearrangementsuggeststhatthedecompositionpathwayissimplebondrupturetoNO2andtheCH3radical.Huetal.[ 59 ]performedcalculationswiththeG2MP2approximationandfoundthattheNMTdissociationoccursviadirectCNbondrupturewithanenergyof61.9kcal/mol.TheNMTdissociationislowerthanboththeNMTMNTrearrangementandtheNMTacinitromethanerearrangementby2.7and2.1kcal/mol,respectively.Nguyenetal.[ 58 ]reportedthattheNMTMNTrearrangementbarrieris69kcal/molandNMTdirectdissociationtoCH3andNO2radicalsis63kcal/molcalculatedattheCCSD(T)/ccpVTZlevelusingCCSD(T)/ccpVDZgeometries.MostrecentlyTaubeandBartlett[ 60 ]reporttheNMTMNTrearrangementat70kcal/molusingCCSD(T). Huetal.[ 59 ]alsoperformedadetailedstudyofunimolecularNMTdecompositionandisomerizationchannelsusingDFTwithG2MP2//B3LYPleveloftheoryina6311++G(2d,2p)basis.SinglepointswerebasissetextrapolatedalongtheminimumenergypathstocorrectforincompletenessusingtheG2MP2method.Theirarticleshowsthatfromaparticularintermediatestate(IS2b)onepossiblereactionpathistheruptureoftheNObondtogivethemethoxyradicalandnitrousoxide.Huetal.indicatethaton 39 PAGE 40 58 61 62 ]. However,thereliabilityofboththeB3LYPand,toalesserextent,CCSD(T)atnonequilibriumgeometries[ 59 ]islimitedandthesameaccuracyshouldnotbeexpectedinregionsotherthanequilibrium. TheworkofLietal.[ 51 ]suggeststhatthethreebodyinteractionenergyinbulkNMTisacriticalcomponentinaccuratelydescribingthepotentialenergysurface.WeinvestigatedtheuseofourThinpredictinghowanaccuratetreatmentofthreebodyeectsdictateschemicalreactivity. 21 theagreementbetweenCCSDandTHCCSDistowithinadegreefortheanglesandlessthanonehundredthofanAforbondlengths. THCCSDisshowninFigure 21 toreproducetheoriginalforcecurveforCNbondrupturetowithinanaccuracyof0.005Hartree/Bohr.TheAM1parametersweredesignedtoreproduceequilibriumstructures,andthustheyhavethecorrectforceattheequilibriumbondlength.However,at0.1Aawayfromequilibriumtheforcesareinconsiderableerror.ArecurringfeatureforAM1forcecurves,whichisproblematicforMD,istheunphysicalrepulsivebehavioratlargeinternucleardistances.ThetransferHamiltonianremovesthisarticialrepulsionandgivesthequalitativelyandquantitativelycorrectdissociation. Onewayofdeterminingenergydierencesinasinglereferencetheory,otherthanHartreeFock,isbyintegrationoftheforcecurve.ForAM1theerrorintheforcesisreectedintheenergyasshowninFigure 22 .AlsoshowninFigure 22 theDFTsolution 40 PAGE 41 63 ]havestudiedtheinteractionenergyofdimersasafunctionofrigidmonomerseparation,wecalculatetherelativeenergiesoftheNMTdimerasafunctionofhydrogenbondingdistance.TheminimumchoseninourNMTdimercalculationsfrom[ 51 ]isstableduetotheformationoftwohydrogenbondsofequallengthandantiparallelarrangementofthemoleculardipoles.InFigure 23 ,weplottheenergiesrelativetotherespectivemethodsminimainROH. Ourdataindicatethatinthecaseoffrozenmonomergeometries,AM1overestimatestherepulsivewallatsmallintermoleculardistances.ItappearsasiftheH,C,N,andOparametersinthecorecorerepulsionfunctionforAM1wereunabletocanceltheeectfromwhichMNDOhastraditionallysuered,namelyoverestimationofenergiesincrowdedmolecularsystems[ 31 64 ].TheTHCCSDreproducestheinteractionforcesquitewellforintermoleculardistancesupto3A,whileAM1isshowntoslightlyoverbind.ThereaftertheTHCCSDunderbindsslightlyindisagreementwithRIMP2butinagreementwiththeOM1Hamiltonian.TheOM1methodhascredenceasamethodwhichmodelsthedipoleandthushydrogenbondinginteractionproperlybycorrectingtheoverlapmatrix[ 32 65 66 ]. AdiabaticCNbonddissociationforsmallNMTclusterswasinvestigatedstartingfrompublishedlocalminima[ 51 ]usingtheTHCCSDandAM1.Instudyingthedimerinteraction,forwhichoneofthemonomersisundergoingunrestrictedCNbondrupture,theTHCCSDforcespredictthatabimolecularprocesstakesplacewhenthestretchedbondisatadistanceof2.42A(1.6Req).Thereactionproductsarenitrosomethane,methoxyradical,andnitrogendioxideradical.Figure 24 showssnapshotsofthisparticulardimerdecompositionchannel. 41 PAGE 42 21 ).ThegeometryoptimizationusingAM1revealsthatthemonomersactuallyseparateby0.45Aandthatnomethylrotationtakesplace. ThebimoleculardimerreactionpredictedbytheTHCCSDissupportedbytheunimolecularrearrangementofintermediateIS2b[ 59 ]forwhichnotransitionstatewaslocated.ThesimilarityisthatthereactionproductsbothinvolvetheformationofmethoxyradicalviathecleavageofanNObond.Unlikethestudyin[ 59 ],ourprocessisbimolecular,wheremethoxyradicalformedisbyNMTundergoingoxygeneliminationorNOrupture.Contraryto[ 59 ],wherethemethylgrouppicksupanoxygenfromitsownmonomer,weobservethattheoxygenispickedupbyamethylgrouponadierentmonomer.Ourresultsforoxygeneliminationareinagreementwithpreviousstudies[ 61 62 ],whichsuggestthisistheprimaryprocessindetonation. 51 ]hasaringstructureinvolvingthreehydrogenbonds.Performingthesamestudydoneforthedimeronthereferencestructuretrimerfrom[ 51 ]showsthatataCNbonddistanceof1.94Aor1.3ReqaconcertedrearrangementofNMTtoMNToccurs.ThedecompositionofNMTintrimerthusoccursataCNdistancethatis0.5Alessthanthatinthedimer.ThisindicatesacooperativereactivityinclustersofNMTinitiatedbythevibrationalexcitationofaCNbond.ThemonomerrearrangementtoMNTiswidelysupported[ 58 59 ]andisbelievedtobeakeystepinthedetonationofNMT.TheinitialintermolecularNNdistancesare4.5A, 42 PAGE 43 24 ).Duringthegeometryoptimization,becauseoneofthemonomershasastretchedCNbond,thesystemisunabletoreachitsringlikehydrogenbondedlocalminimum,leadingtoringstrain.TheTHCCSDpredictsthat,asaresult,theonemonomercleavesataCNbonddistanceof1.8Aand,subsequently,asecondmonomerdissociatesatthesameCNbonddistance.Atthispoint,alltheintermolecularNNbonddistanceshavedecreasedto4.3A,4.5A,and4.6A;theangles,however,remainunchanged.GeometryoptimizationwithAM1showsthatthestretchedmonomershiftsslightlyfromtheremainingtwomonomersallowingoneofthemtorotatetoalowerenergystructurewhereitformsabifurcatedhydrogenbondeachofwhichis2.44A. WehavedemonstratedthattheTHCCSDparametersettrainedforNMTmonomertransferstothedimerandtrimercalculations.ThecriteriafortransferabilityusedhereregardstheextenttowhichtheseforcesyieldproductsthatareinaccordwiththeoreticallysupporteddecompositionmechanismsofNMT[ 58 59 ].NotethatthetrainingofourHamiltonianwasnottailoredtofavoranyparticularrearrangementmechanism,evidentbythefactthatthedimerdecompositiondiersmarkedlyfromthetrimermechanism.AlongtheadiabaticsurfaceAM1doesnotpredictanybondbreakingtooccur.Therefore,theaccuracyofmodelingbondbreakingandformingprocessesinbulksimulationsusingAM1isquestionable,atbest. WehaveshownaparameterizationoftheNDDOHamiltonianthatreproducesproperdissociationofNMTintoradicals,bothfortheenergyandforces.ThisparametersetreproducesthereferenceCCSD/TZPgeometriesandforcesforunrestrictedCNbondbreaking.ThisNDDOparameterizationwouldappeartopotentiallyoerimprovedaccuracyinMDsimulationsforclustersofhighenergymaterials.Wealsoshowthatconcertedreactionsforboththedimer(bimolecularrearrangement)andtrimer(trimolecularrearrangement)arenotgeneratedwithanAM1parameterization,yetarerequiredfortheproperdescriptionofNMTclusters.ThisnewtransferHamiltonianwillallowfor 43 PAGE 44 67 68 ]becauseofitswideuseandpriorsuccesses.However,suchmethodsusingstandardparameterslikeAM1[ 27 ],MNDO[ 30 { 32 ]andPM3[ 28 29 ]arealsoknowntohavesignicantfailingscomparedtoabinitiomethods.OurTH,ontheotherhand,benetsfromthefactthatweonlyrequirespecicparameters[ 45 ]forthesystemsofinterest,forexample,SiO2clustersandH2Oinourhydrolyticweakeningstudies[ 69 ].Theseparameterspermitallappropriateclusteringandbondbreakingforallpossiblepathsforallthecomponents,toallowMDtogowherethephysicsleads.Furthermore,unliketraditionalsemiempiricalmethods,wedonotuseaminimalbasisset,preferringtousepolarizationfunctionsonallatoms,andinsomecases,diusefunctions.TheformoftheHamiltonianemphasizestwocenterinteractions(butisnotlimitedtonearestneighborsasinTB),whichallowsourparameterstosaturateforrelativelysmallclusters.WecheckthisbydoingabinitiocalculationsonlargerclustersinvolvingthesameunitstodocumentthattheTHiseectivelyunchangedwithrespecttosystemsize.Atthatpoint,wehaveaTHthatwecanexpecttodescribeextendedsystemscomposedofSiO2,anditssolutionsproperly,includingallmulticentereectsontheenergiesandforces.OnlytheHamiltonianisshortrange.WealsocheckourTHbyapplyingittorelatedsystems.WewillshowthatourTHdeterminedfromforcesforbondruptureindisiloxane,pyrosilicicacid,anddisilanealsoprovidecorrectstructuresforseveralsmall 44 PAGE 45 NowwewillshowthattheTHCCSDparametersforanNDDOtypeHamiltoniantrainedonasmallsetofreferencemoleculesaretransferableintherstsense.Wechoseourtrainingsettobereasonablerepresentativesofthelocalbehaviorofextendedsystems:disiloxane(H3SiOSiH3),pyrosilicicacid((OH)3SiOSi(OH)3)anddisilane(H3SiSiH3).TheinitialparametersusedwerefromMNDO/d[ 70 ]becauseofavailability,overallperformanceonthetrainingset,andinclusionofdorbitalsonSi. ThetwasbasedonforcesfromtheCCSDapproximationusingDZPbasissets.PrimarilytheTHCCSDtwastotheforcesalongtheintrinsicreactioncoordinateforSiSibondbreakingindisilaneandSiObondbreakinginbothdisiloxaneandpyrosilicicacid.FurtherrenementwasachievedbyttingtotheforcesfortheSiOSibendsofdisiloxaneandpyrosilicicacid.ThemostreliableerrorfunctiontouseinthetypeofttingseemstobethesquareofthedierencebetweenthereferenceCCSDandTHCCSDforces.Thettingprocedureitself,thoughnonlinear,involvediterativelyperformingalinearminimizationofeachparameter,startingwiththosethatmostreducedtheerror. Toverifythetransferabilityoftheparametersoneshouldevaluatepropertiesformoleculesoutsidethereferenceset.VericationforbondbreakingtestswereperformedonSi(OH)4,Si(SiH3)(OH)3,Si(SiH3)2(OH)2,Si(SiH3)3(OH),Si(SiH3)4,H2Si=SiH2, 45 PAGE 46 25 and 26 showtheresultsrelativetoDFT(withtheB3LYPapproximateexchangecorrelationfunctional)values. ItisapparentthattheTHCCSDparametersoutperformmanyofthecommonlyusedNDDOmethodsforavarietyofdierentbondingenvironments.Furthermore,althoughtheparametersweredesignedtoreproduceSiOSi,SiO,andSiSiforces,theytransferwelltoequilibriumOSiO(Figure 27 ). 71 ].UnderstandingsilicaandsiliconbringsusclosertounderstandingdefectsinMOS.Inthem,phenomenasuchascurrentleakage,electricalbreakdownovertime,andthethinningtrendsofthedielectriclayeringateoxidesarestronglyinuencedbydefects[ 72 73 ]. Italsoshouldbenotedthatwefocusontrainingmodelsthatperformwellprimarilyforlocalizedunits(i.e.,clusters)andsecondarilyforbulksystemswithperiodicboundaryconditions.Thisisbecauseitismuchhardertobuildlocalphenomenaintoamodelusedprimarilyforextendedsystemsthanitistoconstructamodelbasedonlocalchemistrythatisapplicableforsmallmolecules,clusters,andbulk.Forexample,periodiccalculationswithplanewavesrequireexceedinglylargenumbersofplanewavestodescribeoxygendefectsinsilica[ 74 ]. 46 PAGE 47 75 ]haveshownexperimentally,basedonthetwofoldcompressionofliquidSiO2,thatthereisessentiallynochangeinnearestneighbordistancesthoughtherearesignicantchangesinmediumrangestructure.Thismediumrangestructureisbestdescribedbyringandclusterstatisticsandclustergeometry[ 76 ],whichwouldbenotbeincorporatedinanytwobodypotential. Toavoidthelimitationsofempiricaltwobodypotentials,onemayconsiderNDDOornonorthogonalTBmethodsdiscussedabove,whichincludesomemanybodyeectsduetoinclusionoftheoverlapmatrix[ 77 ].EspeciallypromisingistheSCCDFTBapproach[ 7 ]whichincludesthelargestcontributionstolongrangeCoulombinteractionswhilehandlingchargeinaselfconsistentway.Thesemethodsoerapracticalbalanceofchemicalaccuracyandcomputationaleciency,asevidencedfromtheirrecentincreaseinpopularity. Figure 28 showsthedeviationsoftheTHCCSDparameterization,AM1,andMNDO/dNDDOcalculatedunitcelldensitiesfromreferenceDFTdensities.ThereferenceDFTcalculationsforallsilicapolymorphsarefromtheVASPcodeusingthegradientcorrectedPW91GGAexchangecorrelationfunctionalwith3Dperiodicity,anultrasoftpseudopotential,andaplanewavebasisset(395eVcuto).A2x2x2kpointgridisemployed,whichconvergestotalenergiestowithin0.05eV/cell,andtheconvergenceofselfconsistentstepsispreciseto104eV/cell.Theoptimizationofcellparametersandatomicpositionsispreciseto<103eVA1. ItisclearthattheTHCCSDparameterstrainedonsmallsilicaclustersaretransferabletoseveralsilicapolymorphs.Itisalsointerestingtonotethatnotonlyis 47 PAGE 48 TheapplicationofthetransferHamiltonianmethodologyhasbeenpresentedfornitrogencontainingenergeticmaterialsandsiliconcontainingcompounds.InbothcaseswehavedemonstratedsignicantimprovementovertheunderlyingsemiempiricalmethodbyreparameterizingwithatrainingsetofhighqualityCCforcesforrepresentativemolecules.ItisapparentthatforfocusedMDstudiesthatusesemiempiricalmodelHamiltoniansthereisanadvantagetocreatingaspecialsetofparametersthatistailoredtothechemicalsystembeingstudied.However,thetransferabilityofthemodelisaectedbywhichmoleculesareselectedforthetrainingsetandthePESreferencepointsused.Thesemodelswilloftenbreakdownwhenappliedtosystemstowhichtheyhavenotbeenparameterized.Oneroutetoachievethisleveloftransferabilityistodesignasimpliedmethodthatdoesnotrequireparameterization,onlycontrolledapproximations,whichleadstotheAATapproachtodeningasuperiorTh. 48 PAGE 49 EquilibriumgeometryforNMT,inAanddegrees MethodRCNRCHRNOAHCNAONC 49 PAGE 50 TheNMTforcecurveforCNrupture 50 PAGE 51 TheNMTenergycurveforCNrupture 51 PAGE 52 TheNMTdimerenergiesrelativeminimumintheintermolecularhydrogenbondingdistance 52 PAGE 53 B C D Figure24. SnapshotsofthegeometryoptimizationsforNMTdimerandtrimerperformedwiththeTHCCSDandAM1.A)NMTdimertreatedwithTHCCSD.B)NMTdimertreatedwithAM1.C)NMTtrimertreatedwithTHCCSD.D)NMTtrimertreatedwithAM1 53 PAGE 54 AveragedeviationofSiOstretch 54 PAGE 55 AveragedeviationofSiSistretch 55 PAGE 56 AveragedeviationofOSiOangle 56 PAGE 57 DensitydeviationfromPW91DFT 57 PAGE 58 Semiempiricalmethoddevelopmenthasalonghistorydatingbacktothe1930'sstartingwithHuckeltheoryfordescribingevenalternanthydrocarbons.Inthecurrentstudywehavebeenguidedbythedevelopmentsoverthelastseventyyearsinregardtobothwhatshouldandshouldnotbeapproximatedinourmodel.WefocushereontheNDDOandSCCDFTBmodels.Bothmodelshavebeenwidelyusedwithreasonablesuccess.Therelativenumericalperformanceofthemostrecentversionsofeachmodelhasbeenreportedrecently[ 2 ],andextensivecomparisonsoftheformalsimilaritiesofeachmodelhavebeenmade[ 78 79 ].ItisimportanttounderstandthecontextinwhichthesemodelsweredesignedtofullyappreciatethedierenceinthedesignphilosophyofthecurrentAATmodel.NDDObasedmethodswereoriginallyconstructedtoreproducetheexperimentalheatsofformationforsmallorganicmolecules.Sincetheheatofformationisdenedonlyformoleculesintheirequilibriumstructures,itisnotsurprisingthatNDDOmodelsareparameterizedspecicallytoreproduceequilibriumproperties.FornonequilibriumstructurestheaccuracyofNDDOmethodsvariesgreatly,aswewillshow.Similarly,SCCDFTBhasbeenparameterizedtoreproduceequilibriumstructures. ThegoaloftheAATapproachistoachieveuniformaccuracyofelectronicproperties,structures,andenergeticsovertheentirePEStotheextentpossibleforaminimalbasisseteectiveoneparticletheory.UniformaccuracyisrequiredtoachievemeaningfulresultswhenperformingMDsimulations,inwhichtransitionstatesandothernonequilibriumstructuresplayacriticalrole. 3.1.1DoesNDDOApproximateHF? 58 PAGE 59 2(j)(3{1) (j)ABCD1 2 (j)ADBC(3{2) where,,,andareoncentersA,B,C,andDrespectivelyandtheterms (j)aresubjecttothemultipolemultipoleexpansionapproximation.ItistruethatatypicalNDDOtypetwoelectronintegralencounteredisoftenlargerinmagnitudethanacorrespondingnonNDDOtypetwocenter,threecenter,orfourcentertwoelectronintegral.Formoleculeswithmorethatthreeheavyatomsinavalenceonlyminimalbasissetrepresentation,therearesignicantlymoretermsinvolvingthreeandfourcentertwoelectronintegrals.ThecombinedeectofthislargenumberofmulticentertermshasasignicanteectonthedensitydependentFockmatrixcontributions.Weconsiderasimpleexample.ShowninFigure 31 aretheindividualcontributionstotheNMTdensitydependentFockmatrixelementofthecarbon2stypeorbitalandnitrogen2stypeorbitalinaminimalbasissetforthecombinedthreeandfourcenter,twocenter,abinitioNDDOtype,andfullmulticentercontributions.Thoughonlyasinglematrixelementisplottedoverthisreactioncoordinate,thetrendisthesameforothermatrixelementsandothermolecules.Thethreeandfourcentercontributionsarenearlyequalandoppositetothetwocentercontributions.Thiscancellationisauniversaltrendandholdsatallenergyscales.And,whileitcouldbesaidthattheNDDOtypeclassoftwoelectronintegralsapproximatesthefullsetoftwocentertwoelectronintegrals,itisclearthattheNDDOtypeintegralsalonebearlittleresemblancetothefulltwoelectronintegralcontribution.Instead,itisthedelicatebalanceamongallmulticentertwoelectronintegralsthatiscriticaltoobtainingtheproperFockmatrixcontribution. IthasbeenknownsincestudieswererstperformedontheNDDOapproximationthataparameterizedsetofNDDOintegralsperformbetterthattheabinitioNDDOintegrals.However,whathasbeenlessclearisthelackofcorrespondencebetweentheNDDOapproximationandthequantityitissupposedlyapproximating.Insteaditseems 59 PAGE 60 80 ]showthattheirPDDG/PM3reparameterizationofPM3yieldsheatsofformationwithameanabsoluteerrorof3.2kcal/molforatestsetof622molecules.However,thereislittleevidencetosuggestthatthesameperformanceshouldbeexpectedforchemicalstructuresthataresignicantlydistortedfromequilibriumstructures.Evenforsimpleactivationbarriers,theerroristypically1015kcal/mol[ 80 ].ForthedeterminationofrateconstantsandinmoleculardynamicsstudiesthePESshouldbedescribedwithin12kcal/mol.Forexample,sinceat300Kanerrorof1kcal/molintheactivationbarrierwillresultinanorderofmagnitudechangeintherateconstant,clearlya1015kcal/molerrorwillnotyieldmeaningfulrateconstants.ThesituationisespeciallybadincasesinwhichdirectbondssionistreatedwithNDDObasedmodelHamiltonians.Insuchcasesthereisanarticialrepulsiveregionnear1.2Reqto1.6Req,asisshowninFigure 32 forNMTUHFdirectbondssion.Generally,incaseswhereitisknownthatthereisbarrierlessdissociation,onewouldrecognizesuchfeaturesasanerrorinthemethodandsubsequentresultswouldbequestioned.Formorecomplicatedchemicalprocesses,forexampletransitionstatesor 60 PAGE 61 ThecorecorerepulsionterminNDDOHamiltoniansinvolvesparametersthataredesignedforenergeticsatequilibrium.Forsuchaspecicsetofpropertiesthedierenceoftwolargenumbers,theattractive(negative)electronicenergyandtherepulsive(positive)nuclearnuclearrepulsionenergy,iseasiertomodelthaneachseparately.Inouropinion,combiningsuchunrelatedtermsistheoriginofthelongstandingcontentionthatelectronicpropertiesandtotalenergypropertiescannotbedescribedwithinthesamesetofparameters.Also,oncethecorecoretermsareparameterized,thecorrespondencebetweenrigoroustheoryandsemiempiricaltheoryislost.Tosolvethisproblem,weneedtosearchforabetterformforourmodelHamiltonianthatdoesnotrequirethatVNNbecombinedwithtermsthatarepurelyelectronic. 61 PAGE 62 1{27 Now,ifweexpectthattheNDDOmodelHamiltonianwillgeneratethecorrectelectronicdensity,whichisarequirementifwewantanyoftherstorderelectronicpropertiesofoursystem,thenwecanrelatetheNDDOandKSDFTenergyexpressions.Makingthetemporaryassumptionthat^fNDDO^fKS,andconsequentlyNDDOKS,thenthetotalenergydierenceis 2Zvxc(r)(r)drXA;B>AZAZB NowitisclearthatiftheNDDOoneparticleHamiltonianweretoyieldthecorrectelectronicdensitythenthetwocentercorecoreterm(ENDDOAB)mustaccountfornotonlythesimplenuclearnuclearrepulsionterm(ZAZB 2Rvxc(r)(r)).Themotivationforapproximatingthislasttermisobvious:itdrasticallysimpliesandavoidsacomplicatedandexpensivestepinthedeterminationofthetotalenergy.However,thereisnoreasontoexpectapriorithatthisapproximationwouldbeaccurate.ForourimprovedmodelHamiltonianwemakeaspecialeorttoensurethatallapproximationsarecontrollableandsystematicallyveriable. 3.2.1SelfConsistentChargeDensityFunctionalTightBindinginanAOFramework 3{4 through 3{8 aretheworkingequationsforSCCDFTBwhichcomedirectlyfromElstner'soriginalpaper[ 7 ]: qA=qAq0A(3{4) 2occXiniX2AatomsX(ciciS+ciciS)(3{5) 62 PAGE 63 ~HSCC=hj^H0ji+1 2SatomsX(A+B)q=~H0+~H1 2atomsXA;BABqAqB+XA6=BErep(RAB)(3{8) andthecorrespondinggradientsare wehaveused2Aand2B. TofacilitatecomparisonoftheSCCDFTBmodelHamiltoniantoothermethodsweforceitintoastructurethatiscomprisedofadensityindependentpartHandadensitydependentpartG.UsingMullikenpopulationanalysistodeterminethenetchargeassociatedwithatomAyields qA=q0A+X2A;PS(3{10) wheretheoccupationnumbersni=1andtheMOcoecientsarereal. 2SatomsX(A+B)q0+1 2SatomsXX2;(A+B)PS=(HSCC)+XP(G]SCC) wehaveintroducedthenotationthat (HSCC)=~H01 2SatomsX(A+B)q0(3{12) 63 PAGE 64 2SS(A+B)(3{13) NotethatunlikeHForNDDO,inSCCDFTBitappearsthatitisnotgenerallytruethatG=G,as2andisunrestricted.IfSCCDFTBweretohavethisproperty,thatwouldbeusefulbothforthecalculationofforcesandforcomparisonamongalleectiveoneparticlemethodsweconsider.Toenforcethissimilaritywemakethefollowingobservation: 2SS(AC+BC)]=XP[1 4SS(AC+BC+AD+BD)](3{14) where,,,andareonA,B,C,andDrespectively,andourworkingequationforthedensitydependentcontributioninSCCDFTBis (GSCC)=1 4SS(AC+BC+AD+BD)(3{15) Equation 3{15 bearsasimilaritytotheMullikenapproximationtoCoulombrepulsion,wewilladdressthisinfurtherdetaillater. Nowwerevisitthetotalenergyexpression.UsingournewdenitionofthedensitydependentandindependenttermsintheSCCDFTBFocklikematrix. 2XP[2HSCC+XP(GSCC)]+XA6=B(1 2ABq0Aq0B+ESCCrep(RAB))(3{16) ThisenergyexpressioncloselyresemblesotheroneparticleHamiltonians,yetisequaltoEquation 3{8 .Similarly,aftersignicantmanipulationthegradientexpressioncanbeshowntobe 2XP[2@HSCC @XA6=B(1 2ABq0Aq0B+ESCCrep(RAB)) (3{17) whichagainisequaltotheoriginalexpression,Equation 3{9 ,butallowsamoresystematiccomparisonwithothermethods. 64 PAGE 65 3{13 ,whichhasbeenshowntobeequivalenttoEquation 3{15 ,thoughthelatterhasaformthatsomemayrecognizeastheMullikenexpansionfortwoelectronintegrals, (j)1 4SS(AC+BC+AD+BD)(3{18) TheoriginoftheMullikenexpansionandthefurtherextensionofthebasicapproximationtotheRudenbergapproximationwillbeexploredlater.Fornowitissucientthatwediscussthenumericalpropertiesofthisapproximationandwhatadvantagesandlimitationsitentails. ItisperhapseasiesttoconsiderthenatureoftheMullikenapproximationbywayofasimpleexample.ForanNDDOtypetwocentertwoelectronintegralwithnormalizedorbitalsitfollowsdirectlyfromEquation 3{18 Therehavebeenmanystudiesonthevariousformsfor,mostnotablybyKlopman[ 39 ],Ohno[ 40 ],andNishimotoMataga[ 41 ].TheparticularformwithinSCCDFTBmodelistheCoulombrepulsionbetweentwo1sSlatertypeorbitalswithanaddedconditionthatastheseparationbetweenthetwocentersapproacheszeroAB!UAwhereUAistheHubbardparameter,whichisrelatedtochemicalhardness,whichinturnisrelatedtothedierencebetweentherstIPandEAofthatatom.TheparticularformoftheABisonlyassignicantastheunderlyingapproximationinwhichitisused.Inthis 65 PAGE 66 IfweapplytheMullikenapproximationinstandardHFandusetheabinitio(j)twocenterintegrals,thebehaviorisquiteunphysical.ShowninFigure 34 areerrorsfromthefullHFwhentheMullikenandRudenbergapproximationsareappliedtothetwoclassesofthreecenterterms,(AAjBC)and(ABjAC),andthefourcentertermsinNMToverthecarbonnitrogenreactioncoordinate.TheerrorintroducedbytheMullikenapproximationisdrastic,particularlyforthe(AAjBC)threecenterterm,forwhichtheerrorisover1000kcal/mol.GiventhattheonlyexplicitdensitydependenceoftheFocklikematrixinSCCDFTBisintroducedbytheMullikenterm,andyettherelativeenergeticsofSCCDFTBinpracticedonotintroducesucherrors,thensomedensityindependentpartofthemodelmustbeaccountingforthedierence.WewillavoidsuchuncontrolledapproximationsintheAATapproach. 2Zvxc(r)(r)drXA;B>AZAZB ThenumericalbehaviorofthiserrorisshowninFigure 33 .Clearly,thisrepulsivepotentialiscarryingtheburdenofaportionoftheelectronicenergy.Thisisnotbyaccident;thevalueofafunction(inthiscasethetotalenergy)thatresultsfromthe 66 PAGE 67 67 PAGE 68 Thisdependenceallowsustouseaselfconsistenteld(SCF)procedure:determinetheFocklikematrixfromaguessdensityanddetermineanewsetofcoecientsandthecorrespondingdensity,repeatuntilthedensityoutputiswithinsomethresholdoftheinputdensity.Alternatively,FPS=PSFandhajheffjii UponconvergenceoftheSCFequationswehaveasetofAOcoecientsthatdenethesetofMOs areferencesingledeterminantwavefunctionisthengeneratedasanantisymmetrizedproductoftheMOs, 0=Af1(1)2(2):::n(n)g(3{23) Itisusefulatthispointtodenetherstorderelectronicenergygivenazerothorderwavefunction,j0i.TheexactelectronicHamiltonianis ^H=Xih(i)+Xi;j>i1 wherethecoreHamiltonianis 2r21XAZA andinitsmatrixrepresentationthecoreHamiltonianis 68 PAGE 69 2Xi;jhijjjiji NotethattherstorderelectronicenergyisexactlythesameastheHFenergyexpressionandfurthermorehasnoexplicitdependenceontheparticularformoftheoneparticleoperator^f. However,thegoalisnottherstorderelectronicenergybuttheexactelectronicenergy.WecanconsidertheCCroute,whereweknowthatforasingledeterminanttheexactelectronicenergyis(withanychoiceoforbitals) TheCCenergyexpressionprovidesgreatinsight,buttheexplicitcalculationoftabijandtaiamplitudesrequiresmuchmorecomputationaleortthanweareabletoexpendforafastsemiempiricalmodel.AmorenaturalenvironmentforasemiempiricalmethodsisgivenbyKSDFT.Assumingthatwehavetheexactexchangecorrelationenergyfunctional,theexactKSDFTenergyis(forKSDFTorbitals) 2Xi;jhijjjii+Exc[]=Tr[PHcore]+1 2Xi;jhijjiji+Exc[] (3{29) ThereisanotherformoftheKSDFTenergyexpressionthatismoreusefulforthecurrentdevelopment,thoughitrequiresthatweuseaslightlymorespecicformforouroneparticleoperator.IfthetargetisaoneparticleHamiltonianthatreturnsthecorrectelectronicdensity,andhencethecorrectrstorderelectronicproperties,thentheidealoneparticleoperatorisindeedtheKSDFToneparticleoperator.And,ifonlythetotaloperatorandnotitsindividualcomponents,i.e.,1 2r2,ext(r),J(r),andxc(r),is 69 PAGE 70 3{29 wouldbesuitable.Byseparatingtermsthataredependentonthedensity(G)wecanwrite 2Tr[P(2HKS+GKS)]1 2Tr[PGKSxc]+Exc[](3{30) where (GKSxc)=hjvxc(r)ji(3{31) andF=H+G.Equation 3{30 isacentralthemeofourwork,asitgivestheunambiguousdenitionoftheexactelectronicenergy.TiedtoitaretheoneparticleKSequationsthatdenetheorbitals,thedensityandassociatedrstorderproperties.Inthisway,theproblemforthetotalenergyandgradientsiscompletelyspeciedby 2Tr[P(2H+G)]+Eresidual[P]+XA;B>AZAZB WewillbringHF,NDDO,DFTB,SCCDFTBandKSDFTintothisframework,tomakeitclearwhatthethreedistinctcomponentsare: Ifamethodistogivethecorrectelectronicdensity,andwecanknowtheexactKSDFTHKSandGKS,thenthemodelcomponentsforHandGabsolutelymust PAGE 71 2r21XAZA 2S(+)if6=(3{34) 2H0(3{35) 2H01 2SNX(A+B)q0(3{36) 2r2XAZA TheonlydensityindependentcontributiontotheidealFocklikematrix(theFocklikematrixwhichreturnsthecorrectoneparticlegroundstatedensity)istheexactcoreHamiltonian,Hcore.ThebiggestbottleneckforcomputingcontributionsforHcoreisthethreecenteroneelectronintegrals,whichscaleasn2orbitalsNatoms.Still,theyrepresentasmallcomputationaleortforlargesystemswhencomparedtothen3orbitalsterms 71 PAGE 72 2(j))(3{38) (j)ABCD1 2 (j)ADBC)(3{39) 4SS(AC+BC+AD+BD)](3{41) ThedensitydependenttermoftheidealFocklikematrixisGKS.Thischoicefollowsdirectlyfromthechosenformforthedensityindependentcomponentandtherequirementthatwereproducetheexactelectronicdensity.ThedensitydependenceoftheFocklikematrixaccountsforthecomplexinteractionsamongelectrons.ThistermisdominatedbyCoulombicrepulsion.Thesecondmostimportanttermisgenerallycalledexchange.Thethirdcomponentiscorrelation,andcanbeconsideredtobethehighlycomplexwayinwhichelectronsinteractwithoneanother,i.e.,theirmotionis 72 PAGE 73 2Tr[PGKSxc]+Exc[](3{47) 73 PAGE 74 3{30 ,thegradientis[ 81 ] d=2@S 2XP@G whichisbasedonPiCi(FiS)=0andorthonormality.Notethatthereisnoexplicitdependenceon@P 3{48 onlyholdsiftermsinGarelinearinthedensity ThederivativesofEresidual,whicharenotpresentinHF,aresimpleforthosemethodsthatonlydependonasumofatompairterms,suchasNDDO,DFTB,andSCCDFTB.InKSDFTadditionalnumericalintegrationisrequiredfortheperturbeddensity.PotentialapproachesforreducingthecomputationalcostassociatedwiththederivativesinKSDFTwillbediscussedinfurtherdetaillater. 74 PAGE 75 75 PAGE 76 MulticentercontributiontotheC2sN2smatrixelementofNMT. 76 PAGE 77 TheAM1articialrepulsivebumpforNMTdirectbondssion. 77 PAGE 78 TheSCCDFTBtotalenergybreakdown. 78 PAGE 79 B C Figure34. EnergiesfromRudenbergandMullikenapproximations.A)Threecenter(AABC).B)Threecenter(ABAC).C)Fourcenter(ABCD). 79 PAGE 80 Wehavespeciedthetargetframeworkofoursimpliedquantummechanicalmethod.Thetasknowistodemonstratesimplewellcontrolledapproximationsthatsignicantlyreducecomputationaleortwhilemaintainingadesiredlevelofaccuracy.Tothisend,aprimaryfocusisonthedensitydependentcomponentoftheFocklikematrix,asitisthemostcomputationallyintensivepartoftheSCFprocedure.Furthermore,theresidualelectronicenergycontributionisexaminedbecausewehaveshownthatitspropertreatmentisessentialforaSMtoachieveuniformaccuracyofelectronicstructureandenergetics. Ourstrategyistostartwithareferencetheoreticalmethodthatisreasonablywellbehaved,inthiscaseKSDFT.ThoughthereismuchdebateaboutthebestDFTenergyfunctional,wechooseahybridDFTfunctional,specicallythewidelyusedB3LYP,becauseitincludesaportionofexactexchangewhichaddsasubstantialdegreeofexibilitytooursimpliedmethod.Tosimplifythecomputationalprocedurewethenexplorevariousapproximationsandverifytheiraccuracy.Itmaybetoomuchtodemandasimpliedmethodthatissystematicallyimprovable(intheCCwavefunctiontheorysense),however,everyapproximationshouldbesystematicallyveriable.WewillshowthatbymakingacoupleofwellcontrolledapproximationsweareabletoreproducethetopographicalfeaturesofahybridDFTPESwithasubstantialreductionincomputationaleort. 82 ],butthiscannotbedonefortheexactexchange,leavingn4orbitalsdependence.Thoughmostquantumchemistryprogramsemployscreeningtechniquestoreducethenumberoftwoelectronintegralscalculated,theyremainasignicantcost,especiallywhenexact 80 PAGE 81 83 { 85 ],thelinearexactexchange[ 86 ],andquantumchemicaltreecode(QCTC)[ 87 ].ThesemethodsstillrelyontheexplicitcalculationofintegralsusingstandardtechniquesforCartesianGaussiantypefunctions.Forlargesystemsthemostnumerousclassofintegralsarethefourcenterintegrals.Furthermore,fourcenterintegralsarenotonlythelargestclassofintegrals,butalsothemostcomputationalinvolvedintegrals.AccordingtoGilletal.[ 88 ]afourcenterintegralistypicallyanorderofmagnitudemorecomputationallydemandingthanatwocenterintegral.Ourfocushereisonsystematicapproximationsthatcanbemadetotheseabinitiomulticenterintegralsthatwillreducethecomputationaleortassociatedwithcalculatingthesetermswhileretainingacceptableaccuracy. ThecornerstoneofourapproachisbasedonideasoriginallypresentedbyRudenbergin1951[ 5 ].TheunderlyingapproximationistoexpandanorbitaloncenterAintermsofasetofauxiliaryorbitalslocatedoncenterB: where IftheauxiliaryorbitalsoncenterBarecompleteandorthogonalthentheexpansionisexact.ByrepeatedapplicationofEquation 4{1 thenumberofcentersthatareinvolvedinthreeandfourcentertwoelectronintegralscanbereduced.ThepossiblecombinationsofthistypeofexpansionhasbeenexploredandenumeratedbyKochetal.[ 89 ],thoughwithoutcomparingthenumericalaccuracyofsuchapproximations.Wewillshowthe 81 PAGE 82 First,considertheeectoftheorbitalexpansiononatwocenterproductoftwoorbitalswithrespecttothesameelectroncoordinate, 21Xk=1[SAkBkB(1)B(1)+SkABA(1)kA(1)](4{3) Sinceageneralmulticentertwoelectronintegralisgivenas substitutingEquation 4{3 intoEquation 4{4 forboththeABandCDorbitalproductsyieldsthefollowingexpression,whichistheRudenbergtwoelectronintegralform, (ABjCD)=1 41Xk=11Xj=1[SAkBSCjD(kBBjjDD)+SAkBSjCD(kBBjCjC)+SkABSCjD(AkAjjDD)+SkABSjCD(AkAjCjC)](4{5) Inthelimitthatthesumsoverauxiliaryorbitalskandjarecompletethisexpansionisexact.Inpractice,theexpansionisapproximatedbylimitingtheorbitalsoverwhichweareexpandingtothoseincludedintheAObasisfunctions,thoughusingaseparateauxiliarysetoforbitalsmayalsobesuitable. TheRudenbergtwoelectronintegralexpression(Equation 4{5 ),ismoregeneralthanthewellknownMullikentwoelectronintegralapproximation.TheMullikenapproximationisequivalenttorestrictingthesuminEquation 4{1 tobeonlyoverorbitalsofthesametype, 2SAB[A(1)A(1)+B(1)B(1)](4{6) 82 PAGE 83 (ABjCD)=1 4SABSCD[(AAjCC)+(AAjDD)+(BBjCC)+(BBjDD)]:(4{7) NotethesimilarityofthisexpressiontothatwhichappearsinSCCDFTB,Equation 3{18 ,inwhich thoughinSCCDFTBtheorbitalsandarerestrictedtobestypefunctions. WeconsiderherealleightdierentapproximationsthatresultfromtheapplicationofEquation 4{1 tothreeandfourcentertwoelectronintegrals.Thetwoclassesofthreecentertermsareoftheform(AAjBC)and(ABjAC),whichwillbereferredtoastypesAandBrespectively.TofacilitateourdiscussionweintroducethenumberingconventioninTable 41 TheRudenbergandMullikenformulascanbeappliedtoeitherclassofthreecenterorthefourcenterintegrals(approximationsIVI)usingEquations 4{5 and 4{7 ,respectively.Forthepartialapproximationstofourcenterterms(VIIandVIII),partialmeansthatthecorrespondingapproximationhasbeenappliedonlyonce.Theexpansionforthesepartialapproximationsisthenintermsofthreecentertwoelectronintegrals(insteadofonlytwocenterterms).FortheMullikenpartial(VII)wehave (ABjCD)=1 4fSAB[(AAjCD)+(BBjCD)]+SCD[(ABjCC)+(ABjDD)]g:(4{9) andtheRudenbergpartial(VIII)is (ABjCD)=1 41Xk=1[SAkB(kBBjCD)+SkAB(AkAjCD)+SCkD(ABjkDD)+SkCD(ABjCkC)](4{10) 83 PAGE 84 Notethatthetotalelectronicenergyisthesumofalltheone,two,three,andfourcentercomponents.And,eachofthesecomponentsindividuallybearlittlesimilaritytothetotalelectronicenergy,againconsidereachcontributiontoanyparticularmatrixelementshowninFigure 31 ,insteadthetotalenergyistheresultofthedelicatecancellationamongtheseterms.Inaddition,sinceonlythenonparallelityerrorsareofimportanceinreconstructingthePES,theenergeticcontributionarisingfromeachindividualclassoftermsmaybeshiftedbyaconstantvalueindependentlyofoneanotherwithoutaectingthefeaturesofaPES.Thesesubtleties,compoundedbythe3N6degreesoffreedomofthePES,makeascribingsignicancetoanysingleclassoftermsdicult.BecausewewanttoassesstheaccuracyofaparticularapproximationasitappliestoasingleclassoftermswewillinsteadfocusontheFockmatrixelements,whichcontainalltheinformationneededforthezerothorderelectronicenergy,toovercomesomeofthesediculties. GivenanysetofpointsonthePES(forexample,adissociationcurve)andtheFockmatrixforeachpoint,eachclassofmulticentertwoelectronintegralswillcontributeaparticularpercentagetothefulldensitydependentFockmatrix.ThebestapproximationforaparticularclassofintegralsisthengivenbyhowwellitreproducestheaveragepercentageoftheFockmatrixforthecorrespondingabinitiosetofintegrals.Consider 84 PAGE 85 41 arethescatterplotsforeachmulticentercontributiontotheCoulombandexchangecontributionsoftheFockmatrixfromHFtheoryinaminimalbasis.TheCNreactioncoordinaterangesfrom1.0Ato5.0Ainstepsof0.2A.Theslopesofthecorrespondingbesttlinesoftheone,two,threeandfourcentercontributionsare0.0063,0.5199,0.4522,and0.0216respectively,notethattheseslopesrepresentthefractionofthefullcontribution,andthattheirsumisone. Table 42 showsthestatisticalaveragesofthedensitydependentcontributionstotheFockmatrixforseveralmoleculesalongthespeciedreactioncoordinate(Rxn)attheHFleveloftheoryinaminimalbasisset.Thedissociationinallcasesisfrom1.0Ato5.0Ainstepsof0.2A.Onaverage,thetwocentertwoelectronintegralcomponentaccountsforroughlyhalfofthefullmatrixelementsandthethreecentertypeAintegralsrecoveranotherthird.Again,thecorrespondingenergeticcontributionswillnotexhibitthesameratios.Thepurposehereistovalidateindividualapproximationsmadetovariousmulticentertwoelectronintegrals.Tothisendweconsidertheaveragecontribution:the 85 PAGE 86 TheaprioribestapproximationtothetwotypesofthreecentercontributionsisthatofRudenberg.ShowninFigure 34 aretheerrorsintroducedbymakingapproximationsIVI,forNMTdirectbondssionalongthecarbonnitrogenreactioncoordinate.FortypeAandBthreecenterintegrals(3CAand3CB)theaveragematrixcontributionsare35.67%and4.98%respectively.TheRudenbergapproximationsIIandIVdierbyabouthalfapercentat35.12%and5.50%.ThoughthepercenterrorforbothapproximationsIIandIVareroughlythesame,thecorrespondingerrorintroducedenergeticallyistwoordersofmagnitudelargerforthe(AAjBC)threecentertermsthanforthe(ABjAC)threecenterterms.Incontrast,theRudenbergapproximationappliedtofourcenterterms(VI)is2.47%whichdiersfromthefourcenterabinitio(4C)valueof2.39%byonly0.08%.TheerrorintroducedbythefullRudenbergfourcenterapproximation(VI)issurprisinglysmall.OnewouldexpectthatthepartialRudenbergapproximationwouldhavethebestagreementwiththeexactvalue,instead,atleastempirically,thefullRudenbergapproximationperformsbest,andhasthesmallesterrorrelativetotheaveragecontributiontotheFocklikematrix.Indeed,theeectthis0.08%dierencehasontheenergeticsofadissociationcurveiswithinanacceptablerangeandreproducestheenergeticswellforavarietyofdierentdissociationcurves.WewillfocusonthemoleculesinTable 42 thatinvolvethedissociationofaCNbond.PlottedinFigure 42 arethedissociationcurvesforNMT,nitroethane(NET),COHNO2,andCH3NH2attheB3LYPleveloftheoryinaminimalbasisset.ThereisacleartrendforthefourcenterRudenbergapproximationtooverestimatetheaveragecontributiontotheFocklikematrix.Theenergetics,ontheotherhand,donotshowthisclearcuttrendsincetheenergyisloweredforNET,raisedforCH3NH2,andessentiallyunchangedforNMT.However,eventhoughthePESarenotpreciselyreproducedtheystillexhibitallthebasicfeaturesoftheab PAGE 87 WehavedemonstratedthattheMullikenandRudenbergapproximationsdonotapproximatethethreecenteroneortwoelectronintegralswell.Theseclassesofintegralsposeachallengeastheyaretheonlyremainingtermsthatrequireexplicitintegrationsinceallonecentertermscanbetabulatedandtwocentertermscanbeinterpolatedbycubicsplines,aswillbedemonstratedlater.However,basedontheGaussianproductrulethesetermscan,withoutapproximation,bereducedtotwocenterterms.Inprinciplethencubicsplineinterpolationcouldbeusedtostoretheintermediatesandcompletelyavoidtheneedtodoanyintegralsexplicitly.TheGaussianproductruleis +(RARB)2e(+)(rRA+RB +)2(4{11) WecangenerateanauxiliarysetoforbitalsforeachGaussiantypeorbitalpairthathavetheorbitalexponent+.Also,becauseproductsofl=1(ptype)orbitalsareinvolvedthisauxiliarysetmustalsonowspanl=2(dtype).Theprefactorisnotincludedinthestoredorbitals,insteaditiscomputedontheytoretainthetwocentercharacterofthetermsinwhichitwillbeused. ThereisasignicantdicultywithusingtheGaussianproductinasimpleway.BecauseweworkincontractedGaussiantypeorbitalsthenumberoforbitalpairsincreasesasthenumberofcontractedfunctionssquared.ThecentersoftheseresultingproductGaussiansarenotnecessarilythesame,RA+RB 87 PAGE 88 IfweconsiderthebasissetSTOnG,asumofcontractedGaussiansapproximateaSlatertypeorbital.ThereisnoSlaterproductrulethatreducesaproductoftwoSlaterstoasingleSlater.Insteadwehave ConsiderasimpleexampleoftheproductoftwoGaussianfunctionsversustwoSlaterfunctionsseparatedbyoneBohr,whichisrepresentedbythebluecurveinFigure 43 .Thechallengeisnowclearer,theunderlyingbasisfunctionsweusearettoSlaterfunctions,andtheproductoftwoSlaterfunctionsdoesnotresembleasimpleGaussian,orevenasumofcontractedGaussians.TheSlaterproductisdenedexactlybytheintersectionoftwofunctions: where,for PAGE 89 Wehavedevelopedthestrategyandframeworkforeventuallyincludingthethreecenterterms.FornowwearesatisedwithreducingthecomputationalcostthroughonlythefourcenterRudenbergapproximationbecausethefourcentertermsaccountformostofthecomputationaleortassociatedwithconstructionoftheFockmatrixforlargesystems. Equation 4{18 isexact.Itonlyrequiresanarbitrarypartitionofthedensityintoatomicparts,thoughconvergenceoftheseriesisaectedbythechoiceofpartition. OurtargetistoincludecontributionsfromtheexchangecorrelationpotentialtotheKSFocklikematrixinawaythatonlyrequiresthestorageofalimitedsetoftwocenter 89 PAGE 90 4{18 thatcanberetainedarev1centerxcandv2centerxc.Furthermore,suchcontributionsoughttobeindependentofchemicalenvironmentandthuscompletelytransferable.Therefore,weinsistondensitiesthataretransferable.Thenaturalchoiceissphericallysymmetricatomcentereddensitiesofneutralatoms,A0A,recognizingwecanhaveandspincomponents.Basedontheseconsiderationswehaveanapproximateexchangecorrelationpotential, ~vxc[XAA](r)XAvxc[0A](r)+XA>B[vxc[0A+0B](r)vxc[0A](r)vxc[0B](r)](4{19) Wedonotwanttomapandstoretherealspaceexchangecorrelationpotential,becauseoftheobvioushighdegreeofcomplexity.Instead,allthatiseverneededarethecontributionsfromthepotentialprojectedinthebasisset.TheFocklikematrixcontributionfromthisapproximateexchangecorrelationpotentialis (eGxc)AB[]=hAj~vxc[]jBi(4{20) Restrictingthistocontributionsthatonlydependontwocentersyields,forthediagonalblock(A=B) (eGxc)AA=hAj~vxc[XA0A]jAi=hAj~vxc[0A]jAi+1 2XB6=AhAj~vxc[0A+0B]~vxc[0A]jAi(4{21) andfortheodiagonalblock(A6=B) (eGxc)AB=hAj~vxc[XA0A]jBi=hAj~vxc[0A+0B]jBi(4{22) 90 PAGE 91 2XC6=A;B1Xk[SAkBhkBj~vxc[0C]jBi+SkABhAj~vxc[0C]jkAi](4{23) Theleadingordertermsare (eG0xc)AB=8><>:hAj~vxc[0A]jAiifA=BhAj~vxc[0A+0B]jBiifA6=B(4{24) andaslightimprovementcouldbeachievedwiththeadditionof (eG1xc)AB=8><>:1 2PB6=AhAj~vxc[0A+0B]~vxc[0A]jAiifA=B1 2PC6=A;BP1k[SAkBhkBj~vxc[0C]jBi+SkABhAjvxc[0C]jkAi]ifA6=B(4{25) ShowninFigure 44 areresultsfortheGZERO(B3LYP)approximation:eG0xcisusedwiththeB3LYPenergyfunctionalinterpolatedoveralltheatompairsCNOH.TheerrorsintroducedarerelativelysmallandGZERO(B3LYP)appearstobegenerallytransferable. IntermsofcomputationalfeasibilityforGZEROwehaveintroducedasingleonecentermatrixperatomtypeandasingletwocentermatrixperatompair.ForC,N,OandHatomtypesthenumberofonecentermatricesistrivialanddoesnotpresentachallenge,astheyarethedimensionofthenumberoforbitalsperatomsquared.ThetwocentermatricesarestoredinthelocalcoordinateframeworkoftheatompairandareroughlythedimensionofthenumberoforbitalsonatomAtimesthenumberoforbitalsonatomB.Theactualdimensioncanbereducedbecausemostofthematrixelementsarezeroduetoorthogonality.Tointerpolatethevaluesofthesetwocentermatrixelementscubicsplinesareused.Cubicsplinesrequireveparameterspernumberinterpolated.Weusedadaptivegrids(nonuniformlyspaced)thataredenedforeverypairtermfrom0.5 91 PAGE 92 InadditiontoapproximationsforintegralswealsodemonstratedapproximationsthatcanbeappliedtotheexchangecorrelationpotentialcontributionoftheFocklikematrix,andtherebyavoidtheneedtoperformnumericalintegrationforeverySCFiteration.Extensionstoincludethreecentercontributionsfromthetelescopingseriesapproximationtothepotentialwerealsopresented. Thelastcomponentwewillconsiderinvolvestheresidualtotalenergy.Thoughitisworthwhiletoavoidnumericalintegrationtogeneratethepotentialandsubsequentmatrixcontribution,analonetimeevaluationofExc[]uponSCFconvergenceshouldnotbeabottleneck.Furthermore,asdemonstratedbyBartlettandOliphant[ 90 ]theexchangecorrelationenergyfunctionalisrelativelyinsensitivetotheinputdensity.Indeed,aswehaveseenalreadyinFigures 42 and 44 theB3LYPenergyevaluationwasusedandintroducednosignicanterror. 92 PAGE 93 2r21+XAZA alloneelectrontermsarecalculatedexactlyastheyrepresentasmallportionofthetotalcomputationalcostforlargesystems.GiventhatweareusingahybridKSDFTfunctionalwehaveafractionofexactexchangeinthedensitydependentFocklikematrix,forB3LYP,=0:2. (GJK1;2;3center)=XP((j)(j))(4{27) (GJK4center)=XP(^(j)^(j))(4{28) where 4Xk=1Xj=1[SAkBSCjD(kBBjjDD)+SAkBSjCD(kBBjCjC)+SkABSCjD(AkAjjDD)+SkABSjCD(AkAjCjC)](4{29) (G0xc)=8><>:hAjvxc[0A]jAiifA=BhAjvxc[0A+0B]jBiifA6=B(4{30) andtheresidualelectronicenergyis 2Tr[PG0xc]+Exc[P](4{32) Finally,showninFigure 45 areUHFdissociationcurvesforAATModelZero(AATM0),whereboththeRudenbergandGZERO(B3LYP)approximationshavebeenappliedsimultaneously.ThismodelreproducestheunderlyingB3LYPdissociationcurveswellespeciallywhencomparedtoAM1. ThecorrespondinggradientexpressionforAATM0isgivenbyEquation 4{34 .TheonlydierencebetweenstandardKStypederivativesarisefromthemodiedfourcenter 93 PAGE 94 4Xk=1Xj=1[@SAkB d=2@S 2XP@f(G)1;2;3center+(G)4center]g 46 aretherelativetimingsfortheevaluationoftheabinitiofourcenterintegralsversustheRudenbergapproximation.Thesetimingsarepreliminaryandmoreworkcouldbedonetooptimizebothmethodsforevaluatingthisfourcentercontributions.However,itisclearthatthereisasignicantadvantagetousingtheRudenbergapproximationand 94 PAGE 95 47 Screeningofsmalltwoelectronintegralscanleadtolinearscalingforlargesystems.SuchscreeningmethodscanbeappliedtotheRudenbergapproximationveryecientlybecauseoftheexplicitdependenceonoverlap. InaneorttominimizecomputationaleortwithbuildingtheFocklikematrixwehavemadeextensiveuseofstoringtransferableonecenterandtwocenterterms.Inourcurrentimplementationallintegralsupthroughthreecentersarecalculatedexactly.Therearetworeasonsforexplicitlyevaluatingtheseintegrals:rst,theyareasmallcontributionwhencomparedtotheexplicitcalculationofthefourcenterterms;second,theapplicationoftheRudenbergorMullikenapproximationstothreecenterterms(decomposingthemintotwocenterterms)introducessignicanterrors.Ifanaccurate,simpleandtransferableapproximationforthethreecentertermswereavailablethenalltermsinthemodelcouldbetwocenter.Itisastraightforwardmattertointerpolatefunctionsoftwocentersusingcubicsplines.Wehaveimplementedsuchtechniquesfortheevaluationofthematrixelementsoftheexchangecorrelation:(G0xc)AB(RAB). Cubicsplinesareapplicablefornumerical2centerfunctions,wheretheexactfunctionalformisnotknown,andforcomplicatedanalytic2centerfunctions.Intherstcase,thebenetofsplinesistoprovideaninterpolationschemeforasetofempiricalreferencepoints.Inthesecondcasesplinesprovideacomputationallyecientframeworkfortheevaluationofotherwiseprohibitivelyexpensiveanalyticfunctions.Weareprimarilyinterestedinthelatter. Theabilitytostorealargenumberofparametersincomputermemory(RAM)isarequirementfortheecientimplementationofcubicsplines.Thememoryrequiredisdirectlydeterminedbythenumberofcubicsplineparameters.Apertinentexample 95 PAGE 96 48 amodestspeedupisachievedoverthealreadyfastNDDOFocklikematrixconstruction.Thesetimingsarefor16proteinsystemsupto20000orbitals.Usingcubicsplinesisupto15%fasterthanthemultipolemultipoleexpansion.Thereisanerrorintroduced,thoughitissmallandcontrollable.ForexampletheaverageerrorperatomisgiveninFigure 49 .Furthermore,whilegoingtohighermultipolesintheexpansionofintegralsisincreasinglyexpensive,cubicsplinescostthesameregardlessoftheunderlyingfunctionbeinginterpolated.However,thenumberofparametersinvolvedinatypicalNDDOmultipolemultipoleexpansionisrelativelysmallversusthenumberofparametersneededforcubicsplines.ThepremiumplacedonRAM(orcore)memorywhenNDDOmethodswererstdevelopedthirtyyearsagowouldhavemadecubicsplinesimpractical.Now,theamountofmemorytypicallyavailableonmostdesktopmachinesismorethansucienttomakecubicsplinesthemethodofchoicewhenevaluatingnearlyanytwocenterintegralinquantumchemistry. Ageneralcubicsplinehastheform 96 PAGE 97 91 ].Theonlydegreesofexibilitywhenusingcubicsplineinterpolationarethenumberofnodesusedandiftheyareevenlyspaced(discrete)orbasedonanadaptivegrid.Inourcase,forbothintegralsandmatrixelements,thefunctionsinallcasesgotozerointheseparatedlimitandatsmallerlengthscalesthefunctionrequiresadensegridtocompletelycapturethesubtlechangesthatoccurinthatregion.Adensenodespacingovertheentirelengthscaleisunnecessarybecausethefunctionschangeverylittleatlonglengthscales. Todetermineasuitablenodemappingfunctionwetriedseveraldierentfunctionalformsandconcludedthatasimpleexponentialfunctionprovidesagoodbalancebetweenhavingadensegridinthemorecriticalbondingregion(fromaround0.5Ato5A)andasparsegridatlongbondlengths(hundredsofA).Theactualnodemappingfunctionusedis Allofthefunctionswehavecurrentlyimplemented,thoseinvolvingtwocenterexchangecorrelationapproximations,arezerobeyond10A.Usingthenodemappingfunctionthenrequiresapproximately600referencepoints.Sincethereare5variablesperreferencepointweneed3000parametersthatwouldformalookuptableforasinglefunction.Ifdoubleprecisionnumbersareusedtostorethesplineparametersthestoragerequirementisabout0.023MBperfunction.FurtherdetailsofthecurrentimplementationoftheAATprocedureintheparallelACESIIIprogramsystemcanbefoundinAppendixA. 97 PAGE 98 92 ].TheG2setismadeupof148neutralmoleculeswithsinglet,doublet,andtripletspinmultiplicities,thoughwefocusonasubsetofthesethatcontainonlyC,N,O,andHbecausethecurrentimplementationofAATM0includesonlythesplineparametersgeneratedforthoseatomatompairs.Thereare71moleculesinthistruncatedsetrangingfrom2to14atoms. TheseparationofelectronicandnuclearcomponentsintheenergyexpressionofourAATapproachsetsitapartfromcommonlyusedSEmethods.Soitfollowsthatweexpectmoreconsistentresultsforpurelyelectronicproperties.TothisendwehaveconsideredtheverticalionizationpotentialsoftherelevantmoleculesintheG2set.Theverticalionizationpotentialispurelyelectronicasitdoesnotdependonthenuclearnuclearrepulsion.WehaveusedthesimpledenitionEIP=E(N)E(N1).Thisrequiresthetotalenergyofthemoleculeand,sinceallmoleculesintheG2setareneutral,thecorrespondingcationwiththeproperlymodiedmultiplicity.WehavecomparedAATM0againsttheunderlyingB3LYPminimalbasissetresultsandndaaverageRMSDfromB3LYPforourentiretestsetof0.057Hartree(1.6eV).TheanalogouscomparisontoAM1yieldsanRMSDfromB3LYPof0.331Hartree(9.0eV).ThisclearlydemonstratestheabilityofAATM0toreproducetheelectronicstructureoftheunderlyingB3LYP.TheerrorofAM1fromB3LYPissubstantial,thoughamorevalidcomparisonforAM1wouldbeagainstexperimentalIPs,theimportantissuehereistheabilityourAATM0toreproducethephysicsofthemethoditisapproximating. Wenowconsiderthedipolemoment.ThedipoleRMSDdeviationofAATM0andAM1fromminimalbasissetB3LYPis0.219a.u.(0.56Debye)and0.595a.u.(1.51Debye)respectively.TherewerenosignerrorsobservedfortheorientationofthedipolemomentforeitherAATM0orAM1.ThoughtheerrorintroducedbyAATM0issignicantlysmallerthanthatofAM1,itislargerthanwouldbeexpectedgiventheperformanceofotherpropertiesgeneratedbytheAATM0procedure. 98 PAGE 99 InadditiontotheaverageRMSDofthedensityofAATM0andfullHFweconsiderthestandarddeviationoftheRMSDofthedensityforeachmoleculeinthetestset.Thestandarddeviationsare0.042and0.057forAATM0andfullHFrespectively.ThissuggeststhatwhiletheoverallaverageRMSDissomewhathigherfromAATM0theerrorisslightlymoresystematic. TurningourattentiontoenergeticpropertiesweconsidertheatomizationenergiesoftherelevantportionoftheG2set.Wetakeastheatomicreferenceenergiestheneutral 99 PAGE 100 SinceonegoalofourapproachistoreproduceacomplexPESoflargesystemswehaveappliedtheAATM0approximationtodescribeasomewhatunphysicalreactionmechanismforC20inordertosubjecttheapproximationtoanextremecaseofbondbreaking.ThepseudoreactionpathcorrespondstosplittingtheC20moleculeinhalf,simultaneouslybreakingtenCCbonds.ShowninFigure 410 arethetotalenergiesofB3LYP,AATM0,AM1,andSCCDFTB.Thereactioncoordinateisthedistancebetweentheendcapsandtheequilibriumvalueis3.223A,therangeisfrom2.8Athrough6A.AsexpectedweseethatAM1exhibitsunphysicalbehavior,andsomeconvergenceproblemsbeyond4.5A.TheSCCDFTBactuallyperformsquitewellforC20,thisresultmaynotbetoosurprisingsincetherepulsiveenergyterminSCCDFTBisparameterizedfromB3LYPtotalenergiesforCCbondbreakingandtheenergyscaleisverysimilartothesumofCCbondenergies.TheAATM0overbindsthissystemsignicantlyandmayinfactbeconvergingtoadierentelectronicstate,astheeigenvaluespectrumdeviatessignicantlyfromtheB3LYPresult.Furtherworkisneededtodeterminethepreciseoriginofthisdiscrepancyandtocorrectit. ThenalcomparisonwemakeistotheRMSDofthegradientsfortheG2set.Wewillconsidertwosetsofdata,onefortheMP2equilibriumstructuresandtheotherfortheB3LYPminimalbasissetequilibriumstructures.ShowninFigure 411 aretheaverageRMSDoftheCartesiangradientsfromtheMP2referencestructurescalculated 100 PAGE 101 AbettercomparisonforhowwelltheAATM0approximatestheunderlyingB3LYPstructureisshowninFigure 412 .AttheB3LYPgeometriestheRMSDis0.05Hartree/Bohrand0.09Hartree/BohrforAM1andSCCDFTB,respectively.WeseetheAATM0hasanaverageRMSDof0.006Hartree/Bohrwithastandarddeviationofonly0.003Hartree/Bohr.ThisclearlyshowsthatAATM0isreproducingtheB3LYPgradientswithexcellentaccuracy. FastandaccurateapproximationstofourcentertwoelectronintegralsviaRudenberg(forCoulomborexactexchange)enablesubstantialsavingsforverylargesystems.SincetheRudenbergapproximationmerelyreducesthecostassociatedwithcalculatingeachfourcentertwoelectronintegral,methodsforscreeningtwoelectronintegralstoachieve 101 PAGE 102 TheFockbuildrequiresnonumericalintegration,buttheAATapproachdoesexplicitlyincludeexchangeandcorrelation.InthecurrentapproachwehaveusedtheminimalbasissetB3LYPexchangecorrelationpotentialsasareference,furtherimprovementispossiblebyevaluatinglargebasissetabinitioDFTexchangecorrelationpotentialsandsubsequentlyprojectingtheminaminimalbasissetrepresentation.ExtensionsoftheAATM0thatincludeamorecompleteexpansionofthetelescopingseriesfortheexchangecorrelationpotentialhavealsobeendescribedandtheirimplementationwouldbeafurtherimprovementontheresultspresented. 102 PAGE 103 Adaptedintegrals (AAjBC)IMulliken3centerTypeAIIRudenberg (ABjAC)IIIMulliken3centerTypeBIVRudenberg (ABjCD)VMulliken4centerVIRudenbergVIIMullikenPartialVIIIRudenbergPartial 103 PAGE 104 AveragepercentageofFockmatrixcontributionbyadaptedapproximation Rxn1C2C3CAIII3CBIIIIV4CVVIVIIVIII COHNO2CN0.6155.1140.5941.8540.722.282.522.361.401.411.401.381.39NMTCN0.6351.9941.0742.3040.854.154.874.472.162.252.202.172.16CH2OCO1.4962.5331.2731.7730.614.846.185.440.130.140.140.130.13CH3NH2CN0.5255.5633.1333.1632.198.0710.048.962.712.982.832.752.73HOOHOO2.2176.3419.8520.1519.220.891.191.100.710.850.790.760.74NETCN0.3541.8148.2848.9947.845.155.985.574.404.654.524.434.41C2H6CC0.1949.4035.4735.1834.399.4911.9510.605.455.935.685.555.50Average0.8656.1135.6736.2035.124.986.105.502.392.562.472.412.40 PAGE 105 Scatterplotsofagreementbetweenapproximateandexacttwoelectronintegrals.A)Onecenter.B)Twocenter.C)Threecenter.D)Fourcenter 105 PAGE 106 DissociationcurveofCNwith4centerRudenbergapproximation.A)NMT.B)NET.C)COHNO2.D)CH3NH2 PAGE 107 Orbitalproducts.A)Gaussianorbitalproduct.B)Slaterorbitalproduct 107 PAGE 108 DissociationcurveofCNwithAATapproximation.A)NMT.B)NET.C)COHNO2.D)CH3NH2 PAGE 109 DissociationcurveofCNwithAATM0approximation.A)NMT.B)NET.C)COHNO2.D)CH3NH2 PAGE 110 TimingofabinitioversusRudenbergfourcenterintegrals 110 PAGE 111 Percentageofmulticentertermsversussystemsize 111 PAGE 112 TimingofFockbuildusingtraditionalNDDOandNDDOwithcubicsplinesasafunctionofsystemsize 112 PAGE 113 Errorintroducedperatomfromcubicsplinesasafunctionofsystemsize 113 PAGE 114 PseudoreactionpathsplittingC20 PAGE 115 ForceRMSDfromMP2 115 PAGE 116 ForceRMSDfromB3LYP 116 PAGE 117 TheACESIIIprogramsystemenvironmenthasbeendesignedfortheecientparallelimplementationofwavefunctiontheorymethodsincomputationalchemistry.TheACESIIIprogramissuitedtotreatinglargesystems,typically5001000basisfunctionsandupto300electronswithpostHFabinitiomethods.Thesuperinstructionprocessor(SIP)managesthecommunicationandprocessingofblocksofdata,suchblocksareintendedtobesomewhatlarge(ontheorderof500,000oatingpointnumbers).ToprovideamoredirectrouteforimplementingnewmethodstheSIPreadsthecodewritteninthehighlevelsymboliclanguage,thesuperinstructionassemblylanguage(SIAL)(pronounced\sail").TheSIPhidesmanyofthemorecomplicatedaspectsofparallelprogramingallowingtheSIALprogrammertofocusonalgorithmandmethoddevelopment.WithinSIALeachindexofanarrayisdividedintosegments,thesegmentsmapelementsofthearrayintoblocks.Analgorithmtheninvolvesoperationsforwhichthesegmentisthefundamentalunitofindexing.Operationsthatinvolvethecontractionoftwoarrayscanthenbebrokenintoseveralsmallercontractionsovertwoblocks.Eachblockpaircontractioncanthenbedistributedovermanyprocessors,allowingfortheparallelimplementationofthecontractionoftwolargearrays. Ageneralissueinparallelcodesisensuringthatlatencydoesnotbecomearatelimitingstep.Thiscanbedonebymakingsurethattheamountofcomputationdoneisonaparwiththelatencyofanyoperation.Suchbalancingismachinedependent,soitrequiressomeamountoftestingtodetermineanoptimalsolutionforaparticularmachine.FortheimplementationoftheAATmethodologythisbalancingiscriticalandrequiresadierenttreatmentthanisusuallyusedforpostHFmethods. ThestructureoftheAATequationsarequitedierentthanpostHFequations.IntheformertheatomicindicesareusedastheprimaryvariableintheloopsusedtoconstructtheFocklikematrix.ForpostHFmethodstheorbitalindicesaretheimportantvariable.Moreover,thenumberoforbitalsperatominAATissmallrelativetothe 117 PAGE 118 ThespeedoftheAATapproachisachievedbyavoidingtheexplicitcalculationofthenumerousandcostlyfourcenterintegrals.Instead,thesetermsareintroducedbyasumofcontractedtwocenterterms.Initially,weusedsmallsegmentblocksizes.TheblocksizeultimatelydeterminesthecomputationaleciencybyaectinghowtheIO,orfetchingofdata,isdone.Intheseinitialstudieseachblockspannedonlyoneatom.Inthiscase,whenthetwoelectronintegralswereevaluateditwasstraightforwardtosimplyneglecttheintegralsthatinvolvedfouratoms(fourcenters.)Thisprocedureworkedforsystemswithfewatoms(<20),but,asthesystemsizeincreased,thenumberofblocksrapidlyincreased,andlatencybecameanproblem. Animprovedapproachthatwasimplemented,whichismoreconsistentwiththeoriginalACESIIIdesignphilosophy,involvedremovingtheatomspanningblocksizerestriction.Thiswasachievedbyrewritingtheoriginalroutinesthatreturnedalltwoelectronintegralsspanningfourblockunitstoexplicitlyexcludeintegralsinvolvingfourcenters.Inaddition,theRudenbergapproximation,whichreliesoncontractionsbetweentheoverlapmatrixandNDDOtypetwocentertwoelectronintegrals((AAjBB)),requiredmodicationtoavoidovercountingsometerms.Similartotheroutinethatexcludesfourcenterintegrals,themodiedroutinetogeneratetheintegralsusedbytheRudenbergapproximationexcludedallintegralsthatwerenotofthetype(AAjBB).Thesetwomodiedroutinesinvolvesomeoverheadthatcouldbereducedwithfurtherrenement.ThesemodiedroutinesarespecictotheAATapproachandarenotincludedinthestandardreleaseversionofACESIII. Sincetwoelectronintegralspossessaneightfoldspatialsymmetrytheoverallnumberofintegralsthatneedtobecalculatedcanbereducedbyroughlyafactorofeight.Todosorequiresrecognizingthesymmetryofeveryintegraltype.InACESIIItheproblemisslightlydierent.InsteadofhavingthiseightfoldsymmetrymanifestitselfbyAO 118 PAGE 119 Do EndIf(=) If6= Do (j)!(j);(j);(j);(j);(j);(j);(j) EndIf(6=and>) EndDo() EndIf(6=) If> [ApplyRudenbergapproximationtotheterms(j)andSj;k](j)!(j);(j);(j) Do EndIf(>and6=and6=) EndDo() EndIf(>and6=) EndDo() EndIf(>) EndDo() EndDo() InadditiontotheRudenbergapproximation,wehavealsoimplementedaninterfacetothestoredcubicsplineparametersneededtoconstructtheapproximationtotheexchangecorrelationcontributiontotheFocklikematrix.Thiscomponentismuchfasterthanthetwoelectronintegralcomponent,thereforeitisnotaratelimitingstep. 119 PAGE 120 Anotheraspectthatneedstobeconsideredistherotationofthetwocentermatrixelementsinalocalatompaircoordinatesystemtotheglobalcoordinatesystemofthemolecule.Thisrotationisunambiguousoncetheglobalcoordinatesystemisspecied,andalthoughorbitalsofthemoleculecanberotatedsimultaneouslywithoutaectingtheenergetics,thecorrespondingmatrixelementswillchangeonrotation.ThisdegreeofexibilitymustbeincorporatedtoensurethattheFocklikematrixcontributionthatarisesfromeveryatompaircorrespondstotheorbitalorientationoftheglobalcoordinatesystemofthemolecule.Toachievethisweuseasimpleprocedurethatgeneratesthematrixthatrotatesapairofatomiccoordinatessothattheyarealignedtotheaxisofthethestoredcubicspline,andthenperformthereverserotationonthatsubblockoftheFocklikematrix. Finally,thethirdcomponentimplementedwasthenumericalintegrationthatisperformedoncetheSCFprocedurehasconverged.TheACESIIIprogramsystemgeneratesajobarchivele(JOBARC)thatcanbereadbyACESIIexecutables,assumingithasbeengeneratedonthesamecomputerarchitecture.TotesttheecacyoftheAATapproximations,thequickestroutewastomodifytheACESIIxintgrtexecutable.ThexintgrtmoduleperformsnumericalintegrationonadensitytodeterminetheexchangecorrelationenergythatisneededfortheresidualelectronicenergycontributionofAAT.Eectively,xintgrtpullstheAOeigenvectormatrixfromtheJOBARC,aswellasotherinformationaboutthesystemthatwouldtypicallybegeneratedbythexvscf ksexecutableinACESIIwhenrunningaKSDFTcalculation(e.g.occupationnumbers).SinceACESIIIdoesnotcreateallofthenecessaryrecordsintheJOBARC,somesmalladjustmentstothexintgrtroutineledtoanewmodiedexecutablethatperformedthe 120 PAGE 121 121 PAGE 122 [1] R.J.Bartlett,V.F.Lotrich,andI.V.Schweigert,J.Chem.Phys.123,062205(2005). [2] K.W.Sattelmeyer,J.TiradoRives,andW.L.Jorgensen,J.Phys.Chem.110,13551(2006). [3] J.McClellan,T.Hughes,andR.J.Bartlett,Int.J.QuantumChem.105,914(2005). [4] R.S.Mulliken,J.Chim.Phys.46,497(1949). [5] K.Rudenberg,J.Chem.Phys.34,1433(1951). [6] W.WeberandW.Thiel,Theor.Chem.Acc.103,495(2000). [7] M.Elstner,D.Porezag,G.Jungnickel,J.Elsner,M.Haugk,T.Frauenheim,S.Suhai,andG.Seifert,Phys.Rev.B58,7260(1998). [8] Q.Zhao,R.C.Morrison,andR.G.Parr,Phys.Rev.A50,2138(1994). [9] R.vanLeeuwenandE.J.Baerends,Phys.Rev.A49,2421(1994). [10] A.BesteandR.J.Bartlett,J.Chem.Phys.120,8395(2004). [11] A.BesteandR.J.Bartlett,J.Chem.Phys.123,154103(2005). [12] W.KohnandL.J.Sham,Phys.Rev.140,A1133(1965). [13] R.J.BartlettandG.D.PurvisIII,Int.J.QuantumChem.XIV,561(1978). [14] J.^C^zekandJ.Paldus,J.Chem.Phys.47,3976(1960). [15] J.PaldusandJ.^C^zek,J.Chem.Phys.54,2293(1971). [16] R.J.Bartlett,Annu.Rev.Phys.Chem.32,359(1981). [17] R.J.BartlettandM.Musial,Rev.Mod.Phys.79,291(2007). [18] T.Helgaker,P.Jrgensen,andJ.Olsen,MolecularElectronicStructureTheory(Wiley,NewYork,2000). [19] R.J.Bartlett,TheoryandApplicationsofComputationalChemistry(Elsevier,2005),p.1191. [20] M.Levy,Proc.oftheNat.Acad.ofSci.76,6062(1979). [21] Q.Zhao,R.C.Morrison,andR.G.Parr,Phys.Rev.A50,2138(1994). [22] Q.ZhaoandR.G.Parr,Phys.Rev.A46,2337(1992). [23] Q.ZhaoandR.Parr,J.Chem.Phys.98,543(1993). 122 PAGE 123 D.J.Tozer,V.E.Ingamells,andN.C.Handy,J.Chem.Phys.105,9200(1996). [25] K.Peirs,D.VanNeck,andM.Waroquier,Phys.Rev.A67,012505(2003). [26] R.Homann,J.Chem.Phys.39,1397(1963). [27] M.J.S.Dewar,J.Friedheim,G.Grady,E.F.Healy,andJ.Stewart,Organometallics5,375(1986). [28] J.J.P.Stewart,J.Computat.Chem.10,209(1989). [29] J.J.P.Stewart,J.Computat.Chem.10,221(1989). [30] M.J.S.DewarandW.Thiel,J.Am.Chem.Soc.99,4899(1977). [31] M.J.S.Dewar,E.G.Zoebisch,E.F.Healy,andJ.J.P.Stewart,J.Am.Chem.Soc.107,3902(1985). [32] M.KolbandW.Thiel,J.Computat.Chem.14,775(1993). [33] M.J.S.DewarandW.Thiel,Theor.Chem.Acc.46,89(1977). [34] J.C.SlaterandG.F.Koster,Phys.Rev.94,1498(1954). [35] W.M.C.FoulkesandR.Haydock,Phys.Rev.B39,12520(1989). [36] R.J.Bartlett,I.Grabowski,S.Hirata,andS.Ivanov,J.Chem.Phys.122,034104(2005). [37] M.WolfsbergandL.Helmholz,J.Chem.Phys.20,837(1952). [38] C.J.BallhausenandH.B.Gray,Inorg.Chem.1,111(1962). [39] G.Klopman,J.Am.Chem.Soc.86,4550(1964). [40] K.Ohno,Theor.Chem.Acc.2,219(1964). [41] K.NishimotoandN.Mataga,Z.Physik.Chem(Frankfurt)12,335(1957). [42] G.Zheng,S.Irle,M.Elstner,andK.Morokuma,J.Phys.Chem.A108,3182(2004). [43] C.E.Taylor,M.G.Cory,R.J.Bartlett,andW.Thiel,Comput.Mat.Sci.27,204(2003). [44] R.J.Bartlett,J.Phys.Chem.93,1697(1989). [45] I.RossiandD.Truhlar,Chem.Phys.Lett.233,231(1995). [46] P.Charbonneau,Astrophys.J.(Suppl.)101,309(1995). [47] R.H.Byrd,LBFGSBversion2.1,UniversityofColorado.(1997). 123 PAGE 124 ACESIIisaprogramproductoftheQuantumTheoryProject,UniversityofFlorida.J.F.Stanton,J.Gauss,J.D.Watts,M.Nooijen,N.Oliphant,S.A.Perera,P.G.Szalay,W.J.Lauderdale,S.R.Gwaltney,S.Beck,A.BalkovaD.E.Bernholdt,K.KBaeck,P.Rozyczko,H.Sekino,C.Hober,andR.J.Bartlett.IntegralpackagesincludedareVMOL(J.AlmlofandP.R.Taylor);VPROPS(P.Taylor)ABACUS;(T.Helgaker,H.J.Aa.Jensen,P.Jrgensen,J.Olsen,andP.R.Taylor)(2006). [49] W.Thiel,MNDO97,OrganischchemischesInstitut,UniversitatZurich,Winterhurerstrasse190,CH8057Zurich,Switzerland(1997). [50] M.Dupuis,A.Marquez,andE.Davidson,HONDO99.6,basedonHONDO95.3,QuantumChemistryProgramExchange(QCPE),IndianaUniversity,Bloomington,IN47405(1999). [51] J.Li,F.Zhao,andF.Jing,J.Comput.Chem.24,345(2003). [52] R.Ahlrichs,M.Bar,M.Haser,H.Horn,M.Hser,andC.Klemel,Chem.Phys.Lett.162,165(1989). [53] M.HaserandR.Ahlrichs,J.Comput.Chem.10,104(1989). [54] F.Weigend,M.Haser,H.Patzelt,andR.Ahlrichs,Chem.Phys.Lett.294,143(1998). [55] P.PolitzerandJ.S.M.Eds.,TheoreticalandComputationalChemistryEnergeticMaterialsPart1.Decomposition,CrystalandMolecularProperties,vol.12(Elsevier,Amsterdam,TheNetherlands,2003). [56] M.J.S.Dewar,J.P.Ritchie,andJ.Alster,J.Org.Chem.50,1031(1985). [57] M.L.J.McKee,J.Am.Chem.Soc.108,5784(1986). [58] M.T.Nguyen,H.T.Le,B.Hajgato,T.Veszpremi,andM.C.Lin,J.Phys.Chem.A107,4286(2003). [59] W.F.Hu,T.J.He,D.M.Chen,andF.C.Liu,J.Phys.Chem.A106,7294(2002). [60] A.G.TaubeandR.J.Bartlett,J.Chem.Phys.accepted(2007). [61] W.D.Taylor,T.D.Allston,M.J.Moscato,G.B.Fazekas,R.Kozlowski,andG.A.Takacs,Int.J.Chem.Kin.12,231(1980). [62] H.A.TaylorandV.V.Vesselovsky,J.Phys.Chem.39,1095(1935). [63] R.Bukowski,K.Szalewicz,andC.F.Chabalowski,J.Phys.Chem.A103,7322(1999). [64] M.J.S.DewarandW.Thiel,J.Am.Chem.Soc.99,4907(1977). [65] W.Thiel,J.Mol.Struc.(Theochem)16,398(1997). 124 PAGE 125 W.WeberandW.Thiel,Theor.Chem.Acc.103,495(2000). [67] J.Stewart,ReviewsinComputationalChemistry(VCHPublishers,1990),vol.1,p.45. [68] W.Thiel,Adv.Chem.Phys.93,703(1996). [69] D.E.Taylor,K.Runge,andR.J.Bartlett,Mol.Phys.103,2019(2005). [70] W.ThielandA.Voityuk,Theor.Chem.Acc.81,391(1992). [71] A.Ourmazd,D.W.Taylor,J.A.Rentschler,andJ.Bevk,Phys.Rev.Lett.59,213(1987). [72] D.Muller,T.Sorsch,S.Moccio,F.Baumann,K.EvansLutterdodt,andG.Timp,Nature399,758(1999). [73] X.Y.Liu,D.Jovanovic,andR.Stumpf,AppliedPhysicsLetters86,82104(2005). [74] A.A.Demkov,J.Orteg,O.F.Sankey,andM.P.Grumback,Phys.Rev.B52,1618(1995). [75] L.StixrudeandM.S.T.Bukowinski,Science250,541(1990). [76] L.StixrudeandM.S.T.Bukowinski,Phys.Rev.B44,2523(1991). [77] M.MenonandK.R.Subbaswamy,Phys.Rev.B55,9231(1997). [78] M.Elstner,J.Phys.Chem.A111,5614(2007). [79] N.Otte,M.Scholten,andW.Thiel,J.Phys.Chem.A111,5751(2007). [80] M.P.Repasky,J.Chandrasekhar,andW.L.Jorgensen,J.Comput.Chem.23,1601(2002). [81] P.Pulay,Mol.Phys.17,197(1969). [82] B.Dunlap,J.Connolly,andJ.Sabin,J.Chem.Phys.71,3396(1979). [83] C.A.White,B.G.Johnson,P.M.W.Gill,andM.HeadGordon,Chem.Phys.Lett.230,8(1994). [84] C.A.White,B.G.Johnson,P.M.W.Gill,andM.HeadGordon,Chem.Phys.Lett.253,268(1996). [85] M.C.Strain,G.E.Scuseria,andM.J.Frisch,Science271,51(1996). [86] J.C.Burant,G.E.Scuseria,andM.J.Frisch,J.Chem.Phys.105,8969(1996). [87] M.ChallacombeandE.Schwegler,J.Chem.Phys.106,5526(1997). 125 PAGE 126 P.M.W.Gill,A.T.B.Gilbert,andT.R.Adams,J.Computat.Chem.21,1505(2000). [89] W.Koch,B.Frey,J.F.S.Ruiz,andT.Scior,Z.Naturforsch58a,756(2003). [90] N.OliphantandR.J.Bartlett,J.Chem.Phys.100,6550(1994). [91] C.S.Duris,ACMTOMS6,92(1980). [92] L.A.Curtiss,K.Raghavachari,P.Redfern,andJ.A.Pople,J.Chem.Phys.106,1063(1997). 126 PAGE 127 IwasborninMichiganbuthavedonetimeinOhio,Washington,Oregon,andnowFlorida.Mytumultuouscareerinsciencebeganinsixthgradewhenmyscienceteachercalledmeadreamer(IknowIamnottheonlyone)forproposingtheconstructionofanelevatortothemoontoavoidthecostsassociatedwithspaceight.Inhindsight,Ishouldnothavespeciedmyoriginaldesignasarigidstructure.ThatsameyearIwasawardedacerticateastheclassroom'sstrangestkid.Aftermakingthetransitiontohighschool,IdoubleduponmathandscienceandnishedallavailableclassesinthoseareasbymyJunioryear,resultinginamisspentSenioryear.SomehowImanagedtogetintotheonlycollegetowhichIapplied.IarrivedaFreshmanlledwithyouthfuloptimismandpersuadedmyacademicadvisortoenrollmeinachemistrycourse2yearsbeyondmylevel.IcontinuedmySophomoreyearwithouttheyouthfuloptimism.IparticipatedintwoResearchExperienceforUndergraduatesprogramssponsoredbytheNationalScienceFoundationatCornellandtheUniversityofChicago.InChicagoIworkedforProf.KarlFreed,essentiallyspendingthesummerbuildingoverlycomplicatedinternalcoordinateinputles.ImusthaveenjoyeditbecausethatsamesummerItattooedmyselfwiththemarkofaquantumchemist,literally.Eventually,ItookaninterdisciplinaryBAinchemistryandphysicsfromReedCollegeinbeautifulPortland,Oregon,in2001. 127 