Theory and implementation of coupled-cluster methods for direct calculation of ionization energies


Material Information

Theory and implementation of coupled-cluster methods for direct calculation of ionization energies
Physical Description:
xi, 166 leaves : ill. ; 29 cm.
Mattie, Renée Peloquin, 1964-
Publication Date:


Subjects / Keywords:
Cluster theory (Nuclear physics)   ( lcsh )
Quantum chemistry   ( lcsh )
Chemistry thesis, Ph. D
Dissertations, Academic -- Chemistry -- UF
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1995.
Includes bibliographical references (leaves 161-165).
Statement of Responsibility:
by Renée Peloquin Mattie.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 002059459
oclc - 33857533
notis - AKP7522
System ID:

Full Text








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.


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 late-night rides back from campus; Ajith Perera and all the students

in the group for taking care of details long-distance; 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.



ACKNOWLEDGEMENTS .................... ........ iii

LIST OF TABLES ................................ viii

LIST OF FIGURES ..... ..... ...... ... .. ..... ... ix

ABSTRACT .......... ............ ............. .. .. x


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 Coupled-Cluster Theory ......................... 13
2.2.1 The Coupled-Cluster Operators .......... .... 13
2.2.2 The H Operator ......................... 14
2.2.3 The Coupled-Cluster Equations .. ........... .. 16
2.3 Configuration Interaction ......................... 17
2.3.1 The CI Operators ......................... 17
2.3.2 The CI Equations ....... ... .............. 18


3.1 The EOM Equations ........................... 19
3.2 Size-Extensivity 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.1 Fock-Space 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 Non-Dyson Equation Approaches ................ .... 41


6.1 State Methods ........ ........................ 44
6.2 Direct Approaches ................... ........ 46
6.2.1 Two-Hole One-Particle Configuration Interaction ....... 46
6.2.2 Spin-Adapted Cluster Methods ................. 47
6.3 Eclectic Approaches ............................ 48


7.1 Perturbation Applied to the CC Hamiltonian ............. 49
7.2 coupled-cluster Response Theory ................ .... 53
7.2.1 Time Dependent Perturbation on a Coupled-Cluster Reference 53
7.2.2 The Linear Response Function . ... 55
7.2.3 Energies and Transition Intensities ............ ... 57


8.1 The Monopole Approximation ................. .... 61
8.2 Transition Moment Expression for Non-Symmetric Eigenstates .. 63
8.3 Transition Moment Expression for the EOMIP problem ........ 63

9 IMPLEMENTATION ................... ........... 67


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 Single-Determinant and Multideterminant References ... 97
10.4.2 The Spectator Triples Problem . ... 102


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 Odd-Length 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 Even-Length 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


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

THEORY .... ................................ .. 156

B.0 The CC EOMIP-SDT Equations ................ 156

B.O A review of the Rayleigh-Schr6dinger Perturbation Theory 158
B.O EOMIP Perturbation ........................... 159

REFERENCES ................... ................ 161

BIOGRAPHICAL SKETCH ............................ 166


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, 10-30 eV .
Effect of orbital choice on ethylene 3ag ionization
Ionization of 1,3 trans-butadiene, 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 closed-shell molecules .















2.1 The HN operator ............................
2.2 HN elements required for CC EOMIP-SD 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. . .



N2 ionizations in three MO sets . ...... 81
EOMIP and SAC-CI ionization spectra of ethylene .. 85
Ionization spectrum of 1,3 trans-butadiene . .... 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, two-determinant 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



Renee Peloquin Mattie

August 1995

Chairman: Rodney J. Bartlett
Major department: Chemistry

The coupled-cluster equation-of-motion (CC-EOM) method, formally equivalent

to the Fock-space coupled-cluster (FSCC) and coupled-cluster linear-response (CCLRT)

methods for electron detachment energies, offers distinct advantages over the single-

reference coupled cluster (SRCC) approach when several electron-detached 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 equation-of-motion ionization potential (EOMIP) method is implemented

using a Davidson-type algorithm for extracting selected roots from a large symmetric

