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 
THEORY AND IMPLEMENTATION OF COUPLEDCLUSTER METHODS FOR DIRECT CALCULATION OF IONIZATION ENERGIES By RENEE PELOQUIN MATTIE 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 1995 DEDICATION To Robert, You encouraged, supported, and pushed me throughout the process. If anyone is happier than I am to see this in its final form, it must be you. ACKNOWLEDGEMENTS My thanks to: Professor Bartlett for sharing his knowledge and wisdom and for bringing me into his group; my entire committee for their patience; John Stanton and Jiirgen Gauss for teaching me the joy of scientific collaboration; Anna Balkova and Hideo Sekino for much practical advice, abstract dialogues, tea, philosophical discussions, and latenight rides back from campus; Ajith Perera and all the students in the group for taking care of details longdistance; Karen Yanke for her help above and beyond the call of duty; Jan Musfeldt and David Bernholdt, Johna Till and Ted Johnson, for all kinds of help and encouragement, both scientific and nonscientific; and everyone else who gave me a kick when I needed it. TABLE OF CONTENTS page ACKNOWLEDGEMENTS .................... ........ iii LIST OF TABLES ................................ viii LIST OF FIGURES ..... ..... ...... ... .. ..... ... ix ABSTRACT .......... ............ ............. .. .. x CHAPTERS 1 INTRODUCTION ............................... 1 1.1 Correlated Methods in Quantum Chemistry .......... ... 1 1.2 The Electron Detachment Problem ............. ... 4 2 ELECTRON CORRELATION ........................ 6 2.1 The Schr6dinger Equation ........................ 6 2.1.1 The Independent Electron Approximation and Fock space 7 2.1.2 Second Quantization ....................... 9 2.2 CoupledCluster Theory ......................... 13 2.2.1 The CoupledCluster Operators .......... .... 13 2.2.2 The H Operator ......................... 14 2.2.3 The CoupledCluster Equations .. ........... .. 16 2.3 Configuration Interaction ......................... 17 2.3.1 The CI Operators ......................... 17 2.3.2 The CI Equations ....... ... .............. 18 3 THE COUPLEDCLUSTER EQUATIONOFMOTION METHOD .... 19 3.1 The EOM Equations ........................... 19 3.2 SizeExtensivity and the EOM problem ................. 21 3.3 The EOM Wavefunction ................... ...... 21 3.4 The Ionization Energy Equations ............. .. 24 3.5 The EOMIP Wavefunction .... .................... .. 26 4 THE FOCKSPACE COUPLEDCLUSTER METHOD ........... 27 4.1 FockSpace Theory ............................ 27 4.2 Comparison of EOM to FSCC ...................... 33 4.2.1 Definitions ................... ......... 33 4.2.2 FSCC Method ........................... 34 4.2.3 EOM Method ................... ....... 35 4.2.4 Comparison ................... ........ 36 5 THE ELECTRON PROPAGATOR ....................... 37 5.1 Superoperator Expression of the Electron Propagator ......... 38 5.2 The Dyson Equation ........................... 40 5.3 NonDyson Equation Approaches ................ .... 41 6 OTHER METHODS FOR IONIZATION CALCULATIONS ........ 43 6.1 State Methods ........ ........................ 44 6.2 Direct Approaches ................... ........ 46 6.2.1 TwoHole OneParticle Configuration Interaction ....... 46 6.2.2 SpinAdapted Cluster Methods ................. 47 6.3 Eclectic Approaches ............................ 48 7 COUPLEDCLUSTER LINEARRESPONSE THEORY .......... 49 7.1 Perturbation Applied to the CC Hamiltonian ............. 49 7.2 coupledcluster Response Theory ................ .... 53 7.2.1 Time Dependent Perturbation on a CoupledCluster Reference 53 7.2.2 The Linear Response Function . ... 55 7.2.3 Energies and Transition Intensities ............ ... 57 8 METHODS FOR TRANSITION MOMENTS IN THE EOMIP APPROACH 60 8.1 The Monopole Approximation ................. .... 61 8.2 Transition Moment Expression for NonSymmetric Eigenstates .. 63 8.3 Transition Moment Expression for the EOMIP problem ........ 63 9 IMPLEMENTATION ................... ........... 67 10 APPLICATION OF THE EOMIP METHOD TO SMALL MOLECULES 72 10.1 A Survey of FSCCSD and EOMCCSD Calculations ......... 72 10.2 Unsaturated Molecules ................... ...... 76 10.2.1 Nitrogen .............................. 76 10.2.2 Ethylene .................... .......... 82 10.2.3 Butadiene ............................. 89 10.3 Saturated Molecules ................... ...... 93 10.3.1 Ammonia ................... ......... 93 10.3.2 W ater ... .. .. .. .. .. .. .. ... 94 10.4 Nitrogen and C2: A study in contrasts . ... 97 10.4.1 SingleDeterminant and Multideterminant References ... 97 10.4.2 The Spectator Triples Problem . ... 102 11 LINEAR CARBON CLUSTER ANION PHOTOELECTRON SPECTRA 105 11.1 C2, C+, and C .................... ....... 108 11.1.1 Ionization of C2 ................. .. ........ 108 11.1.2 Ionization of C ...................... ... .. 109 11.2 Ionization of C3 ... .......... ... ........... 110 11.3 First Ionization Energies ........... .... ..... ..... 112 11.4 Interpretation of Spectra .. ........... .. ... ......... 112 11.4.1 OddLength Chains .. .. .......... ...... .. 115 11.4.2 C ...... ..... ... ... ... .. ... ....... 115 11.4.3 C .... ....... ....... .... .. ...... 120 11.4.4 C7 ....... .. ....... ....... .. 120 11.4.5 C9 ............... .. ...... ....... ..... 123 11.4.6 EvenLength Chains ... ................. 125 11.4.7 C ...... .. ......... .. ... ......... 126 11.4.8 C ....... ........................ 129 11.4.9 C .............................. 133 11.4.10C8 .... .. .... ........ ... 135 11.4.11Summary ..... . .. 135 12 SUMMARY .............. ...... .. ........ ..... 137 APPENDICES A INVARIANCE OF FOCKSPACE ENERGIES TO CHOICE OF ACTIVE SPACE ............. ........ ... ...... .......... 144 A.0 Invariance of the Ionization Energies . .... 144 A.0 A Modification of the Bloch Equation to Preserve Invariance in Any Sector ....... ............. ... .......... 148 A.0 The Bloch Equation and Modified Bloch Equation . ... 153 A.0.11 Equations in the Original Model Space . ... 154 A.0.11 Equations in the Reduced Model Space . 154 B THE IMPORTANCE OF TRIPLES CC EOMIPSDT PERTURBATION THEORY .... ................................ .. 156 B.0 The CC EOMIPSDT Equations ................ 156 B.O A review of the RayleighSchr6dinger Perturbation Theory 158 B.O EOMIP Perturbation ........................... 159 REFERENCES ................... ................ 161 BIOGRAPHICAL SKETCH ............................ 166 LIST OF TABLES A survey of FSCCSD results . . N2 in 5s4pld basis at R = 1.09898 A ...... N2, pVTZ+ basis ................. Summary of N2 EOMIP ionization calculations . N2 EOMIP error ................. Ionization spectrum of ethylene, 1030 eV . Effect of orbital choice on ethylene 3ag ionization Ionization of 1,3 transbutadiene, DZP+sp basis 10.9 NH3 EOMIP error vs. experiment NHa, pVTZ+ basis ......... H20 EOMIP error vs. experiment . Ionization of H20 ......... N2 X'l9 and N+ X2Eg in pVTZ+sp C2 XI+ and C+ a2II in pVTZ+sp basis . basis . Ionization of C2 .............. Ionization of C .............. C3X2,II + C3X 1E+ in pVTZ+sp basis. C,/C; .. . . . C,/Cn (Continued) ............ Ionization of C3 X2IIg in pVTZ+sp basis, Ionization of X2II, C3 in pVTZ+sp basis, C5 Ionization, ROHF orbitals . C7 Ionization, ROHF orbitals . Cg Ionization, UHF orbitals . C2 in pVTZ+sp basis, QRHF orbitals . Ionization of C4 in pVTZ+sp basis, ROH Ionization of C4 in pVTZ+sp basis, QRO Ionization of Cg, ROHF orbitals . C8 Ionization, QROHF orbitals . QRHF reference UHF basis . F orbitals .. HF orbitals . Computational times for EOMIP and FSCC calculations . EOMIP errors, pVTZ+ basis, small closedshell molecules . Table 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 page 73 77 78 80 80 84 87 91 93 95 96 96 99 100 10.10 10.11 10.12 10.13 10.14 11.1 11.2 11.3 11.4 11.4 11.5 11.6 11.7 11.8 11.9 11.10 11.11 11.12 11.13 11.14 12.1 12.2 109 110 111 113 114 116 117 121 122 124 128 130 131 132 134 139 142 . . LIST OF FIGURES page Figure 2.1 The HN operator ............................ 2.2 HN elements required for CC EOMIPSD and FSIP equations . 3.1 The EOMIP equation diagrams . ... .. 8.1 Diagrammatic expression of ro+k terms . . 8.2 Diagrammatic expression of r~fko terms. . . 10.1 10.2 10.3 10.4 10.5 10.6 11.1 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 11.10 N2 ionizations in three MO sets . ...... 81 EOMIP and SACCI ionization spectra of ethylene .. 85 Ionization spectrum of 1,3 transbutadiene . .... 90 NH3 Ionization Spectrum ....... ................ 94 H20 EOMIP ionization spectrum . . ... 97 N2 and C2 potentials and ionization . ..... 100 Ionization spectrum of C3, pVTZ+sp basis . ... 118 Ionization of C5, ROHF basis . .... 121 Ionization of C7, ROHF basis . .. 122 Ionization spectrum of Cg ....................... 124 Ct Ionization spectrum, EOMIP calculation. . 127 C2 ionization spectrum, twodeterminant calculation . 127 C4 ionization spectrum, ROHF orbitals. . .. 130 C4 ionization spectrum, QROHF orbitals . ... 131 Ionization spectrum of Cg ...... .......... ....... 132 Ionization spectrum of Ce, QROHF orbitals. . 134 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 THEORY AND IMPLEMENTATION OF COUPLEDCLUSTER METHODS FOR DIRECT CALCULATION OF IONIZATION ENERGIES By Renee Peloquin Mattie August 1995 Chairman: Rodney J. Bartlett Major department: Chemistry The coupledcluster equationofmotion (CCEOM) method, formally equivalent to the Fockspace coupledcluster (FSCC) and coupledcluster linearresponse (CCLRT) methods for electron detachment energies, offers distinct advantages over the single reference coupled cluster (SRCC) approach when several electrondetached states are to be calculated. The method is capable of providing transition energies for several states at a computational cost less than that of an SRCC calculation for one of them. The equationofmotion ionization potential (EOMIP) method is implemented using a Davidsontype algorithm for extracting selected roots from a large symmetric matrix, making it possible to selectively search for the difficulttocharacterize "shake up" satellite peaks often found in inner valence structure, and several approximations to the ionization intensities are available, making it possible to simulate electron detachment spectra. The EOMIP method is calibrated against experimental results for several small molecules, where the outer and innervalence regions are explored, with special at tention to shakeup transitions. Results compare favorably with the coupledcluster singles and doubles (CCSD) approach in most cases. The method is also applied to the interpretation of the electron detachment spectra of the linear carbon cluster anions, C2 through Cg. CHAPTER 1 INTRODUCTION 1.1 Correlated Methods in Quantum Chemistry In the early days of quantum chemistry, the only ab initio method available was the selfconsistent field (SCF) approximation. This approximation to the electronic Hamiltonian neglects the interaction between individual electrons (electron correla tion), placing each one in an average field of all the other electrons. The approach yields a simple molecular orbital picture of atoms and molecules and, in many cases, yields qualitatively good results for the energy differences between states of a molecule or reactive system. The molecular orbital (MO) picture can provide simple, useful concepts of physical processes, such as reaction, excitation, and ionization. But in many cases this simple picture is not even qualitatively correct. And there are many physical processes for which it is not capable of providing any explanatory tools. The case in point is the photoelectron process. When electromagnetic radiation in the proper frequency range strikes an atom or molecule, an electron will be ex cited into the freeelectron continuum. Photoelectron spectroscopists then observe the absorption spectrum or the spectrum of the free electron kinetic energies. Koop mans' theorem[l] is an SCFbased approach which states that the detached electron is removed from a specific MO, and that the energy required is approximately the energy of that orbital. This approximation gives, in many cases, a qualitatively useful prediction of the principal peaks of the photoelectron spectrum, but it is incapable of making any prediction for the satellite peaks which are close in energy to these principal peaks. These socalled shakeup peaks are, in the MO picture, the result of a mixing of the primary holestate process (ejection of an electron from a previ ously occupied orbital) with a twohole, oneparticle process (simultaneous ejection from an occupied orbital and excitation of an electron from an occupied to a virtual orbital). In order to predict these states, it is necessary to improve upon the SCF approximation  to introduce electron correlation. The configuration interaction (CI) method has been widely used to introduce correlation. In the singlereference formulation, it is simple in concept. Moving an electron from an "occupied" orbital (creating a "hole") to an unoccupied or "virtual" orbital (creating a "particle") defines a single excitation. Doing this twice defines a double excitation, etc. The correlated wavefunction is defined as a linear combina tion of the reference function and all possible single, double, etc., excitations from it. The coefficients are the eigenvectors of the matrix representation of the full electronic Hamiltonian in the basis of all these functions. The lowestenergy solution yields the correlated energy and wavefunction corresponding to the SCF reference, and the higher roots correspond to excited electronic states. The "full CI" solution, as it is called, is the best possible wavefunction for a given basis set and satisfies the varia tional principle. In practice, however, the full CI problem can rarely be solved. As the number of basis functions (M) increases, the number of possible excited determi nants increases as a function approaching M!. Standard implementations therefore truncate the expansion space, typically to only a subset of single and double excita tions from the reference. But this truncation introduces a particularly troublesome error  the energy of the system no longer scales properly with the size of the system being studied. The percentage error in the calculated correlation energy increases rapidly as the size of the system increases, more than doubling each time the size of the system doubles. The error in the energy is reflective of errors in the wavefunction, and calculations of observable properties, which depend on the wavefunction, will not scale properly with size. Typically, these truncated CI methods cannot be applied reliably to systems larger than benzene. The manybody perturbation theory and coupledcluster (CC)[2] methods have been developed specifically to avoid the deficiencies of the CI method. These methods are "sizeextensive"; they scale properly with the size of the system. Implementations of these methods have become increasingly available in the past 10 years and have been shown to be very powerful, capable of predicting experimental results for many types of molecular systems to a high degree of accuracy. Secondorder manybody perturbation theory (MBPT2) calculations, which generally provide 80% or more of the correlation energy, require no more computational resources than do SCF calculations. Coupledcluster singles and doubles (CCSD) calculations are routinely performed on systems up to the size of benzene. With the increasing availability of supercomputer time, it is easy to consider such calculations on molecules up to twice that large. 1.2 The Electron Detachment Problem Since the mid60s, "Molecular Photoelectron Spectroscopy has been recognized as an especially unambiguous method for the study of the molecular electronic structures of substances in the vapor state" [3]. In these experiments, light from an intense source is directed into a molecular gas or beam. The resulting detached electrons are collected and their kinetic energies recorded. At one time, it was common to record the absorption spectrum as well. In some cases fluorescence studies of the gas are performed as well. From the beginning, theoretical methods have been an important tool in inter preting the photoelectron spectra. Koopmans' theorem predictions of ionization en ergies made it possible to identify principal peaks, and later CI singles and doubles (CISD) and electron propagator approaches were used with some success to iden tify the satellites. In the late seventies, several approaches were used to predict the intensities of these principal peaks and their satellites[4, 5, 6]. Nakatsuji and coworkers[7, 8, 9] introduced the symmetryadaptedcluster (SAC) and symmetry adaptedcluster CI (SACCI) approaches to calculating the ionization energies and intensities of molecules. The Fockspace coupledcluster approach, introduced by Mukherjee and Lindgren[10, 11] and developed by others[12, 13, 14, 15] has recently provided a convenient formalism for a computationally advantageous approach to the direct calculation of energies. The CC EOM approach, closely related to CC linearresponse theory (LRT)[16], is based upon a coupledcluster reference. For the calculation of principal ionization energies, this theory can be shown to be formally equivalent to the Fockspace coupledcluster (FSCC) method. In this dissertation, derivations of the CC EOM and FS EOM methods for cal culation of ionization energies will be given, and the formal equivalence of these two methods will be demonstrated. The CC linearresponsetheory will be derived for the ionization problem, and expressions for calculation of transition intensities will be developed as well. Finally, the CC EOMIP method will be applied to several small molecules and to the linear carbon cluster anions, C2 through C9. CHAPTER 2 ELECTRON CORRELATION 2.1 The Schr6dinger Equation All electronic structure methods seek solutions or approximate solutions to the BornOppenheimer approximation to the Schr6dinger equation, HT = ET (2.1) E*Hq! E= T*T where H is expressed in Hartree units as H 2 E ZZ r + a r1, (2.2) Sia i>j which is the sum of a oneelectron and a twoelectron portion, h=H(one) = V ZrQ 2i i 2 = H(two) = r 1. i>j where i, j are indices over all electrons in the system and a over all nuclei, and Z" is the charge of nucleus a. Unfortunately, the kinetic energy term of the Hamiltonian _(( )1+ + taken with the twoelectron term, gives a secondorder differential equation in many variables, for which there is no general analytic solution. It is possible to solve certain approximations to the Schrddinger Equation. 2.1.1 The Independent Electron Approximation and Fock space All the methods to be discussed are based on a finite set of spin orbitals, ,((), oneelectron functions of the combined space and spin coordinate C = {r; w}. For molecular calculations, the spin orbitals are referred to as molecular orbitals and generated using an independent particle method in which the V2 portion of the Hamiltonian is replaced by an independent particle potential VF. The ultimate form of the 0 orbitals is generally arrived at by a selfconsistent field method, such as the HartreeFock theory. In the HartreeFock approximation, the approximate Hamilto nian is referred to as the Fock operator. F = h+VF (2.4) H = F+(V2VF) The basic component of an Nelectron wavefunction is a single Slater determinant 1 S 1~ A (H OP (P) = 7 [l1v(,1)2v(2) ...N,(CN)1, where N can be any number from 0 to the total number of molecular orbitals available, and the index v distinguishes among the (nM possible choices of N orbitals from a set of M. G is the trivial case, with no electrons chosen, and is referred to as the vacuum state. The Slater determinant is an antisymmetrized product of N oneelectron func tions, ensuring that fermion statistics are obeyed the sign of the wavefunction will change if the C coordinates of any two electrons are exchanged. Any manyelectron wavefunction within this truncated space can be represented as a linear combination of Slater determinants. The finite space spanned by the Slater determinants is Nv which is generally a subset of the space spanned by the eigenfunctions of the Hamil tonian. The Slater determinants span the eigenspace of the truncated Hamiltonian  the Hamiltonian projected into the space of the Slater determinants, D NMD (J f HA) M* Nv Mjp Nv Mp. 2.1.2 Second Quantization The second quantized notation is a convenient one which will be used throughout this work. The Hamiltonian can be written in secondquantized form as 1 H = h,,ptq + E (pq? rs)ptrqts (2.5) pq 4 pqrs where hpq= (1) 1V EZar) ,Qd dai (2.6) (pql rs) = / (1)0,(2)(1 P12)r (1,(2)dTidr2dld 2. (2.7) The creation (pt) and annihilation (p) operators represent the creation and annihila tion of orbital p. The fundamental relations between these operators are [p qt] = ptqt +qtp 0 (2.8) [p, q] = pq + qp O (2.9) [pt q]+ = (p Iq). (2.10) When the orbitals are orthonormal, equation 2.10 becomes [pt, q] + py Any Slater determinant can be constructed through the action of creation opera tors on the vacuum determinant (I, = )). qt...pt )j p()...q((N) VN! The properties of these operators ensure the Pauli principle is obeyed tpt = 0 ptpt = 0 ptp ) = 0; a given spin orbital cannot be occupied more than once in the same determinant. If one Slater determinant is chosen as a reference, the operators can be separated into occupied or hole (it, i) and virtual or particle (at, a) creation and annihilation operators. Any Slater determinant can created from a given reference determinant 10) by first annihilating zero or several occupied orbitals (creating zero or several holes), then creating 0 or several virtual orbitals (creating zero or several particles), for example ati 10) represents the excitation of a single electron from occupied orbital i to virtual orbital 0a, while the operator string atji operating on the same reference excites one electron while annihilating another, reducing the number of electrons in the wavefunction by one. A thorough discussion of the nature of secondquantized operators can be found in the first chapter of Jorgensen and Simons' work[17]. In order to simplify further discussions and diagrammatics, it is useful to write operators in normal ordered (or normal product) form[2, 18], represented by placing "{ }" around the product of operators. The normal order of a product of operators is that ordering of the operators whose expectation value with respect to the reference, 10), is zero. In practice, this means "moving" all particle annihilation operators (a) and hole annihilation operators (it) to the right, using the fundamental anti commutation relations given in equations 2.8 through 2.10. Interchanging any two singleparticle operators within a string of operators results in a change in sign of the normalordered product, {pqsrt} = {pqrts} = {prtqs} = {rtpqs} {jti} = {it} = ift {ati} = ati. The last of these, ati, composed only of hole and particlecreation operators, is the only type of normalordered string which may act on the reference 10) without destroying it. With normal order defined, the contraction between any two operators, designated as (pq)c or pq is defined as the difference between the product and the normal order of the product pq =_ {pq} + . Using the fundamental anti commutation relations, the contractions can be evaluated as pq = ptqt = 0 ijt = atb = 0 itj = 6ij abt = ob. Another common normal ordering is that with respect to the vacuum state I). In this ordering, all orbitals are "virtual", and the second quantized operator in equation 2.5 is already in normal order. This is not, however, the case with the normal ordering used here. The normalordered Hamiltonian, HN is that part of the Hamiltonian whose expectation with respect to the reference 10) is nonzero HN = H(OIHIO) 0 (2.11) = H ESCF = fN +WN = HN(one) + HN(two) fN = FpqPtq pq WN = (pq rs)ptqtsr. pqrs (2.12) (2.13) Writing HN in terms of normal ordered operators makes this partitioning into F and ......b i /\.X b  4 c/ "i k  x aV X ++ ++ Figure 2.1: The HN operator nonF pieces possible, and is convenient when it comes to developing diagrammatics. A diagrammatic representation of these HN elements is given in figure 2.1. 2.2 CoupledCluster Theory 2.2.1 The CoupledCluster Operators The coupledcluster theory defines a wavefunction ICC) which is an eigenfunction of the full Hamiltonian H CC) = E CC) (H (01 H 0)) ICC) (2.14) (2.15) = HN ICC) = Ecc ICC), where Ecc is the correlation energy obtained from the application of coupledcluster theory. The CC wavefunction is defined in terms of excitation operators T as ICC) = eT 10) (2.16) The operator T is defined as T1 = tati Ti at' 1 i3 ab T3 tiatbtctijk, "T ijkabc etc. The inverse of this exponential is simply exp(T). Therefore eTH CC) = eTEeT 0) = E O) E = (01eTHeT 0). The correlation energy, Ecc = E ESCF, is Ecc = (0 eTHeT 0) (0 H 0) = (O eTHNeT 0) (2.17) 2.2.2 The H Operator The quantity exp(T)H exp(T), referred to as H, shows up again and again in coupledcluster based methods, so it is worth looking at it in detail. H, expanded in terms of T in the BakerCampbellHausdorff expansion, is S = H+HTTH+ HT2l T2HTHT+... (2.18) 2 2 1 1 1 = H+[H,T]+ [[H,T], T] + [[[H, T], T], T]+ [[[[H, T], T], T], T]. The commutators can be written in terms of their connected and disconnected parts, and it can be shown that the disconnected parts cancel, allowing H to be expressed as a sum of only connected pieces, in which H is contracted to any T operators in the expression H = H + (HT) + (HT2)C + (HT3)C + (HT4)c (2.19) 2 3! c(4! + = (HeT) (2.20) The connectedness of the H elements is what establishes the sizeextensivity of the coupledcluster method. KQ+_ . 4K7= \A v ..\A. + V+ V+ 0~~. + = H +/v + .V + + w /VA/ + F ... Figure 2.2: HN elements required for CC EOMIPSD and FSIP equations     A \V~   v MR = I# = 2.2.3 The CoupledCluster Equations Similar to H is HN = eTHNeT = (HNeT)c, which appears in the coupledcluster energy expression in equation 2.17. Using this formulation, we have Ecc = (0 (HNeT)c 10) (2.21) Since HN 10) = E 0), it is clear that (@INlO0) = 0 (2.22) for any state V orthogonal to the 10) state. This fact is used to define the T amplitudes. If a) = ati 0), (2.23) ab = atibtj 10), (2.24) then the singles and doubles equations are defined by SHN 10) = 0 (2.25) (ab HN 0) = 0 (2.26) Figure 2.2 contains the diagrammatic expressions for the H terms which are necessary for the methods to be discussed in this work. For example, the algebraic expression for the term R (= (il f a)) is (il fl a) + Et (iml lae). me 2.3 Configuration Interaction 2.3.1 The CI Operators In configuration interaction, electron correlation is introduced by constructing a wavefunction from all possible Slater determinants with the same number of occupied orbitals as the reference. This means each additional determinant is related to the reference by an excitation operation ati represents a single excitation, while atbtji represents a double excitation. The CI wavefunction is represented as ICI> = (1 + C) 0), (2.27) where C is defined as C1 + C2 + 3 + . C1 = E cati ia ij ab Sijkabc 2.3.2 The CI Equations The CI energy is calculated from H ICI) = E CI) E = (0 H CI) = E((0I 10) + (0 C10)) EcI+ESCF = (0 (HN + (0H O))(1 +C) 0) Ec = (0 HN(1+ C) 0) and the C operator satisfies (a H(1 +C)l0) = EcICO ab HN(1+C)0) =EclCf \ij etc., (2.30) which can be solved via a matrix eigenvalue approach. In the CI approach, there is no BakerCampbellHausdorff expansion to force the connectedness of H to the wavefunction operator C. (2.28) (2.29) CHAPTER 3 THE COUPLEDCLUSTER EQUATIONOFMOTION METHOD 3.1 The EOM Equations Beginning with the exact coupledcluster ground state ICC), an EOM wavefunc tion is created using a linear operator Q, the full set of which spans the space of excitations from the SCF reference, and therefore the space of eigenfunctions of the secondquantized HN of the appropriate number of electrons. A particular EOM state will be referred to by the index IC. HN JCC) = Ecc CC) = EcceT 0) HN IKcc) = (EK ESCF) lCcc) = (E ESCF)' ICC) = (E' ESCF)QKT 10). (3.1) (3.2) When Q1 is limited to contain excitationtype operators only, it commutes with T, and we derive the expression HN \KCcc) QHN ICC) wKe7Tn 0) eW Te TQ 10) eT W )ceTQ 0) = (E E) ICcc) = HN, K] CC) [= HN, Q eT 0) [HNeT, K] 10) = [eTHNe, ] 0) w OK 10) = [HN,QC] 10) (3.3) where w = Ec E and HN eTHNeT = (HNeT)c. Projection of equation 3.3 on the left by Ih), the orthogonal complement to the reference function 0), and by 10) itself yield Wc (hlQ 10) = (hi [IN, Q] 10) (3.4) This equation defines an eigenvalue problem which is solved to yield the energies (eigenvectors) w (Q). The second term of the commutator, QfHN, contains only uncontracted terms. Dividing the commutator into contracted and uncontracted terms gives [HtN,] = (HQK)c + (IN QC), (fIN K), S(HNQ) + (NC)U (QKHN)U WO 10) = (Hf )C 10). (3.5) The commutator in the EOM equations has the effect of removing all the uncontracted terms from the equations. The EOMIP equations become w: (hl R 10) = (hi (RHNQ)C 10), (3.6) an eigenvalue problem in terms of the transformed Hamiltonian HN. The expansion of the HN terms is limited to a maximum HN excitation level of (Ln 1), where Ln is the highest level of excitation included in QK, and again by the maximum number of upgoing lines in f available for contraction. This will be 2Ln 1 for the IP problem and 2Ln for excitation energies. 3.2 SizeExtensivity and the EOM problem The EOM problem is similar, from a computational standpoint, to an ordinary CI problem except that HN is not symmetric and is not limited to one and twobody terms. From a theoretical standpoint however, it differs considerably from the CI problem, in that all the terms of the HNQ product are contracted. Within this framework, a perturbation theory which is connected to every order, and at any level of truncation of the theory, can be constructed[10] for the principal ionization energy and electron affinity cases (see appendix B). Since the HN operator, being based on the CC solution, is also connected, the EOMIP and EOMEA princi pal energies are connected, and the solutions are therefore sizeextensive. CIbased methods, being unconnected, do not offer this advantage. 3.3 The EOM Wavefunction For the EOM problem, is a linear operator, S = C TA (3.7) where 7, 10) = I,). (3.8) The operator 7, consists of a normalordered product of creation and annihilation operators which, acting on the HartreeFock reference P0), create new determinants. For excited wavefunctions, the set {r,} contains the identity operator, as well as excitation operators ati. atibtj, .... For the electrondetached case, {r,} contains i, atij *, operator products which will create N1 electron determinants. Because the equationofmotion problem is treated as an eigen problem of the HN operator, the solution wavefunctions, IK) = 10) are eigenfunctions of HN within the truncated space defined by the highest level of ex citation in {T1}. But in general, we are interested in wavefunctions of the Hamiltonian H, and so must go back to the expression IICcc) = eR 0) ), (3.9) an approximation to an eigenfunction of H. If we define a matrix C with vectors C' corresponding to the C' which appear in the definition of Q", then Q ( the entire set of QK solution operators) can be written as S= rC / C2 C2 C' \ S T2 ) 2 2(3.10) C1 C2 CI C O the solution space for a given EOM problem can be represented as (01HTC 10) = wC (pj HN I C = wC C 1(IH NIj)C = W Dt(AlHNI>p)C = w. (3.11) The matrix Dt C1, defines a new operator Y = rD. D and Y can also be used to represent the solution space of the EOM equations: Dt (pI HN I) = Dw. (3.12) It is clear that these D vectors, the lefthand eigenvectors of the non Hermitian HN, yield lefthand eigenfunctions of HN, Dt (01 THNT0) = Dtw (3.13) (01 Y'HN = (01 lEtw (3.14) and therefore the projection of the fullspace lefthand eigenfunction into the reduced space, (0I YteTHN c (Ol Yte TE, (3.15) just as the righthand eigenvectors, C, yield the projection of the fullspace righthand eigenfunction into the reduced space, HNQy'eT 10) E cKeT 0). (3.16) It should be noted that, since the untransformed H is a Hermitian operator, its right and left hand eigenfunctions are of course merely the transpose of one other: H R)R = Ei Ri) (LiIH = (LilE, (LI = {Ri)}t. Up to this point, there has not been much attention paid to the T and f l/Y: operators actually used in practical calculations. Although it is easy to show that, in the full expansion of the operators, the lefthand side ground state coupled cluster wavefunction and the transpose of the righthand side wavefunction are identical up to a normalization constant, (0 (1+A)e' = 0(0eTteTlO) (OleT, qualitative differences between the two wavefunctions are apparent in the usual CCSD approximation (John Stanton, private communication), because of the truncation in the operator space. 3.4 The Ionization Energy Equations In an ionization potential problem, the reference state 0) is an Nparticle function, whereas the solution Q' 10) is an (N 1) particle function. The operator f: must therefore contain terms which annihilate on electron in addition to any excitation they might also perform: Q0 = c + cO +... (3.17) c=i + 1Eciaj +. (3.18) i ija I>)=Qa IO) =1 ) ci)+ E c ). (3.19) i 1i 2 ija i 3 I The IP coefficients cl, ccf are indexed by C for convenience. The similarity of this notation to that for the coupledcluster amplitudes tq and ti implies that an electron is being excited out of orbital i and into IC, but this CK does not represent one of the reference orbitals. It may be taken to indicate that the electron has been excited into an "unbound orbital," but this interpretation is not important to the solution of the EOMIP equations. It is important to remember that, when constructing EOMIP equation terms, no contractions may be made to the "freeelectron" index. Sl S2 S3 D1 D2 D3 D4 D5 D6 Figure 3.1: The EOMIP equation diagrams. The EOMIP singles and doubles (EOMIPSD) equations, for which the terms are shown diagrammatically in figure 3.1, can be written as (SINGLES) KC = (n f i) + (n i/a) c" n na 2 cE nmfia) (3.20) 2 nma (DOUBLES)W C = (naizj)c E (<( J f1i) ca (n f j) c@) n n + (al B) ci~; + nmijy) b nm + ((na ibj)c (na bi)c:) nb 1 (nmllbe) t)Kb (3.21) 2 nmb e ) Where the notations f and ( [ ... ) indicated, respectively, one and twobody H elements. The final term of the doubles equation involves the 3electron HN term (nmalfllijb) (figure 2.2). These equations form an eigenvalue problem which may be solved by any number of diagonalization procedures. 3.5 The EOMIP Wavefunction The EOMIPSD wavefunction expressed in equation 3.19 is an eigenfunction of the HN operator projected into the reduced EOMIPSD space. The full wavefunction, eTQ" I0), is an approximation to an eigenfunction of H, and can be written IIcc) = c i) + ( cE +P(ij) ) i g 2\ iia ija I + (Pi/jk) E c(t + (tt tat)) + (ij/k)P(ab) b ab 2 ijka ijkabc l+ (P(i/j/kl)P(a/bc) c + P(ijkl)P(abc) 1 cit btc 2 ijklabc ijklabc +P(ij/kl)P(a/bc) E Ci(t(tI kt tkJ)) abci + (3.22) ijklabc / CHAPTER 4 THE FOCKSPACE COUPLEDCLUSTER METHOD 4.1 FockSpace Theory The Fockspace coupledcluster method (FSCC) uses an exponential operator to describe electron attached, electrondetached, and excited states relative to one single determinant reference. Like EOMIP, it gives direct differences between the correlated reference energy and the energy of the calculated state. The FSCC approaches to the calculation of single electron ionizations and attachments are formally equivalent to the EOMIP approaches to the same, as will be demonstrated for the IP case in section 4.2. The Fock space coupledcluster method is based upon a valence universal operator Q that, when acting upon the appropriate model space, yields the desired wavefunction(s). If we begin with a single determinant function, such as a Hartree Fock wavefunction, a model space can be defined by annihilating n electrons from occupied orbitals (also referred to as creating holes) and creating m electrons in virtual (unoccupied) orbitals (also referred to as creating particles). This is called a rank (m, n) model space, one with particlehole rank (m, n). The rank (m, n) part of the problem is also referred to as the (m, n) sector. For example, the usual coupled cluster problem has particlehole rank (0, 0), while the single electron detachment problem is of rank (0, 1). An active space of occupied and virtual orbitals is chosen in which electrons can be created and annihilated to produce model space determinants. The operator Q is defined by {e)}, where the braces, {}, denote normal ordering. Normal ordering of a product of operators prevents the operators from contracting with, or operating on, each other. For this reason, model spaces of different hole particle ranks are decoupled, and the solution of the equations is hierarchical, start ing with the usual coupledcluster equations, and building solutions to sectors with successively higher rank (the sum of m and n) upon the solutions of lower sectors. The valence universal operator f2 can be written 9tot= {ex} (4.1) where X = S(o'O) + S(o) + S'10) + .. Because of the properties of normal ordered products, Qtot can be written as a product of terms 9tot = {(0,0) (01)* ...} (4.2) where Q(mn) = {exp(S(mn))} = 1 + S(mn) + {(S(mn)2 + S(mn) is always chosen so that, when acting on the (m, n) model space p(m,~), it does not cause scattering between model determinants. That is, P(m",)S(m,")p(mn) = 0, or S(m,n)p(m,n) = Q(m,n)S(m,n)p(mn), where Q(m,n) is the complement to p(m,,). The (0, 0) sector, as was noted above, corresponds to the usual singlereference coupled cluster problem. The operator exp(S(O',)) is identical to the usual eT operator encoun tered in the coupledcluster equations. Because TtoT contractions are impossible, imposing normal order does not change the operator. The operator 2tot has an inverse, tot1 : (0o,)1 (01)1 ...}. (4.3) In the (0, 0) case this is f(0,0) = eT but in the higher sectors, the form becomes more complicated because of the normalordering constraint. These terms can easily be derived using a Taylor series expansion of the inverse of Q(mn), or by other approaches, but are not important to the following discussion. To begin, an Nelectron single determinant reference 10) is chosen, which may be an SCF wavefunction. The reference wavefunction defines a complement, Q, which contains all Nelectron functions orthogonal to the reference. The singlereference coupledcluster problem is solved, yielding HN, which is an effective Hamiltoimal for the (0, 0) sector. This is a canonical (eigenvaluepreserving) transformation of .e Hamiltonian, chosen so that it is the eigenoperator of the reference function, yield ng the desired coupledcluster energy (Ecc) as its eigenvalue. The coupledcluster equations can now be expressed in a compact form. (Ql HN 10) =0, (4.4) and the equations are solved for the T amplitudes. In further discussions, the following quantity will be useful: of(m) = 11 Q n(,k)(m,n)l1 (4.5) old While Q(mn) contains all the contributions from all the Q's of different sectors, the old projection of equations by p(m,n) and Q(m,") limits the contributions to the energy and wavefunction expressions to only certain lowerranked sectors The (0, 0) sector contributes to all higherranked sectors, and is the only one contributing to the (0, 1) and (1, 0) sectors. The (1, 1) sector has contributions from all lowerranked sectors, while the (0, 2) has contributions only from (0, 1) and (0, 0). For a given sector, an effective Hamiltonian is also defined, H(mn), defined H(m,n) = ~(m,n)H Nf o (ldQ(,n), (4.6) eff Nol (46) which will prove useful in the following discussion. Now, with the model space and the Q operator defined, we can begin to define the (m, n) problem. While H(OO)PO,o = EccP(o,O) works fine for the onedimensional problem, it is eff clear that something more is needed to deal with the multideterminant model space. It is necessary to introduce C(m,n), or C for short. The matrix C causes a linear transformation of the modelspace functions. If p(m,n) includes a set of functions pi), then a linear transformation of these functions results in another set ycm,n) = p(m,n)c(m'fn) or IpIC)(m'n) = E piP)m,n)C (4.7) The Cth wavefunction of the sector then is defined to be QtotC mKn)P(mn) and should satisfy HftotP(mn)C_"(n) Q ttP(m,n) cm,n)w(mn) If the full C matrix is used, the result is H totP(m'n)C(m'n) = totP(m'n) C(m,n)W(m,n), (4.8) where W(m,") is a diagonal matrix containing the eigenvalues. This equation, like the analogous (0, 0) equation, can be put into an "Hef form, H(m n)p(m,n) = ,tot HNtottp(mn) (4.9) eff where H( n) p(m,n)C(m,n) = ntotl~tot P.(mn)C(m,n) .n) Heff has the property Q(m,n)H(mn) p(m,n)C(m,n) = Q(m,n) p(m,n)C(m,n)(m,n) eff = 0 (4.10) This fundamental result, unfortunately, is not nearly as useful as the similar equation from the (0, 0) sector is in computing the solution to the coupledcluster equations. The complicated form of Q1 makes for a large number of algebraically messy and computationally difficult terms in the expression Q(mn)HefCP(m,n). Upon conver gence of Q, (m,n) p(m,n),C(m,n) = (p(m,n) + Q(m,n) n)H p(m n)m,) p(m,n (mn)H(m,n) C(m,n) eff eff (4.11) Since the matrix C only mixes the P(m,n) functions among themselves, and cannot promote p(m,n) functions to the Q(m,n) space, it is equally valid to say H(m,n)p(m,n) = p(m,n)H(m,'n)p(m,n). (4.12) eff eff Combining this with the identity HNQ = QHeff, we arrive at an algebraically and computationally more tractable form, Q(mn)(HNQ f P(mn)Heff)P(m,n). Once Q has been determined from the solution of this equation, P(m'")H(m ,)p(m,n) can be diagonalized to yield the C and the energies eff diagonalized to yield the Cr and the energies wk Hf = (ii f j) + S (mi j) m S; (m f Ie) + E Sn mnllje) (4.13) me mne S=(i lj) S(ml fj)+ Sj (m le) m me SiE S mnije) + E SPHie (4.14) 2 mne m 0 = (iajk) E S majk P(jk) S.m( I fk) m m + SE k(al le) P(k) S malej+ E S (mn jk) e me mn 1 E Sm, t' (mnlbc) SE 'H e!f (4.15) mnbc m An important feature of these equations is that the ionization energies are invariant to the choice of active space. This is because the S1 amplitudes and the C coefficients are serving essentially the same purpose; they mix the model space functions. If all occupied orbitals are made active, there will not be any S, amplitudes, and an n x n matrix of C coefficients. If the number of active orbitals is decreased, there will be a decrease in the number of C coefficients and a corresponding increase in the number of S1 amplitudes. In fact, it is easy to show that there is a simple relationship between C and S for the two cases. The proof of the invariance of the energy to the choice of active space can be found in appendix A 4.2 Comparison of EOM to FSCC When FSCCIP and EOMIP calculations are performed using the same active spacethe complete active space, consisting of ionizations from all occupied orbitals certain correspondences can be shown. 4.2.1 Definitions Q2(o) is the usual CC wave reaction operator eT. p(O,1), written P for convenience, is a projector over the model space, which contains all (N 1) electron determinants obtained by annihilating a single occupied electron. Q(',1), or Q, is the complement to P (1 P) in the (N 1) electron space. (0,1) = {exp(S)} =1 + S + {S2} + .* is the solution to the FSCC equation for the (0,1) (IP) sector. (o01)1' satisfies Q(o,1)1Q(o,1)p = .Q(O,1)Q2(,1)1P = P. H = (,1y)'Q((0,0)'HQ(O,O)Q(O,1) = Q(0,1)1ftQ(0,1) 4.2.2 FSCC Method The Fockspace method defines a model space consisting of several (N 1) elec tron determinants and a projector P = p(,1) over this space. An electrondetached wavefunction is defined as Q(0,') IK) E D(o,')CK 0) = Q(0,1) C~i 10). i P can be represented as a sum over these IC) states P = 1K) (KC) c. K The Q(o,1) and C' operators are defined so that RQ(O0,) /K) = EQ(o,') C) (4.16) which leads to a useful definition for Hetf, H P = f(or'1(O)P; (4.17) H ff K) = (o1)l (H Ecc) Q(O01) K) = (E Ecc) 1K) = wI/>), (4.18) which has the property QHff /K) = QEK K) = 0 (4.19) Heff = 2(o0,)lHQ(O,')P = P2(o,1)H (oP. (4.20) Heff acting on any representation of the model space can be diagonalized to yield the entire spectrum of ionization energies wn CHeff C 0) = w 0)} (4.21) This operator eigenvalue equation is easily replaced by the matrix eigenvalue equation, C'HefC = 4.2.3 EOM Method The diagonalization of H yields coefficients C and eigenvalues w. The space of eigenvalues and eigenvectors, and hence the H matrix itself, can be partitioned into principal and shakeup portions. The principal IPs correspond to the solutions of the Fock Space method and can be labeled as the P portion of the solution space. ( PP HpQ ( Cp CPQ ( CPP Q up 0 HQP HQQ Q QQ CQP CQQ 0 WQ or ( HPPCPP + HpQCQP HPPCpQ + HPQCQQ C_ PPWP CPQWQ HQpCpp + HQQCQp HQpCpQ + HQQCQQ CQPWP CQQWQ " The PP portion of this equation can be rewritten as Cpp (THpp + HpQCQpCpp) Cpp = Wp. Since up has been chosen to match the solutions to the FS equations, there is quite naturally an identification: C'Hf C = Wp = Cpp (Hpp + PQCQpCp) Cpp, which will lead to Hgff = (Hpp + HPQCQpCp) as long as the model functions chosen for both methods are identical. 4.2.4 Comparison The eigenvalues and eigenvectors obtained from the diagonalization of the effective Hamiltonian Hef are the (PP) blocks of the respective matrices obtained from the diagonalization of the full matrix H Hef = Hpp+ (HS)p (4.22) = Hpp + (HpQS)c (4.23) = RPP + HPQCQPCpp (4.24) SQP = CQpCpJ (4.25) CQP = SQPCPP (4.26) This means that the eigenvectors of principal IP states obtained from an EOM calculation can be directly compared to those from the FSCC calculation of the same states. The Cpp vectors obtained in both examples are identical. CHAPTER 5 THE ELECTRON PROPAGATOR The electron propagator[19], also called the one electron Green's function, was one of the first improvements made on Koopmans' theorem for the calculation of IPs and EAs [20], and is still in use. Propagator methods are direct methods, yielding energy differences rather than state energies. Unlike the methods described so far, which are based upon correlated wavefunctions for the reference and final states, most Green's function methods provide transition properties rather than wavefunctions for the final states. And it is feasible to begin a correlated calculation of ionization en ergies with an uncorrelated, single determinant reference, if a suitably large manifold of annihilation/excitation operators, such as those used in the the construction of the EOMIP wavefunction, is chosen. However, the symmetric form of the electron propagator requires that the manifold include, in addition to IPtype operators, i and atij, the EAtype operators at and atbti be included as well. Energy differences will correspond to the difference between a correlated reference and a correlated final state, neither of which wavefunction is, in general, directly constructed. Despite these apparent differences, there are similarities between the results of electron propagator and CI calculations. More recent formal developments have bridged the theoretical gap between Green's function and coupledcluster methods, and propagator researchers have independently arrived at equations which are iden tical to those derived from coupledcluster EOM theory. The general form of Green's function in the energy representation is + (0 A k) (k) lB0) (0 B 0 k) (k\A 0) ((A; B))E = lim 1 + EklA0 (5.1) E 77+0+ Eo Ek + E + ir k EkEO E where the states k) are eigenfunctions of the full Hamiltonian, H. The electron propagator is a matrix G of elements GP = (ONpEkN ) k Ni) (Ok1 q 0O) (0 q ikN+) (kN+llPt .0) Gq = lim + (5.2) + k E EkN + E + i77 kN EkN+ Eo + E i7I where the IkN1) and IkN+i) functions are, respectively, electrondetached and electron attached functions. The poles of this expression occur at Eo EkN_ + E = 0 and Ek~,+ Eo + E = 0; the residue at a particular pole yields the transition moment (0NI P IkN1) (kNi q 0) or (01 q kN+) (kN+ I pt 0). This form of the Green's func tion is not directly useful for obtaining energies or wavefunctions, so another approach is used which calculates the energy differences and transition moments directly, with out first determining 10), jk), or any of the state energies. 5.1 Superoperator Expression of the Electron Propagator Propagator equations are most compactly written in terms of superoperator no tation. For any operator C, the superoperator C acts on other operators, and a new "binary superoperator product" is defined between any two operators: CD = [C,D]_=CD DC (5.3) ID D (5.4) (CID) (0l CtD 0}) (0 DCt 0), (5.5) where the sign in equation 5.5 depends on the number of creation and annihilation operators in C and D. In this notation, the Green's function can be written as ((A; B))E = (Bt(Ei + H)'A). (5.6) Inserting the spectral resolution of the identity, i = Ft (FtlFt) (Ft (5.7) = E Ft) (FtlFt)k' (Ft (5.8) kl into the expression for the propagator gives ((A;B))E = (Bt (Ei + H) A) = (BtIFt) (FtFt)1 (Ft (Ei + H)' Ft) (Ft Ft) l (FA) = (Bt Ft) (Ft (Ei + ft) F') 1 (Ft A), (5.9) which are the working equations from which propagator approximations are derived. The expression (E + Hft) is called the superoperator resolvent. The poles and residues of the resolvent give the energies and transition moments. For the oneelectron propagator, propagator matrix elements become Gp,(E) = ((pt; q)) = (pl (Ei + I)l qt) = (pIF)(FI (E + H) F)(FIq). (5.10) The manifold F can be separated into a, which contains all the oneelectron operators, p and q, and its orthogonal complement f, which contains all higher oddelectron op erators. Because the f portion of the manifold is orthogonal to the p and q operators, only the a a projected portion of (F (El + ^) F) contributes. This portion can be written as G(E) = (al (El + H) a) (al (El + i) f)(f (Ei +H) f) l'(f (Ei + /) a) = El (aI\\a) (alHtf) (fI (Ei + Ht) If) (f lHla) (5.11) 5.2 The Dvson Equation The Dyson equation is an expression of the difference between the inverses of the propagator G1(E) and the zeroth order propagator Go (E), the Green's function in terms of the zeroth order Hamiltonian, H0. G'(E) = Go'(E) E(E) (5.12) where Go'(E)p, = (E c)6pq (5.13) E(E) = (aljHa)cr + (alff) (f (El + R) f) l(flfa) (5.14) (pHq)corr = (pIIHq) ep6pq. In this formulation, E(E), the selfenergy, is the correction to the Koopmans' theorem approximation to the propagator given by Go(E). The self energy is itself divided into two parts, the energy dependent part E'(E) and the constant part E(oo) E(E) = E(oo) + '(E) (5.15) i(oo) = (alfa)corr E'(E) = (afl^t(fif(Ei+ )'(fI^a). (5.16) A great many approaches to the electron propagator are concerned with successive orders of approximation to the energydependent part of the self energy[21],[22], in combination with appropriate choices for the reference state 10)[23], and the most important portions of the manifold F. 5.3 NonDvson Equation Approaches The development, in the world of CI theory, of the Davidson procedure for con vergence of selected roots of a large matrix[24] spurred Baker and Pickup to develop a method for extracting the roots of the one electron propagator matrix rather than us ing selfenergy expressions[25]. More recently, Ortiz has experimented with the use of a highlycorrelated reference based on the coupled cluster wavefunction [26], a process referred to as renormalization of the ground state, in a procedure which maintains the symmetry of the electron propagator matrix. Snijders and Nooijen [27], using an RSPT analysis of the diagrams contained in the one electron propagator, showed that they can be represented using the coupledcluster wavefunction as the reference, with two major consequences. The symmetry of the resulting matrices is broken, and the IPtype manifold is decoupled from the EAtype manifold, decreasing the dimension of the manifold space to less than half for the IP case when typical modern basis sets are used. The resulting equations are identical to the CCEOM equations. One other Green's functionbased approach which gives results identical to the EOM approach is the coupledcluster response function approach developed by Koch and Jorgensen[16]. In this approach, a time dependent formulation is used to derive expressions for the linear and quadratic response functions of the coupledcluster equations. The coupledcluster Jacobian (the first derivative of the coupledcluster equations with respect to the cluster amplitudes) used in their work is identical to the coupledcluster EOM matrix. Transition matrix elements between the ground and excited states and between two excited states are derived from residuals of the quadratic response function. CHAPTER 6 OTHER METHODS FOR IONIZATION CALCULATIONS The many approaches to calculating ionization energies can be broken down into two main categories statebystate approaches which involve a separate SCF refer ence and correlated calculation for each state of interest, and direct methods which use a single reference (SCF or correlated) as the starting point for calculation of the energy difference between the reference state and several states of interest. Stateby state approaches can use any SCF or correlated (CI, perturbation, or coupledcluster theory) method in the calculation of differences. Among the direct approaches in clude two holeone particle CI[28, 6], the symmetry adapted clusterCI approach of Nakatsuji[29], the FockSpace Coupled Cluster approach[10, 13, 30, 14], and one elec tron propagator approaches[31, 32, 23, 33], as well as the coupledcluster equationof motion approach, which has been derived independently several times, and is known variously as the coupled cluster linear response function[16] and the coupledcluster Green's function[27] methods. 6.1 State Methods The conceptually simplest way to consider calculating energy differences is to calculate the individual energies of the states of interest and take the differences. The problems inherent in comparing SCF energies of systems with different spins (as in the common case of the closedshell singlet neutral reference and the doublet ionized state) are well known and tend to make the SCF energy differences not even qualitatively useful. Configuration interaction, perturbation theory, and coupled cluster methods may all be applied, and for small molecules which are welldescribed by a single SCF determinant, any of these methods might give a good result for valence ionizations. There are a few difficulties which can make the problem less straightforward, however. Modern SCF programs usually can be induced to converge on valence ionized states, and the fact that core molecular orbitals are generally well described as a linear combination of core atomic orbitals and do not change much when the occupation changes means that core ionized states can be converged as well. But when an innershell ionization is desired, there is no good way to characterize the excited SCF state. SCF programs will quickly collapse back to the valenceionized state. One remedy for this problem is to use a quasirestricted HartreeFock (QRHF) reference, in which an SCF calculation is performed on suitable reference such as the closedshell neutral, following which the orbitals are reoccupied and the Fock matrix recalculated. This nonHartreeFock reference is then used as the basis for the correlated (typically coupledcluster) calculation[34]. Because each state is handled as a separate case, there is considerable freedom in the choice of the reference. For approximate correlation methods, poorly chosen orbitals result in a substantial shift in the energy of the state considered. And if the orbitals are chosen well for one state but poorly for another, the lack of balance means the accuracy of the difference will suffer. UHF references lack spin symmetry. ROHF and QRHF references for highspin states will always be pure spin states, but the orbitals themselves may be unphysical. In some difficult cases, it is a matter of trial and error to determine which set of reference orbitals is best suited to the problem  the size of the correlation corrections for each reference is often the best guide. One more factor to consider is molecular geometry. Experimentalists may report vertical or adiabatic ionizations, or may be unsure. Proper vertical ionizations can only be approximated from the optimized geometry of the Nelectron species, while adiabatic ionizations involve finding the potential energy minimum for each of the N 1electron states. In most cases, a singlereference approach like CCSD(T), an approximation to the full CCSDT method, is significantly flexible to provide adequate correlation to all of the states in question. For calculations on certain openshell systems, however, it becomes necessary to consider a multideterminant reference. For example, the excited states of C2 accessible from ionization of the C2 ion[35, 36] include open shell singlet states. Attempting singlereference calculations on these states results in significant spin contamination, and the splitting between these and the related triplets (which can be calculated from a single highspin determinant calculation) are too large to ignore. For such states, the multireference coupledcluster (MRCC) method[37, 38] can be used. In this approach, two (or more) determinants serve as the reference function, and an exponential operator (analogous to the singlereference coupledcluster theory described in chapter 2) is applied to each reference in turn to create the correlated wavefunction. The number of amplitudes required grows as the number of unique determinants used in the reference; a typical application of the theory might require two or four reference determinants. Since higher and higher levels of SRCC theory are required to obtain better spinprojection, approaching a purespin result, the MRCC method gives substantial savings in cases where a full triples or even CCSDTQ approach might be required to gain the desired degree of projection. 6.2 Direct Approaches 6.2.1 TwoHole OneParticle Configuration Interaction Configuration interaction approaches can also be applied to the direct calcula tion of a spectrum of ionized states based on a single SCF reference. Martin and Davidson[6, 28] used a 2hlp CI approach in which the CI expansion for each of the cationic states is based on a closedshell SCF reference. The excitation operators used in creating the ionized wavefunction are identical to those used in the EOM method, but the reference is the SCF wavefunction. The ionized wavefunctions thus may be very different from the corresponding EOMIP wavefunction. Also, the hpCI en ergy will be the difference, between the correlated energy of the ionized state and the SCF energy of the reference state. This means that only the differences between the ionization energies are considered "good", and it is common practice to offset the entire calculated spectrum by a constant energy in order to match one theoretical peak (typically the lowest ionization energy) to its experimental value. 6.2.2 SpinAdapted Cluster Methods The Symmetry Adapted Clusters (SAC) method of Nakatsuji[7, 39] is similar in most respects to the CC method presented in this work. The important features of this approach are the explicit spin adaptation of the excitation operators and the projection of the final wavefunction by a spin projection operator. The result, though requiring fewer independent variables, has a much more complicated algorithmic ex pression, a necessity to customizing the program for each type of spin function, and a proliferation of terms to be coded. As a practical matter, the developer found it nec essary to truncate the exponential expansion to no more than quadratic terms. For open shell references, this is not quite equivalent to truncation of the CC equations to T2 terms as it includes some higher spinflip terms as a consequence of the spin adaptation. The SAC configuration interaction (SACCI) method[7, 8] is very much akin to the EOM approach, but using explicitly symmetryadapted EOMtype operators, and based on a SAC rather than a CC reference, resulting in a slightly different H, which is further altered by the neglect of certain contributions to H elements. In addition, implementations of SACCI have also employed configuration selection procedures similar to those used by CI practitioners to limit the size of the matrices that must be diagonalized, unlike the EOMCC implementation presented here. 6.3 Eclectic Approaches When the most straightforward methods fail, it is possible to approach a difficult state from an unusual direction. For example, inner valence ionization states can be accessible as excitations from an outer valence ionized reference, or as electron attached states from a diionized reference. Bernholdt [40], finding that the structures of the lowlying E states of the C+ were difficult to impossible to converge using QRHF at the CCSD level of theory as well as evidence of discontinuities in the 2E state, was able to use the FSCC(0,1) approach with the wellbehaved, closedshell C5 as the starting point, in order to examine portions of the potential energy surfaces of all three ionized states. CHAPTER 7 COUPLEDCLUSTER LINEARRESPONSE THEORY The formalism developed by Koch and Jorgensen[16] for coupledcluster linear response functions may be as easily applied to ionization potentials as to excitation energies. All that is required is to construct a perturbation operator capable of causing an "excitation" from the occupied to "freeelectron" space, and a new wavefunction of the appropriate form. In the following sections, a linear response theory formalism, based on the timedependent perurbation theory (TDPT) approach of Koch and Jorgensen, will be developed. 7.1 Perturbation Applied to the CC Hamiltonian An RSPT approach is used, with the usual coupledcluster quantities defined as the zerothorder solutions, and the full wavefunction expanded in terms of corrections to the coupledcluster T operator. Although there is a single zeroth order expression, a spectrum of solutions to the full Hamiltonian will be sought, the maximum number of which can be determined from the form of the corrections to T. The definitions to be used are: 1. HSCF, the zeroth order approximation to the Hamiltonian, is the usual Fock operator. The perturbation is a firstorder correction to the Hamiltonian, V1P. The full Hamiltonian becomes HSCF + 3VIP (7.1) 2. The zeroth order wavefunctions and energies are those from the coupledcluster and A solutions; K(o) ) E(o) HSCF JC(O) (LC (Lo HSCF 1 0) = CCO) = eT 10) SEcc = Eo(O) () E Ecc C(0)) (01 (1 + Ao)e = (o Ecc. 3. The T operator for the state labeled KI is written as the zeroth order T plus a sum of perturbative corrections: TIC = T(O) +Tl') +T(2 +"'" = T(O) + T( + T2) +* (7.7) (7.8) The full wavefunction can be written as IK) eT HF) (7.9) = e ()*e ((l+ 2) e+...)+) HF) r (7.2) (7.3) (7.4) (7.5) (7.6) which works out the zeroth order wavefunction plus a sum of perturbative cor reactions = 1K0)0+ K'C) + +c2)) = CC) + K~ (2) + ..j. = T)IK(O))  (T(2 (7.10) (7.11) (7.12) + T12") :K(O) 2", > etc. The A operator can also be expanded in a perturbation series Al = A() +A + ..., so that the full wavefunction is ( I = (01 (1 + A + A() + ..)e(T()+T a sum of the zeroth order wavefunction and perturbative corrections: C(l= (01 [A)+ ( +(+ A ')) T ]e ( = (0 [A ) A) T,) + (1+ A.))( ' 2))] ec c l 1 J (7.13) (7.14) (7.15) 4. The perturbation operator, VIP, must allow firstorder corrections to T and A to be defined as follows: TF)) tati + j tatibtj + ... ai ijab (7.16) ) T t(1)KiXi + t(1)"Xiatj +.*. (7.17) i ija S (1) AO (0) (, t (7.18) A( ') E(') t. (7.19) Here, as expected, is where the treatment of the IP problem is different from that of the EE problem. Although the T(0) amplitudes have nonzero contri butions only to the usual excitations from occupied to virtual space, a single excitation to a "freeelectron" function XK must be introduced in the first order to describe an ionization process while preserving the number of electrons. In the limit where the "escaping" electron has zero kinetic energy, this form of T() allows the the firstorder corrections to the wavefunctions, IC()) and (L() to be represented as (N l)electron eigenfunctions of the zerothorder Hamilto nian. Keeping in mind that in the zerothorder wavefunction coefficients to the r\(r')t operators are zero, while in the firstorder correction the coefficients to the \(tau)t operators are zeror e o, it is possible to treat the perturbation problem in a manner analogous to that of Koch and Jorgensen. So far, the only assumption made in this development of a perturbation theory is that the only processes initiated by the perturbing operator, VIP, are electron detachment type processes. This is justified in describing, for example, a photode tachment process where there is a sizable gap between the frequencies of light which will cause electron detachment and those which cause excitations. One does not expect excitations and ionizations to occur at the same frequency. 7.2 coupledcluster Response Theory 7.2.1 Time Dependent Perturbation on a CoupledCluster Reference With these definitions, the coupledcluster response theory for electron detach ment processes can be developed in a manner analogous to that of Koch and Jorgensen. Although the determination of the energies requires only the TF') amplitudes, the formation of the linear response function requires the calculation of the firstorder corrections to the wavefunction. Once the CC solution is found, an operator HN = eTx7 HSCFeT: is defined, which has right and left eigenfunctions 10) and (01 (1 + A)) HN 10) = (01 (1 + AO) HN = HSCFeT I0) = (0l (1 + A4)eT HSCF = The next step is to find the response of Projecting the Schr6dinger equation Ecc 0) (01 (1 + A) Ecc, Ecce' 10) (01 (1 + A)e Ecc. (7.20) (7.21) (7.22) (7.23) both the left and right eigenfunctions. d i ICC(t)) = H ICC(t)) dt (7.24) by (pI eTK and collecting terms by order in 3 leads to the following expressions for the time evolution of TC: dT = i(H + ) CC T iCC) = i(HscF + V'P) ICC); dt dt(o dt tCCO) dt Si (LO HSCF CCO) = i (Ol HlN 10)= 0; (7.25) (7.26) = [T dT() Ty CC) = dt HSCF T VIP CC dt (AIl Tk1) HSCF JCCO) + (fil IHSCF TQ) JCCO) + (f1AI VIP ICCO) = i [(L [Hl, T1)] I0) + (Fl VI CC)] ; (7.27) where (i = (1pi eTk HIN = eTr HSCFe The response of the A amplitudes, i(L(t) d eT H F ( dT1 dA)+ (HFi (1 + A, Ac 1 j+ dt dt) = (C(t) HeTK (7.28) = i (HF (1 + AK) eT(HSCF + VIP)eT, can be derived with the help of equation 7.25 = i (CI (HSCF + VIP) T, ICC) i (01 (1 + AK) eTc (HSCF + VIP) ICC) = i(01 (1 + Ar) eT [(HSCF + VIP), ,] eT' 10), d = i (01 (1 + A() [HN 7] 10) = 0, dt) dt = i(0 (1+ A() [[ftN,)r+] ,Tk], ]0)+i(0IA [HN,7 ] 10) +i(01 (1 + A()) eTc [VIP, ,] eTx) 10). dAC dt (7.29) (7.30) (7.31) 7.2.2 The Linear Response Function Following the development of Olsen and Jorgensen[41], as applied by Koch and Jorgensen to the coupledcluster reference[16], we define the operator perturbation operator, VI'. Assuming the perturbation vanishes at t = oo, it can be written as VIP = d VWe(iw+a)t, (7.32) 00 where a real, positive, and small a ensures that VP is zero at t = oo, and Vwe"t is assumed to reach a finite value at noninfinite times. If V" contains one constant frequency component Wb, this becomes VI = f dw F [6(w Wb) + ( + b)] V(w)e(i+a)t (7.33) = (V(wb) e( +")t + V() e(b+)) (7.34 which could describe the amplitude of a periodic electric or magnetic field of frequency Wb. The expectation value of an operator O(t) is (C(t) O(t) ICC(t)) and can be written as a perturbation expansion of C and \CC) (IOlCC) = (C o CCO) + (Lc 1o CCo)+ (co o CC1) + (C2 0 CCO + (C1 0 CC1) + (C0o 0 CC2) +., (7.35) (7.36) Or, in response function notation, (C 0 CCo) + dw, e(iw i ((0; VW )), fooD + rO dwl dw2 e( Q')t e(i2+)t ((0; ; V^;V))L;, + (7.37) 2 oo 00 eL So that (v1 oCCO ) + (C OO Cc1) = J & e(r'l+a' ((O; V')),,v (7.38) The t and A amplitudes can also be written in terms of their Fourier transforms: (t = dw X (wL + ic) e(i''x+)t A( = dl Y (w1 + ia) e(iW1+a)t. K OO (7.39) (7.40) Substituting these expressions into equations 7.27 and 7.31 gives X )(w + ia) Y() (w + ia) h\U(T "U+ = (A + (wi + ia)I);'J = + ('(wi + ) I X(w + (A + (w, + ia) I)2 AP, = (I [fH ,,] 0), '.1 = (vjeT )V'ei 10), F. = (O(1+A 0) [[fN,T,], ,7] 0). Now the linear response function can be written in terms of the Fourier transform of the firstorder perturbed wavefunction: ((0; V)), = (wi) (I o C) + C (L [O, r,] CC) X (wi) (7.43) = Z ([VWTp] CCO) + F,+EZ(A+Wl),l x ,u v where (7.41) ia)) x (7.42) E(A wI) (a O CC) + E (K0 [0, 7,] cco) E(A + wlI)U (7.44) v 7.2.3 Energies and Transition Intensities The transition energies are found at the poles of the linear response function, which are found at +/ the eigenvalues of the matrix A. This matrix, A, = (AI [N, ,] 10), is exactly the matrix derived in the coupledcluster Equation of Motion method, by a completely different route. The Tk operators are the EOMIP wavefunctions. There is no quantity analogous to the T() amplitudes. What the LRT provides, from the EOM practitioner's point of view, is an unambiguous definition of the transition intensities. I = lim (1 k)((O;)) (7.45) where wk is a simple pole of ((O; Vt)), an eigenvalue of the matrix A. If (S'AS)m = Wnm = 6nm,n, (7.46) the bra and ket functions can be transformed: (kl = ESk2 (II Tk = 7TSLk. A The X(') and Y'() vectors can also be transformed: X1)(wi) = S'X()(wi) = S'(A + wI)IS Sl '(w1) = (W+wlI)1 V' (7.47) V L"i X (w)k (7.48) () 1 Wk) y ) = Y(1)(wi)S= rf'SS'(A+wlI)'S X(1)tS't SFtS S1 (A + wlI)S = V,"(W +wil) VWt ( + w,)t G (W +wi)1 (7.49) y ) = Gnk (7.50) (wk + w,) 1 n (W1 Wk)(W + W,) where 1/, = ECs' = < v1, cc), VW' = (EkjS, = (k CC ) Gmn = F,,SamSn = 0 (1 + AO) [[iN, rn], m] 0). The linear response function can now be written 0((O; V"') = Y () ( O CCo) (0 [H, I 0}X n m S' (n1 o c ) V Gmn (i 1O CCo) n (nJrn + 1) m (cl m) (l1 + n) V"i mZK 0 r CO O) (W1 Wm) and the transition intensity as VW, 0 cc, lim (w1 k) [ V' (IO CC) W k  (w, + wl) +v ((CO 1 [0,i \ C ) G OCCVO} 0 Gmn (+t W)O Co) ro X fEkV (7.51) where o oo (0ol [[O, 7n], Tk kccO)( i O CCO) 0k (( [ O, rk] CCO) ( = (0 (1 + AO)(Ork) 0) S (0 (1 + AO) (OrnTk)c 10)(n 0 I0), (7.52) n (wk (7+ 52) and rF k = (k V CCO) k+O I are referred to as the transition matrix elements and o = e Oe = (Oe1 )c V k = eTIC) k eT X = (V1W eT )c. In the usual propagator approach to transition expectation values, an expression of the form (01 0 Ik)(kl VWk 10) would be expected. The expression given above does not have the same symmetric appearance because of the difference between the right and lefthand solutions. The calculation of the transition expectation value is more complicated because of this asymmetry. CHAPTER 8 METHODS FOR TRANSITION MOMENTS IN THE EOMIP APPROACH For an experiment in which it is proposed to observe an absorption spectrum, the dipole operator F is the appropriate choice for 0. However, in the case of photoelec tron experiments, it is the electron momentum which is observed. The proper choice of O must be f, the momentum operator for the outgoing electron. In either of these cases, O is a onebody operator. The second term of FOOk contains the integral (0 (1 + AO) (Onrk)c)O) , which includes a product of two onebody excitation operators (m and Tk) contracted with a onebody operator (0). If the outgoing electron functions are orthogonal to the occupied and virtual orbitals, there is no contraction of these three operators which can survive the left projection by (01 (1 + A), so that term becomes zero and the transition matrix element becomes ro = (01(1 + Ao) (OTk) 10) (8.1) No matter what is used for 0, the nature of the outgoing oneelectron function must be known in order to calculate F_ k and FkV Freeelectron functions are dependent on the experimental method and not included in the basis sets for standard electronic structure methods; integrals of F and p involving these functions are not readily available. Therefore, it is usual to make some kind of approximation for both Fr0 k and Fk1O. The most common approximation is the socalled monopole approximation as developed by, for example, Martin and Shirley in 1976[4]. 8.1 The Monopole Approximation A symmetric Hamiltonian is assumed, and the dipole approximation to the cross section for photoemission by an incident photon field is employed: UDA O( (~1)1 [r. ( /(N)Ip11 (N))]2 (Ef), (8.2) where p(Ef) is the density of states in the continuum at the final energy. The authors proposed to determine the intensities of the shakeup satellite peaks relative to the re lated principal core ionization potentials, introducing two important approximations. 1. The principal ionizations are characterized by annihilation of an electron from a single orbital, and the satellites also have a nonnegligible contribution from this annihilation. Therefore the momentum operator, acting on an initial state dominated by an ionization from orbital i, is approximated as: S(ml P Ix) mtX 1f(N)) } itX 1f(N)) (i p Ix). (8.3) m 2. Both p and the momentum integrals (i i 'X) are approximately constant over the range of energies spanned by the principal and its satellites. This means that the ratio of the intensity of the nth satellite peak, I(n), to its prin cipal, I(0), is I(n) (Xnl jp i) 12 1 (ln itXn 0) 2 p(E,) I(0) I (xol I i) 2 1 (olit XO 0) 12 p(Eo) ^l~(8.4) .(oI itXO 107 (8.4) ( ol itZxo 10)2 And the calculation of the intensities is reduced to the calculation of transition density matrix elements for the principal and its satellites. Further simplifications are possible. In 1977, Martin and Davidson[6] employed a simplification which included only final state correlation and only one term from the resulting transition density matrix. For the ionized states, they used a (lh, 2hlp) CI calculation and argued that, for the principal and its satellites, only one contribution was important  that hole state which dominates both. So for a CI wavefunction of the form ,, = Cn, + ECna( (8.5) i ija they approximated the ratio of intensities as (I(n)/( oliW 10)2 (Cp)2 (8.6) ) / (0) 2 (20) (co) In a 1992 paper by Murray and Davidson[28], the same features were calculated with a multireference CI approach, and the same expression was used for the relative intensities. The authors report good agreement with experiment. 8.2 Transition Moment Expression for NonSymmetric Eigenstates For the non symmetric eigenvalue problem, there are two different approximations to the wavefunction in the full space. Equation 8.2 must be modified to recognize that the left and right handed wavefunctions are not identical. This results in UDA C( (hi) )1 [r. (L[N) p f(N))] [r. (0,(N)I )RP(N))] p(Ef), (8.7) which yields, under the same set of approximations as was made above, I(n) I (Xn.P I) I2 (OL ItXn ? n) (VI itXn OR) p(En) I(0) I(X0o i) 2 (OL ito IXoR) ( 0L ito 0OR) p(Eo) o(o ifxn 9")( itR JoR) \UI Xn Vn (V)'Xn (8.8) (OL itXO iVoR) (LOLI itxo OR) In this case, the Murray and Davidson approximation gives OL it i o) ( l )OR) C) I(n)/I(O) Z (OL I i (D ) (8.9) (0 iit o) (oL 'I O ) (8.9) <0L it {V\o> In cases where HN is close to being symmetric, there is some justification, consid ering the number of approximations already made, in reducing this back to equation 8.6. This approximation can, of course, only be applied to those cases where one of the states of interest is a principal ionization state. 8.3 Transition Moment Expression for the EOMIP problem For the EOMIP case, the expressions for the transition matrix elements are: 0rk = (0 (1 + A)eTOeTXty 0) 13 0 0  ^  0 ao ^  0 Figure 8.1: Diagrammatic expression of r0_, terms = (i O x) Dix + (a0 I) Dax i a + (al x) a (i i Ob) Fibra + E (ilO j)ra + 1 (bl Ic) Fbcax Sib ij be Dix c ce + C tm)em 1 c.tef Aef zt E c2mi M WiAm 1: Ci in mn me menf S )C r\a +Z Ke ae Dax 1Z m mn mn Damm Cmn mn i mne r C~ b irt 1 E Aae rCl e S"ibax AZ m mi 2i M mn m m mne rijax = A ae Aae Ss jmCim me r ax e = ebc ic Fbcax = 2 ^ mnCmn" T2 mn The other matrix element is: k0o = (01 YK~xTOeT I0) where (8.10) 0 V S0(T0 a0 Oo .... k 0 Figure 8.2: Diagrammatic expression of P terms = Z(xOli)Dxi+ Z(x Oa)Dxa (8.11) + a) (o Ob) Fiba + (i i Xa + 0(b c) bcx a bc / where Dxi = Fn2  V0 :*, 1 ta m mne m ( +mvm m = 2im  j ri me  0 mn mn abcXa e I n m mn In cases where the outgoing electron function is orthogonal to the virtual orbitals, all F terms, which involve (x p) or (p x) will disappear, leaving the much simpler D expressions in terms of (p 0 lx) and (x 0 p). Using the approximations of Martin and Shirley, these expressions reduce to 'k (i O z\L Dix) (8.12) c i r o (Xh i) Di i (XI D (8.13) where i is the occupied orbital which dominates the ionization event. where i is the occupied orbital which dominates the ionization event. That means that the ratio of intensities of satellite state A" and its principal state C is approximated as. I(AP) DAT1Y ____ ~ xt Zx I(K) D D' ^Y (c + E (cine cm t t) A E cA me mirenf (8 me menf Because YK is, to first order, identical to CK, one can make the further approxi mation I () me menf ci c + E c 2 cim zJ Am ~ c tff A2m n \ me 2 men f/ The simplest possible approximation, then, to EOMIP relative intensities is I(nA) Ci2 S(8.16) I(K) Ci2 ' and the next simplest approximation is: () ((C + C (c$ ctCmti) / me One more approach is simply to use the approximations to the transition matrix elements directly, as electron propagator practitioners often use pole strengths to approximate intensities. This procedure is most defensible when the pole strengths are compared to lowincident energy photoelectron spectra. I(A) ( c. + E ( cA c A). (8.18) me CHAPTER 9 IMPLEMENTATION The EOMIP equations have been implemented in the ACESII framework, us ing the direct product decomposition approach of Stanton, Gauss, and Bartlett. Rather than explicitly forming the matrix defined by equations 3.20 and 3.21, the nonsymmetric generalization of the Davidson[24] method developed by Hirao and Nakatsuji[42] for eigenvector extraction has been used. In this approach, matrix vector products are formed and stored on disk, and only a smaller, projected "mini matrix" is kept in memory. It is computationally convenient to form a small matrix in memory and diagonalize it to recover all the eigenvectors and eigenvalues. But in many cases, only a few of these solutions are of interest, and the diagonalization becomes increasingly expensive. Therefore, it is advantageous to use methods in which the matrix to be diagonalized is never formed explicitly, or at least never held in memory, and the matrix is never diagonalized. All of these large matrix methods depend on the fact that the Raleigh quotient of a nonsymmetric matrix A ytAx p(y, x) = x (9.1) ytx is stationary with respect to variation of the elements of x and y if and only if these vectors are corresponding left and right eigenvectors of A. Davidson's method, like many other large matrix methods, uses a variation on the Lanczos method. In both these methods, the nonsymmetric eigenvalue problem AXi = XiAi (9.2) YiA = AYit (9.3) is solved by projecting the nonsymmetric matrix A into a smaller space spanned by a set of trial vectors B. Since the matrix is nonsymmetric, it is formally correct to use two different biorthogonal trial spaces, B and L. The projected "minimatrix" M = LtAB (9.4) LtB = I (9.5) is diagonalized to yield approximate eigenvalues and eigenvectors which are expansion coefficients to approximate eigenvectors of A DtMC = d (9.6) X, = CCikBk (9.7) k Yi = DikLk. (9.8) k Residual vectors X and y are a measure of the error in the approximate eigenvector. These error vectors are defined to satisfy (A di)(X, Xi) 0 (9.9) (A di)(Yi Yi) 0. (9.10) It is the method of approximating the residual vector which characterizes these meth ods. The Davidson method approximates these as Xii = (di Aj)' [(AX)ji diX i] (9.11) y, = (d, Aj)1 [(Aty)j dY] (9.12) These new correction vectors are then biorthonormalized and added to the B and L spaces, after which the entire procedure is repeated. If only the righthand eigenvectors, X, are required, there is another approach which can be used to save computational time. The righthand projection space B can be used in place of the L space, and the method becomes almost identical to the symmetric Davidson method except that the algorithm for diagonalization of the minimatrix M must be suitable for nonsymmetric matrices. By the time convergence is reached, B generally spans a space sufficient to represent, to a good approximation, the lefthand eigenvectors of A. In test calculations, these two approaches have been shown to have very similar convergence properties[42]. At each iteration, it is necessary to choose some subset of roots of the minimatrix to be improved. There are several approaches. If one is interested in the N lowest roots of the matrix A, one can simply always choose the lowest roots of the minimatrix. If, however, it is desired to converge directly to higher roots, some other approach must be used. Because the principal IP eigenvectors are, in almost all cases, dominated by the annihilation of an electron from a single occupied molecular orbital, it is useful to "preload" the B space with guess vectors having Ck = 1 for each root K desired corresponding to the ionization of the electron in orbital k. Then, at each iteration, those minimatrix roots which are most similar to the original guess vector are chosen for improvement until convergence is reached. Alternatively, those roots which are most similar to roots saved from the previous iteration can be chosen, but the au thor has found that the use of this approach can lead to choosing the wrong root, as the desired root may be poorly represented in the first few iterations, and the algo rithm may follow several different roots successively, finally settling on an unexpected solution. Shakeup roots which are characterizable as satellites of a particular principal ion ization are still dominated by the annihilation of a single electron, but also contain large contributions from one or several excitations of the ionized species. These so lutions are more difficult to characterize for two reasons: 1) it is not always known ahead of time which excitations are most important to the desired shakeups, and 2) even in cases where the important excitations are known, the unoccupied orbitals available from a typical HartreeFock calculation in a larger than minimal basis set are not useful in a "molecular orbital" analysis of the system being studied they have very little correspondence to the orbitals which would be produced in an excitedstate or N+1 (or more) electron calculation which occupied portions of the desired orbital space. In this case, it is useful, once the principal root is found, to specify an energy range which removes the principal root from consideration, and then converge all roots found within that range, or impose the additional constraint that chosen roots should have a significant overlap with the principal root of which they are satellites. In the author's experience that first converging the principal roots in question results in a B space containing significant contributions to the space of shakeup roots. Sev eral of these satellites near in energy to the principal may be selected within the first few shakeup iterations, especially in cases where these shakeups are large compared to the "principal", and the two wavefunctions are similar. In other cases, the author uses other procedures to load the B space with relevant vectors, such as relaxing energy range restrictions for a few iterations, even choosing random vectors. CHAPTER 10 APPLICATION OF THE EOMIP METHOD TO SMALL MOLECULES 10.1 A Survey of FSCCSD and EOMCCSD Calculations Table 10.1 is a survey of FSCCSD for the ionization of a range of small molecules. Bernholdt [40] reported results differing by almost zero to approximately .6 eV from experiment for valence ionizations of a range of molecules, mostly in DZP basis sets, mostly at experimental geometries. There are a few general trends. Energies in a larger basis at the same geometry are larger. DZP ionization energies are generally smaller than experimental energies. For wellbehaved closedshell molecules, the first ionization energy, even at the DZP level is generally within .2 eV or so of experiment. A notable exception to this rule is the stetrazine example. In the larger basis sets, there is a trend for higherenergy EOMIP ionizations to be increasingly high relative to their corresponding experimental energies. Increasing the basis set size, as in the cases of the Stanton, Bartlett, and Rittby 02 study and the Pal et al. study of formaldehyde, improves agreement with experiment. Table 10.1: A survey of FSCCSD results Final Energy (eV) Molecule Basis State Theory Expt. CH2NH DZP' A 10.34 10.52, 10.70, 10.56 2A" 12.29 12.43, 12.48, 12.44 2A' 15.08 15.13, 15.11, 15.00 2A' 17.19 17.04, 17.07, 17.00 stetrazine DZPh 2B3g 9.20 9.7 2BI. 11.87 11.9 2B2g 12.12 12.1 2A, 12.69 12.8 2B2. 13.01 13.3 2Big 13.52 13.5 H2CO [5s3pld/2slp]a 2B2 10.51 10.88 2B1 14.36 14.38 2A1 15.75 16.00 H2CO [5s3p/2s]d/ 2B2 10.07d, 10.43h 10.88 [5s3p1d/2s1p]e 2B1 14.14d, 14.13h 14.38 2A1 15.20d, 15.67h 16.00 CH2 DZPc 2A, 10.20 10.26 (DZPFCI) 2B2 14.91 14.85 (DZPFCI) 2A 22.55 22.14 (DZPFCI) CH2PH DZP 2A" 10.08 10.30 2A' 10.28 10.70 2A' 13.17 13.20 3g 02 DZPa/PVQZ++g 2IIg 11.76a, 12.399 12.35 4II 16.40a, 16.76g 16.85 4Eg 17.75", 18.269 18.33 4E 24.58a, 24.829 24.66 N2 [5s4pldja 2Eg 15.44 15.6b 2Iu 17.14 16.88b S, 18.64 18.6b ketene DZP" B 9.40 9.64(a) 2B2 14.26 13.84(a) H20 DZPc 2Bi 11.97 12.61 2A 14.28 14.73 2B2 18.75 18.55 diazomethane DZP 2B1 8.60 9.00(a) 2B2 13.79 14.13(v) HF [5s3p/3s]/ 2I 15.08 16.1 2 19.43 19.9 SE 39.75 39.7 a Bernholdt [40]; b Experimental results [34] c Vaval, Ghose, Pal, and Mukherjee [43] d Pal, Rittby, Bartlett, Sinha, Mukherjee [44], [5s3p/2s] basis; e ibid, [5s3pld/2slp] basis (a) and (v) indicate adiabatic/vertical ionizations of ketene and diazomethane f Sinha et al. [45] 9 Stanton, Bartlett, and Rittby [46] h Rittby and Bartlett [13] Watts, Rittby, and Bartlett [47] These results are encouraging; it is easy to surmise that the EOMIP approach can be applied to a variety of molecules with very good results. In addition, the EOMIP method can be used to examine shakeup states, and is not plagued (as is FSCC) by intruderstate problems in the innervalence region. The EOMIP approach is conceptually well suited to the study of ionic spectra, including first, principal, and shakeup ionization energies and intensities. For N2, for example, there are six molecular orbitals (five a and one 7r) to be explored. Using a wellcorrelated singlestate method such as CCSD or CCSD(T) would require one calculation for the neutral N2, and a separate calculation for each of the N' states. All the N' calculations would have to be openshell, doubling the effective size of the orbital space involved in the SCF and correlated calculations. In addition, it is not easy to obtain a suitable set of molecular orbitals for every openshell state. While it is a simple matter to remove an electron from the highestenergy orbital of any symmetry, it requires a bit of trickery and a few more SCF iterations to converge on a corehole state. And it is generally impossible to converge an SCF state with an electron removed from an innervalence orbital. The QRHF reference, a nonSCF state can be used in this case. There is no reliable way to describe a state corresponding to a shakeup peak, and no way to specify which of the two 2ag nitrogen peaks (table 10.4 and figure 10.1) will be found in this approach. When the EOMIP approach is used, only one CCSD calculation, on the closedshell state, is required. The EOMIP calculations which follow are considerably less expensive in terms of computer time. And shakeup peaks can be specified in terms of energy ranges and the molecular orbitals with with they are associated. The equationofmotion excitation energy (EOMEE) approach can also be used as a route to finding ionization states in the excitation spectrum of N' species. There are two disadvantages to this approach. First, not all of the desired ionization states can be described primarily as a single excitation from one reference. Multiple references, and thus multiple CCSD calculations must be employed. Secondly, although the linear form of the EOMEE equations makes iteration less computationally demanding than CCSD iterations, EOMEE is still more demanding than EOMIP there are more amplitudes and more and bigger terms to calculate in EOMEE. What follows in this chapter is an EOMIP study of the ionization spectra of several wellunderstood molecules, using large basis sets, with attention to prominent shakeup features. All calculations were performed with the ACES II program system1 and the author's ACES IIbased CC EOMIP program2. The EOMIP intensities which appear in this and the following chapter were calcu lated using the approximation of equation 8.18, except where principal peak intensities have been scaled to experimental intensities, as noted. For the figures, Lorentzian functions were applied to the theoretical transition energy/intensity information in order to simulated the degree of experimentallyobserved broadening. In some cases, the experimental curves were also created from available energy/intensity data by 'ACES II is a quantum chemical program package especially designed for CC and MBPT en ergy and gradient calculations. Elements of this package are: the SCF, integral transformation, correlation energy, and gradient codes written by J. F. Stanton, J. Gauss, J. D. Watts, W. J. Laud erdale, and R. J. Bartlett; the VMOL integral and VPROPS property integral programs written by P. R. Taylor and J. Almlif; a modified version of the integral derivative program ABACUS written by T. Helgaker, H. J. Aa. Jensen, P. Jorgensen, J. Olsen, and P. R. Taylor[48] 2CCEOMIP is a quantum chemical program designed to calculate vertical electron detachment energies and intensities, written by Renee Peloquin Mattie. It requires ACESII CCSD T and H lists [49] application of Lorentzian lineshapes. In the carbon cluster plots, the experimental curves were created by digitization of experimental plots provided by Yang et al[35]. 10.2 Unsaturated Molecules 10.2.1 Nitrogen Nitrogen presented one of the earliest challenges to theoreticians. The ionization energies obtained from Koopmans' theorem are not predictive the observed spectrum, and are in fact incorrectly ordered (table 10.3). Selfconsistent field calculations do not order these ionizations properly. Correlated calculations are necessary, and highly correlated calculations are needed to describe the breakdown of the molecular orbital picture in the innervalence region correctly. Two basis sets were used in the nitrogen calculations. The smaller basis was based on the triplezeta plus polarization basis set of Dunning[50], with the innermost p or bital uncontracted, as was used by Rittby and Bartlett[34]. The larger is the Dunning pVTZ+ basis set [51], a polarizedvalence triplezeta basis set which includes diffuse s,p,d and f functions Here, however, the CCSD optimized internuclear separation, R = 2.0778 was used rather than the experimentallyobtained geometry R = 2.0693 used by Rittby and Bartlett. In the smaller basis set (table 10.2), results show agreement with experimental results to within .2 eV for the outervalence region. The 2a, inner valence hole state, however, shows an error of 1.3 eV and the lo, core hole state is in error by .7 eV. Table 10.2: N2 in 5s4pld basis at R = 1.09898 A N2 N' Reference EOMIP (Itheo) Expa X E+ RHF 109.358913 au 9 X2Eg 3ag 15.376 eV 15.5 eV 2II u lr 17.022 eV 16.8 eV 2E 2a, 18.607 eV 18.6 eV 25.0 eV (shakeup) 28.2 eV (shakeup) 31.5 eV (shakeup) 2E+ 2ag, 1~tt17r,2a, 32.153 eV (.5b) E+ 2ag 38.616 eV (1.06) 37.3 eV 2Eg 2a, 5ug3a,3ag 42.543 eV (.2b) 2Eg lag 410.623 eV 409.9 eV a Banna and Shirley [52] b Simple approximation relative to 2ag principal Moving the the larger pVTZ+ basis seems to afford no improvement of the valence and core peaks, in fact shifting each peak toward higher energies by an amount which increases with the transition energy. As shown in table 10.3, the valence a, and a, holes are still within .2 eV of experiment, but the energy of the 17r~ hole seems to have gotten worse it is not quite within .5 eV of experiment. However, the inner valence region is now much more similar to the experimental spectrum. The inner valence 2ag hole ionization is characterized experimentally by a fairly broad peak, relative to the outer valence or core ionization. In this region, the simple MO picture breaks down, and it is difficult to speak of a single line corresponding to the ejection of the 2a, electron. The calculations show two lines of similarly large intensity, one near 33 and one near 39 eV. Vibrational and vibronic effects can be expected to be quite strong in this region [31] as well. Table 10.3: N2, pVTZ+ basis Ionization energy eV (Intensity) Koop. EOMIP Other" Exp Orbital approx CLOSED CORE Assignment 17.277 16.793 16.793 21.147 15.683 17.332 17.332 18.888 29.887 32.679 33.319 40.189 38.695 39.307 42.400 42.519 43.030 43.461 44.566 45.452 46.919 426.723 410.587 426.623 410.768 (.931) (.958) (.961) (.895) (.02) (.320) (0.004) (.518) (.002) (.034) (.001) (.045) (.067) (.067) (.001) (.005) (.806) (.807) 15.499 (.880) 17.159 (.918) 17.159 (.964) 18.674 (.858) 28.917 31.983 (.026) (.284) 15.666" 15.5 17.227 16.8 17.227 16.8 18.908 18.6 25.0 29.0734 38.106 (.470) 409.609 409.799 (.753) (.755) 3"9 3erg 1 7r2 2o7 28.2 2a,, 27r3aglru, 31.5 2og, 27rT2a~,17 17r, 37r]lau3ag 37.3 2ag 2ao, 4art3a 3a 2og, 47! 37u37ru 2a9, 5 3og3og 29g, 2TT37g27ru 27g, 4~3 u37r, 27g + 37, 27T17u,3g 2au, 6 l17fl7iru 2au, 53o3g337g 409.9 la, la a EOMEE, except where otherwise noted b ACCSD In the plots in figure 10.1, Lorentzian lineshapes have been applied to the EOMIP peaks to simulate the line broadening found in lowresolution experimental spectra. Each EOMIP plot is superimposed on a plot of an experimentallyderived photoe mission spectrum based on MgKa [52] and YM( [53] experiments. Peak positions were taken from the YM( experiment (at 1253.6 eV), which provided much better resolution in the region between 20 and 45 eV, while relative peak intensities were taken from the MgKa experiment (at 132.3 eV) which is more suitable for compar ison to EOMIP intensities, as they should most resemble threshold photodetachment intensities. The simulated spectra in figure 10.1 indicate a good match between the experi mental peaks at 28.2 and 31.5 eV and the 2cr, peaks at 29.9 and 32.7 eV with pVTZ+ closed shell orbitals and 28.9 and 31.9 eV with pVTZ+ core hole orbitals. The core hole orbitals give somewhat better agreement over the entire spectrum, especially at higher energies. In the innervalence region and for the core peaks, improvement becomes notable, as seen in table 10.5 No peaks were found which can reasonably be assigned to the 25.0 eV experimental peak. Schirmer, Cederbaum, Domke, and von Niessen, in their 1977 paper, saw a peak right at 25 eV, but they were using a much smaller basis set (11s 7p/5s 3p), which is lacking in polarization functions com pared to the pVTZ+ (11s 6p 3d 2f/5s 4p 3d 2f). The author's experience and that of others indicates that this means some larger shakeup peaks in the smaller basis set, and Schirmer et al. show more shakeup peaks than are perhaps indicated by the experimental spectra, making a doubletwithashoulder of the 37 eV 2Ug ionization, and giving them a significant 2og peak at about 40 eV, similar to the appearance of Table 10.4: Summary of N2 EOMIP ionization calculations Energies eV pVTZ+ C.S. C.H. Exp 15.683 15.499 15.5 17.332 17.159 16.8 18.888 18.674 18.6 25.0 29.887 28.917 28.2 32.153 32.679 33.319 38.616 38.695 39.307 42.400 42.519 43.030 43.461 44.566 45.452 410.587 410.623 410.768 31.983 31.5 38.106 37.3 409.609 409.799 409.9 Table 10.5: N2 EOMIP error, 5s4pld and pVTZ+ bases. Energies and errors are in eV. Expt. 5s4pld pVTZ+ closed core 15.5 0.1 0.2 0.0 16.8 0.2 0.5 0.4 18.6 0.0 0.3 0.1 28.2 1.0 0.7 31.5 0.6 1.2 0.4 37.3 1.3 1.4 0.8 409.9 0.7 0.8 0.2 5s4P1D 15.376 17.002 18.607 N2 lonizations, EOMIP 0 12 1'4 116 1iB 2! 212 24 26 2 30 312 314 Ionization Energy (eV) N2 Ionizations. EOMIP (6 3b 4b d2 4 I I' II ii ii II ii II 30o 1C 10; 20u 20a A Expt  EOMIP  I II Il 'I II ; I11 I ) 12 114 1b 1b 2O 22 4 26 2A 8 3b 312 34 36 38 kb 02 4 Figure 10.1: N2 ionizations in three MO sets 3oa 1u, Io; Expt  EOMIP I I I  I ' 1 2 N2 Ionizations, EO MIP N2 lonizations, EOMIP Expt 30g 1"u a.* 20; OMIP  Ii i i I p I' II III II\ J & ~'~ '   ~ J. I , r\ r Ij\ significant EOMIP peak at 38.6 eV when the 5s4pld basis is used. This region shows much less intensity in the pVTZ+ closedshell or corehole orbital calculations, in agreement with experiment. In the corehole orbital set, the largest of the 2a, peaks has moved .8 eV closer to the experimental peak, while its companion peak has moved .3 eV closer to experiment and the ratio of the two peaks is more similar to that in the experimental spectrum. In the core region, the la, and lao peaks have both moved slightly closer to the experimental value. In general, the energies of all the major peaks are reduced in the corehole basis, moving closer to experiment. Interestingly, corehole orbitals do not improve the performance of EOMIP by decreasing the size of correlation corrections to the wavefunction. The sizes of the EOMIP doubles contri butions in the two cases are similar, although there is some difference in the ordering of these contributions because the largest T1 amplitudes are approximately twice as large (.075 versus .015) in the corehole as in the closedshell basis. 10.2.2 Ethylene Ethylene has been studied in two different basis sets Nakatsuji's Basis I[54] and the Dunning pVTZ+ basis[51]. In the smaller basis, the Herzberg[55] ground state geometry was used and relative intensities were calculated using the simple approximation of equation 8.16. For the larger basis set study, the geometry was optimized at the CCSD level in the pVTZ basis set, and the improved approximation to relative intensities from equation 8.17 was used. The inner and outervalence region of the ethylene ionization spectrum has been studied before, with an eye toward the assignment of the strong satellite peak at 27.4 eV[6, 28, 56, 54]. The valence ionizations, as well as the strongest of the related shake ups, are given in table 10.6 along with the experimental numbers and the theoretical results of Nakatsuji in the same basis set. The outervalence results show excellent agreement with experiment, while the innervalence principle peak at 24 eV, like the SACCI, appears at a slightly higher energy. The shakeup at 27.4 eV appears in the experimental spectrum as a strong, somewhat broad satellite. Cederbaum[56] implies that this peak should be assigned to a 2ag shakeup, and in this Davidson[28] concurs. Nakatsuji attributes this peak to a sum of the nearby 1b2z and 2b1i shake ups. The difference in the ordering of the 1b2a and 2a, peaks just below the the experimental peak is puzzling for two methods which are apparently so similar, but there are several other theoretical shakeup peaks with very small intensity in the same region. It would seem that the nearby 2ag shakeup is another likely candidate. It is interesting to note that both Cederbaum and Davidson report significantly greater relative intensities for the shakeup peaks than have been obtained from either the EOMIP or SACCI calculations. Davidson has noted that the satellite intensities decreased as the virtual space size was increased, so that the addition of Rydberg orbitals to relatively small basis sets such as the Nakatsuji basis I could be a significant effect. To clear up some of the questions raised, the study was repeated using the larger, Dunning pVTZ+ basis[51], which includes more diffuse s,p,d and f func tions. The pVTZ, CCSD optimized geometry is Rcc = 1.327A, RCH = 1.077A, Table 10.6: Ionization spectrum of ethylene, 1030 eV Ionization energy eV (Intensity) EOMIP SACCI" Experimentb Orbital Assignment pVTZ+ BASIS I 10.72 (.95) 10.39 (.95) 13.14 (.93) 12.85 (.93) 14.97(.92) 14.59 (.92) 16.37(.88) 16.03 (.88) 19.65 (.85) 19.36 (.85) 24.52 (.007) 23.55 (.008) 24.51(.74) 24.25(.74) 10.25 (.95)c 12.78 (.93)c 14.50 (.92)c 15.93 (.88)c 19.32 (.85)c 23.45(.01) 24.37(.74)c 25.86 (.04) 10.51 (.95) 12.85 (.93) 14.66 (.92) 15.87 (.88) 19.1 (.85) 23.7 (.74) 26.13 (.015) 26.16 (.052) 26.73(.052) 26.81 (.006) 28.18 (.001) 28.6 (.017) 29.26 (.044) 30.03(.012) 30.77 (.075) 26.97(.01) 27.4 (.29) 29.76 (.010) 30.06 (.094) 29.28 (.02) 30.15(.09) 1b3, lb3g 3ag lb2u 2bi, 2bij, 2ag 1b2u, 2ag, lb2u, lb2u, 2ag, 2bij, 2ag, 2biu, lb39g, 1b3g, 2ag, 2ag, 2a., 31.2 (.1)d 2b2g1b3ula, 2bt lb3glb3u 7at 1b3ulb3 3b2glb3glb3u 2bag lb3g lb3u 5ag lb3ulb3u 3btf lb3u lb3u 7at 1b3ulb3u 2b2g 3ag lb3u 2bt2 lb2u lb3u 3bt3 lba3u lb3 lbtg lb3u 2blu 6at lb3ulb3u, 2b2g 2b lb3u 2b 2bu lb 2g, 2b1, 1b3, a Nakatsuji [54] b Banna and Shirley [52], Gelius [57] c Theoretical intensities of principal peaks scaled d Energy and intensity from Weigold [58] to experimental data Ethylene Ioniations, EOMIP, Nakatsuji basis EO 0.9 0.8 0.7 III II 0. II I II 0.4I III i I 0.3 I I I I I 1 0. I I I I Figure 10.2: EOMIP and SACCI ionization spectra of ethylene Ethylene lonizaions. EOMP, pVTZ+ beasi LHCH = 121.4590. The virtual space is significantly larger, but the diffuse functions less so than those in the Rydbergsupplemented basis used by Nakatsuji. It was an ticipated that this basis set would represent an improvement over the smaller basis and, since ethylene is well represented at the CCSD level of theory, the end result would be improved ionization energies. Instead, what is observed is a spectrum of ionized functions which are generally similar to those obtained in the smaller basis (allowing for the presence of diffuse orbitals in the lowlying virtual space), but with ionization energies which grow progressively larger than the experimental energies and those calculated in the smaller basis. A small study of the effects of basis set choice and orbital relaxation on the 3a, ionization (table 10.7) reveals that removing the diffuse functions or replacing them with Nakatsuji's moleculecentered Rydberg set has little discernible effect on the ionization energies calculated by the CCSD, CCSD(T), or EOMIP methods. It is apparent, however, that orbital relaxation is an important effect. The ACC (QRHF calculations for the cation) and EOMIP results in the C2H4 orbitals overestimate the energy by about .19 to .34 eV, while in the C2Hf ROHF orbitals, ACC (QRHF calculations for the neutral) methods actually underestimate the energy by over 1 eV, and the EOMIP overestimates it by a similar amount. In the more natural calculation  both neutral and cation correlated calculations are based on the appropriate SCF calculation both of the CC methods perform well, within .2 eV of experiment in either basis set. The problem for the CC calculations is apparently in large part one of orbital relaxation. Neither set of orbitals provides a reasonable starting point for both the neutral and the cation. The coupledcluster triples correction can hardly be said to be much help in improving either of the semiQRHF calculations. In the pVTZ closedshell orbitals, triples lower the energy of the neutral by .42 and the cation by .38 eV. In the pVTZ cation orbitals, triples lower the neutral by .47 and the cation by .44 eV. Table 10.7: Effect of orbital choice on ethylene 3ag ionization basis MOs Energy eV ACCSD ACCSD(T) EOMIP pVTZ C2H4 14.846 14.877 14.935 C2H+ 13.458 13.486 15.902 MIXED 14.879 14.855 pVTZ C2H4 14.847 14.880 14.936 +Rydberg C2H+ 13.669 13.752 15.886 MIXED 14.881 14.858 Basis I C2H4 14.610 14.627 14.686 C2H+ 13.491 13.533 15.660 MIXED 14.640 14.604 "ROHF performed with 3ag orbital half filled. bNeutral in Neutral orbitals, cation in ROHF orbitals The difference between the ACC and EOMIP calculations in either set of orbitals are largely attributable to spectator triples excitations which are CC doubles but EOMIP triple excitations. In the CCSD calculation on the neutral in these orbitals, there is a large T2 amplitude (.14272) describing the excitation of both Ibl~ (7rcc) HOMO electrons to a lb2, (7ric) orbital, followed by another double excitation r to 7r*(.05468). The EOMIP wavefunction is dominated by the annihilation of the 3ag electron, with other contributions from the ir * 7r* double excitation, followed by a single excitation from the lba3 HOMO to the 2b3u virtual orbital and one from the 1b2u to 2b2,, giving a wavefunction dominated by the action of .143 2b 2g21bglb3, 1b3ag + .0953 2b, lb3a3a + .0839 2b2lb2u3ag. For the cation in these orbitals, the largest amplitude is again the T2 providing the double excitation from the Ib3 lrcc bonding HOMO to the 2b3u, rc anti bonding orbital, followed by a double excitation from the bonding to the anti bonding orbital, combined with an excitation from the halffull 3ag orbital to the empty 4ag (.127 2btg2bglb3u,1b .0701 3a2blu,2bg1b3u .0547 24g1,3 lb31 )3a, It is the second term here that is missing from the EOMIP wavefunction expression above, where it contributes as a T2 C1 term, and has an amplitude of less than .001 relative to the principal ionization. EOMIPSD has no way to recover this determi nant. Strangely enough, the same phenomenon is observed in the C2H4 orbitals, but the effect on the energy is not nearly as large. The ACC results are much closer to the "mixed" results (around .2 eV) and the EOMIP is closer as well (within .1 eV of the ACCSD, and within .3 eV of either "mixed' result). In Nakatsuji's Basis I, the very same kind of behavior is again observed. In the closedshell orbitals, the largest contributors to EOMIP wavefunction are 3,g+.098 3b3l1b3u3a,+.0531 4bglb3u3b2glb3,3ag.0525 3bglb33a,.0489 4bu1b2u3ag, while the largest contributors to C2H4 CCSD wavefunction are 3ag+.143 32g91 33bg1lb3u3g .103 3bulb3a3~g .078 2btlb33ul g .076 lbtg3a 2blu. The first three determinants are the same in either case, but with very different relative weights. Still, the EOMIP energy is very slightly closer to the ACCSD energy than is the case in the pVTZ basis, with or without Rydberg orbitals. Going from the small BASIS I to the larger pVTZ basis lowers the energy of the closedshell ethylene more than it does the cation, which means that all the calculated ionization energies are lower in the Nakatsuji basis. The ACC energies in the closedshell orbitals are actually below the experimental value of 14.66 eV, and the EOMIP value is a very good match to experiment. Because Nakatsuji's basis I lowers C2H4 ionization energies relative to more com plete basis sets it manages, without improving the agreement between EOMIP and ACC methods, to provide a very good match between SACCI or CC EOMIP calcu lations and the experiment. 10.2.3 Butadiene As the experimental spectrum in figure 10.3 indicates, there is noticeably more complexity in the butadiene than in the ethylene valence spectrum. The broadening of the peaks in the region above 20 eV indicates the possibility of shakeup peaks as well as the onset of shakeoff effects double ionizations. Because the PS (photoemission spectroscopy) experiment was performed in the xray region, at 1487 eV, the experi mental intensities are not directly comparable to those obtained using the formula in equation 8.18. However, there appears to be good correlation between the positions of experimental and theoretical peaks. The tendency of higher energy peaks to be increasingly higher in energy than experiment is consistent with the trends observed at this level of theory for large basis sets with several polarization functions. 