UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
TRIPLE EXCITATION EFFECTS IN THE FOCKSPACE COUPLED CLUSTER METHOD By DAVID EDWARD BERNHOLDT 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 1993 To Jan, for so many reasons. ACKNOWLEDGMENTS Graduate school has been an educational experience in many ways, and a great many people have made important contributions to that education. I would especially like to thank my advisor, Rod Bartlett, for the singular role he has played in my life. His research group members, especially John Stanton, have also shared knowledge and offered both help and challenges throughout my work with them. My colleagues, past and present, in the Quantum Theory Project have created a wonderful environment in which to learn and do research, and it has been my great pleasure to experience it. I would like to thank Drs. Vala and Weltner for numerous discussions and suggestions; I regret that I did not have time to examine more of their interesting molecules. My research in this field began while I was an undergraduate at the University of Illinois, and I am deeply indebted to Cliff Dykstra and Shiyi Liu for getting me started and for conveying to me their enthusiasm for the subject. The daily existence of a graduate student often revolves around research, but there is much more to life. I thank my friends and family, who have tried to tell me this and who have shared other gifts with me over the years: humor, diversion, perspective, encouragement, sympathy, and love. Finally, words cannot express my gratitude to Jan, who rode two different roller coasters at the same time and lived to tell the tale. TABLE OF CONTENTS ACKNOWLEDGMENTS ......... ABSTRACT ................. CHAPTERS 1 INTRODUCTION ......... 2 THEORETICAL BACKGROUND ............ . 5 2.1 Overview ................ 2.2 The Schrodinger Equation ...... 2.3 Potential Energy Surfaces ...... 2.4 Calculations on Potential Energy Sul 2.5 Vibronic Interactions ......... 2.6 Basis Set Expansion ......... 2.7 Oneelectron Approximations ... 2.8 The Correlation Problem ....... 2.9 Occupation Number Representation. 3 COUPLED CLUSTER METHODS . 3.1 Introduction .............. 3.2 Basic Equations ............ 3.3 Reduction to Spin Orbital Equations faces . . . 29 . iii . vii . . . . ........................ 3.4 Numerical Solution of the Coupled Cluster Equations 3.5 The Coupled Cluster Effective Hamiltonian . 3.6 Truncation of the Cluster Operator . . 4 FOCK SPACE COUPLED CLUSTER THEORY ...... 4.1 Introduction ........................ 4.2 Structure of FSCC .................... 4.3 Comparison with Single Reference CC Theory . 4.4 Basic Equations ..................... 4.5 Reduction to Spin Orbital Equations . . 4.6 Approximations to the Full FSCCSDT Model ..... 5 A CRITICAL ASSESSMENT OF FSCCSD AND THE TI TRIPLES APPROXIMATIONS ............... 5.1 Introduction ........................ 5.2 Effect of Changes in the Active Space . . 5.3 Ionization Potentials ................... 5.4 Potential Energy Surfaces and Related Properties . 5.5 Conclusions ........................ 6 THE STRUCTURE OF THE C,+ MOLECULE ...... 6.1 Background ........................ 6.2 Computational Methods ................. 6.3 The Heart of the Problem ................ V . . 33 HIRDORDER ........, . 33 . 38 . 42 S. 42 S. 43 . 50 S. 54 . 55 . 56 S. 64 S. 64 . 66 S. .71 . 79 S. 87 . 90 S. 90 . .91 . 93 6.4 Geometries of the Linear Stationary Points of C5 . ... .95 6.5 Relative Energies ........... ........ ... ................97 6.6 Stability of the Stationary Point Structures . ... 102 6.7 Implications for Other Carbon Clusters . .... 105 6.8 Recapitulation ................................. .. 108 7 SUM M ARY ..................................... 109 APPENDIX IMPLEMENTATION OF THE PERTURBATIVE TRIPLES MODELS ...... 115 REFERENCES ........................ ................. 125 BIOGRAPHICAL SKETCH ........................ ......... 131 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 TRIPLE EXCITATION EFFECTS IN THE FOCKSPACE COUPLED CLUSTER METHOD By DAVID EDWARD BERNHOLDT May 1993 Chairman: Rodney J. Bartlett Major Department: Chemistry The Fockspace coupled cluster (FSCC) approach to the electron correlation problem offers several advantages over the more familiar single reference coupled cluster (CC) methods. It is capable of treating states with a multireference character, and it is capable of providing statetostate energy differences for several states in a single calculation. As explained in this work, basic differences between the two theories result in higherexcited cluster operators being more important in FSCC methods than CC methods. In order to improve the accuracy of FSCC calculations, this work involves the addition of triple excitation effects to the widely applied FSCC model restricted to single and double excitations (FSCCSD) for ionization potentials. Equations for the full model, including single, double, and triple excitations (FSCCSDT), are derived and various approximations are proposed as a compromise between the improvement in the wavefunction and energy due to the inclusion of triples and the computational costs, which are significantly larger than for the FSCCSD model. Two approximations to the full FSCCSDT model have been implemented and their performance is examined for a variety of molecules. Results are found to be generally good, providing results comparable to single reference coupled cluster methods with approximate inclusion of triple excitations. Finally, both CC and FSCC methods are used to study the structures of several small carbon clusters, a class of molecules which has been shown to require highquality theoretical methods and which can make use of both FSCC's multireference and direct energy difference capabilities. Detailed calculations show that the ground state of C5+ is a 2E+, with Coov symmetry, but with low frequency bending vibrational modes. Additional calculations for other odd carbon cluster cations suggest that the "linear" structures of C7' and C9' observed experimentally may in fact be distorted from true linearity due to the RennerTeller effect. CHAPTER 1 INTRODUCTION In recent years, advances in both experimental and theoretical chemistry have fueled an exciting synergy between the two fields [1]. Experimental innovations have allowed study of molecules in conditions similar to those commonly employed in theoretical calculations (low temperature, low density), producing data of increasing accuracy and specificity. Likewise, through increases in computational power and methodological sophistication, theoretical chemists have been able to tackle increasingly complex systems with high accuracy. Where experimental and theoretical data are directly comparable, critical evaluations of both approaches are possible. In other cases, the data are often complementary, allowing one or the other researcher to bring the problem into finer focus. This has been especially apparent in research on atomic clusters, a field which has expanded rapidly in the past decade. Clusters are interesting for a variety of rea sons. Metal clusters are used as models for the study of catalytic processes. Clus ters of lighter atoms, particularly carbon, have been identified in circumstellar space and are expected also to be found in interstellar molecular clouds, where they may be important in the creation of new stars. They have also been found in combustion products on Earth. Clusters may have useful and unique physical properties, mak ing them interesting from a materials viewpoint, and they are studied as models of how properties of bulk materials are approached by increasingly large collections of atoms. 2 Although the literature on clusters is expanding rapidly, the sought after understanding has not yet emerged. One reason for this is that clusters, even relatively small ones, are rather unlike "traditional" molecules which scientists have studied for centuries, and they often defy the "rules" and intuitions which researchers try to bring from traditional chemistry to cluster chemistry. Consequently, workers in the area are having to develop new ideas and new "rules" about the structure and properties of clusters. An important component of creating this new understanding is having reliable data on the structure and properties of many clusters. Since there are few experiments which give direct, unambiguous structural data, theoretical approaches are particularly useful in this context. When the capabilities of experiment and theory do overlap, it is possible to compare results critically and perhaps improve both. To be truly useful, however, the theoretical methods employed must be accurate and reliable. The subject of this dissertation, the inclusion of triple excitation effects in the Fockspace coupled cluster (FSCC) method, was undertaken in response to exactly this need. Coupled cluster (CC) and related manybody perturbation theory (MBPT) methods (see Chapter 3) have proven their ability to deliver extremely highquality properties and energetic for a wide range of complex systems, and because of their size extensivity (proper scaling of the energy with the size of the system under study) [2,3], are preferred over configuration interaction (CI) methods for theoretical studies of clusters. Not every situation, however, is ideally suited to the use of traditional single reference CC methods. For example, closelying electronic states may cause problems for even the most sophisticated CC methods, thus requiring a multireference description of the wavefunction 3 to obtain meaningful results. Ionization potentials (IPs) and electron affinities (EAs) are often obtained using experimental studies of clusters in photoionization and photoelectron spectroscopies. In order to provide comparable theoretical results for interpretation or confirmation of the spectra, separate calculations for the ground state and each excited state of interest are required when using familiar CC methods. It would be useful to have a method capable of giving IPs and EAs for several states of interest in a single calculation. Fockspace coupled cluster methods are one way to approach both of these problems. Often referred to as Fockspace multireferenme coupled cluster theory, FSCC (see Chapter 4) can provide a multireference description of chosen electronic states of a molecule. Considered in another context, it can calculate the energies of several electron detached or attached states at the same time, thus simplifying the calculation of IPs and EAs. Because of formal differences between single reference CC and FSCC theories, however, a given truncation of the two methods does not necessarily give results of equivalent quality. For example, limiting both methods to single and double excitations, we find that an FSCCSD description of a particular state may leave out contributions which could be present in the CCSD description of the same state. Consequently, in order to improve the quality of the FSCC description, this work examines the inclusion of triple excitation effects in the FSCC method. After setting the stage with a review of theory and terminology leading up to the correlation problem, the traditional coupled cluster equations are presented, emphasizing their use as an "input" to the FSCC method. Then, in some detail, equations for the FSCC 4 method are derived, focusing on the FSCCSDT model for IPs. The importance of triple excitations is examined, in comparison with single reference CCSD and CCSDT methods, and several approximations are proposed as a compromise between the computational cost of evaluating the complete triples contribution and quality required of the wavefunction and energies. Although this development was spurred by an interest in carbon clusters, FSCC methods, with or without triple excitations, are by no means limited to this application. In order to assess the performance of the FSCC methods, calculations are presented for a variety of molecules. Finally, Fockspace and single reference coupled cluster methods are used to study several carbon cluster cations. CHAPTER 2 THEORETICAL BACKGROUND 2.1 Overview Much of the theoretical chemistry literature assumes the reader is a practitioner in the field. This chapter is intended to provide a common ground for the concepts and notation used throughout this dissertation. It is by no means exhaustive, but should help non specialists to understand the material that follows. References are given to the literature for more detailed treatments of various topics. 2.2 The Schrodinger Equation We seek to solve the timeindependent Schr6dinger equation, HI}) = El\ ) (2.1) for a molecular system (one or more molecules of interest isolated from interaction with anything else). This equation tells us that a system in an eigenstate I ) of the Hamiltonian operator H has an energy E. The usual nonrelativistic Hamiltonian can be divided into five different terms, H = Te + Tn + Ven + Ve + Vnn (2.2) representing, respectively, the electronic and nuclear kinetic energies and the potentials due to interaction of electrons and nuclei, electrons with other electrons, and nuclei with 6 other nuclei. The individual terms, in Hartree atomic units, are defined as T(r) = 1 2 T.(R) = AV A Ven(r, R) = ZA (2.3) AIri RAI(2.3) iA 1 A Vee ( ) = 1 1 (R 1 ZAZB Vn2(R) = E IRA RBI AOB In these expressions, i and j refer to electrons, A and B refer to nuclei. RA and ri refer to the spatial coordinates of nucleus A and electron i, respectively, and Z, and MA are the charge and mass of nucleus A. The unsubscripted symbols r and R refer to the complete set of position vectors of the electrons and nuclei, respectively. 2.3 Potential Energy Surfaces The terms of the full Hamiltonian, Eq. 2.2, can be grouped according to their dependence on the positions of the nuclei into a purely electronic component, He = Te + Vee, (2.4) a purely nuclear term, Hn = Tn, (2.5) and a potential which couples them, Vnt(r, R) = Ven(r, R) + Vn(R). (2.6) If we choose one set of positions R as the origin, denoted by 0, we can expand the electronnuclear interaction potential near this point in a Taylor series, Vi (r,nR)= V t(r,o0)+ :(RA A )A B RA+t ) AR +. (2.7) A ARA, 0RAdRB o The first term in the expansion can be thought of as the potential felt by the electrons due to a fixed set of nuclei, leading to the electronic Schrodinger equation, (He(r) + V.t(r,0))lPk(r)) = E'Ik(r)), (2.8) for a particular state Ik(r)) with energy Ek. In order to see how varying the nuclear positions effects these solutions, we must return to the full Schrodinger equation 2.1. By expanding the full wavefunction in terms of the solutions to the electronic wavefunction and a nuclear component, I (r, R)) = 1 Dk(r))Xk(R), (2.9) k and recognizing that the Hamiltonian of Eq. 2.8 is a subset of the full Hamiltonian, we obtain {He + Vnt,o + H, + W} IE Ik)Xk = EE kk)Xk, (2.10) k k where V,,t,o denotes the leading term in the Vnt expansion, and Wdenotes the remainder, W(r, R) = Vint(r, R) Vnt(r,0), known as the vibronic interaction. He + Vit,o is often 8 referred to as the electrostatic Hamiltonian, and is dependent on the nuclear coordinates only parametrically. Projecting on the left with an electronic eigenfunction, (4m(r)l, we obtain a coupled set of equations, (note that H, is independent of the electronic coordinates r), (E" + H,)Xm + E( Im Wlk)Xk = EXm, (2.11) k or using the more compact matrix element notation for W, (E E' H,)Xm = 1 WmkXk. (2.12) kfm If there is no vibronic interaction between electronic states (Wink = 0), then the equations decouple, giving (H.(R) + Em)Xm(R) = EXm(R) V m, (2.13) which is the Schridinger equation for the nuclei moving in the mean field of electronic state I\4m(r)), and the solution to the full Schridinger equation is simply I(r, R)) = I m,, (r))Xm(R). This is the origin of the idea of a potential energy surface (PES) as the set of nuclear configurations which satisfy Eq. 2.13 for a particular electronic state. Taking Wmk = 0 is known as the simple adiabatic or BornOppenheimer approximation [4,5]. In practice, although Wnk is not in general zero, the BornOppenheimer approxima tion is usually quite good because of the fact that the nuclei of a molecule are about 1836 times more massive than the electrons, so we can usually think of the nuclei moving slowly in the average field of the electrons, which are able to adapt almost instanta neously to the nuclear motions. This approximation is used in the majority of electronic 9 structure calculations to date, which simply evaluate the electronic wavefunction and energy, Eq. 2.8, at various points on the PES (nuclear configurations) and ignore the nuclear dynamics. 2.4 Calculations on Potential Energy Surfaces The BornOppenheimer potential energy surface for a particular state of a molecule with N nuclei is (neglecting invariance to rigid translations and rotations of the entire nuclear framework) a 3Ndimensional space which can be characterized according to its local shape. A stationary point is one at which the gradient of the energy with respect to nuclear displacements be zero. Stationary points may be minima, transition states, higher order saddle points, or maxima. If the curvature (second derivative) at a stationary point is positive, the point is a local minimum, which is a locally stable geometric conformation for the state in question. If the curvature in one of the 3N directions is negative, with all others positive, the point is a transition state, which is a local maximum on a path connecting two local minima; transition states are very important in the theory of chemical reactions. If the stationary point has more than one direction of negative curvature, it is a saddle point (or local maximum), which is not chemically meaningful. Stationary points on potential energy surfaces can be located using numerical opti mization techniques, such as NewtonRaphson, quasiNewtonRaphson, etc., but because individual calculations of the energy, gradient, and hessian (matrix of second derivatives of the energy), as required by the optimization method, are extremely expensive relative to evaluation of the optimization step, it is usually necessary to use chemical intuition 10 to choose regions of the PES on which to focus and sometimes to restrict the search to a subset of the full 3Ndimensional space. The hessian, which allows a stationary point to be characterized as a minimum, transition state, or saddle point, is also the primary input for evaluation of the vibrational frequencies at that structure (in the harmonic approximation). Since different molecules and different structures will generally have different vibrational spectra, it serves as a useful tool for the identification and characterization of unknown species. If the electronic state of the molecule of interest is not known from other information, it may be necessary to consider potential energy surfaces for several states in order to locate the lowest energy global minimum of all the states, which is known as the ground state of the system. 2.5 Vibronic Interactions Although the BornOppenheimer approximation and the idea of potential energy surfaces have served quantum chemistry well, there are certainly cases in which the W, matrix elements connecting different electronic states are significantly different from zero, and the BornOppenheimer picture is no longer adequate. Consider the vibronic interaction term, W, truncated after the quadratic term, which is adequate for purposes of explanation. Making a linear transformation from cartesian coordinates R to symmetrized normal coordinates Q (obtained by diagonalizing the hessian matrix), we have W(r,1 Q) = v (ai.tQ + Q 1 Z 0 t QaQ6. (2.14) : \d_2 EQ ) + 0 Each normal coordinate, Qc, belongs to an irreducible representation (irrep) Fa of the molecular point. Likewise, the electronic states Im) belong to an irrep, Fm. Thus from group theory, we know that the linear vibronic constant, F(kFm) = (2.15) F\ Q\9aQ I is nonzero if and only if rFo [Fk x Fm]. This is a general result for the interaction of any two states, but of particular importance is the case in which the two electronic states belong to the same irreducible representation. If the irrep is not degenerate, then [Fk x rk] = Fk x Fk = FA1, so that the coupling is via totally symmetric motions, which by definition cannot change the symmetry of the molecule. If, on the other hand, the two electronic states (or components of states) belong to a degenerate irrep, the direct product results in nontotally symmetric representations, and thus the states are coupled by motions which break (lower) the symmetry of the molecular framework. This is the basis of the JahnTeller (JT) theorem [6], which says that for a nonlinear polyatomic molecule in an nfold degenerate electronic state, the high symmetry point at which the n potential surfaces intersect cannot be a minimum on all n surfaces. In other words, there will always be at least one direction for which the energy is lowered when the degeneracy is lifted (by moving away from the high symmetry point). The tendency of a molecule in a degenerate state to lower its symmetry and thereby lift the degeneracy and lower its energy is known as the JahnTeller effect. 12 In the more general case, interaction to two nearlying states via a suitable nuclear motion causes the surfaces to repel each other. Looking at a slice of the two surfaces along the symmetrybreaking motion that couples the two states, the upper surface becomes more curved as the interaction increases, while the lower state is flattened. If the interaction is strong enough, the lower state may display a double minimum shape, with the original high symmetry point a transition state between the two minima. This is usually referred to as the pseudoJahnTeller effect (PJTE) [7]. The magnitude of the PJT effect depends both on the separation of the two states and the magnitude of the linear vibronic constant. The exclusion of linear molecules from the JT effect is also on the basis of symmetry arguments. The wavefunction of a linear molecule is either even or odd with respect to reflections containing the axis of the molecule or perpendicular to it, so of course the product of such functions will always be even. The nuclear motions that bend the molecule are odd in this respect, and so can never be contained in the direct product of the wavefunction representations the linear vibronic constant F(k k) is zero for linear molecules. For linear molecules, and other cases for which the linear vibronic constant is zero or negligibly small, the quadratic term of Eq. 2.14 comes into play. For linear molecules in degenerate states, the quadratic vibronic coupling leads to a condition similar to the JT effect, called the RennerTeller effect (RTE) [8,9]. Lifting of the degeneracy can be easily visualized in this case. Consider a linear molecule oriented along the z axis, in a 2II state, so that one component of the 7r orbital, say xrx contains an electron. Bending motions 13 in the orthogonal directions x and y would normally be equivalent in a linear molecule, but if one component of the 7r is occupied, the electronic state of the bent structure will be different depending on whether the electron is now in an orbital in the plane of the molecule or perpendicular to it. Unlike the JT case, the RennerTeller effect does not guarantee that lifting the degeneracy by bending will result in a lower energy. The two states will not be degenerate (except the linear geometry), but neither, one, or both may still be minima at the linear geometry, or the linear geometry may be a transition state between two equivalent bent minima depending on the magnitude of the interaction. In the more general case, with states belonging to different irreps interacting, there is little to be gained from a symmetry analysis of the quadratic vibronic constant. There are usually many combinations of vibrational motions that can couple two electronic states for either degenerate or nondegenerate point groups, however these effects are usually relatively small. Clearly vibronic coupling effects can cause major problems for computational meth ods based on the BornOppenheimer approximation, but if we are looking for minima on potential surfaces, it is clear that such points are usually not of interest when the interactions are strong enough, they cause the symmetry of the molecule to be broken, thus moving it away from the degenerate or pseudodegenerate point. As a result, it is often enough to recognize points on the PES at which strong vibronic interactions occur and then locate and follow the symmetry breaking to lower energies. Of course this approach will not always work, and even calculations near a (pseudo) degenerate point can be quite taxing for the method, but it has been used in a number of studies 14 (see Bersuker's bibliography [10] for examples). In 1984, Lee et al. [11] demonstrated that analytically evaluated hessians could be used to compute bending frequencies for RennerTeller molecules. In fact if due care is taken to stay on the correct PES during the calculation, hessians evaluated from finite differences of gradients, or even energies, can also give RT bending frequencies. 2.6 Basis Set Expansion Several approaches may be taken to the solution of the electronic Schrodinger equation, but the most widely used and most generally applicable is expansion in a basis of Gaussian functions. These functions are usually (but not necessarily) chosen to mimic the atomic orbital (AO) description used in chemistry and are referred to by that name. A complete (infinite) basis would span all of space and thus allow an exact description of the wavefunction. In practice, however, computational resources place limits on the size of basis that may be employed, and it is necessary to compromise between the cost of the calculation and the accuracy required. It is now possible to evaluate matrix elements of the Hamiltonian operator in terms of the AO basis set. For a set of functions {x,}, we have, for example, the one and twoelectron integrals (or matrix elements) [12] (lhlv) = h, = J (xi)hx;(xl)dxi (2.16) (Pv" Ar) = X,(xl)X*(x2) X(xA)x(x2)dxd I r12 where r{ a(wi), or xi = () (2.17) is a combination of the position and spin of electron i, and r12 = Irl r2l. The lefthand sides of Eq. 2.16 are shown the braket notation, originally due to Dirac [13], which is in common use and for the oneelectron integral, a matrix element notation is also shown. Both the braket and matrix element notations are used in the literature, and both will be used in this work. 2.7 Oneelectron Approximations Given a basis of atomic orbitals, which are (very rough) approximations to the eigenstates of the atom, the first step in an electronic structure calculation is usually to find a set of oneelectron functions for the molecule as a whole, called molecular orbitals (MOs). These orbitals correspond to the qualitative ideas about molecular orbitals often used by chemists.The most common prescription for calculating MOs is the HartreeFock (HF) SelfConsistent Field (SCF) procedure [12]. The result of an SCF calculation is a set of MOs {(,}, which are represented as linear combinations of the AO basis functions, p = CSPXA (2.18) with the matrix C determined by the SCF procedure. Of these MOs, some will be occupied by electrons and some will be unoccupied, or virtual orbitals. The 16 antisymmetrized product of the occupied orbitals is the Slater determinant for that state. Depending on the particular formulation (primarily the treatment of electrons with different spins) used in the HF procedure, it may be referred to as RHF, UHF, or ROHF, among other possibilities. The most straightforward is restricted HartreeFock (RHF), applicable only to closed shell molecules (where all electrons are paired), where the alpha and beta spin orbitals in each pair are restricted to have the same spa tial function. The restricted openshell HartreeFock (ROHF) formalism extends the RHF idea to open shell systems by requiring that all paired electrons have share the same spatial function, but allowing additional orbitals which may be occupied by sin gle unpaired electrons. In the unrestricted HartreeFock (UHF) method (also known as DODS different orbitals for different spins), the alpha and beta spin electrons have distinct spatial parts, allowing the possibility that they might be quite differ ent. Only when pairs of electrons have the same spatial component (RHF, ROHF) is it guaranteed that the HF wavefunction is an eigenfunction of the S2 operator, satisfying S21HF) = s(s + 1)IHF). UHF wavefunctions may suffer from spin contamination, in which the wavefunction is an admixture of states of different S9 values. Because the true wavefunction is always an eigenfunction of 52, it is thus possible that a UHF wavefunction may be to some extent unphysical. Nevertheless, the UHF method has proven quite useful, particularly as the basis for formulating general (i.e. applicable to both closed and open shell systems) correlated methods. 2.8 The Correlation Problem Although the SCF procedure provides a very useful qualitative description of molecules, it is generally inadequate for quantitative uses requiring high ac curacy. The method considers each electron the average field of all others, which ignores the fact that the motion of each electron is instantaneously cor related with all others because of the Pauli Principle. For the higher ac curacy required of most interesting applications, it is necessary to go be yond the oneelectron approximation and treat correlation effects in the sys tem. The correlation problem is usually formulated in terms of the interaction between different configurations (occupations) of a set of oneelectron functions. The molecu lar orbitals are usually used for this purpose because the idea of occupied and unoc cupied orbitals is well defined and the derivation of equations for correlated methods are thus greatly simplified. The SCF determinant serves as a reference for the defi nition of the excited state configurations in which the correlated wavefunction is ex panded. 2.9 Occupation Number Representation The occupation number, or second quantized, formalism [14] is a convenient representation for the development of correlated methods because it represents in a compact way the changes in occupation of the various determinants. In this de 18 scription and throughout this dissertation, the convention will be used that subscripts i, j, ., n refer to orbitals which are occupied in the reference determinant, a, b, ., f refer to virtual orbitals, and p, q, ., u refer to orbitals of either type. Creation and annihilation operators are defined that act to create or destroy an electron in a particular oneelectron function (MO). Thus the creation operator at would create an electron in orbital p, and the annihilation operator aq would remove an electron from orbital q. The Pauli Principle is satisfied by the fact that aa = aq = 0, which is a result of the anticommutation relations among the creation and annihilation operators. Strings of creation and annihilation operators act on the reference determinant, 10), to produce excited state determinants. For example, an electron can be promoted from an occupied orbital i to a virtual orbital a, aail0) = I0) = ') (2.19) where the resulting excited state has been shown in two common notations. Similarly, higher excitations may be obtained ata"aajiO) = aa0 = 0) a aka aaa l) = ~i ) = i?) (2.20) Operators such as the Hamiltonian are represented in the second quantized formalism by the matrix element of the operator in the basis of oneelectron functions multiplied 19 by a string of creation and annihilation operators: H hpqaaq + (pql rs)ataasar (2.21) p,q p,q,r,s where we have used the antisymmetrized twoelectron integral, (pqllrs) = (pqlrs)  (pqlsr). Strings of creation and annihilation operators may be evaluated using Wick's theorem [15] and the concept of a normal order of operators. Normal ordered operator strings are those which have all at and all aa to the right of the string, which can be accomplished using the anticommutation relations for second quantized operators. This is done to exploit the fact that a!10) = aalO) = 0. The normal order of an operator string A is denoted N[A], and has the property (OIN[A]I0) = 0. (2.22) This result is obvious if all a! and aa operators are brought to the right, where they act on the reference giving zero immediately. If there are no at or aa, then the operators act on the reference to create an excited state on the left hand side, which is orthogonal to the reference function on the right, (/O \ =0. (2.23) Note too that in the absence of at or aa, A=N[A]. Wick's theorem tells us that the contraction of a pair of operators is defined as aq = N t aaq] + N a aq or (2.24) apaq = N apart ] + N apa , 20 with the following definitions according to whether p and q are occupied or virtual orbitals (they must be the same type since 6ai = 0) N [aaj]= 6ij N[aia = 0 (2.25) N aab =0 N [aaaa = ab Longer strings of operators are evaluated by taking all possible single contractions, then all double contractions (two simultaneous contractions), up to full contractions, with the sign in each case determined by the permutations necessary to make pairs of contracted operators adjacent. As an example, consider the Hamiltonian, Eq. 2.21. The oneelectron term is cast in normal ordered form by recognizing that ataq = N a4aq + N [aq] (2.26) = N a[aq] + bij so that Shpqaaqy = hpqN [aaq] + hii. (2.27) p,q p,q i The twoelectron term has a string of four operators, so all possible single and double contractions must be considered. There are four single contractions all of which can be 21 manipulated into the same form, and similarly for the two double contractions, 1 1 (pqjrs)N[aataa] 4 p,g,r,s p q 4 q q p,q,r,s (2.28) + (pil qi)N [aaq] + ( ij) p,q,i i,j When evaluating a product of several normal ordered operators, the same rules hold, but all contractions which are internal to any given operator have already been accounted for, so only those contractions connecting different normal ordered operators need be considered. When developing correlated methods, it is convenient to employ the normal ordered product form, ON, of a second quantized operator, 0, defined by ON 0 (01010). (2.29) This form has the property, (0oON10) = (010 (01o00)10) = (01010) (01010) = 0 (2.30) which allows one to focus on the correlation correction to the result rather than carrying the reference expectation value (a constant) through all of the equations. For example, consider the Hamiltonian. The normal ordered operator is the sum of Eqs. 2.27 and 2.28, or N[H] h= hN[aq + ap,qi ]qi)N[a] + (pjjpqr, (Pqlrs)N[ap ,ar] p,q + hi + Y (Iijlij) i i,j (2.31) and the expectation value is (01HI0) = 0 h hpqaaq + P ,,r, (pqrs)aptaaar 0 \ P pqrs q (2.32) = hii + i (ijllij) The normal ordered product form of the Hamiltonian is thus HN = H (OIHIO) = hpN [ataq] + ,,(piqi)N [a + pq,,s(pqllrs)N aataa p,q = fpN [aaq] + (pqllrs)N aaasaT] S4 p,,r,s Q p,q (2.33) with the Fock operator, fN = fN [ataq] hqN [aa]+ (pi [qi) N[aaq] (2.34) p,q p,q p,q,s used in the final form. The normal ordered product Hamiltonian is also frequently written in terms of its one and twoelectron components, HN = fN + WN. (2.35) CHAPTER 3 COUPLED CLUSTER METHODS 3.1 Introduction Two of the most common approaches currently in use for obtaining correlated energies and properties are configuration interaction (CI) and coupled cluster (CC) methods. Both methods use an expansion space of excited Slater determinants to represent the wavefunction; the difference between them originates with the particular form chosen for the expansion. In the CI case, a linear expansion is used, 'Icl) = (1 + C)0) = (1 + C1 + C2 + .)10) (3.1) where C, generates all possible ifold excitations from the reference determinant weighted by coefficients giving the importance of each determinant in the total CI wavefunction. In most cases, it is necessary for computational reasons, to truncate the expansion space somewhere short of including all possible excitations (up to nfold for an n electron system). This results in a variationally determined wavefunction (and energy), but sacrifices the size extensivity property. Size extensivity [2,3] is the requirement that the energy of a system scale properly with the number of electrons. It is most easily understood by considering an assembly of N identical noninteracting molecules. If each molecule, p, has an energy and reference function satisfying HIO0) = EP)O()) (3.2) 24 (the oneelectron functions of O1(p)) are assumed to be localized on that molecule), then we can write the reference energy and wavefunction of the entire assembly as N E(total) _S EP) p=l (3.3) N 10(total)) = 1 (p)) p=1 Since the molecules are noninteracting, we expect the correlated energy and wavefunc tions to be representable as similar sums and products, respectively. When a truncated CI expansion is applied to this system, however, this is not so. For an individual molecule, p, the coefficients C C) ,... C'P), with m equation, H(1 + C ) + CP + .) I0() = E(' (1 + C' + C2' + )0 O()} (3.4) and give a product wavefunction which includes excitations up to Nm due to products of excitations on different molecules. If the assembly were treated as a whole, however, only up to mfold excitations would appear in the wavefunction. Thus, they can never be equivalent. The energy of the assembly suffers similar problems. For example, a lattice of N noninteracting hydrogen molecules treated with CI restricted to double excitations (CID) has an energy which scales as /N instead of N, resulting in an error of 12.3% for 10 molecules and 21.1% for 20 [16]. Calculations by Ahlrichs and coworkers at the CID level show that the difference between supermolecule (the two molecules separated by a very large distance calculated together) and the sum of individual molecular energies are 7 kcal mol1 for BH3 + BH3 [17], and 15 kcal mol1 for CH3F + F [18]. 25 The well known property of exponentials, eAeB = eA+B (where [A,B]=0), suggests an alternative solution to the problem. If the individual molecule wavefunctions are of the form I(p)) = eTP)0(P)) (3.5) with P" being an expansion operator analogous to C, in Eq. 3.1, and satisfying He') O()) = E( e P) (p)) (3.6) then the product wavefunction is N p=1 o>P) N (3.7) .(ep ) = IO(p)) = e 0(=1 1 10(P)) p=l = ei' 'otal)o0(total)) which is identical to the total wavefunction resulting from a calculation of the entire assembly. This is the form of the wavefunction employed in coupled cluster theory. This simple physical picture demonstrates the need for the exponential wave operator, but it is the BruecknerGoldstone linked diagram theorem [1923], which CC methods satisfy [24], that guarantees size extensivity when this form is used. Thus when the system is composed of interacting fragments, the simple physical picture becomes clouded, but the size extensivity property is still guaranteed. Having motivated the use of coupled cluster methods and the exponential wave operator, in the remainder of this chapter, we will derive the working equations of 26 coupled cluster theory and other results which will be needed for the Fockspace CC methods discussed in the next chapter. 3.2 Basic Equations In order to expand the coupled cluster wavefunction in the space of Slater deter minants, we define a set of cluster operators which are capable of representing the contributions to the wavefunction due to determinants of particular levels of excitation, j= = tatai i,a T2 = a t attaj .>b (3.8) T3 = tibkaaabacaiajak t>j>k a>b>c where the cluster amplitudes, t:I are the scalar coefficients for each excited Slater determinant in the wavefunction. The cluster operators for each level of excitation are collected together into a single operator, T = Ti + T2 + T3 + (3.9) for convenience. To obtain the multiplicatively separable coupled cluster wavefunction, the wave operator (which produces the correlated wavefunction from the reference) is exponential in the cluster operators,  Qcc) = QOcl0) = eTO). (3.10) This wavefunction must satisfy the Schridinger equation He T0) = EcceO) (3.11) where E,, is the coupled cluster energy and H is the second quantized electrostatic Hamiltonian. In order to simplify later discussion, we change to the normal product form of H by subtracting Eo = (0HI0) from both sides of Eq. 3.11, (H (OIHO))eT0O) = (Ec Eo)eTIO) (3.12) HNeTIO) = EcorreTIO) where the correlation contribution to the energy has been designated E.,,. Finally, we multiply both sides of Eq. 3.12 by eT to obtain eTHNeTI) = Ecorr0) (3.13) Notice that the operator eTHNeT acts as an effective Hamiltonian operator, acting on the reference state to give the correlation energy; this quantity is often referred to as HN in postCC methods (i.e. those for which the derivation relies on having first solved for a CC energy and wavefunction). The effective Hamiltonian can be simplified by the use of a BakerCampbell Hausdorff expansion eTHeT = H + [H,T] + I[[H,T],T] + (3.14) 28 The commutators of strings of second quantized operators may be evaluated using the contraction theorem of Harris, Jeziorski and Monkhorst [25], with the result that the expansion terminates after the quadruplynested commutator term because H (Eq. 2.21) contains only four operators. This result is often denoted eTHeT= (HeT) (3.15) where the subscript c indicates that only socalled connected terms remain. The ter minology arises from the diagrammatic formulation of the equations often used in CC theory: surviving terms are those in which their diagrammatic representation forms a connected figure, from which no component may be removed except by breaking lines. The algebraic analog is that the term cannot be separated into a product of distinct sum mations. Disconnected terms are those which can be separated into two components without breaking a line, or algebraically the term may be cast as the product of two distinct summations. Use of Eq. 3.15 greatly simplifies reduction of the operator form of the coupled cluster equations to a spin orbital form suitable for numerical work. The equations which determine the CC energy and amplitudes are obtained from Eq. 3.13 by projecting from the left with the reference determinant and each class of excited determinants (single, double, triple, etc.). These projections provide equations which determine the energy, T,, Ts, T3, etc., respectively. (O1(HNeT) 10) = Ecorr (3.16) 'I (HeT) o) = 0 K (He T) > = 0 I("I 0) 0 (3.17) abc (He T)J = 0 abij 3.3 Reduction to Spin Orbital Equations The operator expression for the energy, Eq. 3.16, can be expanded to Ecorr= 0o(fN+ WN) + T +T2 + T + + T3 + T T2+ T + ) 0 (3.18) by replacing H, with its one and twobody components (see Eq. 2.21) and expanding the exponential. Since the Hamiltonian is in normal product form, the lead term in the exponential expansion is (0 HN 0) = (01(H (0 HIO))I0) = 0. The remaining cluster operators act on the reference to give an excited determinant and a product of amplitudes. As a onebody operator, fN can create at most a singly excited state acting to the left on the bra reference determinant, consequently nonzero contributions will result only from a singly excited state on the right. Similarly, the twobody WN operator requires a doubly excited state on the right. Thus, the equation simplifies to Ecorr = ( fNTI + (T+ WN( 2 +T ) O). (3.19) \ ~ ~ 2 \ *//c 30 It should be observed that the CC energy includes contributions due only to T, and T2 regardless of higher level cluster operators in T; higher clusters contribute only indirectly via their effects on T1 and Ts. Reduction of the amplitude equations (Eq. 3.17) to a more explicit operator form proceeds in a similar fashion. For example, the T, equation can be expanded to o=( (fn+W ) (l+Tl+T2+ T+ T3+TIT2+ T+*) ) (3.20) A nontrivial equation is obtained when the combined action of H, and the cluster expansion produces a single excitation I1) on the right. The onebody operator fN is capable of changing the excitation level by 1, 0, or +1, and W, by 2, 1, 0, +1, or +2. The fN +1 term and the WN +2 term cannot form connected contractions with T, so they can only contribute to the T1 and T2 equations, respectively, via the "1+" term in the exponential expansion. Because of the nature of the Hamiltonian, mentioned above, it should be evident that the equation for a particular cluster amplitude Tm will involve contractions of fN with products of cluster operators resulting in m and m+l excitations, and WN with cluster operator products yielding m1 to m+2 excitations. Equations for the first three cluster operators are thus o=( fI + T 2(+T2+T) +W( T+T+ ) +W0 T I + T2 + WT{ l + TT + TIT2 + i+T+T 2! T)3! 0= iab.f T2+ T3 + TT21 + T WN + TI + T2 + Ti + Ta + T1T2 1 4 ) \ + T3 +T4+ TiTa + T + TIT2T2 + T 0 (3.21) 3! 2!2 2! 4! C K= .. T Ta+bcTIT3 +1T2 + WN T2 + T3 + T T2 + T4 + T T3 + TIT2+ TiTT 2 S2 T3 + T2 + I T23 + Ts5 + TIT4 + T2T3 + T 3 T2 2! 2! 2! 3T C where the operator determined by the equation has been underlined. Conversion of the operator equations into the spin orbital equations useful for numerical work is done by considering in turn each of the operator products in the foregoing equations and using Wick's Theorem [15] to evaluate the possible contractions among the strings of second quantized operators contained in the Hamiltonian and cluster operators. This analysis can be carried out using a purely algebraic approach or a diagrammatic approach [26], but rather than adding greatly to the length and complexity of this overview we prefer to omit the detailed procedure here, and simply report the result. 32 Consider, for example, the first term of the CC energy as given by Eq. 3.19: (OIfNTlj0) = fpqt(O{ aq}{aia}O) p,q,i,a = fpqtfpisqa(OlO) (3.22) p,q,i,a = fiat? i,a Similarly for the second term, ( WN (T2 + 7) = E (pqlrs) ( It+ i p,q,r,s,i,j,a,b S0 aa{4 asar}({aaaaiaj} + {4ai}{a aj}) O) ~1 t4( 1 ab = (pq rs) 1 t + tit)Spi6jqi6s6ra(010) p,q,r,s,i,j,a,b = 2 (ij la b) t it+ j i,j,a,b (3.23) which gives the spin orbital equations for the coupled cluster energy: Ecorr = E fiat? + E (ij ab) (t + tit) (3.24) i,a i,j,a,b Spin orbital forms of the amplitude equations are obtained in the same way. Complete spin orbital equations for the CCSDTQ model, the most complete CC model implemented at the present time, may be found in the recent paper of Kucharski and Bartlett [27]. 33 3.4 Numerical Solution of the Coupled Cluster Equations The overall result, as can be seen in Eq. 3.21 is a set of coupled nonlinear equations for the amplitudes tjb... based on the MO integrals fpq and (pqllrs). These equations are usually solved by fixedpoint iteration starting with a guess for ta obtained from secondorder perturbation theory, and it has been found that acceleration methods based on conjugate direction techniques [2832] are very useful in obtaining faster convergence, and thus reducing the computational effort. Once the amplitudes are converged, Eq. 3.24 may be evaluated to obtain the CC energy. In many cases, the amplitudes and energy will be used, along with the MO integrals, as inputs to further calculations to determine the CC density matrix, gradients of the energy, or other results, for which it is often useful to work in terms of the CC effective Hamiltonian. 3.5 The Coupled Cluster Effective Hamiltonian When the coupled cluster calculation is carried out as the first step in a Fockspace CC calculation, it is useful to form the elements of HN= eTHNeT = (HNT) (3.25) explicitly from the MO integrals and CC amplitudes because the FSCC equations can be expressed in a compact and computationally efficient form using HN. Note that this operator is nonHermitian, so the braket symmetry of the MO integrals (i.e. (pqllrs) = (rsllpq)) is not carried over to the elements of HN:(pqllrs) 5 (rslpq). 34 Because HN contains at most four second quantized operators, HN will include contractions of HN with as many as four T clusters. Like HN, the effective Hamiltonian is described as having one and twobody components. In fact, while HN terminates at twobody terms, HN has up to n+2body terms in an n electron system. Just as the Hamiltonian is often written in terms of fN and WN, the components of the effective Hamiltonian are usually denoted fN, WN, IIIN,... The HN elements required for the Fockspace coupled cluster methods considered in the next chapter include all of the one and twobody terms and a few of the threebody terms. They are obtained in the same manner as the spin orbital equations for the CC energy and amplitudes, above, and the results are shown below. The first stage is to isolate the components of the expansion (H eT)= ( + +WN( I+Ti+T2+T + T3 + TT2 + T +)) (3.26) from which contractions can form. At the operator level, the equations are closely related to the CC amplitude equations (Eq. 3.21), but with two important differences. First, HN is defined based on converged cluster amplitudes, so Eq. 3.21 must be satisfied prior to forming the effective Hamiltonian elements. Secondly, like the Hamiltonian itself, the effective Hamiltonian can have n, ... 0, n excitation character for an nbody component. So while ( HN 10) = 0 is identical to the CC amplitude equations given above, the effective Hamiltonian itself contains numerous terms besides those required 35 to satisfy the equation. Thus the onebody component is fN =f I I +T +T2+ T12 (3.27) +WN TI ++ T2 T + T3 1 TT2 + 3 and the two and threebody components of HN, are WN =(fN {T2 + T3 + TT2 + W 1N{ + Ti + T2 + T + T + TiT2 (3.28) +T +T4+ TT + T + I TTr2 + TTP 3! 2! 2! 4 c and MIN =(IN( T+3 T+4 TT+ TT3 + T2 1 1 1 12 + WN T2 + T3 + TT2 + T4 + TiT3 + 7TT2 + T T2 (3.29) +TS5 + TIT4 + T2 T3 + Tf T3 + T1T + IT T2) Reduction to spin orbital form follows the same procedure used for the CC equations. In general, the nbody component of the effective Hamiltonian will have (n+1)2 different terms according to the distinct labelings of the indices as occupied or virtual orbitals. The one corresponding to the nfold excitation will, of course, be zero and for higher n, an increasing number of other terms will be zero. For example in IIIN it is impossible to form terms with 3, 2, or 1 excitation character, as well as half of the possible 0 excitation terms (the all occupied and all virtual labelings). Thus while the number and complexity of higherbody effective Hamiltonian terms does increase, it is not as explosive as that seen in manybody perturbation theory (MBPT) [26]. The spin orbital equations for HN are give in Tables 3.13.3. Onebody coupled cluster effective Hamiltonian elements in spin orbital form. level of operator is given in the first column. See text for explanation of fia + (imjlae)t fii + fifth (imllje)t + (imllfe)tf + (imllfe)tft fab fmbta + (amllbe)tm (nmllbe)te (nmllbe)tat, fai frmit + faet m + fneti fmett (mallie)te }(nml ie)te (nmllie)tte + t (amllf e)te + (amllfe)tft' +l(mnlef)t~ef + (mnllef)t?,eti (mnllef)tmt 1(mnllef)t tit (mnllef)tltatf Table 3.2. Twobody coupled cluster effective Hamiltonian elements in spin orbital form. Excitation level of operator is given in the first column. See text for explanation of notation. 2 (ijflab) = 1 (ijIlka) = 1 (ai lbc) = o (ijl kl)= 0 (abllcd) = 0 (ia ljb) = (ijllab) (ijilka) + (ijilea)tk (aillbc) (millbc)it (ijilkl) + ()PP(kl)(ijllke)ty + (ijllef)tP +()PP(kl)(ijlle f)tet (abllcd) ()PP(ab)(amI cd)tb + (mn Icd)ta +()PP(ab)(mnl cd)tatb (ial\jb) (imrnjb)ta + (ia[[eb)t (imlleb)tf (im\\eb tt Table 3.1. Excitation notation. 1 fai = 0 fij = 0 fab = +1 fia 37 Table 3.2. (Continued) Twobody coupled cluster effective Hamiltonian element, spin orbital form. Excitation level of operator is given in the first column. See text explanation of notation. +1 (ialjk) = (iaIIjk) + fiet, (imlljk)ta + ()'P(jk)(iallje)t' +()Pp(jk)(imlje)ta ()PP(jk)(im Ije)tekt +()PP(jk)(iaIIef)ttf + (ialIef)tIf + mi( ief)\tefa 3)kk22m 2 JI mjk +(mizllef) tftf l(iml fe)t ta ()PP(jk)(imI fe)tL ()PP(jk)(imI fe)tjt tm +1 (ac(ab(l) abllci) ftci + (abIce)t ()PP(ab)(amllci)t +(mn lci)ta + ()PP(ab)(mnl Ici)tt b +()"P(ab)(amIce)t ()PP(ab)(am Ice)titb (mnI ec)t (mnl ec)tert + ()PP(ab)(nm Ice)tattb +(nmlce)t mt + ()pP(ab)(nm lce)ttatqt +2 (abllij)= (ab lij) ()PP(ij)f .. tab + ()PP(ab)fbeta +fmetj )P(i)fme ()PP(ab)fmetiftb ()PP(ab)(mbl ij)ta + ()PP(ij)(ab lej)t + (mnllij)ta + (ab ef)tf ()PP(ij)(mb lie)ta, + ()PP(ab)(mnjjij)tatb +()PP(ij)(ablef)t/tf ()PP(ij)(mbllie)t at ()PP(ij)(nmllje)tb + ()PP(ab)(bmr fe)ltf ()PP(ij)(mnei)tet + ()PP(ab)(ma ef)te;t ()PP(ij)(ab)(nmIIie)ttebi + ()PP(ij)(ab)(amllfe)tfteb + l()PP(ij)(mnllej)t1tab ()PP(ab)(mbllef) t/ Inn 1()PP(ab)(mbIIef)tit j +(mnlei f)abefn + (mnlle f)taqbtl ()P(ij)(mnlle f) ttafb ()PP(ab)(mn ef)tn + ()PP(ij)(mnllef)tift? + (mniIe f)t~t ()PP(ab)(mnlle f)tetn s in for 38 Table 3.2. (Continued) Twobody coupled cluster effective Hamiltonian elements in spin orbital form. Excitation level of operator is given in the first column. See text for explanation of notation. ()'P(ij) mn \ ef)t? ()PP(ij)(ab)(mni efp t tfb + ()PP(ij)(mn e f)t \tt + ()PP(ab)(mnlef)t ttt ()'P(ij)(nm Ife)tt'tf ()'P(ab)(nml fe)ttfato +()PP(ij)(mn[ ef)t t tb The summation convention (all repeated indices are summed) has been used to save space and simplify presentation of these formulae. The permutation symbols, ()PP(), are used in the same fashion as Lee et al. [33] and Kucharski and Bartlett [27]. A permutation of two indices, e.g. ()PP(ij), denotes two terms. One is the identity permutation, with signature p=+l, while the second is the term with indices i and j interchanged. A permutation of three indices, e.g. ()PP(ijk), represents six terms with signature p=+l for cyclic permutations and p=l for acyclic permutations. The symbol ()PP(i/jk) indicates that permutations i j and i <+ k should be taken, but not j < k; the signature is the same as for a binary permutation. Finally, a composite symbol, such as ()PP(i/jk,ab), indicates that the direct product of the two sets of permutations (i/jk and ab) should be taken. 3.6 Truncation of the Cluster Operator It is worth noting that we have made no approximations in deriving the equations given above, and the energy, wavefunction, and effective Hamiltonian thus determined would be exact (for the chosen atomic orbital basis). Unfortunately, the number of 39 Table 3.3. Threebody coupled cluster effective Hamiltonian elements needed for Fockspace (1,0), (0,1), and (1,1) sectors in spin orbital form. (The full IIIN includes an element at excitation levels +2, and +3 in addition to those shown here.) Excitation level of operator is given in the first column. See text for explanation of notation. 0 (aijllklc) = (ijlec)ta 0 (abilljcd) =(milcd)t! +1 (ijajjklm) ()PP(ml/k)(jiIlek)t" + (jilef)ti, +()PP(k/lm)a(jilef)tae tf +1 abclldei) = ()pp(b/ac)(mblIde)ta + (mnllde)tab +()PP(b/ac)(mnl de)tia tb +1 (iablljck) = (miiec)tb + () P(jk)(mi\lec)tb t +()P(ab)(mizlec)tq tb ()PP(jk)(milljc)ti +()PP(ab)(iallec)tb +2 (iab jkl) = ()PP(jk/1)(iml Ijk)t~t + ()P(j/kl,ab)(ia re)t prn/jl ab)P 1,mai) (a Ia eb e ab () P(k/jl, a)(mi Ike)tmt ()PP(jkl)(im lek)tjtfi +() P(jk/l, ab)(im lfe)tte I ()PP(jk/l)(im\ ef )teta +()PP(j/kl,ab)(iaIef)tetf + f, +()pp (j/kl)(irlmje)tea + ()PP(ab)(ia Ief)te +()PP(j/kl)(im Ife)tfiea 1()P(ab)(ZmIIef)taj i 2\( eftkjltm +()PP(jk/l, ab)(i lfe)tate ()PP(jk/l)(im llefteft ()PP(k/jl,ab)(miI ef)}tta tb ()PP(jkl)(im\ fe)tf t t +_(mlIIefefab +e fab 1(rnllef)tmikL + (millef)tit ikl amplitudes to be determined is O(N!) for N AO basis functions. This rapidly becomes too expensive to compute using current and foreseeable computer facilities, so it is necessary to make approximations. The obvious choice is to limit the cluster operators included in T to lower levels of excitation because the number of amplitudes grows by a factor 40 of n(Nn) for each additional excitation. Because of the exponential nature of the CC wave operator, some higher excitation effects are incorporated by products of lowerlevel excitations (i.e. T, T2 is a triple excitation). In fact in most cases, such products dominate the wavefunction at a given level of excitation [34]. The single excitations alone cannot have any contribution to the energy, so the first reasonable approximation is T= T2, known as the CCD (coupled clusterdoubles) model. With the doubles present, the singles can also contribute, so T=TI+T2, the CCSD model, would be next on the hierarchy. Increasingly sophisticated methods are obtained by adding the next level cluster operator: CCSDT, CCSDTQ, .. . To obtain approximate versions of the equations given in this chapter, simply consider all cluster operators not present in the desired model to be zero and eliminate them from the equation. For example, CCSD equations would have Ts=T4=T5=0, and all terms involving these operators would be eliminated. Observe, however, that products of TI and T2 giving overall triple, quadruple, and higher excitations would still appear. This is a consequence of the exponential expansion and is the reason why a given truncation of coupled cluster theory recovers more of the correlation energy than the same truncation of operator manifold in the linear CI expansion (Eq. 3.1). In practice, the CCSD model has proven widely applicable, and the more expensive CCSDT model, and approximations which involve perturbative inclusion of the triple excitation effects, provide higher accuracy and have been shown capable of handling some problems that would otherwise require multireference techniques [34]. The CCSDT model, or perturbative approximations, are currently the practical limit for the majority 41 of interesting systems. Further, the computer code used in this work does not presently include the full CCSDT model, so calculation of HN is limited to the CCSD level. CHAPTER 4 FOCK SPACE COUPLED CLUSTER THEORY 4.1 Introduction In the preceding chapter, the CC equations were derived in a spin orbital framework. This approach allows the use of wide range of reference functions (RHF, UHF, ROHF, QRHF, etc.) thus allowing application to both closed and open shell systems with relative ease. For open shell molecules, these spin orbital based CC methods may be characterized as eigenfunctions of the S2 operator in a projected sense, (0 S2 cc,) = s(s + 1), for ROHF and QRHF references, or not at all, for UHF references. Although they have been very successful [34], failure to be a rigorous eigenfunction, S2 cc) = s(s + 1) cc), makes this approach less desirable to some workers in the field. The Fock space coupled cluster method grew out of a search for an approach to open shell systems which would be a rigorous eigenfunction of S' and be capable of handling the multiple reference determinants needed to represent many open shell states. There are two primary approaches to the treatment of multireference systems. The Hilbert space approach, also referred to as statespecific, or valanceuniversal, treats a manifold of states with a constant number of electrons, and each Hilbert space CC calculation focuses on a particular state. The Fock space (FS) approach treats states with different numbers of electrons (mathematically, a Fock space is a sum of Hilbert spaces with different numbers of electrons), simultaneously calculating several different states. 43 The FSCC method is, at present, the more developed of the two, and offers the useful feature, mentioned in Chapter 1, that it can provide data for more than one state in a single calculation. In this chapter, we will develop the FSCC equations in more detail, obtaining spin orbital equations for IPs at the FSCCSDT level. The FSCC approach will be compared with the familiar, and (reasonably) well understood single reference CC approach to point up the inequivalence between FSCC and CC calculations at the same level of truncation and note the increased importance of higher excitations (triples in particular) in the FSCC approach. Finally, a number of approximations to the full inclusion of triple excitations are proposed in order to include some triple excitation effects but reduce the computational effort required. 4.2 Structure of FSCC A multireference description of a correlated wavefunction means that the determi nantal expansion is partitioned into two components. The model space includes the determinants desired as references, and the complementary space, which is orthogonal to the model space, contains excited determinants resulting from each of the model space functions. The same division holds trivially for single reference CC, with the one reference determinant being the model space and all excited determinants being in the complementary space. These two spaces may be represented by projection operators, P, and Q, for the model and complementary spaces respectively, defined from the functions they include: P= E I0p)(OPI p (4.1) Q = )(Cq = 1 P q which are orthogonal PQ = QP = 0 (4.2) and idempotent PP = P, QQ = Q. (4.3) For convenience, we choose a single determinant oneelectron reference function, 10), to provide an unambiguous identification of occupied and virtual orbitals. The idea of a particlehole rank, denoted (m,n), is used to indicate the states related to 10) by the addition of m electrons (particles in the otherwise empty virtual orbitals) and the removal of n electrons (holes in the occupied orbitals). Thus singly ionized states are in the (0,1) sector, electronattached states are in the (1,0) sector, and the (1,1) sector includes states obtained by exciting a single electron from 10). Higher particlehole ranks give multiply ionized or multiply attached states, excitations involving more than one electron, and combinations such as "shakeups" (combination of detachment of an electron and excitation of another). The idea of the model space and its orthogonal complement can also be specialized to a particular particlehole rank, giving projection operators p(m,n) and Q(m,), respectively. 45 virtual inactive orbitals (particles) f, active occupied orbitals T1 (holes) T1 inactive Figure 4.1. Division of orbital space into active and inactive groups. In order to define the model space, the orbital space of the reference determinant, 10), is divided into two classes. Active orbitals are those which change occupation in the various model functions, while inactive orbitals maintain the same occupancy. Since the reference function has already been divided into virtual and occupied orbitals, the result is a four different categories of orbitals: active occupied, active virtual, inactive occupied, and inactive virtual. These divisions are shown in Figure 4.1. Once the active orbitals are chosen, the model space for each particlehole rank is determined all appropriate occupations of the active space are included. For example in Figure 4.1, the (0,1) sector (one electron removed) model space includes six determinants and the (1,0) sector has four because there are four ways to add a single electron to the active virtual space; there are 12 spinconserving single excitations within the active space, hence 12 model functions in the (1,1) sector. It should be noted that this scheme provides for exactly one model function, 10), in the (0,0) sector, and the calculation in this sector is exactly the single reference CC method discussed in the preceding chapter. 46 Model spaces can be divided into two types, complete or incomplete [35]. A complete model space is one that contains all possible occupations of the active orbitals, while in an incomplete model space (IMS), some occupations of the active orbitals belong to the Q space. The pure ionization and pure electron addition sectors, (0,n) and (m,0) sectors are always complete, but any mixed sector is, in general, incomplete. For example, the (1,1) sector has all single excitations within the model space in P, but places all double, triple, and higher excitations into the Q space. The (1,1) sector is, however, part of a special class of IMS known as the quasicomplete model space [35,36]. In this case, the active orbitals are divided into disjoint sets, each with a fixed number of electrons. If each of these subsets is complete within itself, the entire model space is quasicomplete. The FSCC wavefunction for a state is obtained from a wave operator, Q, operating on a linear combination of model space functions: I'FSCC) = 'lk). (4.4) where k Ckii). (4.5) The wave operator used in FSCC is an exponential of cluster operators, just as in the single reference theory, but with a normal product used to avoid contractions within the cluster operators itself (this is not necessary in the single reference formalism since it is impossible for those cluster operators to contract among themselves), S= {eS. (4.6) 47 The cluster operator itself, however, has a somewhat different structure from the single reference case: S = S(mn) = T(kl) (4.7) k=0 ...m l=0...n where (m, n) is the highest sector relevant to the problem, and the T(k'l) for a given sector includes single, double, and higher excitation cluster operators for that sector: T(k,') = T(k,') + T (kk) + T ..l) (4.8) Such a wave operator is valanceuniversal because it applies to all sectors from (0,0) to (m,n). As mentioned above, the (0,0) sector is a CC calculation with 10) as the reference, so T(o,0) is nothing more than the familiar CC cluster operator described in the previous chapter. Cluster operators for other sectors are different from the (0,0) ones because they apply to spaces with different numbers of electrons, or at least electrons in different places than in 10) (i.e. active virtual may contain electrons and active occupied may receive electrons). In the (0,1) sector, for example, one of the virtual orbital labels on a cluster amplitude becomes an active occupied, corresponding to the fact that in this sector, an electron has been removed from that group, making it possible for the empty orbital to receive an electron. The differences by sector in the cluster operators is depicted in graphical form by Rittby and Bartlett [37] for a number of sectors, and some implications of these differences are discussed in Section 4.3, below. To derive the FSCC equations themselves, we begin from the Schridinger equation for the FSCC wavefunction, (4.9) HS2IkIk) = Ekf~lx~k) 48 This equation can be rewritten in terms of the model space functions by defining a matrix e' having the E, as eigenvalues and C as eigenvectors, and using the model space projector, P, to effectively encompass the entire model space, HfQP = fl'P. (4.10) Using the fact that the (0,0) sector cluster operator commutes with the others [38], we can decouple the (0,0) sector calculation from the rest of the FSCC calculation. First, we extract T(O,0) from the wave operator, S= eToofi = eo'o) { eS(m.n), (4.11) where g(m,n) = T(k,), (4.12) k=O...m,1=0..,n k+14#0 and multiply both sides of the Schr6dinger equation by exp (T(,0)) to obtain eT~oa) HNeoo) iP = i'P (4.13) Recalling that the coupled cluster effective Hamiltonian, derived in the preceding chapter, is HN = e oo He 0 (0HI0) (4.14) = eo'O)He T(O Ecc,o) where Ecc,0o) is the CC energy of the reference function 10), and defining a new matrix, e = e' Ecc,j), which, when diagonalized, represents energy differences between Ecc,jo) and the final state, we can write the FSCC Schridinger equation in the simplified form: fHNIP = eiP. (4.15) Projecting Eq. 4.15 on the left with PO1, we can isolate e as the FSCC effective Hamiltonian, PHN,effP = Pf HNP = PeP, (4.16) which, upon diagonalization, gives the statetostate energy differences we seek. Note that this result could have been obtained by requiring intermediate normalization, PQP = 1, but Mukherjee [39] has shown that for the general case of an incomplete model space, intermediate normalization is incompatible with size extensivity. Although this will not cause problems for the (0,1) sector, our ultimate focus, we prefer to derive generally applicable equations as far as possible. To determine the cluster operators, we project Eq. 4.15 on the left with the complementary space determinants, giving QHNQP QQPHN,effP = 0, (4.17) which is often referred to as the Bloch equation. Mukherjee and coworkers have shown that for complete and quasicomplete model spaces, HN,eff in Eq. 4.17 is connected and the energies obtained by diagonalizing it are size extensive [3944]. Haque and Mukherjee [38,45] have also shown that there is a subsystem embedding condition (SEC) in the FSCC equations which allows Eq. 4.17 to be solved separately for each sector in a hierarchical fashion. As already observed, the (0,0) sector reference problem decouples from the rest of the Fockspace problem. The (1,0) and (0,1) sectors can be solved with only the (0,0) results they do not depend on any higher sectors. Given the (0,0), (1,0), and (0,1) sector amplitudes and HNeff, the (2,0), (1,1), and (0,2) 50 sectors can be solved once again, no information from higher sectors is required. In this fashion, the solution of the desired (m,n) sector is built up from lower sectors, rather than requiring the solution to a set of coupled equations for all sectors. The SEC is also very useful in the development of computer programs for FSCC because the implementation of a new sector should not require any changes to the lower sectors. 4.3 Comparison with Single Reference CC Theory Before going on to derive the more specific equations for Fockspace methods, it is worth taking a look at how the FSCC approach differs from the more familiar single reference formulation. As has been noted, the Fock space coupled cluster amplitudes differ from (0,0) sector amplitudes in that they involve excitations between the occupied or virtual orbitals and orbitals of the same type in the active space. In the reference function, 10), such amplitudes would be zero because active occupied orbitals would be fully occupied, so no excitation could deposit an electron there, and active virtual would be completely unoccupied, providing no electron to excite to another virtual. In the model space, however, there will be determinants with holes in the active occupied space and particles in the active virtual space, so such amplitudes may be nonzero. This fact changes the character of FSCC single, double, etc., excitation amplitudes from their familiar (0,0) sector interpretations and means that a given truncation of FSCC will not provide the same results for a given state as the same truncation of the (0,0) sector equations applied to the model space function. T I T I T T I I T T J T I T I T I 10) i 0,1)) ( 0,01)) 00'1 ) i 0,1) Figure 4.2. Prototypical reference functions and model space determinants. Active orbitals are shaded. T I T I T I T I T I Figure 4.3. Determinants appearing both in the single excitation space of CCSD calcu lations on the four model functions and in the FSCCSD singles space. To make this idea more concrete, consider a FSCCSD IP calculation compared to CCSD calculations based on the individual determinants in the model space. Using an RHF reference function with three occupied orbitals and one virtual orbital, if the two highest occupied orbitals are made ac tive, we have four determinants in the model space, as shown in Figure 4.2. A CCSD calculation on any of the model space functions has seven determinants in the single excitation space, and the union of the four single excitation spaces contains 20 distinct determinants. FSCCSD, on the other hand, restricts single excitations to be from an inactive occupied orbital to an active occupied, which allows only four amplitudes and yields two distinct determinants, shown in Figure 4.3. The remaining 18 T T T T T T T TI I T I t T I T I T T I t24 34 t24 t24 24 23 32 21 2 t34 t34 t34 31 32 31 Figure 4.4. Examples of determinants arising from spectator double excitations in FSCC, corresponding to single excitations in the CC model. Shown are all determinants generated from the 10(o1)) model function with the relevant amplitudes from T(o') indicated below each one. Orbitals are numbered from lowest energy (bottom) to highest (top); beta spin electrons are indicated by a bar over the orbital index. CCSD determinants are generated in the Fock space calculation by a class of excitations referred to as spectator doubles. FS double amplitudes include one excitation from an occupied orbital to an active occupied. There is no restriction on the initial orbital, so it may be active as well in fact it may be the same as the final orbital. This null excitation, which is necessary to distinguish between excitations connecting different model functions to the same final determinant [46,45], makes spectator doubles effectively single excitations in the usual coupled cluster sense. Figure 4.4 shows the spectator doubles arising from model function 1(o1')), along with the corresponding amplitudes. When considered at truncation to only double excitations (CCD and FSCCD), the Fock space approach is going to include some (but not all) effects of single excitations  a potential advantage. When higher cluster operators are included, however, FSCC is at an unequivocal disadvantage. Compare, for the same model problem, the double ex citations obtained from the CC and FSCC approaches. FSCCSD generates 22 "true" tT T T 1 T 1 T TT T I T I T I T I T T 1T 1 T T t244 t344 t244 t244 3344 t344 232 322 212 231 321 311 t344 t244 t244 t244 312 231 221 211 t244 211 Figure 4.5. Examples of determinants arising from spectator triple excitations in FSCC, corresponding to double excitations in the CC model. Shown are all determinants generated from the 4(0,1)) model function with the relevant amplitudes from T(o') indicated below each one. Orbitals are numbered from lowest energy (bottom) to highest (top); beta spin electrons are indicated by a bar over the orbital index. doubly excited determinants, while CCSD calculations on the individual model func tions yield 18 additional determinants. The "extra" CCSD determinants are created in FSCC by spectator triple excitations, a few of which are displayed in Figure 4.5. From these examples, it should be evident that for a FSCC IP calculation to match the expansion space of a single reference CC of a given truncation on one of the model functions, the FSCC calculation must include spectator contributions from the next higher cluster operator. The same would be true of EA calculations, and in higher sectors, more indices of a given amplitude terminate in the active space (p=m+n for the (m,n) sector), so it would be necessary to include spectators from the ph higher cluster operator to obtain an expansion space equivalent to that of the single reference CC. 4.4 Basic Equations From this point, we will focus on the (0,1) sector, which describes oneelectron ionizations with respect to the reference state and are of particular interest in this work. From Eq. 4.17, the detailed operator equations which determine each cluster operator are easily determined by expanding the exponential wave operator, as in the previous chapter. Note that in FSCC, products of cluster operators in the exponential expansion contribute to other sectors, making the (0,1) equations a good deal simpler than their CC counterparts. P(o,) fN,effP(O'l) = P(o1) f N(1 + T0',1) + T~1)) (4.18) + WNT (TO,) + T('1) + IIINT(,'1)}P(o,1) o = Q( ') (1 + T,(O) + T('1)) + WN (T(O1) T(0,1)) (4.19) +I7INTo01,) T(o'l)f1) cP(ol,) 0 = QO')(fN (TO'1) T(O"'))+ N ( + T 1 T ) T() + T3(,)) (4.20) +IIN (T(0,o1) + T(O',)) T(Ol)f p(o,1) (llNU2 +3 ) 2 NJVeff 1, (0,1) 0,1) 0,1) O = Io3 l)fNT3l) + WN (T2~1) )+ T(0'1)) + IIN (i + T(l) + T( ) + T3(01)) (4.21) T(o,1) fr(O,) p(o,) 3 N,effJc 4.5 Reduction to Spin Orbital Equations Equations 4.184.21 are transformed to spin orbital form following the same proce dures used in the previous chapter for the CC amplitude equations and HN. The resulting equations for the FSCCSDT method for ionization potentials are given in Eq. 4.224.25, below. The notation is the same as that used in the previous chapter with a few additions. The Fockspace cluster amplitudes are called r to avoid any possibility of confusion with the (0,0) sector amplitudes, t. A single dot on an orbital index, ri, indicates restriction to inactive orbitals, and a double dot, ih, indicates restriction to active orbitals (this choice of notation is meant to remind readers familiar with the diagrammatic notation of the "circled arrow" used to denote inactive lines, and the double arrow for active lines). (HNeff)k = f r + fme mn \e /r + mn\ ef (4.22) 0= fj i+f mnke)rm + mnlef)r +{(HN,eff)iH (4.23) fkrm + k 4h+ fm' m 2 4kmn k 0 = (iaIki) + fea (1)P(ki)fmirT (rhaIki)rm + 2(l),P(ki)(mnLki)rC 7 (ma rke)e i 2(manfikie)re. + fmer + (1)PP(ki)(mnf]ke)ri + I(mnalef)rkmi + 4 (HN,eff)i (4.24) 0 = ilabI\kij) ()PP(k/ij)(mb ij))rka + ()'P(ki/j)(ablej})T, + ()PP(ab)(man[Ikij)rTn ()PP(ki/j)(mibi)kie)rEn (rhabIkij)r} + () P(ab)fberkij ()PP(ki/j)fm,7jTki + () P(ki/j)(mnIrki)r, ()'P(k/ij,ab)(ma\lke)rL +ki(HN,eff)l (4.25) 4.6 Approximations to the Full FSCCSDT Model In the single reference CC world, the full CCSDT model has seen limited use for a number of reasons * Complexity of the method the triples equations have many more terms than single and double excitations, thus requiring substantially more effort to code them into a computer program. * Primary and secondary memory (i.e. physical memory and disk space) requirements for storage and manipulation of the cluster amplitudes, which are significantly more 57 numerous than the double excitation amplitudes a factor of n(Nn) in the language of section 3.6. Time constraints, since many more amplitudes must be determined, the time required increases in a similar fashion. As a result, a number of approximations to the full CCSDT method have been proposed and have proven quite effective at incorporating important effects of triple excitations while at the same time minimizing the above concerns [33,4752]. These concerns might also be raised with respect to the full FSCCSDT method, so in the same spirit, we propose several approximations of FSCCSDT as a compromise between the cost and performance of the full model. There are two different approaches to approximating FSCC triples which may be considered. One, in analogy to the approximate triples treatments of single reference CC, involves truncations and perturbative approximations to the triples equations. The second method is to restrict the triples amplitudes which are evaluated to something less than the full set. This approach, which has not been employed in CC theory, is suggested by the earlier observation that spectator triples must be included to, in some sense, duplicate the expansion space of a single reference CCSD calculation. The two approaches are not exactly orthogonal, but can be considered, by and large, independently. In single reference CC theory, there are no selfevident subsets of the triples amplitudes, so this approach has not been tried (or at least not reported in the literature). In FSCC, however, there are two obvious subsets which might be considered. The smallest set would be the spectator amplitudes, as suggested by the discussion in the 58 previous section. A larger set would be to include all activeactive amplitudes. This set, which includes the spectators, allows excitations from any active orbital into any other (spectators involve an active orbital scattering into itself, and the full set of triples allows scattering from inactive orbitals to occupied orbitals as well as active to active). As we saw in the previous section, it is necessary to introduce spectator triples in the Fock space approach in order to match the expansion space of a CCSD treatment of each of the model functions in turn. On the other hand, since FSCC describes all states of interest with a single set of amplitudes (with diagonalization of the effective Hamiltonian distinguishing among individual states), it can be argued that it is impossible for it to be equivalent to a separate CC calculation for each model function, and the differences in the expansion space between the two approaches should not be a concern. Indeed, a number of calculations in the literature [53,54,43,55, for example] support the idea that FSCCSD is wellbalanced and the results obtained are numerically on a par with other methods. The CCSD method has been widely applied, and its behavior is reasonably well understood certainly there is far more experience with CCSD than with any FSCC method and this, more than any special quality of CCSD underlies the desire for a FSCC method that can mimic it. Only a critical examination of the proposed spectator and activeactive approximations for a wide range of systems can determine if they will improve the results over FSCCSD or destroy the balance of the FSCC method. The second approach to reducing the cost of a FSCCSDT calculation is to approximate or ignore the most expensive terms to evaluate. In this sense, the simplest possible approximation would be to evaluate the triple excitation contribution perturbatively from 59 converged FSCCSD amplitudes. The single reference analog to this approach is the CCSD+T(CCSD) method, in which the triple excitation contribution is correct through fourth order the first order at which it appears. For the FSCC perturbation expansion, we follow the choice of Pal, Rittby, and Bartlett [56] which was also implicit in the work of Haque and Kaldor [57], of counting orders in terms of H, the reference state Hamiltonian. Alternatives, such as using H would lead to inconsistencies between the (0,0) sector definitions and those of higher sectors. For example, f and W are first order in H, while the higher components, III,... are at least second order. Counting orders in H would place all of these terms on the same level, making a perturbation expansion of Eq. 4.21 useless. Although in both CC and FSCC the triple excitation amplitudes first appear in second order of perturbation theory, in the CC formalism excitations beyond singles and doubles contribute only indirectly to the energy via their contributions to the singles and doubles amplitudes, and thus triples first appear in the energy in fourth order. In FSCC, all clusters contribute directly to HNeff and thus the energy, so the triples first manifest themselves in third order. Inclusion of the third order triples effects based on converged FSCCSD amplitudes was first proposed by Haque and Kaldor [57], and called CCSD+T. Later, Pal, Rittby, and Bartlett [56] discuss the third order model using the name MRCCSD+T(3) and propose several other approximations. The MRCCSD+T*(3) model attempts to include some higher order terms by replacing the f and W required in T(3) with their counterparts from the effective Hamiltonian, f and W. The nature of this method, which will be referred to as FSCCSD+T(3) in this work, can be seen from the spin orbital 60 FSCCSDT equations given above. The terms which contribute to the T(0'1) amplitudes in second order are those on the first row of Eq. 4.25, with appropriate restrictions on the HN elements: the WN elements of the second and third terms are simply the bare twoelectron integrals in first order, and only two terms of (lab lkij) appear at this level (the first two terms of (iabfljkl) in Table 3.3). The second order triples amplitudes in turn contribute to HN,eff (and thus the energy) in third order thought the final term of HN,eff in Eq. 4.22. The extension of this approach to higher orders of perturbation theory is obvious. A fourthorder method, FSCCSD+T(4), would be obtained by the following steps: 1. Secondorder T o1) amplitudes, as evaluated above, give thirdorder contributions when inserted in the T,(o') and T(O'1) amplitude equations. This correction to the T"'1) and To'1) amplitudes can be directly contracted into HN,eff, giving a fourth order correction to the energy. 2. The thirdorder T(o1) must be evaluated. This requires the fourth, fifth and sixth terms from Eq. 4.25, and additional terms from (labllkij) (terms 310 of (iab\jjkl) in Table 3.3). Note that even if the (0,0) reference calculation does not include triple excitations, approximate T 0'0) amplitudes can be formed at second order in pertur bation theory, and therefore must be included in the general case. Computational experiments will be required to determine the effect of approximate T(0'0) in Ny. Of course such a perturbation expansion could go on for a very long time, but in practice, better results are likely to be obtained (based on SRCC experience) by iterating an approximate triples equation as part of the FSCCSD calculation. The 61 simplest approximation for which this makes any sense is T(4), described in the preceding paragraph it is the first perturbative method for which T,("') T(o'1) and T3(0')  T0",) contributions occur. The procedure is simple: In each iteration of the FSCC calculation, after the single and double excitation amplitudes have been evaluated, the T(4) procedure is used to calculate the contribution of triples excitations to HN,eff and the singles and doubles amplitudes. Instead of this being the end of the calculation, the new ( 1) and T'1) amplitudes, and HN,eff are fed back into the next cycle of iteration. Note that there are no T(o'1) T3(o') contributions in this method, so the triple excitation amplitudes do not need to be kept from one iteration to the next, and also avoids the most expensive part of the T(''1) equation to evaluate. This approach is the same as that taken by Urban et al. [48] in approximating single reference CCSDT. In the hopes of maintaining some consistency in nomenclature between SRCC and Fockspace approaches, this method should be called FSCCSDT1. In the single reference case, the idea has been extended in a series of additional iterative approximations, CCSDT2 [48,50], CCSDT3 [48,50], and CCSDT4 [51], which include more and more of the terms due to triples, but continuing to exclude T3(') To'O), thus avoiding the most expensive terms both computationally and for storage space. In the FSCC case, the T(o'1) equation has a simpler structure, and only T3o'1) T(o'1) terms remain after the FSCCSDT1 approximation. Over the years, experience with various SRCC approximate triples methods has led to refinements in several of them. One of the main observations is that in cluding selected higher order terms in a perturbative or iterative method can pro 62 duce better results. For example, Watts, Gauss, and Bartlett [58] have recently intro duced the CCSD[T] method, which is a modification of CCSD+T(CCSD), the sim plest SRCC perturbative triples method. They differ in the inclusion of a fifth order term and in the requirement that nonHartreeFock (i.e. QRHF or ROHF) ref erences be treated properly. This requires that terms involving offdiagonal f ele ments be considered since they are zero for a HartreeFock reference, but nonzero in the nonHF case. Transforming the reference functions to semicanonical orbitals [59] has the effect of diagonalizing the occupiedoccupied and virtualvirtual blocks of the Fock matrix, so that terms involving fij or fab can be avoided. This leaves occupiedvirtual Fock matrix elements, fai to be accounted for, which involves shift ing these terms down one order of perturbation theory, from first to zeroth or der. In the Fockspace coupled cluster framework, however, the situation is somewhat different. The transformation to semicanonical orbitals diagonalizes blocks of f, not fN. Even with such a transformation, the leading terms in fij and fab are first order in the non HF case. This would result in FT~'1) T(o' 1) contributions, one of the most expensive terms in FSCCSDT, beginning with the thirdorder amplitudes (FSCCSD+T(4)). Thus as far as approximations to FSCCSDT go, such a requirement would be reasonable only for the most basic thirdorder method, FSCCSD+T(3). This would add the FT(0'0) term of (labllkij) (see Table 3.3) and the FT(o'O) terms of (mbhij) and (abllej) (see Table 3.2), and would require evaluation of the FT(o') T(o1 ")H H term of Eq. 4.24 since fme is now zeroth order. 63 There is one further matter regarding approximation to the FSCCSDT equations. Mattie [60] has proven that for the (0,1) sector of Fockspace, the values of the roots obtained by diagonalizing HN,eff are independent of the active space used for the calculation. In other words, if two orbitals, a and b are taken as active, the resulting ionization potentials would be identical to those obtained from separate calculations with either a or b alone active. This is a very useful result because it means that we do not have to worry about choosing the "right" active space for a given calculation in order to get good results. However the proof of this invariance rests on the FSCC amplitudes and HN,eff satisfying the Fockspace Bloch equation 4.17. None of the approximate FSCCSDT methods described above satisfy Eq. 4.17, so the invariance is lost. It remains to be seen how large this effect will be, or if there is an alternative route to the invariance condition that will work with approximate solutions. As a start towards understanding the FSCC triple excitation correction as well as we now understand the single reference triples correction, the remainder of this work will focus on the thirdorder triples approximations, FSCCSD+T(3) and FSCCSD+T*(3). The implementation of these methods for both closed and openshell reference functions is described in some detail in the Appendix, including the final, spinintegrated equations from which the computer code was written. In the following chapters, these methods are applied to a variety of molecules in order to assess their performance. CHAPTER 5 A CRITICAL ASSESSMENT OF FSCCSD AND THE THIRDORDER TRIPLES APPROXIMATIONS 5.1 Introduction Before applying a new theoretical method to otherwise unknown molecular systems, it is important to get a feel for its behavior strengths, weaknesses, and overall quality  by comparing the results obtained with new method to other theoretical and experimental data for similar types of problems. This "calibration" of a method is a continuous and dynamic process each time the method is used for a new problem, more information about its performance is gained, thus modifying the assessment of its accuracy and applicability. Choosing the appropriate standard of comparison for a new method is in most cases not a simple matter. Experimental results are an obvious choice, but there some potential problems. Experimental corresponding to readily computed results are not always available or may not be known accurately. For example, FSCC most readily provides vertical ionization potentials, but only adiabatic may be available from experiment. For potential energy surfaces geometries, relative energies, and the like there is usually very limited experimental data. Structures may be known for minima on the PES, but experiments which can probe transition states are quite new, and have not yet produced a sizeable body of data on different systems. The vibrational spectra of a given species 65 may be only partially known, and in many experiments, assignment of spectral lines to particular species can be a very hard problem. It is also possible for a calculation to obtain the right (experimental) answer but for the wrong reasons a fortuitous cancellation of errors. If such problems are given due attention, experimental results can be a very useful measure the performance of a theoretical method, but it is worth considering other standards as well. Comparison of a new method with experimental results often introduces more vari ables than one would like. For example, the choice of basis set may not be adequate to obtain highquality theoretical results. Full CI provides the exact answer for a given basis set and Hamiltonian, and therefore offers an excellent standard of comparison for new methods. Unfortunately, FCI calculations are very expensive, so the results available for comparison are somewhat limited. Another option is to compare with methods that, from a welldefined theoretical standpoint, are better than the method in question. An example of this would be comparing methods which approximates triple excitation effects to the full triples model such as CCSD(T) vs. CCSDT, or FSCCSD+T(3) vs. FSCCSDT. In this work, however, as is often the case, the approximate method has been implemented before the full method because of the additional complexity of the full model, forcing such comparisons to wait until the more complete method is implemented. A final al ternative is comparison of the new method with one for which the behavior is relatively well understood, but which cannot be clearly designated as a "better" method. Although these comparisons are the least definite, they are, of necessity, the most common in the calibration of any new method. 66 This chapter is a beginning at the calibration of the FSCCSD+T(3) and FSCCSD+T*(3), models. The molecules chosen for these calculations are drawn largely from applications of the FSCCSD method already in the literature, and in this work extended to include triple excitation contributions. The intent is to get an overview of the performance of the new methods for a variety of systems, most of which have also been the subject of single reference CC calculations, thus providing a useful point of comparison. 5.2 Effect of Changes in the Active Space The triples approximations, as was noted in the previous chapter, do not have the same invariance with respect to the choice of active space as occurs for com plete methods. In order to simplify later comparisons, we first examine the effect of the change in active space on the thirdorder triples models. The "error bars" ob tained in this section may then be kept in mind later on, where a single choice of active space will be used in comparisons of the new FSCC methods with other re sults. Because HN,eff is block diagonal due to molecular symmetry, the active space within a particular symmetry block can be changed without affecting the others. If, in a series of calculations, the active space for a particular irrep remains unchanged, the IPs of that symmetry will not change. This is the case regardless of changes to the active space in other irreps. Conversely, changing the active space within a given irrep will result in slightly different IPs. 67 In order to assess the magnitude of this effect, calculations have been carried out on a number of molecules using different active spaces. The molecules, listed in Table 5.1, include a number for which the IPs themselves are studied below, and several additional systems with fairly low symmetry actual or computational symmetry. The idea behind these additional molecules was to see if having more occupied orbitals in each irrep (and therefore a larger set of possible active spaces) would result in larger variations in the calculated IPs. FSCCSD+triples calculations were carried out on each molecule using a DZP basis at the DZP/SCF optimized geometry for a set of active spaces consisting of 1. the largest active space for which the underlying FSCCSD could be converged (referred to as the "reference" active space), 2. single orbital active spaces for each of the orbitals in the reference active space, and 3. a series of active spaces obtained by starting with the lowest IP and successively adding the next highest IP. This choice provides a slightly larger set of data on the variation of lower IPs, but those are usually the roots of greatest interest anyway. Analysis of these calculations involved first reducing the data to obtain a set of unique active spaces within each symmetry block, then computing for each IP, the mean, range, and mean difference from the available data. The mean was chosen as the reference point because there is no theoretical basis for preferring one active space over another as the "correct" one, except perhaps for having all occupied orbitals active, which proved impossible to converge in most cases. Table 5.1 presents the largest range and mean difference observed among the various IPs 68 Table 5.1. Effect of changes in active space in the T(3) and T*(3) models. for a variety of molecules. Ranges and mean differences of IPs are given in mH. Point Active FSCCSD+T(3) FSCCSD+T*(3) Point Active Molecule nC) Mean Mean group) orbitalsb) Ranged) Ranged) Diff.e) Diff.e) C3H60 Cs 3 3 0.023 0.010 0.023 0.010 CHg) C2v 3 3 0.078 0.035 0.065 0.029 CH2NHg) Cs 4 4 0.744 0.220 0.545 0.165 CH2PHe) Cs 4 2/4 0.103 0.052 0.074 0.028 CH3CIf) C (C3v) 4 4 0.008 0.004 0.006 0.003 CH3CNf Cs (C3v) 4 4 0.354 0.133 0.284 0.107 CH3NCf) Cs (C3v) 4 4 0.044 0.022 0.044 0.022 CH3NH20 Cs 4 4 0.226 0.083 0.175 0.066 H2COg) C2v 2 2 1.118 0.559 0.931 0.466 HOC10 C, 3 2 1.844 0.921 1.581 0.790 HOFO Cs 3 2 3.467 1.734 2.757 1.379 N2g) D2h (Dooh) 3 3 1.297 0.576 0.861 0.383 a) Computational point group (full point group). b) Largest active space used in the irrep containing the root showing the largest variation. c) Number of different active spaces sampled for the root showing the largest variation. Same for T(3)/T*(3) unless two values are given. d) Largest range (difference between largest and smallest values) for a single root. (in mH). e) Mean difference (I =1 IP, IP) for the root showing the largest variation (in mH). 0 Calculated with a DZP basis at the DZP/SCF optimized geometry. g) Calculated using the geometry and basis used in the ionization potential calculations in the next section. for each molecule, along with an indication of the number of distinct active spaces in which that IP was computed. 69 The largest variation is observed for HOF, with HOC1, N2, and H2CO close behind. For a number of other species, however, the variation is quite small. Interestingly, the largest variations are seen among molecules with higher symmetry and smaller active spaces rather than large ones. This obviously suggests that the number of occupied or bitals is not the most important determinant of the variation, however examination of the FSCCSD and reference CCSD results provided no clues to distinguish molecules like HOF or H2CO from ones with a small variation, such as CH3Cl. Comparison of the T(3) and T*(3) results shows that the range and mean difference for T*(3), which includes some higherorder terms, are at least as good as, and generally better than, those for T(3). Although there are a few cases where going to the T*(3) approximation changes which of the IPs has the largest variation, generally it is the same ionization for both. In order to verify that the variation is not due to an inadequacy of the DZP basis sets used for most results in Table 5.1, a second set of calculations was carried out on HOF using Dunning's correlation consistent PVQZ basis [61], which is 5s4p3d2flg on O and F, and 4s3p2dlfon H for a total of 125 basis functions (using spherical harmonic functions), compared to 4s2pld and 2slp for the DZP basis, totalling 37 basis functions (using cartesian functions). Table 5.2 presents the mean ionization potential and its variation with active space for the lowest five IPs of HOF calculated using both basis sets. There are changes in the range of variation from one basis to the other some roots have a larger range with the PVQZ basis, others have a smaller range. In general, however, the they are roughly the same, suggesting that the variation is largely independent of the basis set. 70 Table 5.2. Effect of changes in the active space on the calculated ionization potentials (in mH) of HOF for two basis sets. DZP basis PVQZ basis Sym. n Mean Mean Mean Mean Range Range IP Diff. IP Diff. FSCCSD+T(3) a' 3 545.7 1.158 0.608 559.0 1.176 0.480 a' 3 672.8 3.186 1.064 686.0 3.529 1.281 a' 2 775.7 3.467 1.734 765.3 3.714 1.856 a" 2 466.5 0.171 0.085 482.0 0.118 0.059 a" 2 657.6 0.138 0.069 674.2 0.093 0.047 FSCCSD+T*(3) a' 3 538.0 1.391 0.547 550.5 1.076 0.431 a' 3 665.9 2.465 0.853 678.1 2.872 1.033 a' 2 749.8 2.757 1.379 758.4 3.080 1.540 a" 2 459.7 0.073 0.037 474.3 0.033 0.017 a" 2 646.9 0.044 0.022 662.5 0.012 0.006 Given the data above, a reasonable estimate of the error bars due to changes in the active space is 2 mH for both methods. (It must be pointed out that there are certainly to be systems with larger variations, but more experience will be required to understand the origin of the large variations and therefore be able to identify the worst offenders.) These correspond to about 0.05 eV or 1.3 kcal mol1. Although this may seem like a fairly large uncertainty, particularly for comparing energies at different points on a potential energy surface, there are several important factors worth noting. First, this is a conservative estimate for entirely unknown molecules. Many of the molecules in Table 5.1 are more than an order of magnitude better. Also, the table presents only 71 the largest variation many others were significantly smaller. Secondly, it must be remembered that for molecules with symmetry, many of the IPs exhibit no variation because there is only one possible choice of active space in that symmetry. Finally, it is a fairly simple matter, and not computationally expensive, to obtain a feel for the variation of any molecule of interest. Only one (0,0) sector calculation is required to try out several FSCC active spaces, and the FSCC calculations themselves are very inexpensive compared to the (0,0) reference. As a result, while the uncertainties associated changing the active space in T(3) and T*(3) calculations for the general case might be considered large for some applications much tighter limits may often be obtained with minimal effort. For the remainder of this chapter, the active spaces used will correspond directly to those results which are reported. Uncertainties due to the active space will be introduced into discussions as appropriate. 5.3 Ionization Potentials Although it can be used in many other ways, the "natural" result of a (0,1) sector FSCC calculation is the energy of the vertical ionization process extracting an electron from a molecule without allowing the geometry to relax to that of the resulting ion. Cal culation of ionization potentials is, thus far, the predominant application of the FSCCSD method and such calculations will provide the bulk of the data for characterization of the approximate triples methods. These results will be compared with SRCC calculations, experimental results, and a benchmark full CI calculation for the methylene molecule by Bauschlicher and Taylor [62]. 72 Table 5.3. Vertical ionization potentials of 'A, and 3B1 CH2 (in eV) calculated with a DZP basis. Method A 3B 12A, 22A1 2B2 2A1 2B1 FSCCSD 10.21 22.54 14.92 10.14 11.30 FSCCSD+T(3) 10.32 22.74 15.00 10.26 11.43 FSCCSD+T*(3) 10.28 22.68 14.96 10.23 11.39 Full CI 10.26 22.14 14.85 10.16 11.31 QRHFCCSD 10.20 22.25 14.84 10.14 11.30 QRHFCCSD[T] 10.25 22.28 14.85 10.16 11.31 QRHFCCSDT1 10.25 22.28 14.85 QRHFCCSDT 10.25 22.15 14.85 Note: Geometry, basis set, and full CI results taken from Bauschlicher and Taylor [62]. QRHFCC results for 'A1 ionizations from Watts, Gauss, and Bartlett [58]. 3B calculations used a UHFCCSD reference function. Table 5.3 presents FSCC calculations on ionization potentials of methylene in com parison with both the full CI data from Bauschlicher and Taylor, and several QRHFCC calculations. At the FSCCSD level, four of the five calculated IPs are within 0.1 eV of the full CI number, with the 'A,+22A result 0.4 eV off. The effect of the T(3) correction to FSCCSD is to move all of the results away from full CI, but four of the roots are still within 0.1 eV, while the A1i+22A, number has changed to 0.6 eV away from full CI. Introducing some higher order effects via the T*(3) model shifts all of the results back toward full CI, but not as close as the FSCCSD result. These differences from full CI cannot result from uncertainties due to the choice of active space, which from Table 5.1 are about 0.001 eV (0.04 mH) for this particular molecule. By contrast, 73 the single reference QRHFCC calculations start out with all five ionizations within 0.11 eV of full CI, and the inclusion of triple excitations at various levels improves the results to match full CI within 0.01 eV. This behavior suggests that the single reference CCSD provides a good approximation to the various states in question, and that triple excitations are relatively unimportant, which would seem to conflict with the FSCC results. It is important to recall, however, that there is a fundamental difference between FSCC and SRCC with respect to the nature of their excitation operators. As discussed in the previous chapter, spectator triple excitations must be added to the FSCCSD method to achieve an expansion space that is in some sense equivalent to that of the SRCC. But even then, the FSCC amplitudes are restricted by the necessity to represent several ion ized states simultaneously. Thus, we should not expect the FSCCSD to produce results identical to a single reference CCSD, nor should we expect triple excitation corrections to behave in the same way. The differences between FSCC and SRCC shown in Table 5.3 and others, below, should be interpreted as a manifestation of these differences. These differences between single reference and Fockspace methods can also be seen in the other calculations presented here (Tables 5.45.9). The sign of the FSCC triples correction is frequently opposite that of the SRCC triples corrections, and their magni tudes do not appear to be related a large SRCC triples contribution does not necessarily imply a large FSCC triples correction, and vice versa. The size of the SRCC and FSCC triples contributions are rather different as well. The largest SRCC effect is about 0.2 eV, while the FSCCSD+T(3) correction is as large as 0.8 eV, and 0.5 eV for T*(3). All this is further support for the idea that the nature of the excitation operator hierarchy is 74 Table 5.4. Vertical ionization potentials of CH2NH (in eV) calculated at the experimental geometry with a DZP basis. Orbital Basis/Method FSCCSD FSCCSD+T(3) FSCCSD+T*(3) QRHFCCSD QRHFCCSDT1 5a' 10.34 10.71 10.55 10.40 10.43 10.52 10.70 10.56 Experiment [63] Experiment [64] Experiment [65] la" 12.29 12.35 12.29 12.14 12.28 12.43 12.48 12.44 4a' 15.08 15.30 15.20 15.17 15.11 15.13 15.11 15.00 3a' 17.19 17.46 17.34 17.32 17.14 17.04 17.07 17.00 Notes: Geometries, basis sets, and QRHFCC results taken from [66]. Table 5.5. Vertical ionization potentials of CH2PH (in eV) calculated at the experimental geometry with a DZP basis. Orbital Basis/Method la" 5a' 4a 3a FSCCSD 10.08 10.28 13.17 15.14 FSCCSD+T(3) 10.01 10.33 13.25 15.41 FSCCSD+T*(3) 9.98 10.28 13.25 15.27 QRHFCCSD 9.90 10.22 13.15 15.15 QRHFCCSDT1 10.05 10.24 13.15 Experiment [67] 10.30 10.70 13.20 15.00 Notes: Geometries, basis sets, and QRHFCC results taken from [66]. different in single reference and Fockspace coupled cluster methods. 75 Table 5.6. Vertical ionization potentials for 3EEg 02 (in eV) at the experimental geometry. Method Basis 02 Final State Method Basis 2rIg 4E 4iu 4Eu FSCCSD DZP 11.76 17.75 16.40 24.16 FSCCSD+T(3) DZP 12.13 17.99 16.30 24.98 FSCCSD+T*(3) DZP 12.08 17.87 16.29 24.66 UHFCCSD DZP 12.13 17.80 16.15 24.64 UHFCCSD(T) DZP 11.95 17.77 16.26 24.40 FSCCSD PVQZ++ 12.37 18.34 16.94 24.74 FSCCSD+T(3) PVQZ++ 12.82 18.48 16.74 25.43 FSCCSD+T*(3) PVQZ++ 12.57 18.34 16.71 25.11 UHFCCSD(T) PVQZ++ 12.39 18.26 16.76 24.82 Expt. (Adiabatic) 12.07 18.17 16.11 24.58 Expt. (Vertical) 12.35 18.33 16.85 24.66 Notes: Geometry taken from Herzberg [68]. Basis sets, SRCC, and experimental result taken from [69]. Examining the various tables, we see that the SRCC triples correction usually improves agreement with experiment, while the FSCC triples corrections appear to shift the calculated result away from the experimental one as often as towards it. It is also worth noting that the T*(3) correction is generally smaller than T(3) and gives better agreement with experimental results. Overall, FSCCSD and FSCCSD+T*(3) seem to do about the same in comparison with experiment, with an average absolute difference from the experimental results of roughly 0.2 eV. The single reference CCSD results have an average absolute difference of closer to 0.3 eV, which is improved to 0.2 eV 76 Table 5.7. Vertical ionization potentials of H2CO (in eV) evaluated at the experimental geometry with a [5s3pld/2slp] basis. Orbital of ionized electron Method 2b2 lbl 5al FSCCSD 10.51 14.36 15.75 FSCCSD+T(3) 10.96 14.49 16.30 FSCCSD+T*(3) 10.77 14.40 16.06 QRHFCCSD 10.59 14.19 15.77 QRHFCCSD[T] 10.63 14.37 15.83 QRHFCCSDTlb 10.63 14.39 15.83 QRHFCCSDT3 10.58 14.31 15.76 Experiment [70] 10.88 14.38 16.00 Notes: Geometry taken from Herzberg [68]. The basis set consisted of Dunning's (9s5p)[5s3p] set on C and 0, and his (4s)[2s] set with exponents scaled by 1.44 on H [71], with polarization functions exponents of 0.654 (C), 1.211 (0), and 0.7 (H) [72]. by the addition of triple excitation effects. Thus, it would appear that FSCCSD and FSCCSD+T*(3) perform about as well as single reference CCSD(T) in comparison with experiment, but the single reference calculations (so far) have an edge in the predictability. In order to give proper weight to the performance of these methods relative to experiment, we must consider the potential problems inherent in experimental data, as mentioned in the Introduction. Most of these problems are displayed by one or more of the experimental results in this chapter. The natural result of a FSCC calculation is a vertical ionization potential, while experiment may provide vertical, adiabatic, or both. Adiabatic IPs are not directly comparable with calculated vertical ionization potentials, since they depend on the change Table 5.8. Vertical ionization potentials of N2 basis. (in eV) at R=2.0693 au with a [5s4pld] Orbital of ionized electron Method 3aCg l 2au FSCCSD 15.44 17.14 18.64 FSCCSD+T(3) 15.58 16.86 19.04 FSCCSD+T*(3) 15.44 16.88 18.81 QRHFCCSD 15.39 16.69 18.71 QRHFCCSDT1 15.36 16.87 18.57 Experiment 15.60 16.98 18.78 Notes: Geometry, basis, and QRHFCC and experimental results taken from Rittby and Bartlett [73]. Table 5.9. Vertical ionization potentials of ketene and diazomethane (in eV) calculated at experimental geometries calculated with a DZP basis. Ketene Diazomethane Method 2bl 2b2 2bl 2b2 FSCCSD 9.40 14.25 8.60 13.79 FSCCSD+T(3) 9.39 14.49 8.52 13.85 FSCCSD+T*(3) 9.35 14.38 8.47 13.77 QRHFCCSD 9.29 14.31 8.43 13.78 QRHFCCSD(T) 9.38 14.24 8.54 13.73 Experiment 9.64 (a) 13.84 (a) 9.00 (a) 14.13 (v) Notes: Geometries, basis sets, and experimental results taken from Rittby, Pal, and Bartlett [55]. Experimental IPs are adiabatic or vertical, as indicated in the table. in geometry from the neutral to the ion. If the ion and neutral geometries are very similar, the vertical and adiabatic IPs will also be nearly the same, but in general, the adiabatic 78 IP is only a lower bound for the vertical one. In cases where the experimental value is adiabatic or unspecified, a larger difference from the computed vertical IP must be expected. A second problem with experimental results is that, as with theoretical methods, different experimental techniques can give somewhat different results. This is illustrated in Table 5.4, which lists three sets of experimental data for methyleneamine, differing among themselves by as much as 0.2 eV. A third complication of using experimental data is that details of the calculation other than the correlation method become important, particularly the basis set. Most of the calculations presented here use a DZP or slightly better quality basis set. While very useful for comparing one theoretical method to another, these basis sets may not be adequate, however, to obtain good results compared to experiment. Table 5.6 shows the ionization potentials of 02 as calculated using both DZP (32 functions) and an augmented PVQZ basis (126 functions) in order to demonstrate the effect of basis set on both the FSCC and SRCC IPs. It is clear that the PVQZ++ basis provides a much better comparison with the experimental results than the DZP basis, for both FSCC and SRCC calculations. The change in FSCC results due to increasing the basis set is about 0.50.6 eV, and for the SRCC 0.40.5 eV. Although there is no experimental comparison, a similar basis set effect can also be seen in the HOF results in Table 5.2, where the change from a DZP to a PVQZ basis changes the IPs by 0.20.4 eV. It would seem, then, that correlated calculations with DZPquality basis sets should not be expected to reproduce experimental results too well. 79 Given that the basis sets used here are generally too small, coupled with the other problems of the experimental data themselves, it would seem that an error of 0.2 eV for these calculations is not unreasonable, and is not necessarily an indication of intrinsic problems with the method. Indeed, the fact that the SRCC results show a similar margin of error with respect to experiment is a good sign for FSCC. It would seem that the biggest difference in this regard between the single reference and Fockspace results is the relative lack of consistency of the FSCC triples corrections, which can probably be better understood with more calculations. 5.4 Potential Energy Surfaces and Related Properties Another use of FSCC methods is for the study of potential energy surfaces, especially with the ability to look at several surfaces in the same calculation. The basic result of the FSCC calculation is still an energy difference, but when added to the absolute energy of the (0,0) sector reference, it provides an absolute energy for the state of interest. In the (0,1) sector case, the state of interest has one electron removed relative to the reference (0,0) state. While this may seem a round about way of calculating the PES of a molecule, it has practical advantages in many cases. For calculations on molecules with lowlying electronic states or symmetry breaking phenomena, direct calculation of the PES by traditional single reference methods frequently suffers from convergence and other problems due to the proximity of other states to the one of interest. On the other hand, the same molecule with another electron is often wellbehaved at the same geometries. This the same as the ionization process, but from the opposite point of 80 view in this case, we want to talk about the final state, while in IPs, it is traditional to look at it from the point of view of the initial state. This application of FSCC methods is a very recent development, and there are only a few such applications in the literature. Moreover, because they deal with symmetry breaking problems, there is little conclusive experimental evidence to allow concrete conclusions to be drawn about the performance theoretical methods. Comparisons are possible with single reference calculations, though as we shall see, the comparison of highlevel SRCC calculations with Fockspace results does not always lead to an unambiguous determination of which approach is more reliable. Nevertheless, FSCC is very well suited to such applications and will undoubtedly see increasing use in the study of potential energy surfaces (beginning with the next chapter of this work), so it is worth trying to get a feel for its performance. Engelbrecht and Liu, in a 1983 paper [74], used multiconfiguration selfconsistent field (MCSCF) plus CI calculations to characterize the 3A2 state of CO2. This state, which is bent (C2v symmetry), has the SCF occupation (core)3a12b64a13b25a1 b la24b26al and interacts strongly with the configuration (core)3a 2b 4a 3b 5a lb 1a 4b 6a because both the a, and bl orbitals transform as a" in C, symmetry. As a result, SCF calculations can obtain a lower energy by going to a Cs structure. 81 Asymmetric Stretch Potential of 3A2 CO2 12 1  FSCCSD 10 10 S1 FSCCSD+T(3) 8 FSCCSD+T*(3) 0 QRHFCCSD *' 6 a QRHFCCSD(T) 2 1 0.10 0.05 0.00 0.05 0.10 AR (A) Figure 5.1. Potential curves for the asymmetric stretch mode (AR = !(R1 R2)) of 3A2 CO2 at various levels of theory. The C2v structure is the DZP/FOCIB geometry of Engelbrecht and Liu [74]. With suitable methods, such interactions can be treated properly, and reason able potential energy surfaces obtained, but this requires the ability to describe inherently multireference states. Most such methods are not manybody in na ture, and of course we prefer a size extensive approach if possible (Chapter 3). It has been observed that QRHFbased coupled cluster methods are capable of han dling many multireference problems when sufficient triple excitation effects are incor porated [34]. In Figure 5.1, we see two QRHFCC potential curves for the asymmet ric stretch of CO2, which were calculated from closedshell CO22 SCF orbitals. The QRHFCCSD curve has a double minimum shape, while the QRHFCCSD(T) method 82 gives a stable C2v structure. This is consistent with the MCSCF+CI calculations of Engelbrecht and Liu and their conclusion that the equilibrium structure of this state is C2v. The results of three sets of Fockspace coupled cluster calculations using 2A1 CO2 as the (0,0) reference are also displayed in Figure 5.1. FSCCSD gets the proper shape of the PES even without triple excitation effects, and the inclusion of triples makes small adjustments to the shape of the surface. (Note that, as in the IP calculations, the FSCCSD+T*(3) results lie between FSCCSD and FSCCSD+T(3).) This behavior is what we might expect from a method that is designed in a multireference framework. Another symmetrybreaking problem to which a number of single and multiref erence approaches have been applied is the nitrate radical, N03. Different exper iments have been interpreted in terms of either a totally symmetry D3h equilib rium geometry, or a lower symmetry C2v structure with one NO bond length dif ferent from the other two. However, as Kawaguchi and coworkers [75] comment, no one has advanced a model which explains all of the available spectroscopic data. In the most recent theoretical examination of the problem, Stanton, Gauss, and Bartlett [76] located three stationary points on the NO3 surface using the QRHFCCSD method. The D3h structure was determined to be a local minimum, and two C2v structures were found. They are characterized by the relationship between the unique NO bond distance and the other two: one long/two short (1L2S) or one two long/one short (2L1S). The 1L2S structure was also found to be a local minimum, and the 2L1S a transition state 192 cm' 83 Table 5.10. Calculated relative energies (in kcal mol1') of C2v and D3h NO3 structures. DZP basis PVTZ basis Method 1L2SD3h 2L1SD3b 1L2SD3h 2LISD3h FSCCSD 2.50 1.05 2.94 1.05 FSCCSD+T(3) 3.03 1.64 3.72 1.73 FSCCSD+T*(3) 2.74 1.35 3.36 1.41 QRHFCCSD 2.59 2.04 2.96 QRHFCCSD(T) 0.34 0.57 QRHFCCSDT3 0.21 Notes: Geometries, optimized at the DZP/QRHFCCSD level, taken from [76]. The C2v geometries are characterized by one of the NO bonds being longer (1L2S) or shorter (2L1S) than the other two. QRHFCC relative energies taken from [77]. higher in energy than the 1L2S state; both C2v structures being lower in energy than the totally symmetry structure. Their conclusion is that the molecule has a C2v equilibrium geometry, but probably has an extremely low barrier to pseudorotation, which would give spectra with characteristic of a D3h system. In a later paper focusing on different orbital choices for single reference treatments of symmetry breaking problems, Stanton et al. [77] calculated relative energies of the D3h and C2v minima using more sophisticated methods than in their first paper. Their QRHF CC results are shown in Table 5.10 along with FSCC results using the same basis sets and geometries. Using both DZP and PVTZ basis sets, QRHFCCSD calculations favor the C2v minimum by about 3 kcal mol1. As soon as triples corrections are introduced, either using the noniterative CCSD(T) approximation or the iterative CCSDT3 model, the energies shift in favor of the totally symmetric structure, but by only 0.5 kcal mol'1. 84 FSCC, on the other hand, predicts the D3h structure to be lower for all three methods shown, but by roughly 3 kcal mol1. Given previous experience with QRHFCC methods applied to problems of this sort, we should have greater confidence in the QRHFCCSD(T) and T3 results than in results without triple excitations, but with the available data, it is impossible to tell whether the single reference numbers including triple excitations or the FSCC results are more reliable. As Stanton et al. [77] observe, "An extremely sophisticated (and expensive) calculation would be required to firmly establish the nature of the D3h stationary point." C3+ is another molecule for which both experimental and theoretical studies have produced differing opinions on the nature of the ground state. Coulomb explosion data can be interpreted in terms of either a floppy linear molecule or a bent structure, depending on the vibrational temperature of the sample [78,79]. On the theoretical side, the consensus is that the ground state is a bent 2B2 with the linear 2 u' higher in energy [80]. The energy difference between these two states is still not a settled issue. This is important because if the barrier is small enough, the molecule may be quasilinear although the true ground state is bent, a large amplitude bending motion could bring the molecule through the linear structure. Tables 5.11 and 5.12 present the results of FSCC optimizations, vibrational frequency calculations, and relative energies in comparison with single reference CC results from Watts, Stanton, Gauss, and Bartlett [80] which include full CCSDT calculations. All three Fockspace methods give essentially the same geometry for both structures, comparing most closely with the single reference CCSD results. For the vibrational frequencies of the Table 5.11. Optimized geometries and basis. vibrational frequencies of 2B2 C3+, using a DZP Metd R () 0 () wl (a,) w2 (al) w3 (b2) Method R (A) (0) c ( (m' (cm ) (cm ) (cm1) FSCCSD 1.328 71.6 1685 583 1269 FSCCSD+T(3) 1.336 70.2 1643 614 1221 FSCCSD+T*(3) 1.335 70.5 1651 603 1225 CCSD 1.337 68.3 1668 687 1215 CCSD(T) 1.350 68.3 1601 515 1074 CCSDT1 1.350 69.0 1558 569 1078 CCSDT2 1.347 69.0 1612 595 1213 CCSDT 1.350 68.0 1601 638 1021 Notes: Basis set and SRCC results taken from Watts, et al. [80]. Table 5.12. Optimized geometries, and vibrational frequencies of 2Eu+ C3+, and energy relative to 2B2 C3+, using a DZP basis. E(2 +  (Method R) w(ag) w(ru) E(2  Method R (A) (cm1 (cm1) (cm1) 2B2) (cm4) (cm') (cm') (kcal mol') FSCCSD 1.312 88 1227 492i 7.1 FSCCSD+T(3) 1.315 71i 1203 1911 10.9 FSCCSD+T*(3) 1.316 532 1205 1462 8.6 CCSD 1.314 2500i 13.2 CCSD(T) 1.327 1431i 2.9 CCSDT1 1.339 1135 1.7 CCSDT2 1.329 1595 2.2 CCSDT 1.327 451i 4.9 Notes: Basis set and SRCC results taken from Watts, et al. [80]. 2B2 state, the FSCC results are also basically the same, and correspond most closely with 86 the CCSDT2 approximate iterative triples model among the single reference calculations. On the basis of several calculations with a larger basis set, Watts and coworkers suggest that the true harmonic frequencies should lie in the ranges 15901629, 592638, and 10941173 cm1, and compared to these, the FSCC wl and w3 are slightly high. The asymmetric stretching frequency (au) of the 2Eu+ state (Table 5.12) is a more complex problem. An imaginary frequency in this mode indicates that the Doh structure is a transition state, and a lower energy can be obtained by changing the bond lengths asymmetrically. For both single reference and Fockspace methods, the calculated frequency varies widely. Also, the bending mode (Tru) changes dramatically with the choice of FSCC method. Such large changes due to different correlation methods are usually an indication that important aspects of the problem are not being taken into account, and a more detailed examination, using both single reference and Fockspace methods would be in order. Table 5.12 also presents the energy differences between the linear and bent structures. Note that the single reference CCSDT1 method places the linear structure lower in energy than the bent one. Watts et al. attributed this to an apparent tendency for the iterative nature of the CCSDT1 method to overemphasize problems in the closelyrelated CCSD(T) model in difficult cases, thus producing poorer results than the perturbative CCSD(T) method. Their CCSDT AE of 4.9 kcal mol1 is in reasonable accord with CISDTQ calculations of Grev, Alberts, and Schaefer [81,82] indicating a 4.1 kcal mol' difference. Taylor et al. [83] used multireference CI methods to obtain a separation of 1.68 kcal mol' with 87 a DZP basis set, and claim that this should be very close to the full CI limit. These authors also show that multireference effects are important for both the bent and linear structures. Despite the multireference character of the FSCC approach, all three methods presented here give a significantly larger energy gap than Taylor's calculations. The problem with the relative energies of the two states may be due, at least in part, to the apparent problems in properly describing the 2Eu+ state. It would seem, on this minimal evidence, that where single reference methods including triple excitations have problems, Fockspace methods do as well. 5.5 Conclusions A variety of results has been presented using the Fockspace coupled cluster methods FSCCSD, FSCCSD+T(3), and FSCCSD+T*(3), for both ionization potential problems at fixed geometries, and in cases where the potential surfaces must be explored. Overall, the performance of the FSCC methods is good. FSCCSD results usually compare well with single reference CCSD(T) calculations. The FSCC triples corrections seem initially to be less predictable, than the triples correction in single reference methods. The Fockspace triple excitations have a different nature than their SRCC counterparts, and it will require more experience with these methods to understand them and discern the system they follow. Including some higherorder contributions in the T*(3) method seems to have a "dampening" effect compared to the basic T(3) model, which causes the T*(3) results to compare better with experimental results or other calculations. This is probably similar 88 to the single reference coupled cluster case, where the simplest perturbative inclusion of triples, CCSD+T(CCSD) was found to be significantly improved by the inclusion of certain higher order terms, resulting in the CCSD(T) model. For the present, the FSCCSD+T*(3) method seems to be preferable to the purely thirdorder FSCCSD+T(3). In the long run, it is probably worthwhile to investigate more systematic ways of incorporating higherorder terms in the FSCC triples contribution, as discussed in the previous chapter. Numerically, effects of changing the active space in the FSCCSD+triples methods have been demonstrated to introduce uncertainties as large as 2 mH in calculated energies (IPs or state energies), however most of the molecules examined have much lower uncertainties. On the whole, both FSCCSD and FSCCSD+T*(3) seem to do about as well compared to experiment as SRCC methods including perturbative triples, and better than CCSD. For very accurate comparison with experiment for both FSCC and single reference methods, however, a large basis is required. For energy differences on potential energy surfaces, where the requirements are tighter than on IPs, Fockspace results do not seem to compare as well with single reference methods as for IPs. Although the signs of relative energies evaluated using FSCC seem to match higherquality single reference calculations, the Fockspace energy gaps are noticeably larger than in the SRCC case. The cases examined are hard problems, both experimentally and theoretically, and it is not clear that either approach can be claimed to be better. Consequently, great care must be taken with such results, and more study would be most useful. 89 The FSCC methods all seem to do quite well in determining local properties of potential surfaces geometries, vibrational frequencies, etc. In most of these cases, there does not appear to be a substantial difference among the three FSCC methods examined, except for linear C3*, for which both single reference methods including triple excitations and Fockspace methods seem to have problems. In addition to the quality of results of a new method, one should also consider its ease of use and cost in comparison with competitive methods. By this measure, FSCC methods have a decided advantage over the single reference approaches. An FSCC calculation, even including triple excitation effects, amounts to a small fraction of the cost of the (0,0) sector CCSD reference calculation. In contrast, each of the IP calculations required separate SRCC calculations for the neutral, and for each ion state, and at least one of those on an openshell molecule, for which the SRCC is about four times more expensive than a closed shell calculation. In most cases, the reference CCSD calculation used in the Fockspace approach is a closed shell, and all of the IPs can be computed in a single calculation. For potential energy surfaces, the ability to explore several surfaces at once is also a distinct advantage over single reference methods. CHAPTER 6 THE STRUCTURE OF THE C5+ MOLECULE 6.1 Background Although there is relatively little experimental data on C3+, discussed in the preceding chapter, there is even less known about Cs+. The experiments which have observed Cs+ have done so with a mass spectrometer in reactivity [84], frag mentation [8587], or mobility [88] studies. The essential conclusion of these pa pers (with respect to C5+) is that there is a single isomer (multiple isomers show up starting with C7) and that its reactivity and mobility are similar to those of its neighbors, C4+ and C6+, and the whole set of clusters from C2z to C6 are commonly labeled as "linear" molecules [84,88]. (C3+, which theoretical results agree is actually bent, fits the trends for "linear" clusters in the experimental data as well as C5+ does, suggesting that the "linear" label cannot be interpreted rigor ously.) Recently, matrixisolated infrared spectroscopy experiments, using isotopically mixed graphite, produced an IR band that could be attributed to a planar pentago nal (Dsh symmetry) C5+ molecule [89]. The neutral Cs cluster has been well es tablished, both theoretically and experimentally, to have a linear ground state [90], with cyclic structures much higher in energy. In this situation, it would seem un likely that the cation would be so different from the neutral that a cyclic cation 91 would be observable, but the existence of isoenergetic linear and cyclic structures of neutral C4, which is now generally accepted [90,91], was once also thought un likely. A series of SCF and MBPT(2) calculations [92] established that D5h C5+ could not be responsible for the observed band, and a reanalysis of the experimental data led to assignment of the band to a C3*H20 adduct, in agreement with a later assignment by Ortman et al. [93]. The lowlevel calculations pointed toward a linear ground state for Cs5, but there were major problems with spin contamination of the UHFbased calculations, stability of the SCF solution, symmetry breaking, and failure to converge calculations on some states. In this chapter we attempt to understand the source of the troubles with lowerlevel calculations and to establish the true ground state structure of Cs5 using highquality manybody methods. 6.2 Computational Methods The majority of the calculations in this work used a DZP basis set comprised of Dunning's (9s5p)/[4s2p] doublezeta set [71] and a polarization function (a0.654) as optimized by Redmon, Purvis, and Bartlett [72]. All six cartesian d components were retained, giving a total of 80 contracted basis functions. At selected points on the various potential energy surfaces, a larger PVTZ basis [61] was used to compute the energy. This is a (10s5p2dlf)/[4s3p2dlJ] generallycontracted basis, which was used with spherical harmonic polarization components (5d and 7J), for a total of 150 contracted basis functions. 92 Both single reference QRHFbased and Fockspace coupled cluster methods have been used and in both types of calculations, triple excitation effects have been con sidered, mostly through perturbative corrections: QRHFCCSD(T), FSCCSD+T(3), and FSCCSD+T*(3). Both QRHFCC and FSCC calculations used the closed shell neutral Cs molecule as a reference. The reasons for using both SRCC and Fockspace approaches are both historical and practical. When work on this project was first begun, no effi cient FSCC program was available and QRHFbased approaches were the best choice to deal with the problems observed in the UHF calculations mentioned earlier. When the newlycompleted ACES II FSCCSD code was applied to this molecule, the nature of the problems became apparent, and it was clear that understanding the Cs5 molecule was a problem wellsuited to this method because of its closelyspaced, strongly interacting electronic states. In some calculations, the carbon Is orbitals and the corresponding virtual orbitals (five occupied and five virtual for the DZP basis with cartesian angular momentum functions, five occupied and four virtual for the spherical harmonic PVTZ basis) were not correlated in order to reduce the computational cost. The effect of dropping these orbitals on the relative energies is negligible. Geometry optimizations have been carried out using QRHFCCSD and all three Fockspace methods. Analytic gradient techniques were used with the QRHFCCSD calculations, but gradients have not been formulated for FSCC yet, so optimizations were performed via numerical differentiation of FSCC energies. Similarly, vibrational frequencies were evaluated from the numerical second derivatives of the FSCC energy, 