matrix, making it possible to selectively search for the difficult-to-characterize "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 inner-valence regions are explored, with special at-

tention to shake-up transitions. Results compare favorably with the coupled-cluster

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.



1.1 Correlated Methods in Quantum Chemistry

In the early days of quantum chemistry, the only ab initio method available was

the self-consistent 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 free-electron continuum. Photoelectron spectroscopists then observe

the absorption spectrum or the spectrum of the free electron kinetic energies. Koop-

mans' theorem[l] is an SCF-based 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 so-called shake-up peaks are, in the MO picture, the result

of a mixing of the primary hole-state process (ejection of an electron from a previ-

ously occupied orbital) with a two-hole, one-particle 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 single-reference 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 lowest-energy 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 many-body perturbation theory and coupled-cluster (CC)[2] methods have

been developed specifically to avoid the deficiencies of the CI method. These methods

are "size-extensive"; 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. Second-order many-body

perturbation theory (MBPT-2) calculations, which generally provide 80% or more

of the correlation energy, require no more computational resources than do SCF

calculations. Coupled-cluster 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 mid-60s, "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 symmetry-adapted-cluster (SAC) and symmetry-

adapted-cluster CI (SAC-CI) approaches to calculating the ionization energies and

intensities of molecules. The Fock-space coupled-cluster 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

linear-response theory (LRT)[16], is based upon a coupled-cluster reference. For the

calculation of principal ionization energies, this theory can be shown to be formally

equivalent to the Fock-space coupled-cluster (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 linear-response-theory 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.



2.1 The Schr6dinger Equation

All electronic structure methods seek solutions or approximate solutions to the

Born-Oppenheimer approximation to the Schr6dinger equation,

HT = ET (2.1)

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 one-electron and a two-electron portion,

h=H(one) = V Z-rQ
2i i
2 = H(two) = r -1.

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 two-electron term, gives a second-order 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, ,((),

one-electron 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 self-consistent field method, such as the

Hartree-Fock theory. In the Hartree-Fock approximation, the approximate Hamilto-

nian is referred to as the Fock operator.

F = h+VF (2.4)

H = F+(V2-VF)

The basic component of an N-electron wavefunction is a single Slater determinant


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 one-electron 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 many-electron

wavefunction within this truncated space can be represented as a linear combination

of Slater determinants.

The finite space spanned by the Slater determinants is


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 second-quantized form as

H = h,,ptq + E (pq? rs)ptrqts (2.5)
pq 4 pqrs

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, = -)). -)j- p()...q((N)|

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 second-quantized

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

single-particle operators within a string of operators results in a change in sign of the

normal-ordered product,

{pqsrt} = -{pqrts} = {prtqs} = -{rtpqs}

{jti} = {it} = -ift

{ati} = ati.

The last of these, ati, composed only of hole- and particle-creation operators, is

the only type of normal-ordered 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


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 normal-ordered Hamiltonian, HN is that part of the

Hamiltonian whose expectation with respect to the reference 10) is non-zero

HN = H-(OIHIO) 0



= fN +WN

= HN(one) + HN(two)

fN = FpqPtq

WN = (pq |rs)ptqtsr.



Writing HN in terms of normal ordered operators makes this partitioning into F and

......b i

a b i c

i jb a

/\-.-X b --

a a+b

k a j i

4 c/ "i k
-- acdb j i k
acdb j~ik

x aV -X

I j a+i

++ ++
k j a i a i c b

j a i b a b j i

Figure 2.1: The HN operator

non-F 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 Coupled-Cluster Theory

2.2.1 The Coupled-Cluster Operators

The coupled-cluster theory defines a wavefunction ICC) which is an eigenfunction

of the full Hamiltonian

H CC) = E CC)

(H (01 H 0)) ICC)



= HN ICC) = Ecc ICC),

where Ecc is the correlation energy obtained from the application of coupled-cluster


The CC wavefunction is defined in terms of excitation operators T as

ICC) = eT 10)


The operator T is defined as

T1 = tati

Ti at'
i3 ab
T3 tiatbtctijk,
"T ijkabc

etc. The inverse of this exponential is simply exp(-T). Therefore

e-TH CC)

= e-TEeT 0)

= E O)

E = (01e-THeT 0).

The correlation energy, Ecc = E ESCF, is

Ecc = (0 e-THeT 0) (0 H 0)

= (O e-THNeT 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

coupled-cluster based methods, so it is worth looking at it in detail. H, expanded in

terms of T in the Baker-Campbell-Hausdorff expansion, is

S = H+HT-TH+ HT2-l T2H-THT+... (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 size-extensivity of the

coupled-cluster method.

KQ+_ .

4K7= \A

v ..\A.--
+ V+ V+
0~~. +

--= H- +/v +

.V +

+ w

/VA/ + F ...-

Figure 2.2: HN elements required for CC EOMIP-SD and FSIP equations

--- -- -- --



---- ------ v

MR =

I# =

2.2.3 The Coupled-Cluster Equations

Similar to H is HN = e-THNeT = (HNeT)c, which appears in the coupled-cluster

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.


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).

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

ij ab


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



which can be solved via a matrix eigenvalue approach.

In the CI approach, there is no Baker-Campbell-Hausdorff expansion to force the

connectedness of H to the wavefunction operator C.





3.1 The EOM Equations

Beginning with the exact coupled-cluster 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

second-quantized 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).



When Q1 is limited to contain excitation-type operators only, it commutes with

T, and we derive the expression


wKe7Tn 0)
e-W Te TQ 10)
e-T W )ceTQ 0)

= (E E) ICcc) = HN, K] CC)

[= HN, Q eT 0) [HNeT, -K] 10)

= [e-THNe, ] 0)

w OK 10) = [HN,QC] 10) (3.3)

where w = Ec E and HN e-THNeT = (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),


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 up-going lines in f available for contraction. This will be

2Ln 1 for the IP problem and 2Ln for excitation energies.

3.2 Size-Extensivity 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 two-body

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 size-extensive. CI-based

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)


7, 10) = I,). (3.8)

The operator 7, consists of a normal-ordered product of creation and annihilation

operators which, acting on the Hartree-Fock reference P0), create new determinants.

For excited wavefunctions, the set {r,} contains the identity operator, as well as

excitation operators ati. atibtj, .... For the electron-detached case, {r,} contains

i, atij *, operator products which will create N-1 electron determinants.

Because the equation-of-motion 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


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 left-hand eigenvectors of the non Hermitian HN,

yield left-hand eigenfunctions of HN,

Dt (01 THNT0) = Dtw (3.13)

(01 Y'-HN = (01 lEtw (3.14)

and therefore the projection of the full-space left-hand eigenfunction into the reduced


(0I Yte-THN c (Ol Yte- TE, (3.15)

just as the right-hand eigenvectors, C, yield the projection of the full-space right-hand

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 left-hand side ground state coupled cluster

wavefunction and the transpose of the right-hand 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 N-particle 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 coupled-cluster 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 "free-electron" index.

Sl S2 S3

D1 D2



D5 D6

Figure 3.1: The EOMIP equation diagrams.

The EOMIP singles and doubles (EOMIP-SD) 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:)

-1 (nmllbe) t)Kb (3.21)
2 nmb e )

Where the notations f and (- [ ... ) indicated, respectively, one- and two-body H

elements. The final term of the doubles equation involves the 3-electron 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 EOMIP-SD wavefunction expressed in equation 3.19 is an eigenfunction of

the HN operator projected into the reduced EOMIP-SD 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 /



4.1 Fock-Space Theory

The Fock-space coupled-cluster method (FSCC) uses an exponential operator to

describe electron attached, electron-detached, 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 coupled-cluster 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 particle-hole 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 particle-hole 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 coupled-cluster 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 single-reference coupled-

cluster problem. The operator exp(S(O',)) is identical to the usual eT operator encoun-

tered in the coupled-cluster equations. Because T-to-T contractions are impossible,

imposing normal order does not change the operator.

The operator 2tot has an inverse,

tot-1 : (0o,)-1 (01)-1 ...}. (4.3)

In the (0, 0) case this is f(0,0) = e-T but in the higher sectors, the form becomes more

complicated because of the normal-ordering 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 N-electron single determinant reference 10) is chosen, which may be

an SCF wavefunction. The reference wavefunction defines a complement, Q, which

contains all N-electron functions orthogonal to the reference. The single-reference

coupled-cluster problem is solved, yielding HN, which is an effective Hamiltoimal

for the (0, 0) sector. This is a canonical (eigenvalue-preserving) transformation of .e

Hamiltonian, chosen so that it is the eigen-operator of the reference function, yield ng

the desired coupled-cluster energy (Ecc) as its eigenvalue.

The coupled-cluster 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)
While Q(mn) contains all the contributions from all the Q's of different sectors, the
projection of equations by p(m,n) and Q(m,") limits the contributions to the energy

and wavefunction expressions to only certain lower-ranked sectors The (0, 0) sector

contributes to all higher-ranked sectors, and is the only one contributing to the (0, 1)

and (1, 0) sectors. The (1, 1) sector has contributions from all lower-ranked 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 one-dimensional problem, it is

clear that something more is needed to deal with the multi-determinant model space.

It is necessary to introduce C(m,n), or C for short. The matrix C causes a linear

transformation of the model-space 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)


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


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)


H( n) p(m,n)C(m,n) = ntot-l~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)

= 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 coupled-cluster equations.

The complicated form of Q-1 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

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)

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 ma|jk -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 FSCC-IP and EOMIP calculations are performed using the same active

space-the 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 Fock-space method defines a model space consisting of several (N 1) elec-

tron determinants and a projector P = p(,1) over this space. An electron-detached

wavefunction is defined as

Q(0,') IK) E D(o,')CK |0) = Q(0,1) C~i 10).
P can be represented as a sum over these IC) states

P = 1K) (KC)- c.
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

C-Heff C 0) = w 0)}


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 shake-up 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


The P-P 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.



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 IP-type operators, i

and atij, the EA-type 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 coupled-cluster methods,

and propagator researchers have independently arrived at equations which are iden-

tical to those derived from coupled-cluster EOM theory.

The general form of Green's function in the energy representation is

+ (0 A k) (k) lB|0) (0 B 0 k) (k\A 0)
((A; B))E = lim 1 + EklA0 (5.1)
E 77-+0+ Eo Ek + E + ir k Ek-EO E--

where the states |k) are eigenfunctions of the full Hamiltonian, H. The electron

propagator is a matrix G of elements

GP =- (ONpEkN- ) k N-i) (Ok-1 q 0O) (0 q ikN+) (kN+llPt .0)
Gq = lim + (5.2)

-+ k E EkN- + E + i77 kN EkN+ Eo + E i7I

where the IkN-1) and IkN+i) functions are, respectively, electron-detached 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 IkN-1) (kN-i 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


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)
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 one-electron 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 one-electron operators,

p and q, and its orthogonal complement f, which contains all higher odd-electron 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 G-1(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)


Go'(E)p, = (E c)6pq (5.13)

E(E) = (aljHa)cr + (alff) (f (El + R) f) l(flfa) (5.14)

(p|Hq)corr = (pIIHq) ep6pq.

In this formulation, E(E), the self-energy, 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)


i(oo) = (alfa)corr

E'(E) = (afl^t(fif(Ei+ )'(fI^a).


A great many approaches to the electron propagator are concerned with successive

orders of approximation to the energy-dependent 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 Non-Dvson 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 self-energy expressions[25]. More recently, Ortiz has experimented with the use of

a highly-correlated 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 coupled-cluster wavefunction as the reference, with

two major consequences. The symmetry of the resulting matrices is broken, and the

IP-type manifold is decoupled from the EA-type 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 function-based approach which gives results identical to the

EOM approach is the coupled-cluster 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 coupled-cluster

equations. The coupled-cluster Jacobian (the first derivative of the coupled-cluster

equations with respect to the cluster amplitudes) used in their work is identical to

the coupled-cluster 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.



The many approaches to calculating ionization energies can be broken down into

two main categories state-by-state 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. State-by-

state approaches can use any SCF or correlated (CI, perturbation, or coupled-cluster

theory) method in the calculation of differences. Among the direct approaches in-

clude two hole-one particle CI[28, 6], the symmetry adapted cluster-CI approach of

Nakatsuji[29], the Fock-Space Coupled Cluster approach[10, 13, 30, 14], and one elec-

tron propagator approaches[31, 32, 23, 33], as well as the coupled-cluster equation-of-

motion approach, which has been derived independently several times, and is known

variously as the coupled cluster linear response function[16] and the coupled-cluster

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 closed-shell 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 well-described

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 inner-shell ionization is desired, there is no good way to characterize the

excited SCF state. SCF programs will quickly collapse back to the valence-ionized

state. One remedy for this problem is to use a quasi-restricted Hartree-Fock (QRHF)

reference, in which an SCF calculation is performed on suitable reference such as

the closed-shell neutral, following which the orbitals are re-occupied and the Fock

matrix recalculated. This non-Hartree-Fock reference is then used as the basis for the

correlated (typically coupled-cluster) 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 high-spin 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 N-electron species, while

adiabatic ionizations involve finding the potential energy minimum for each of the

N 1-electron states.

In most cases, a single-reference 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 open-shell systems, however,

it becomes necessary to consider a multi-determinant reference. For example, the

excited states of C2 accessible from ionization of the C2 ion[35, 36] include open-

shell singlet states. Attempting single-reference 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 high-spin determinant calculation)

are too large to ignore. For such states, the multireference coupled-cluster (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 single-reference

coupled-cluster 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 spin-projection, approaching a

pure-spin 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


6.2 Direct Approaches

6.2.1 Two-Hole One-Particle 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 closed-shell 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 Spin-Adapted 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 spin-flip terms as a consequence of the spin


The SAC- configuration interaction (SAC-CI) method[7, 8] is very much akin to

the EOM approach, but using explicitly symmetry-adapted EOM-type 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 SAC-CI 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 di-ionized reference. Bernholdt [40], finding that the structures

of the low-lying 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 well-behaved, closed-shell C5 as

the starting point, in order to examine portions of the potential energy surfaces of all

three ionized states.



The formalism developed by Koch and Jorgensen[16] for coupled-cluster 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 "free-electron" space, and a new wavefunction of

the appropriate form. In the following sections, a linear response theory formalism,

based on the time-dependent 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 coupled-cluster quantities defined as

the zeroth-order solutions, and the full wavefunction expanded in terms of corrections

to the coupled-cluster 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 first-order correction to the Hamiltonian, V1P.

The full Hamiltonian becomes



2. The zeroth order wavefunctions and energies are those from the coupled-cluster

and A solutions;

K(o) )





1 0) = CCO) = eT 10)


= 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)


The full wavefunction can be written as

IK) eT- HF) (7.9)

= e ()*e ((l+ 2) e+...)+-) HF)






which works out the zeroth order wavefunction plus a sum of perturbative cor-


= 1K0)0+ K'C) + +c2))

= |CC) + K~ (2) + ..j.

= T)IK(O))

- (T(2




+ T12") :K(O)
2", >


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




4. The perturbation operator, VIP, must allow first-order corrections to T and A

to be defined as follows:

TF)) tati + j tatibtj + ...
ai ijab



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 "free-electron" 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 first-order corrections to the wavefunctions, IC()) and (L() to

be represented as (N l)-electron eigenfunctions of the zeroth-order Hamilto-

nian. Keeping in mind that in the zeroth-order wavefunction coefficients to the

r\(r')t operators are zero, while in the first-order 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 coupled-cluster Response Theory

7.2.1 Time Dependent Perturbation on a Coupled-Cluster Reference

With these definitions, the coupled-cluster 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 first-order

corrections to the wavefunction.

Once the CC solution is found, an operator HN -= e-Tx7 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)e-T 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.





both the left and right eigenfunctions.

i ICC(t)) = H ICC(t))


by (pI e-TK 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);




S-i (LO HSCF CCO) = -i (Ol HlN 10)= 0;



= [T dT() Ty CC)


= -i [(L [Hl, T1)] I0) + (Fl VI CC)] ;


where (i -= (1pi e-Tk HIN = e-Tr 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) e-T(HSCF + VIP)eT,

can be derived with the help of equation 7.25

= i (CI (HSCF + VIP) T, ICC) i (01 (1 + AK) e-Tc (HSCF + VIP) ICC)

= i(01 (1 + Ar) e-T [(HSCF + VIP), ,] eT' 10),

d = i (01 (1 + A() [HN 7] 10) = 0,
dt = i(0 (1+ A() [[ftN,)r+] ,Tk], ]0)+i(0I|A [HN,7 ] 10)

+i(01 (1 + A()) e-Tc [VIP, ,] eTx) 10).





7.2.2 The Linear Response Function

Following the development of Olsen and Jorgensen[41], as applied by Koch and

Jorgensen to the coupled-cluster 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)

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 non-infinite 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


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)


Or, in response function notation,

(C 0 CCo) + dw, e(-iw i ((0; VW )),

+ -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


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.



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 = (vje-T )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 first-order 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



ia)) x


E(-A wI)- (a O CC)

+ E (K0 [0, 7,] cco) E(-A -+ wlI)U (7.44)

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 coupled-cluster 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.

The X(') and Y'() vectors can also be transformed:

X1)(wi) = S-'X()(wi) = S-'(-A + wI)-IS S-l '(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 S-1 (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,)


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)
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)

o oo (0ol [[O, 7n], Tk kccO)( i O CCO)
0-k (( [ O, rk] CCO) (

= (0 (1 + AO)(Ork) 0)
S (0 (1 + AO) (OrnTk)c 10)(n 0 I0), (7.52)
n (wk-- (7+ 52)

rF k = (k V CCO)
k---+O I

are referred to as the transition matrix elements and

o = e Oe = (Oe1 )c

V k = e-TIC) 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 left-hand solutions. The calculation of the transition expectation value is more

complicated because of this asymmetry.



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 one-body operator. The second term of FOOk contains the integral

(0 (1 + AO) (Onrk)c)O) ,

which includes a product of two one-body excitation operators (m- and Tk) contracted

with a one-body 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 one-electron function

must be known in order to calculate F_ k and FkV Free-electron 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 so-called 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 non-negligible 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 1|f(N)) } itX 1f(N)) (i p Ix). (8.3)

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)

.(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 multi-reference 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 Non-Symmetric 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)e-TOeTXty 0)



0 ------ ^ -------

0 ao ---^ --- 0

Figure 8.1: Diagrammatic expression of r0_, terms

= (i O x) Dix + (a0 I) Dax
i a

+ (al x)

(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
r ax e = ebc ic
Fbcax = 2 ^ mnCmn"
T2 mn

The other matrix element is:

k-0o = (01 YK~x-TOeT I0)



0--- V


Oo-- ....-

k 0

Figure 8.2: Diagrammatic expression of P terms

= Z(x|Oli)Dxi+ Z(x| Oa)Dxa (8.11)
+ a) (o----- O--b) Fiba + (i i Xa + 0(b- -c) bcx

a bc /


Dxi =

Fn2 -- V0 :*, 1 -------ta
m mne
-m ( +mvm m
= 2im -- j ri
----- 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.

____ ~ 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-


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
I(K) Ci2 '

and the next simplest approximation is:

() ((C + C (c$ ctCmti) /
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 low-incident energy photoelectron spectra.

I(A) ( c. + E ( cA c A). (8.18)



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

non-symmetric 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 non-symmetric matrix A

p(y, x) = x (9.1)

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 bi-orthogonal 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)
Yi = DikLk. (9.8)
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 bi-orthonormalized and added to the B and L

spaces, after which the entire procedure is repeated.

If only the right-hand eigenvectors, X, are required, there is another approach

which can be used to save computational time. The right-hand 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 left-hand 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

"pre-load" 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


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 Hartree-Fock 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 excited-state

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.



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 well-behaved closed-shell 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 s-tetrazine example.

In the larger basis sets, there is a trend for higher-energy 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


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
s-tetrazine 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 (DZP-FCI)
2B2 14.91 14.85 (DZP-FCI)
2A 22.55 22.14 (DZP-FCI)
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 shake-up states, and is not plagued (as is FSCC) by

intruder-state problems in the inner-valence 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 well-correlated single-state 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 open-shell, 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 open-shell state.

While it is a simple matter to remove an electron from the highest-energy orbital

of any symmetry, it requires a bit of trickery and a few more SCF iterations to

converge on a core-hole state. And it is generally impossible to converge an SCF

state with an electron removed from an inner-valence orbital. The QRHF reference,

a non-SCF 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 closed-shell 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 equation-of-motion

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 well-understood 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 II-based 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 experimentally-observed 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 ACES-II 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). Self-consistent 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 inner-valence region correctly.

Two basis sets were used in the nitrogen calculations. The smaller basis was based

on the triple-zeta 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 polarized-valence triple-zeta 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 experimentally-obtained 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 outer-valence 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
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



40.189 38.695
426.723 410.587
426.623 410.768



15.499 (.880)
17.159 (.918)
17.159 (.964)
18.674 (.858)



15.666" 15.5
17.227 16.8
17.227 16.8
18.908 18.6


38.106 (.470)



1 7r2

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,

a EOMEE, except where otherwise noted

In the plots in figure 10.1, Lorentzian lineshapes have been applied to the EOMIP

peaks to simulate the line broadening found in low-resolution experimental spectra.

Each EOMIP plot is superimposed on a plot of an experimentally-derived 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


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 inner-valence 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 doublet-with-a-shoulder 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
C.S. C.H. Exp
15.683 15.499 15.5
17.332 17.159 16.8
18.888 18.674 18.6
29.887 28.917 28.2

32.153 32.679
38.616 38.695
410.623 410.768

31.983 31.5

38.106 37.3



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



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





30o 1C 10; 20u 20a


Expt ---



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 ---

I |

I '

1 2 N2 Ionizations, EO MIP
N2 lonizations, EOMIP

Expt ---30g 1"u a.* 20;

Ii i

i I
p I'


J &


' ---- ---



I ,

significant EOMIP peak at 38.6 eV when the 5s4pld basis is used. This region shows

much less intensity in the pVTZ+ closed-shell or core-hole orbital calculations, in

agreement with experiment. In the core-hole 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 core-hole basis, moving closer to experiment. Interestingly,

core-hole 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 core-hole as in the closed-shell 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 outer-valence 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 outer-valence results show excellent

agreement with experiment, while the inner-valence principle peak at 24 eV, like

the SAC-CI, 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 SAC-CI 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


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, 10-30 eV

Ionization energy eV (Intensity)
EOMIP SAC-CI" Experimentb Orbital Assignment

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

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.81 (.006)

28.18 (.001)
28.6 (.017)
29.26 (.044)


30.77 (.075)


27.4 (.29)

29.76 (.010)

30.06 (.094)

29.28 (.02)





31.2 (.1)d


2bt lb3glb3u
7at 1b3ulb3
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




0. II

0.3- I I I
I I 1
0. I I

Figure 10.2: EOMIP and SAC-CI 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 Rydberg-supplemented 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 low-lying 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 molecule-centered 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 over-estimate

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 coupled-cluster triples correction can hardly be

said to be much help in improving either of the semi-QRHF calculations. In the pVTZ

closed-shell 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
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 half-full 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. EOMIP-SD 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

closed-shell 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 closed-shell 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 closed-shell 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 SAC-CI 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 shake-off effects double ionizations. Because the PS (photoemission

spectroscopy) experiment was performed in the x-ray 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.