Triple excitation effects in the fock-space coupled cluster method


Material Information

Triple excitation effects in the fock-space coupled cluster method
Physical Description:
viii, 131 leaves : ill. ; 29 cm.
Bernholdt, David Edward, 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, 1993.
Includes bibliographical references (leaves 125-130).
Statement of Responsibility:
by David Edward Bernholdt.
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 - 001890698
oclc - 29642259
notis - AJW5922
System ID:

Full Text







To Jan,
for so many reasons.


Graduate school has been an educational experience in many ways, and a great many

people have made important contributions to that education.

I would especially like to thank my advisor, Rod Bartlett, for the singular role he has

played in my life. His research group members, especially John Stanton, have also shared

knowledge and offered both help and challenges throughout my work with them. My

colleagues, past and present, in the Quantum Theory Project have created a wonderful

environment in which to learn and do research, and it has been my great pleasure to

experience it. I would like to thank Drs. Vala and Weltner for numerous discussions

and suggestions; I regret that I did not have time to examine more of their interesting


My research in this field began while I was an undergraduate at the University of

Illinois, and I am deeply indebted to Cliff Dykstra and Shi-yi Liu for getting me started

and for conveying to me their enthusiasm for the subject.

The daily existence of a graduate student often revolves around research, but there

is much more to life. I thank my friends and family, who have tried to tell me this

and who have shared other gifts with me over the years: humor, diversion, perspective,

encouragement, sympathy, and love.

Finally, words cannot express my gratitude to Jan, who rode two different roller

coasters at the same time and lived to tell the tale.



ABSTRACT .................


1 INTRODUCTION .........


. 5

2.1 Overview ................

2.2 The Schrodinger Equation ......

2.3 Potential Energy Surfaces ......

2.4 Calculations on Potential Energy Sul

2.5 Vibronic Interactions .........

2.6 Basis Set Expansion .........

2.7 One-electron Approximations ...

2.8 The Correlation Problem .......

2.9 Occupation Number Representation.


3.1 Introduction ..............

3.2 Basic Equations ............

3.3 Reduction to Spin Orbital Equations


. . . 29

. iii

. vii

. . . .


3.4 Numerical Solution of the Coupled Cluster Equations

3.5 The Coupled Cluster Effective Hamiltonian .

3.6 Truncation of the Cluster Operator . .


4.1 Introduction ........................

4.2 Structure of FSCC ....................

4.3 Comparison with Single Reference CC Theory .

4.4 Basic Equations .....................

4.5 Reduction to Spin Orbital Equations . .

4.6 Approximations to the Full FSCCSDT Model .....


5.1 Introduction ........................

5.2 Effect of Changes in the Active Space . .

5.3 Ionization Potentials ...................

5.4 Potential Energy Surfaces and Related Properties .

5.5 Conclusions ........................


6.1 Background ........................

6.2 Computational Methods .................

6.3 The Heart of the Problem ................


. . 33


. 33

. 38

. 42

S. 42

S. 43

. 50

S. 54

. 55

. 56

S. 64

S. 64

. 66

S. .71

. 79

S. 87

. 90

S. 90

. .91

. 93

6.4 Geometries of the Linear Stationary Points of C5 . ... .95

6.5 Relative Energies ........... ........ ... ................97

6.6 Stability of the Stationary Point Structures . ... 102

6.7 Implications for Other Carbon Clusters . .... 105

6.8 Recapitulation ................................. .. 108

7 SUM M ARY ..................................... 109



REFERENCES ........................ ................. 125

BIOGRAPHICAL SKETCH ........................ ......... 131

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy




May 1993

Chairman: Rodney J. Bartlett
Major Department: Chemistry

The Fock-space coupled cluster (FSCC) approach to the electron correlation problem

offers several advantages over the more familiar single reference coupled cluster (CC)

methods. It is capable of treating states with a multireference character, and it is capable

of providing state-to-state energy differences for several states in a single calculation. As

explained in this work, basic differences between the two theories result in higher-excited

cluster operators being more important in FSCC methods than CC methods.

In order to improve the accuracy of FSCC calculations, this work involves the

addition of triple excitation effects to the widely applied FSCC model restricted to

single and double excitations (FSCCSD) for ionization potentials. Equations for the

full model, including single, double, and triple excitations (FSCCSDT), are derived and

various approximations are proposed as a compromise between the improvement in the

wavefunction and energy due to the inclusion of triples and the computational costs,

which are significantly larger than for the FSCCSD model.

Two approximations to the full FSCCSDT model have been implemented and their

performance is examined for a variety of molecules. Results are found to be generally

good, providing results comparable to single reference coupled cluster methods with

approximate inclusion of triple excitations.

Finally, both CC and FSCC methods are used to study the structures of several

small carbon clusters, a class of molecules which has been shown to require high-quality

theoretical methods and which can make use of both FSCC's multireference and direct

energy difference capabilities. Detailed calculations show that the ground state of C5+ is a

2E+, with Coov symmetry, but with low frequency bending vibrational modes. Additional

calculations for other odd carbon cluster cations suggest that the "linear" structures of

C7' and C9' observed experimentally may in fact be distorted from true linearity due

to the Renner-Teller effect.


In recent years, advances in both experimental and theoretical chemistry have fueled

an exciting synergy between the two fields [1]. Experimental innovations have allowed

study of molecules in conditions similar to those commonly employed in theoretical

calculations (low temperature, low density), producing data of increasing accuracy and

specificity. Likewise, through increases in computational power and methodological

sophistication, theoretical chemists have been able to tackle increasingly complex systems

with high accuracy. Where experimental and theoretical data are directly comparable,

critical evaluations of both approaches are possible. In other cases, the data are often

complementary, allowing one or the other researcher to bring the problem into finer focus.

This has been especially apparent in research on atomic clusters, a field which has

expanded rapidly in the past decade. Clusters are interesting for a variety of rea-

sons. Metal clusters are used as models for the study of catalytic processes. Clus-

ters of lighter atoms, particularly carbon, have been identified in circumstellar space

and are expected also to be found in interstellar molecular clouds, where they may

be important in the creation of new stars. They have also been found in combustion

products on Earth. Clusters may have useful and unique physical properties, mak-

ing them interesting from a materials viewpoint, and they are studied as models of

how properties of bulk materials are approached by increasingly large collections of



Although the literature on clusters is expanding rapidly, the sought after understanding

has not yet emerged. One reason for this is that clusters, even relatively small ones, are

rather unlike "traditional" molecules which scientists have studied for centuries, and

they often defy the "rules" and intuitions which researchers try to bring from traditional

chemistry to cluster chemistry. Consequently, workers in the area are having to develop

new ideas and new "rules" about the structure and properties of clusters.

An important component of creating this new understanding is having reliable data

on the structure and properties of many clusters. Since there are few experiments which

give direct, unambiguous structural data, theoretical approaches are particularly useful in

this context. When the capabilities of experiment and theory do overlap, it is possible

to compare results critically and perhaps improve both. To be truly useful, however,

the theoretical methods employed must be accurate and reliable. The subject of this

dissertation, the inclusion of triple excitation effects in the Fock-space coupled cluster

(FSCC) method, was undertaken in response to exactly this need.

Coupled cluster (CC) and related many-body perturbation theory (MBPT) methods

(see Chapter 3) have proven their ability to deliver extremely high-quality properties and

energetic for a wide range of complex systems, and because of their size extensivity

(proper scaling of the energy with the size of the system under study) [2,3], are preferred

over configuration interaction (CI) methods for theoretical studies of clusters. Not

every situation, however, is ideally suited to the use of traditional single reference CC

methods. For example, close-lying electronic states may cause problems for even the most

sophisticated CC methods, thus requiring a multireference description of the wavefunction


to obtain meaningful results. Ionization potentials (IPs) and electron affinities (EAs) are

often obtained using experimental studies of clusters in photoionization and photoelectron

spectroscopies. In order to provide comparable theoretical results for interpretation or

confirmation of the spectra, separate calculations for the ground state and each excited

state of interest are required when using familiar CC methods. It would be useful to

have a method capable of giving IPs and EAs for several states of interest in a single

calculation. Fock-space coupled cluster methods are one way to approach both of these


Often referred to as Fock-space multireferenme coupled cluster theory, FSCC (see

Chapter 4) can provide a multireference description of chosen electronic states of a

molecule. Considered in another context, it can calculate the energies of several electron-

detached or -attached states at the same time, thus simplifying the calculation of IPs and

EAs. Because of formal differences between single reference CC and FSCC theories,

however, a given truncation of the two methods does not necessarily give results of

equivalent quality. For example, limiting both methods to single and double excitations,

we find that an FSCCSD description of a particular state may leave out contributions

which could be present in the CCSD description of the same state. Consequently, in

order to improve the quality of the FSCC description, this work examines the inclusion

of triple excitation effects in the FSCC method.

After setting the stage with a review of theory and terminology leading up to the

correlation problem, the traditional coupled cluster equations are presented, emphasizing

their use as an "input" to the FSCC method. Then, in some detail, equations for the FSCC


method are derived, focusing on the FSCCSDT model for IPs. The importance of triple

excitations is examined, in comparison with single reference CCSD and CCSDT methods,

and several approximations are proposed as a compromise between the computational cost

of evaluating the complete triples contribution and quality required of the wavefunction

and energies. Although this development was spurred by an interest in carbon clusters,

FSCC methods, with or without triple excitations, are by no means limited to this

application. In order to assess the performance of the FSCC methods, calculations are

presented for a variety of molecules. Finally, Fock-space and single reference coupled

cluster methods are used to study several carbon cluster cations.


2.1 Overview

Much of the theoretical chemistry literature assumes the reader is a practitioner in the

field. This chapter is intended to provide a common ground for the concepts and notation

used throughout this dissertation. It is by no means exhaustive, but should help non-

specialists to understand the material that follows. References are given to the literature

for more detailed treatments of various topics.

2.2 The Schrodinger Equation

We seek to solve the time-independent Schr6dinger equation,

HI}) = El\ ) (2.1)

for a molecular system (one or more molecules of interest isolated from interaction with

anything else). This equation tells us that a system in an eigenstate I ) of the Hamiltonian

operator H has an energy E. The usual nonrelativistic Hamiltonian can be divided into

five different terms,

H = Te + Tn + Ven + Ve + Vnn (2.2)

representing, respectively, the electronic and nuclear kinetic energies and the potentials

due to interaction of electrons and nuclei, electrons with other electrons, and nuclei with


other nuclei. The individual terms, in Hartree atomic units, are defined as

T(r) =-
1 2

T.(R) = AV

Ven(r, R) = ZA (2.3)
AIri RAI(2.3)

1 A

Vee ( ) = 1-- 1

Vn2(R) = E IRA RBI

In these expressions, i and j refer to electrons, A and B refer to nuclei. RA and ri refer to

the spatial coordinates of nucleus A and electron i, respectively, and Z, and MA are the

charge and mass of nucleus A. The unsubscripted symbols r and R refer to the complete

set of position vectors of the electrons and nuclei, respectively.

2.3 Potential Energy Surfaces

The terms of the full Hamiltonian, Eq. 2.2, can be grouped according to their

dependence on the positions of the nuclei into a purely electronic component,

He = Te + Vee, (2.4)

a purely nuclear term,

Hn = Tn,


and a potential which couples them,

Vnt(r, R) = Ven(r, R) + Vn(R). (2.6)

If we choose one set of positions R as the origin, denoted by 0, we can expand the

electron-nuclear interaction potential near this point in a Taylor series,

Vi (r,nR)= V t(r,o0)+ :(RA A )A B RA+t ) AR +-. (2.7)
The first term in the expansion can be thought of as the potential felt by the electrons

due to a fixed set of nuclei, leading to the electronic Schrodinger equation,

(He(r) + V.t(r,0))lPk(r)) = E|'Ik(r)), (2.8)

for a particular state |Ik(r)) with energy Ek.

In order to see how varying the nuclear positions effects these solutions, we must

return to the full Schrodinger equation 2.1. By expanding the full wavefunction in terms

of the solutions to the electronic wavefunction and a nuclear component,

I| (r, R)) = 1 Dk(r))Xk(R), (2.9)
and recognizing that the Hamiltonian of Eq. 2.8 is a subset of the full Hamiltonian,

we obtain

{He + Vnt,o + H, + W} IE Ik)Xk = EE |kk)Xk, (2.10)
k k
where V,,t,o denotes the leading term in the Vnt expansion, and Wdenotes the remainder,

W(r, R) = Vint(r, R) Vnt(r,0), known as the vibronic interaction. He + Vit,o is often


referred to as the electrostatic Hamiltonian, and is dependent on the nuclear coordinates

only parametrically. Projecting on the left with an electronic eigenfunction, (4m(r)l,

we obtain a coupled set of equations, (note that H, is independent of the electronic

coordinates r),

(E" + H,)Xm + E( Im Wlk)Xk = EXm, (2.11)
or using the more compact matrix element notation for W,

(E E' H,)Xm = 1 WmkXk. (2.12)

If there is no vibronic interaction between electronic states (Wink = 0), then the

equations decouple, giving

(H.(R) + Em)Xm(R) = EXm(R) V m, (2.13)

which is the Schridinger equation for the nuclei moving in the mean field of electronic

state I\4m(r)), and the solution to the full Schridinger equation is simply I|(r, R)) =

I m,, (r))Xm(R). This is the origin of the idea of a potential energy surface (PES) as the set

of nuclear configurations which satisfy Eq. 2.13 for a particular electronic state. Taking

Wmk = 0 is known as the simple adiabatic or Born-Oppenheimer approximation [4,5].

In practice, although Wnk is not in general zero, the Born-Oppenheimer approxima-

tion is usually quite good because of the fact that the nuclei of a molecule are about 1836

times more massive than the electrons, so we can usually think of the nuclei moving

slowly in the average field of the electrons, which are able to adapt almost instanta-

neously to the nuclear motions. This approximation is used in the majority of electronic


structure calculations to date, which simply evaluate the electronic wavefunction and

energy, Eq. 2.8, at various points on the PES (nuclear configurations) and ignore the

nuclear dynamics.

2.4 Calculations on Potential Energy Surfaces

The Born-Oppenheimer potential energy surface for a particular state of a molecule

with N nuclei is (neglecting invariance to rigid translations and rotations of the entire

nuclear framework) a 3N-dimensional space which can be characterized according to its

local shape. A stationary point is one at which the gradient of the energy with respect to

nuclear displacements be zero. Stationary points may be minima, transition states, higher-

order saddle points, or maxima. If the curvature (second derivative) at a stationary point

is positive, the point is a local minimum, which is a locally stable geometric conformation

for the state in question. If the curvature in one of the 3N directions is negative, with

all others positive, the point is a transition state, which is a local maximum on a path

connecting two local minima; transition states are very important in the theory of chemical

reactions. If the stationary point has more than one direction of negative curvature, it is

a saddle point (or local maximum), which is not chemically meaningful.

Stationary points on potential energy surfaces can be located using numerical opti-

mization techniques, such as Newton-Raphson, quasi-Newton-Raphson, etc., but because

individual calculations of the energy, gradient, and hessian (matrix of second derivatives

of the energy), as required by the optimization method, are extremely expensive relative

to evaluation of the optimization step, it is usually necessary to use chemical intuition


to choose regions of the PES on which to focus and sometimes to restrict the search to

a subset of the full 3N-dimensional space.

The hessian, which allows a stationary point to be characterized as a minimum,

transition state, or saddle point, is also the primary input for evaluation of the vibrational

frequencies at that structure (in the harmonic approximation). Since different molecules

and different structures will generally have different vibrational spectra, it serves as a

useful tool for the identification and characterization of unknown species.

If the electronic state of the molecule of interest is not known from other information,

it may be necessary to consider potential energy surfaces for several states in order to

locate the lowest energy global minimum of all the states, which is known as the ground

state of the system.

2.5 Vibronic Interactions

Although the Born-Oppenheimer approximation and the idea of potential energy

surfaces have served quantum chemistry well, there are certainly cases in which the W,

matrix elements connecting different electronic states are significantly different from zero,

and the Born-Oppenheimer picture is no longer adequate.

Consider the vibronic interaction term, W, truncated after the quadratic term, which

is adequate for purposes of explanation. Making a linear transformation from cartesian

coordinates R to symmetrized normal coordinates Q (obtained by diagonalizing the

hessian matrix), we have

W(r,1 Q) = v (ai.tQ + Q 1 Z 0 t QaQ6. (2.14)
: \d_2 EQ ) + 0

Each normal coordinate, Qc, belongs to an irreducible representation (irrep) Fa of the

molecular point. Likewise, the electronic states I|m) belong to an irrep, Fm. Thus from

group theory, we know that the linear vibronic constant,

F(kFm) =- (2.15)
F\ Q\9aQ I

is nonzero if and only if rFo [Fk x Fm]. This is a general result for the interaction

of any two states, but of particular importance is the case in which the two electronic

states belong to the same irreducible representation. If the irrep is not degenerate, then

[Fk x rk] = Fk x Fk = FA1, so that the coupling is via totally symmetric motions, which

by definition cannot change the symmetry of the molecule. If, on the other hand, the

two electronic states (or components of states) belong to a degenerate irrep, the direct

product results in non-totally symmetric representations, and thus the states are coupled

by motions which break (lower) the symmetry of the molecular framework. This is the

basis of the Jahn-Teller (JT) theorem [6], which says that for a nonlinear polyatomic

molecule in an n-fold degenerate electronic state, the high symmetry point at which the

n potential surfaces intersect cannot be a minimum on all n surfaces. In other words,

there will always be at least one direction for which the energy is lowered when the

degeneracy is lifted (by moving away from the high symmetry point). The tendency of

a molecule in a degenerate state to lower its symmetry and thereby lift the degeneracy

and lower its energy is known as the Jahn-Teller effect.


In the more general case, interaction to two near-lying states via a suitable nuclear

motion causes the surfaces to repel each other. Looking at a slice of the two surfaces

along the symmetry-breaking motion that couples the two states, the upper surface

becomes more curved as the interaction increases, while the lower state is flattened.

If the interaction is strong enough, the lower state may display a double minimum shape,

with the original high symmetry point a transition state between the two minima. This

is usually referred to as the pseudo-Jahn-Teller effect (PJTE) [7]. The magnitude of the

PJT effect depends both on the separation of the two states and the magnitude of the

linear vibronic constant.

The exclusion of linear molecules from the JT effect is also on the basis of symmetry

arguments. The wavefunction of a linear molecule is either even or odd with respect to

reflections containing the axis of the molecule or perpendicular to it, so of course the

product of such functions will always be even. The nuclear motions that bend the

molecule are odd in this respect, and so can never be contained in the direct product

of the wavefunction representations the linear vibronic constant F(k k) is zero for

linear molecules.

For linear molecules, and other cases for which the linear vibronic constant is zero

or negligibly small, the quadratic term of Eq. 2.14 comes into play. For linear molecules

in degenerate states, the quadratic vibronic coupling leads to a condition similar to the JT

effect, called the Renner-Teller effect (RTE) [8,9]. Lifting of the degeneracy can be easily

visualized in this case. Consider a linear molecule oriented along the z axis, in a 2II state,

so that one component of the 7r orbital, say xrx contains an electron. Bending motions


in the orthogonal directions x and y would normally be equivalent in a linear molecule,

but if one component of the 7r is occupied, the electronic state of the bent structure will

be different depending on whether the electron is now in an orbital in the plane of the

molecule or perpendicular to it. Unlike the JT case, the Renner-Teller effect does not

guarantee that lifting the degeneracy by bending will result in a lower energy. The two

states will not be degenerate (except the linear geometry), but neither, one, or both may

still be minima at the linear geometry, or the linear geometry may be a transition state

between two equivalent bent minima depending on the magnitude of the interaction.

In the more general case, with states belonging to different irreps interacting, there is

little to be gained from a symmetry analysis of the quadratic vibronic constant. There are

usually many combinations of vibrational motions that can couple two electronic states

for either degenerate or nondegenerate point groups, however these effects are usually

relatively small.

Clearly vibronic coupling effects can cause major problems for computational meth-

ods based on the Born-Oppenheimer approximation, but if we are looking for minima

on potential surfaces, it is clear that such points are usually not of interest when the

interactions are strong enough, they cause the symmetry of the molecule to be broken,

thus moving it away from the degenerate or pseudo-degenerate point. As a result, it

is often enough to recognize points on the PES at which strong vibronic interactions

occur and then locate and follow the symmetry breaking to lower energies. Of course

this approach will not always work, and even calculations near a (pseudo-) degenerate

point can be quite taxing for the method, but it has been used in a number of studies

(see Bersuker's bibliography [10] for examples). In 1984, Lee et al. [11] demonstrated

that analytically evaluated hessians could be used to compute bending frequencies for

Renner-Teller molecules. In fact if due care is taken to stay on the correct PES during

the calculation, hessians evaluated from finite differences of gradients, or even energies,

can also give RT bending frequencies.

2.6 Basis Set Expansion

Several approaches may be taken to the solution of the electronic Schrodinger

equation, but the most widely used and most generally applicable is expansion in a

basis of Gaussian functions. These functions are usually (but not necessarily) chosen

to mimic the atomic orbital (AO) description used in chemistry and are referred to by

that name. A complete (infinite) basis would span all of space and thus allow an exact

description of the wavefunction. In practice, however, computational resources place

limits on the size of basis that may be employed, and it is necessary to compromise

between the cost of the calculation and the accuracy required.

It is now possible to evaluate matrix elements of the Hamiltonian operator in terms

of the AO basis set. For a set of functions {x,}, we have, for example, the one- and

two-electron integrals (or matrix elements) [12]

(lhlv) = h, = J (xi)hx;(xl)dxi
(Pv" Ar) = X,(xl)X*(x2)- X(xA)x(x2)dxd
I r12

r{ a(wi), or
xi = () (2.17)

is a combination of the position and spin of electron i, and r12 = Irl r2l. The left-hand

sides of Eq. 2.16 are shown the bra-ket notation, originally due to Dirac [13], which

is in common use and for the one-electron integral, a matrix element notation is also

shown. Both the bra-ket and matrix element notations are used in the literature, and both

will be used in this work.

2.7 One-electron Approximations

Given a basis of atomic orbitals, which are (very rough) approximations to

the eigenstates of the atom, the first step in an electronic structure calculation

is usually to find a set of one-electron functions for the molecule as a whole,

called molecular orbitals (MOs). These orbitals correspond to the qualitative ideas

about molecular orbitals often used by chemists.The most common prescription for

calculating MOs is the Hartree-Fock (HF) Self-Consistent Field (SCF) procedure


The result of an SCF calculation is a set of MOs {(,}, which are represented as

linear combinations of the AO basis functions,

p = CSPXA (2.18)

with the matrix C determined by the SCF procedure. Of these MOs, some will

be occupied by electrons and some will be unoccupied, or virtual orbitals. The


antisymmetrized product of the occupied orbitals is the Slater determinant for that


Depending on the particular formulation (primarily the treatment of electrons with

different spins) used in the HF procedure, it may be referred to as RHF, UHF, or

ROHF, among other possibilities. The most straightforward is restricted Hartree-Fock

(RHF), applicable only to closed shell molecules (where all electrons are paired), where

the alpha and beta spin orbitals in each pair are restricted to have the same spa-

tial function. The restricted open-shell Hartree-Fock (ROHF) formalism extends the

RHF idea to open shell systems by requiring that all paired electrons have share the

same spatial function, but allowing additional orbitals which may be occupied by sin-

gle unpaired electrons. In the unrestricted Hartree-Fock (UHF) method (also known

as DODS different orbitals for different spins), the alpha and beta spin electrons

have distinct spatial parts, allowing the possibility that they might be quite differ-


Only when pairs of electrons have the same spatial component (RHF, ROHF) is it

guaranteed that the HF wavefunction is an eigenfunction of the S2 operator, satisfying

S21HF) = s(s + 1)IHF). UHF wavefunctions may suffer from spin contamination,

in which the wavefunction is an admixture of states of different S9 values. Because

the true wavefunction is always an eigenfunction of 52, it is thus possible that a UHF

wavefunction may be to some extent unphysical. Nevertheless, the UHF method has

proven quite useful, particularly as the basis for formulating general (i.e. applicable to

both closed and open shell systems) correlated methods.

2.8 The Correlation Problem

Although the SCF procedure provides a very useful qualitative description

of molecules, it is generally inadequate for quantitative uses requiring high ac-

curacy. The method considers each electron the average field of all others,

which ignores the fact that the motion of each electron is instantaneously cor-

related with all others because of the Pauli Principle. For the higher ac-

curacy required of most interesting applications, it is necessary to go be-

yond the one-electron approximation and treat correlation effects in the sys-


The correlation problem is usually formulated in terms of the interaction between

different configurations (occupations) of a set of one-electron functions. The molecu-

lar orbitals are usually used for this purpose because the idea of occupied and unoc-

cupied orbitals is well defined and the derivation of equations for correlated methods

are thus greatly simplified. The SCF determinant serves as a reference for the defi-

nition of the excited state configurations in which the correlated wavefunction is ex-


2.9 Occupation Number Representation

The occupation number, or second quantized, formalism [14] is a convenient

representation for the development of correlated methods because it represents in a

compact way the changes in occupation of the various determinants. In this de-


scription and throughout this dissertation, the convention will be used that subscripts

i, j, ., n refer to orbitals which are occupied in the reference determinant, a,

b, ., f refer to virtual orbitals, and p, q, ., u refer to orbitals of either


Creation and annihilation operators are defined that act to create or destroy an electron

in a particular one-electron function (MO). Thus the creation operator at would create

an electron in orbital p, and the annihilation operator aq would remove an electron from

orbital q. The Pauli Principle is satisfied by the fact that aa = aq = 0, which is a

result of the anticommutation relations among the creation and annihilation operators.

Strings of creation and annihilation operators act on the reference determinant, 10),

to produce excited state determinants. For example, an electron can be promoted from

an occupied orbital i to a virtual orbital a,

aai|l0) = |I0) = |') (2.19)

where the resulting excited state has been shown in two common notations. Similarly,

higher excitations may be obtained
ata"aajiO) = |aa0 = |0)

a aka aaa l) = ~i ) = |i?) (2.20)

Operators such as the Hamiltonian are represented in the second quantized formalism

by the matrix element of the operator in the basis of one-electron functions multiplied

by a string of creation and annihilation operators:

H hpqaaq + (pql rs)ataasar (2.21)
p,q p,q,r,s
where we have used the antisymmetrized two-electron integral, (pqllrs) = (pqlrs) -


Strings of creation and annihilation operators may be evaluated using Wick's theorem

[15] and the concept of a normal order of operators. Normal ordered operator strings are

those which have all at and all aa to the right of the string, which can be accomplished

using the anticommutation relations for second quantized operators. This is done to

exploit the fact that a!10) = aalO) = 0. The normal order of an operator string A is

denoted N[A], and has the property

(OIN[A]I0) = 0. (2.22)

This result is obvious if all a! and aa operators are brought to the right, where they act

on the reference giving zero immediately. If there are no at or aa, then the operators

act on the reference to create an excited state on the left hand side, which is orthogonal

to the reference function on the right,

(/O| \ =0. (2.23)

Note too that in the absence of at or aa, A=N[A].

Wick's theorem tells us that the contraction of a pair of operators is defined as

aq = N t aaq] + N a aq or
apaq = N apart ] + N apa ,


with the following definitions according to whether p and q are occupied or virtual orbitals

(they must be the same type since 6ai = 0)

N [aaj]= 6ij

N[aia = 0
N aab =0

N [aaaa = ab

Longer strings of operators are evaluated by taking all possible single contractions, then

all double contractions (two simultaneous contractions), up to full contractions, with the

sign in each case determined by the permutations necessary to make pairs of contracted

operators adjacent. As an example, consider the Hamiltonian, Eq. 2.21. The one-electron

term is cast in normal ordered form by recognizing that

ataq = N a4aq + N [aq]
= N a[aq] + bij

so that

Shpqaaqy = hpqN [aaq] + hii. (2.27)
p,q p,q i

The two-electron term has a string of four operators, so all possible single and double

contractions must be considered. There are four single contractions all of which can be


manipulated into the same form, and similarly for the two double contractions,
1 1 (pqj|rs)N[aataa]
4 -p,g,r,s p q 4 q q
+ (pil qi)N [aaq] + ( ij)
p,q,i i,j
When evaluating a product of several normal ordered operators, the same rules hold, but

all contractions which are internal to any given operator have already been accounted

for, so only those contractions connecting different normal ordered operators need be


When developing correlated methods, it is convenient to employ the normal ordered

product form, ON, of a second quantized operator, 0, defined by

ON 0 (01010). (2.29)

This form has the property,

(0oON10) = (010 (01o00)10) = (01010) (01010) = 0 (2.30)

which allows one to focus on the correlation correction to the result rather than carrying

the reference expectation value (a constant) through all of the equations. For example,

consider the Hamiltonian. The normal ordered operator is the sum of Eqs. 2.27 and

2.28, or

N[H] h= hN[aq + ap,qi ]qi)N[a] + (pjjpqr, (Pqlrs)N[ap ,ar]

+ hi + Y (Iijlij)
i i,j

and the expectation value is

(01HI0) = 0 h hpqaaq + P ,,r, (pqrs)aptaaar 0
\ P pqrs q
= hii + i (ijllij)

The normal ordered product form of the Hamiltonian is thus


= hpN [ataq] + ,,(piqi)N [a + pq,,s(pqllrs)N aataa

= fpN [aaq] + (pqllrs)N aaasaT]
S4 p,,r,s Q

with the Fock operator,

fN = fN [ataq] hqN [aa]+ (pi [qi) N[aaq] (2.34)
p,q p,q p,q,s
used in the final form. The normal ordered product Hamiltonian is also frequently written

in terms of its one- and two-electron components,

HN = fN + WN. (2.35)


3.1 Introduction

Two of the most common approaches currently in use for obtaining correlated energies

and properties are configuration interaction (CI) and coupled cluster (CC) methods.

Both methods use an expansion space of excited Slater determinants to represent the

wavefunction; the difference between them originates with the particular form chosen for

the expansion. In the CI case, a linear expansion is used,

'Icl) = (1 + C)0) = (1 + C1 + C2 + .)10) (3.1)

where C, generates all possible i-fold excitations from the reference determinant weighted

by coefficients giving the importance of each determinant in the total CI wavefunction.

In most cases, it is necessary for computational reasons, to truncate the expansion space

somewhere short of including all possible excitations (up to n-fold for an n electron

system). This results in a variationally determined wavefunction (and energy), but

sacrifices the size extensivity property.

Size extensivity [2,3] is the requirement that the energy of a system scale properly

with the number of electrons. It is most easily understood by considering an assembly of

N identical non-interacting molecules. If each molecule, p, has an energy and reference

function satisfying

HIO0) = EP)O()) (3.2)


(the one-electron functions of O1(p)) are assumed to be localized on that molecule), then

we can write the reference energy and wavefunction of the entire assembly as
E(total) _-S EP)
10(total)) = 1 (p))
Since the molecules are non-interacting, we expect the correlated energy and wavefunc-

tions to be representable as similar sums and products, respectively. When a truncated CI

expansion is applied to this system, however, this is not so. For an individual molecule,

p, the coefficients C C) ,... C'P), with m

H(1 + C- ) + CP + .) I0() = E(' (1 + C' + C2' + )0 |O()} (3.4)

and give a product wavefunction which includes excitations up to Nm due to products

of excitations on different molecules. If the assembly were treated as a whole, however,

only up to m-fold excitations would appear in the wavefunction. Thus, they can never be

equivalent. The energy of the assembly suffers similar problems. For example, a lattice

of N non-interacting hydrogen molecules treated with CI restricted to double excitations

(CID) has an energy which scales as /N instead of N, resulting in an error of 12.3% for

10 molecules and 21.1% for 20 [16]. Calculations by Ahlrichs and coworkers at the CID

level show that the difference between supermolecule (the two molecules separated by a

very large distance calculated together) and the sum of individual molecular energies are

-7 kcal mol-1 for BH3 + BH3 [17], and -15 kcal mol-1 for CH3F + F- [18].


The well known property of exponentials, eAeB = eA+B (where [A,B]=0), suggests

an alternative solution to the problem. If the individual molecule wavefunctions are of

the form

I|(p)) = eTP)0(P)) (3.5)

with P" being an expansion operator analogous to C, in Eq. 3.1, and satisfying

He') O()) = E( e P) (p)) (3.6)

then the product wavefunction is


o>P) N (3.7)
.(ep ) = IO(p))
= e 0(-=1 1 10(P))

= ei' 'otal)o0(total))

which is identical to the total wavefunction resulting from a calculation of the entire

assembly. This is the form of the wavefunction employed in coupled cluster theory.

This simple physical picture demonstrates the need for the exponential wave operator,

but it is the Brueckner-Goldstone linked diagram theorem [19-23], which CC methods

satisfy [24], that guarantees size extensivity when this form is used. Thus when the system

is composed of interacting fragments, the simple physical picture becomes clouded, but

the size extensivity property is still guaranteed.

Having motivated the use of coupled cluster methods and the exponential wave

operator, in the remainder of this chapter, we will derive the working equations of


coupled cluster theory and other results which will be needed for the Fock-space CC

methods discussed in the next chapter.

3.2 Basic Equations

In order to expand the coupled cluster wavefunction in the space of Slater deter-

minants, we define a set of cluster operators which are capable of representing the

contributions to the wavefunction due to determinants of particular levels of excitation,

j= = tatai

T2 = -a t attaj
.>b (3.8)

T3 = tibkaaabacaiajak

where the cluster amplitudes, t:I are the scalar coefficients for each excited Slater

determinant in the wavefunction. The cluster operators for each level of excitation are

collected together into a single operator,

T = Ti + T2 + T3 + (3.9)

for convenience.

To obtain the multiplicatively separable coupled cluster wavefunction, the wave

operator (which produces the correlated wavefunction from the reference) is exponential

in the cluster operators,

| Qcc) = QOcl0) = eTO). (3.10)

This wavefunction must satisfy the Schridinger equation

He T0) = EcceO) (3.11)

where E,, is the coupled cluster energy and H is the second quantized electrostatic

Hamiltonian. In order to simplify later discussion, we change to the normal product form

of H by subtracting Eo = (0|HI0) from both sides of Eq. 3.11,

(H (OIHO))eT0O) = (Ec Eo)eTIO)
HNeTIO) = EcorreTIO)
where the correlation contribution to the energy has been designated E.,,. Finally, we

multiply both sides of Eq. 3.12 by e-T to obtain

e-THNeTI) = Ecorr0) (3.13)

Notice that the operator e-THNeT acts as an effective Hamiltonian operator, acting on

the reference state to give the correlation energy; this quantity is often referred to as HN

in post-CC methods (i.e. those for which the derivation relies on having first solved for

a CC energy and wavefunction).

The effective Hamiltonian can be simplified by the use of a Baker-Campbell-

Hausdorff expansion

eTHeT = H + [H,T] + I[[H,T],T] + (3.14)


The commutators of strings of second quantized operators may be evaluated using the

contraction theorem of Harris, Jeziorski and Monkhorst [25], with the result that the

expansion terminates after the quadruply-nested commutator term because H (Eq. 2.21)

contains only four operators. This result is often denoted

e-THeT= (HeT) (3.15)

where the subscript c indicates that only so-called connected terms remain. The ter-

minology arises from the diagrammatic formulation of the equations often used in CC

theory: surviving terms are those in which their diagrammatic representation forms a

connected figure, from which no component may be removed except by breaking lines.

The algebraic analog is that the term cannot be separated into a product of distinct sum-

mations. Disconnected terms are those which can be separated into two components

without breaking a line, or algebraically the term may be cast as the product of two

distinct summations. Use of Eq. 3.15 greatly simplifies reduction of the operator form

of the coupled cluster equations to a spin orbital form suitable for numerical work.

The equations which determine the CC energy and amplitudes are obtained from Eq.

3.13 by projecting from the left with the reference determinant and each class of excited

determinants (single, double, triple, etc.). These projections provide equations which

determine the energy, T,, Ts, T3, etc., respectively.

(O1(HNeT) 10) = Ecorr (3.16)

'I (HeT) o) = 0
K (He T) > = 0
I("I 0) 0
abc (He T)J = 0

3.3 Reduction to Spin Orbital Equations

The operator expression for the energy, Eq. 3.16, can be expanded to

Ecorr= 0o(fN+ WN) + T +T2 + T + + T3 + T T2+ T + ) 0 (3.18)

by replacing H, with its one- and two-body components (see Eq. 2.21) and expanding
the exponential. Since the Hamiltonian is in normal product form, the lead term in the
exponential expansion is (0 HN 0) = (01(H (0 HIO))I0) = 0. The remaining cluster
operators act on the reference to give an excited determinant and a product of amplitudes.
As a one-body operator, fN can create at most a singly excited state acting to the left on
the bra reference determinant, consequently nonzero contributions will result only from a
singly excited state on the right. Similarly, the two-body WN operator requires a doubly
excited state on the right. Thus, the equation simplifies to

Ecorr = ( fNTI + (T+ WN( 2 +T ) O). (3.19)
\ ~ ~ 2 \ -*//c

It should be observed that the CC energy includes contributions due only to T,

and T2 regardless of higher level cluster operators in T; higher clusters contribute only

indirectly via their effects on T1 and Ts.

Reduction of the amplitude equations (Eq. 3.17) to a more explicit operator form

proceeds in a similar fashion. For example, the T, equation can be expanded to

o=( (fn+W ) (l+Tl+T2+ T+ T3+TIT2+ T+--*) ) (3.20)

A nontrivial equation is obtained when the combined action of H, and the cluster

expansion produces a single excitation I1) on the right. The one-body operator fN

is capable of changing the excitation level by -1, 0, or +1, and W, by -2, -1, 0, +1, or

+2. The fN +1 term and the WN +2 term cannot form connected contractions with T, so

they can only contribute to the T1 and T2 equations, respectively, via the "1+" term in

the exponential expansion. Because of the nature of the Hamiltonian, mentioned above,

it should be evident that the equation for a particular cluster amplitude Tm will involve

contractions of fN with products of cluster operators resulting in m and m+l excitations,

and WN with cluster operator products yielding m-1 to m+2 excitations. Equations for

the first three cluster operators are thus

o=( fI + T 2(+T2+-T)

+W( T+T+ )
+W0 T I + T2 + WT{ l + TT + TIT2 + i+T+T
2! T)3!

0= iab.f T2+ T3 + TT21 + T WN + TI + T2 +- Ti + Ta + T1T2

1 4 ) \
+ T3 +T4+ TiTa + T + TIT2T2 + T 0 (3.21)
3! 2!-2 2! 4! C
K= .. T Ta+bcTIT3 +-1T2

+ WN T2 + T3 + T T2 + T4 + T T3 + TIT2+ TiTT 2

S2 T3 + T2 + I T23
+ Ts5 + TIT4 + T2T3 + T 3 T2
2! 2! 2! 3T C

where the operator determined by the equation has been underlined.

Conversion of the operator equations into the spin orbital equations useful for

numerical work is done by considering in turn each of the operator products in the

foregoing equations and using Wick's Theorem [15] to evaluate the possible contractions

among the strings of second quantized operators contained in the Hamiltonian and cluster

operators. This analysis can be carried out using a purely algebraic approach or a

diagrammatic approach [26], but rather than adding greatly to the length and complexity

of this overview we prefer to omit the detailed procedure here, and simply report the


Consider, for example, the first term of the CC energy as given by Eq. 3.19:

(OIfNTlj0) = fpqt(O{ aq}{aia}O)

= fpqtfpisqa(OlO) (3.22)

= fiat?

Similarly for the second term,

( WN (T2 + 7) = E (pqlrs) ( It+ i

S0 aa{4 asar}({aaaaiaj} + {4ai}{a aj}) O)
~1 t4( 1 ab

= (pq rs) 1 t + tit)Spi6jqi6s6ra(010)

= 2 (ij la b) t it+ j

which gives the spin orbital equations for the coupled cluster energy:

Ecorr = E fiat? + E (ij ab) (t + tit) (3.24)
i,a i,j,a,b

Spin orbital forms of the amplitude equations are obtained in the same way. Complete
spin orbital equations for the CCSDTQ model, the most complete CC model implemented
at the present time, may be found in the recent paper of Kucharski and Bartlett [27].


3.4 Numerical Solution of the Coupled Cluster Equations

The overall result, as can be seen in Eq. 3.21 is a set of coupled nonlinear equations

for the amplitudes tjb... based on the MO integrals fpq and (pqllrs). These equations

are usually solved by fixed-point iteration starting with a guess for ta obtained from

second-order perturbation theory, and it has been found that acceleration methods based

on conjugate direction techniques [28-32] are very useful in obtaining faster convergence,

and thus reducing the computational effort. Once the amplitudes are converged, Eq. 3.24

may be evaluated to obtain the CC energy. In many cases, the amplitudes and energy

will be used, along with the MO integrals, as inputs to further calculations to determine

the CC density matrix, gradients of the energy, or other results, for which it is often

useful to work in terms of the CC effective Hamiltonian.

3.5 The Coupled Cluster Effective Hamiltonian

When the coupled cluster calculation is carried out as the first step in a Fock-space

CC calculation, it is useful to form the elements of

HN= eTHNeT = (HNT) (3.25)

explicitly from the MO integrals and CC amplitudes because the FSCC equations can

be expressed in a compact and computationally efficient form using HN. Note that

this operator is non-Hermitian, so the bra-ket symmetry of the MO integrals (i.e.

(pqllrs) = (rsllpq)) is not carried over to the elements of HN:(pqllrs) 5 (rslpq).

Because HN contains at most four second quantized operators, HN will include

contractions of HN with as many as four T clusters. Like HN, the effective Hamiltonian

is described as having one- and two-body components. In fact, while HN terminates

at two-body terms, HN has up to n+2-body terms in an n electron system. Just as the

Hamiltonian is often written in terms of fN and WN, the components of the effective

Hamiltonian are usually denoted fN, WN, IIIN,...

The HN elements required for the Fock-space coupled cluster methods considered in

the next chapter include all of the one- and two-body terms and a few of the three-body

terms. They are obtained in the same manner as the spin orbital equations for the CC

energy and amplitudes, above, and the results are shown below.

The first stage is to isolate the components of the expansion

(H eT)= ( + +WN( I+Ti+T2+T + T3 + TT2 + T +))

from which contractions can form. At the operator level, the equations are closely related

to the CC amplitude equations (Eq. 3.21), but with two important differences. First, HN

is defined based on converged cluster amplitudes, so Eq. 3.21 must be satisfied prior to

forming the effective Hamiltonian elements. Secondly, like the Hamiltonian itself, the

effective Hamiltonian can have -n, ... 0, n excitation character for an n-body

component. So while ( HN 10) = 0 is identical to the CC amplitude equations given
above, the effective Hamiltonian itself contains numerous terms besides those required

to satisfy the equation. Thus the one-body component is

fN =f I I +T +T2+ T12
+WN TI ++ T2 T + T3 1 TT2 + 3

and the two- and three-body components of HN, are

WN =(fN {T2 + T3 + TT2 + W 1N{ + Ti + T2 + T + T + TiT2
+T- +T4+ TT + T + I TTr2 + TTP
3! 2! 2! 4 c

MIN =(IN( T+3 T+4 TT+ TT3 + T2

1 1 1 12
+ WN T2 + T3 + TT2 + T4 + TiT3 + 7TT2 + T T2 (3.29)

+TS5 + TIT4 + T2 T3 + Tf T3 + T1T + IT T2)

Reduction to spin orbital form follows the same procedure used for the CC equations.

In general, the n-body component of the effective Hamiltonian will have (n+1)2 different

terms according to the distinct labelings of the indices as occupied or virtual orbitals.

The one corresponding to the n-fold excitation will, of course, be zero and for higher n,

an increasing number of other terms will be zero. For example in IIIN it is impossible

to form terms with -3, -2, or -1 excitation character, as well as half of the possible

0 excitation terms (the all occupied and all virtual labelings). Thus while the number

and complexity of higher-body effective Hamiltonian terms does increase, it is not as

explosive as that seen in many-body perturbation theory (MBPT) [26]. The spin orbital
equations for HN are give in Tables 3.1-3.3.

One-body coupled cluster effective Hamiltonian elements in spin orbital form.
level of operator is given in the first column. See text for explanation of

fia + (imjlae)t
fii + fifth (imllje)t + (imllfe)tf + (imllfe)tft
fab fmbta + (amllbe)tm (nmllbe)te (nmllbe)tat,
fai frmit + faet m + fneti fmett (mallie)te

-}(nml ie)te (nmllie)tte + t (amllf e)te + (amllfe)tft'
+l(mnlef)t~ef + (mnllef)t?,eti (mnllef)tmt
-1(mnllef)t tit (mnllef)tltatf

Table 3.2. Two-body coupled cluster effective Hamiltonian elements in spin orbital form.
Excitation level of operator is given in the first column. See text for explanation of

-2 (ijflab) =
-1 (ijIlka) =
-1 (ai lbc) =
o (ijl kl)=

0 (abllcd) =

0 (ia ljb) =

(ijilka) + (ijilea)tk
(aillbc) (millbc)it
(ijilkl) + (-)PP(kl)(ijllke)ty + (ijllef)tP
+(-)PP(kl)(ijlle f)tet
(abllcd) (-)PP(ab)(amI cd)tb + (mn Icd)ta
+(-)PP(ab)(mnl cd)tatb
(ial\jb) (imrnjb)ta + (ia[[eb)t (imlleb)tf (im\\eb tt

Table 3.1.

-1 fai =
0 fij =
0 fab =
+1 fia

Table 3.2. (Continued) Two-body coupled cluster effective Hamiltonian element,
spin orbital form. Excitation level of operator is given in the first column. See text
explanation of notation.
+1 (ia|ljk) = (iaIIjk) + fiet, (imlljk)ta + (-)'P(jk)(iallje)t'
+(-)Pp(jk)(iml||je)ta (-)PP(jk)(im Ije)tekt
+(-)PP(jk)(iaIIef)ttf + (ialIef)tIf + mi( ief)\tefa
3)kk22m 2 JI mjk
+(mizllef) tftf l(iml fe)t ta
-(-)PP(jk)(imI fe)tL (-)PP(jk)(imI fe)tjt tm
+1 (ac(ab(l) abllci) ftci + (abI|ce)t (-)PP(ab)(amllci)t
+(mn lci)ta + (-)PP(ab)(mnl Ici)tt b
+(--)"P(ab)(am||Ice)t (-)PP(ab)(am Ice)titb
(mnI ec)t (mnl ec)tert + (-)PP(ab)(nm Ice)tattb
+(nm|lce)t mt + (-)pP(ab)(nm lce)ttatqt
+2 (abllij)= (ab lij) (-)PP(ij)f .. tab + (-)PP(ab)fbeta
+fmetj )P(i)fme (-)PP(ab)fmetiftb
-(-)PP(ab)(mbl ij)ta + (-)PP(ij)(ab lej)t + (mnllij)ta
+ (ab |ef)tf (-)PP(ij)(mb lie)ta, + (-)PP(ab)(mnjjij)tatb
+(-)PP(ij)(abl|ef)t/tf (-)PP(ij)(mbllie)t at
-(-)PP(ij)(nmllje)tb + (-)PP(ab)(bmr| fe)ltf
-(-)PP(ij)(mn|ei)tet + (-)PP(ab)(ma ef)te;t
-(-)PP(ij)(ab)(nmIIie)ttebi + (-)PP(ij)(ab)(amllfe)tfteb
+ l(-)PP(ij)(mnllej)t1tab (-)PP(ab)(mbllef) t/
Inn 1(-)PP(ab)(mbIIef)tit j

+(mnl|ei f)abefn + (mnlle f)taqbtl (-)P(ij)(mnlle f) ttafb
-(-)PP(ab)(mn ef)tn + (-)PP(ij)(mnllef)tift?
+ (mniIe f)t~t (-)PP(ab)(mnlle f)tetn

s in

Table 3.2. (Continued) Two-body coupled cluster effective Hamiltonian elements in
spin orbital form. Excitation level of operator is given in the first column. See text for
explanation of notation.

-(-)'P(ij) mn \ ef)t? (-)PP(ij)(ab)(mni efp t tfb
+ (-)PP(ij)(mn e f)t \tt + (-)PP(ab)(mn|lef)t ttt
-(-)'P(ij)(nm Ife)tt'tf (-)'P(ab)(nml fe)ttfato
+(-)PP(ij)(mn[ ef)t t tb

The summation convention (all repeated indices are summed) has been used to save

space and simplify presentation of these formulae. The permutation symbols, (-)PP(),

are used in the same fashion as Lee et al. [33] and Kucharski and Bartlett [27]. A

permutation of two indices, e.g. (-)PP(ij), denotes two terms. One is the identity

permutation, with signature p=+l, while the second is the term with indices i and j

interchanged. A permutation of three indices, e.g. (-)PP(ijk), represents six terms with

signature p=+l for cyclic permutations and p=-l for acyclic permutations. The symbol

(-)PP(i/jk) indicates that permutations i j and i <-+ k should be taken, but not

j <-- k; the signature is the same as for a binary permutation. Finally, a composite

symbol, such as (-)PP(i/jk,ab), indicates that the direct product of the two sets of

permutations (i/jk and ab) should be taken.

3.6 Truncation of the Cluster Operator

It is worth noting that we have made no approximations in deriving the equations

given above, and the energy, wavefunction, and effective Hamiltonian thus determined

would be exact (for the chosen atomic orbital basis). Unfortunately, the number of

Table 3.3. Three-body coupled cluster effective Hamiltonian elements needed for
Fock-space (1,0), (0,1), and (1,1) sectors in spin orbital form. (The full IIIN includes
an element at excitation levels +2, and +3 in addition to those shown here.) Excitation
level of operator is given in the first column. See text for explanation of notation.
0 (aijllklc) = (ijlec)ta
0 (abilljcd) =-(milcd)t!
+1 (ijajjklm) (-)PP(ml/k)(jiIlek)t" + (jilef)ti,
+(-)PP(k/lm)a(jilef)tae tf
+1 abclldei) = -(-)pp(b/ac)(mblIde)ta + (mnllde)tab
+(-)PP(b/ac)(mnl de)tia tb
+1 (iablljck) = -(mii|ec)tb + (-) P(jk)(mi\lec)tb t
+(-)P(ab)(mizlec)tq tb (-)PP(jk)(milljc)ti
+2 (iab jkl) = (-)PP(jk/1)(iml Ijk)t~t + (-)P(j/kl,ab)(ia| re)t
--prn/jl ab)P 1,mai) (a Ia eb e ab
-(-) P(k/jl, a)(mi| Ike)tmt (-)PP(jkl)(im lek)tjtfi
+(-) P(jk/l, ab)(im lfe)tte I (-)PP(jk/l)(im\ ef )teta
+(-)PP(j/kl,ab)(iaIef)tetf + f,
+(-)pp (j/kl)(irlmje)tea + (-)PP(ab)(ia Ief)te
+(-)PP(j/kl)(im Ife)tfiea 1(-)P(ab)(ZmIIef)taj
-i 2\( eftkjltm

+(-)PP(jk/l, ab)(i lfe)tate (-)PP(jk/l)(im llefteft
-(-)PP(k/jl,ab)(miI ef)}tta tb (-)PP(jkl)(im\ fe)tf t t
+_(mlIIefefab +e fab
1(rnllef)tmikL + (millef)tit ikl

amplitudes to be determined is O(N!) for N AO basis functions. This rapidly becomes too
expensive to compute using current and foreseeable computer facilities, so it is necessary
to make approximations. The obvious choice is to limit the cluster operators included
in T to lower levels of excitation because the number of amplitudes grows by a factor

of n(N-n) for each additional excitation. Because of the exponential nature of the CC

wave operator, some higher excitation effects are incorporated by products of lower-level

excitations (i.e. T, T2 is a triple excitation). In fact in most cases, such products dominate

the wavefunction at a given level of excitation [34].

The single excitations alone cannot have any contribution to the energy, so the first

reasonable approximation is T= T2, known as the CCD (coupled cluster-doubles) model.

With the doubles present, the singles can also contribute, so T=TI+T2, the CCSD

model, would be next on the hierarchy. Increasingly sophisticated methods are obtained

by adding the next level cluster operator: CCSDT, CCSDTQ, .. .

To obtain approximate versions of the equations given in this chapter, simply consider

all cluster operators not present in the desired model to be zero and eliminate them from

the equation. For example, CCSD equations would have Ts=T4=T5=0, and all terms

involving these operators would be eliminated. Observe, however, that products of TI

and T2 giving overall triple, quadruple, and higher excitations would still appear. This is

a consequence of the exponential expansion and is the reason why a given truncation of

coupled cluster theory recovers more of the correlation energy than the same truncation

of operator manifold in the linear CI expansion (Eq. 3.1).

In practice, the CCSD model has proven widely applicable, and the more expensive

CCSDT model, and approximations which involve perturbative inclusion of the triple

excitation effects, provide higher accuracy and have been shown capable of handling

some problems that would otherwise require multireference techniques [34]. The CCSDT

model, or perturbative approximations, are currently the practical limit for the majority


of interesting systems. Further, the computer code used in this work does not presently

include the full CCSDT model, so calculation of HN is limited to the CCSD level.


4.1 Introduction

In the preceding chapter, the CC equations were derived in a spin orbital framework.

This approach allows the use of wide range of reference functions (RHF, UHF, ROHF,

QRHF, etc.) thus allowing application to both closed and open shell systems with relative

ease. For open shell molecules, these spin orbital based CC methods may be characterized

as eigenfunctions of the S2 operator in a projected sense, (0 S2 cc,) = s(s + 1), for

ROHF and QRHF references, or not at all, for UHF references. Although they have been

very successful [34], failure to be a rigorous eigenfunction, S2 cc) = s(s + 1)| cc),

makes this approach less desirable to some workers in the field. The Fock space coupled

cluster method grew out of a search for an approach to open shell systems which would

be a rigorous eigenfunction of S' and be capable of handling the multiple reference

determinants needed to represent many open shell states.

There are two primary approaches to the treatment of multireference systems. The

Hilbert space approach, also referred to as state-specific, or valance-universal, treats

a manifold of states with a constant number of electrons, and each Hilbert space CC

calculation focuses on a particular state. The Fock space (FS) approach treats states with

different numbers of electrons (mathematically, a Fock space is a sum of Hilbert spaces

with different numbers of electrons), simultaneously calculating several different states.


The FSCC method is, at present, the more developed of the two, and offers the useful

feature, mentioned in Chapter 1, that it can provide data for more than one state in a

single calculation.

In this chapter, we will develop the FSCC equations in more detail, obtaining spin

orbital equations for IPs at the FSCCSDT level. The FSCC approach will be compared

with the familiar, and (reasonably) well understood single reference CC approach to

point up the inequivalence between FSCC and CC calculations at the same level of

truncation and note the increased importance of higher excitations (triples in particular)

in the FSCC approach. Finally, a number of approximations to the full inclusion of triple

excitations are proposed in order to include some triple excitation effects but reduce the

computational effort required.

4.2 Structure of FSCC

A multireference description of a correlated wavefunction means that the determi-

nantal expansion is partitioned into two components. The model space includes the

determinants desired as references, and the complementary space, which is orthogonal

to the model space, contains excited determinants resulting from each of the model

space functions. The same division holds trivially for single reference CC, with the one

reference determinant being the model space and all excited determinants being in the

complementary space. These two spaces may be represented by projection operators, P,

and Q, for the model and complementary spaces respectively, defined from the functions

they include:

P= E I0p)(OPI
Q = )(Cq = 1 P
which are orthogonal

PQ = QP = 0 (4.2)

and idempotent

PP = P, QQ = Q. (4.3)

For convenience, we choose a single determinant one-electron reference function,

10), to provide an unambiguous identification of occupied and virtual orbitals. The idea

of a particle-hole rank, denoted (m,n), is used to indicate the states related to 10) by

the addition of m electrons (particles in the otherwise empty virtual orbitals) and the

removal of n electrons (holes in the occupied orbitals). Thus singly ionized states are

in the (0,1) sector, electron-attached states are in the (1,0) sector, and the (1,1) sector

includes states obtained by exciting a single electron from 10). Higher particle-hole

ranks give multiply ionized or multiply attached states, excitations involving more than

one electron, and combinations such as "shake-ups" (combination of detachment of an

electron and excitation of another). The idea of the model space and its orthogonal

complement can also be specialized to a particular particle-hole rank, giving projection

operators p(m,n) and Q(m,), respectively.




f, active

orbitals T1

(holes) T1 inactive

Figure 4.1. Division of orbital space into active and inactive groups.

In order to define the model space, the orbital space of the reference determinant,

10), is divided into two classes. Active orbitals are those which change occupation in the

various model functions, while inactive orbitals maintain the same occupancy. Since the

reference function has already been divided into virtual and occupied orbitals, the result

is a four different categories of orbitals: active occupied, active virtual, inactive occupied,

and inactive virtual. These divisions are shown in Figure 4.1. Once the active orbitals

are chosen, the model space for each particle-hole rank is determined all appropriate

occupations of the active space are included. For example in Figure 4.1, the (0,1) sector

(one electron removed) model space includes six determinants and the (1,0) sector has

four because there are four ways to add a single electron to the active virtual space;

there are 12 spin-conserving single excitations within the active space, hence 12 model

functions in the (1,1) sector. It should be noted that this scheme provides for exactly one

model function, 10), in the (0,0) sector, and the calculation in this sector is exactly the

single reference CC method discussed in the preceding chapter.


Model spaces can be divided into two types, complete or incomplete [35]. A complete

model space is one that contains all possible occupations of the active orbitals, while in

an incomplete model space (IMS), some occupations of the active orbitals belong to the

Q space. The pure ionization and pure electron addition sectors, (0,n) and (m,0) sectors

are always complete, but any mixed sector is, in general, incomplete. For example, the

(1,1) sector has all single excitations within the model space in P, but places all double,

triple, and higher excitations into the Q space. The (1,1) sector is, however, part of a

special class of IMS known as the quasi-complete model space [35,36]. In this case, the

active orbitals are divided into disjoint sets, each with a fixed number of electrons. If

each of these subsets is complete within itself, the entire model space is quasi-complete.

The FSCC wavefunction for a state is obtained from a wave operator, Q, operating

on a linear combination of model space functions:

I'FSCC) = 'lk). (4.4)


k Ckii). (4.5)

The wave operator used in FSCC is an exponential of cluster operators, just as in the

single reference theory, but with a normal product used to avoid contractions within the

cluster operators itself (this is not necessary in the single reference formalism since it is

impossible for those cluster operators to contract among themselves),

S= {eS. (4.6)


The cluster operator itself, however, has a somewhat different structure from the single

reference case:

S = S(mn) = T(kl) (4.7)
k=0 ...m
where (m, n) is the highest sector relevant to the problem, and the T(k'l) for a given sector

includes single, double, and higher excitation cluster operators for that sector:

T(k,') = T(k,') + T (kk) + T ..l) (4.8)

Such a wave operator is valance-universal because it applies to all sectors from (0,0) to

(m,n). As mentioned above, the (0,0) sector is a CC calculation with 10) as the reference,

so T(o,0) is nothing more than the familiar CC cluster operator described in the previous

chapter. Cluster operators for other sectors are different from the (0,0) ones because they

apply to spaces with different numbers of electrons, or at least electrons in different places

than in 10) (i.e. active virtual may contain electrons and active occupied may receive

electrons). In the (0,1) sector, for example, one of the virtual orbital labels on a cluster

amplitude becomes an active occupied, corresponding to the fact that in this sector, an

electron has been removed from that group, making it possible for the empty orbital

to receive an electron. The differences by sector in the cluster operators is depicted in

graphical form by Rittby and Bartlett [37] for a number of sectors, and some implications

of these differences are discussed in Section 4.3, below.

To derive the FSCC equations themselves, we begin from the Schridinger equation

for the FSCC wavefunction,


HS2IkIk) = Ekf~lx~k)


This equation can be rewritten in terms of the model space functions by defining a

matrix e' having the E, as eigenvalues and C as eigenvectors, and using the model space

projector, P, to effectively encompass the entire model space,

HfQP = fl'P. (4.10)

Using the fact that the (0,0) sector cluster operator commutes with the others [38], we

can decouple the (0,0) sector calculation from the rest of the FSCC calculation. First,

we extract T(O,0) from the wave operator,

S= eToofi = eo'o) { eS(m.n), (4.11)


g(m,n) = T(k,), (4.12)
and multiply both sides of the Schr6dinger equation by exp (-T(,0)) to obtain

e-T~oa) HNeoo) iP = i'P (4.13)

Recalling that the coupled cluster effective Hamiltonian, derived in the preceding chapter,

HN = e- oo He 0 (0HI0)
= e-o'O)He T(O Ecc,o)
where Ecc,0o) is the CC energy of the reference function 10), and defining a new matrix,

e = e' Ecc,j), which, when diagonalized, represents energy differences between Ecc,jo)

and the final state, we can write the FSCC Schridinger equation in the simplified form:

fHNIP = eiP.


Projecting Eq. 4.15 on the left with PO-1, we can isolate e as the FSCC effective


PHN,effP = Pf- HNP = PeP, (4.16)

which, upon diagonalization, gives the state-to-state energy differences we seek. Note that

this result could have been obtained by requiring intermediate normalization, PQP = 1,

but Mukherjee [39] has shown that for the general case of an incomplete model space,

intermediate normalization is incompatible with size extensivity. Although this will not

cause problems for the (0,1) sector, our ultimate focus, we prefer to derive generally

applicable equations as far as possible. To determine the cluster operators, we project

Eq. 4.15 on the left with the complementary space determinants, giving

QHNQP QQPHN,effP = 0, (4.17)

which is often referred to as the Bloch equation. Mukherjee and coworkers have shown

that for complete and quasi-complete model spaces, HN,eff in Eq. 4.17 is connected and

the energies obtained by diagonalizing it are size extensive [39-44].

Haque and Mukherjee [38,45] have also shown that there is a subsystem embedding

condition (SEC) in the FSCC equations which allows Eq. 4.17 to be solved separately

for each sector in a hierarchical fashion. As already observed, the (0,0) sector reference

problem decouples from the rest of the Fock-space problem. The (1,0) and (0,1) sectors

can be solved with only the (0,0) results they do not depend on any higher sectors.

Given the (0,0), (1,0), and (0,1) sector amplitudes and HNeff, the (2,0), (1,1), and (0,2)


sectors can be solved once again, no information from higher sectors is required.

In this fashion, the solution of the desired (m,n) sector is built up from lower sectors,

rather than requiring the solution to a set of coupled equations for all sectors. The SEC

is also very useful in the development of computer programs for FSCC because the

implementation of a new sector should not require any changes to the lower sectors.

4.3 Comparison with Single Reference CC Theory

Before going on to derive the more specific equations for Fock-space methods, it

is worth taking a look at how the FSCC approach differs from the more familiar single

reference formulation.

As has been noted, the Fock space coupled cluster amplitudes differ from (0,0)

sector amplitudes in that they involve excitations between the occupied or virtual orbitals

and orbitals of the same type in the active space. In the reference function, 10), such

amplitudes would be zero because active occupied orbitals would be fully occupied, so

no excitation could deposit an electron there, and active virtual would be completely

unoccupied, providing no electron to excite to another virtual. In the model space,

however, there will be determinants with holes in the active occupied space and particles

in the active virtual space, so such amplitudes may be nonzero. This fact changes the

character of FSCC single, double, etc., excitation amplitudes from their familiar (0,0)

sector interpretations and means that a given truncation of FSCC will not provide the

same results for a given state as the same truncation of the (0,0) sector equations applied

to the model space function.

10) i 0,1)) ( 0,01)) 00'1 ) i 0,1)

Figure 4.2. Prototypical reference functions and model space determinants. Active
orbitals are shaded.

Figure 4.3. Determinants appearing both in the single excitation space of CCSD calcu-
lations on the four model functions and in the FSCCSD singles space.

To make this idea more concrete, consider a FSCCSD IP calculation

compared to CCSD calculations based on the individual determinants in the

model space. Using an RHF reference function with three occupied orbitals

and one virtual orbital, if the two highest occupied orbitals are made ac-

tive, we have four determinants in the model space, as shown in Figure 4.2.

A CCSD calculation on any of the model space functions has seven determinants in

the single excitation space, and the union of the four single excitation spaces contains

20 distinct determinants. FSCCSD, on the other hand, restricts single excitations to

be from an inactive occupied orbital to an active occupied, which allows only four

amplitudes and yields two distinct determinants, shown in Figure 4.3. The remaining 18

TI I T I t T I
t24 34 t24 t24 24
23 32 21 2
t34 t34 t34
31 32 31
Figure 4.4. Examples of determinants arising from spectator double excitations in FSCC,
corresponding to single excitations in the CC model. Shown are all determinants
generated from the 10(o1)) model function with the relevant amplitudes from T(o')
indicated below each one. Orbitals are numbered from lowest energy (bottom) to highest
(top); beta spin electrons are indicated by a bar over the orbital index.

CCSD determinants are generated in the Fock space calculation by a class of excitations

referred to as spectator doubles. FS double amplitudes include one excitation from an

occupied orbital to an active occupied. There is no restriction on the initial orbital, so

it may be active as well in fact it may be the same as the final orbital. This null

excitation, which is necessary to distinguish between excitations connecting different

model functions to the same final determinant [46,45], makes spectator doubles effectively

single excitations in the usual coupled cluster sense. Figure 4.4 shows the spectator

doubles arising from model function 1(o1')), along with the corresponding amplitudes.

When considered at truncation to only double excitations (CCD and FSCCD), the

Fock space approach is going to include some (but not all) effects of single excitations -

a potential advantage. When higher cluster operators are included, however, FSCC is

at an unequivocal disadvantage. Compare, for the same model problem, the double ex-

citations obtained from the CC and FSCC approaches. FSCCSD generates 22 "true"

tT T T 1 T 1 T
T 1T 1 T T
t244 t344 t244 t244 3344 t344
232 322 212 231 321 311
t344 t244 t244 t244
312 231 221 211
Figure 4.5. Examples of determinants arising from spectator triple excitations in FSCC,
corresponding to double excitations in the CC model. Shown are all determinants
generated from the 4(0,1)) model function with the relevant amplitudes from T(o')
indicated below each one. Orbitals are numbered from lowest energy (bottom) to highest
(top); beta spin electrons are indicated by a bar over the orbital index.

doubly excited determinants, while CCSD calculations on the individual model func-

tions yield 18 additional determinants. The "extra" CCSD determinants are created

in FSCC by spectator triple excitations, a few of which are displayed in Figure 4.5.

From these examples, it should be evident that for a FSCC IP calculation to match

the expansion space of a single reference CC of a given truncation on one of the model

functions, the FSCC calculation must include spectator contributions from the next higher

cluster operator. The same would be true of EA calculations, and in higher sectors, more

indices of a given amplitude terminate in the active space (p=m+n for the (m,n) sector),

so it would be necessary to include spectators from the ph higher cluster operator to

obtain an expansion space equivalent to that of the single reference CC.

4.4 Basic Equations

From this point, we will focus on the (0,1) sector, which describes one-electron
ionizations with respect to the reference state and are of particular interest in this work.

From Eq. 4.17, the detailed operator equations which determine each cluster operator
are easily determined by expanding the exponential wave operator, as in the previous
chapter. Note that in FSCC, products of cluster operators in the exponential expansion
contribute to other sectors, making the (0,1) equations a good deal simpler than their
CC counterparts.

P(o,) fN,effP(O'l) = P(o1) f N(1 + T0',1) + T~1))
+ WNT (TO,) + T('1) + IIINT(,'1)}P(o,1)

o = Q( ') (1 + T,(O) + T('1)) + WN (T(O1) T(0,1))
+I7INTo01,) T(o'l)f1) cP(ol,)

0 = QO')(fN (TO'1) T(O"'))+ N ( + T 1 T ) T() + T3(,))
+IIN (T(0,o1) + T(O',)) T(Ol)f p(o,1)
(llNU2 +3 ) 2 NJVeff 1,

(0,1) 0,1) 0,1)
O = Io3 l)fNT3l) + WN (T2~1) )+ T(0'1))

+ IIN (i + T(l) + T( ) + T3(01)) (4.21)

-T(o,1) fr(O,) p(o,)
-3 N,effJc

4.5 Reduction to Spin Orbital Equations

Equations 4.18-4.21 are transformed to spin orbital form following the same proce-
dures used in the previous chapter for the CC amplitude equations and HN. The resulting
equations for the FSCCSDT method for ionization potentials are given in Eq. 4.22-4.25,
below. The notation is the same as that used in the previous chapter with a few additions.
The Fock-space cluster amplitudes are called r to avoid any possibility of confusion with
the (0,0) sector amplitudes, t. A single dot on an orbital index, ri, indicates restriction to
inactive orbitals, and a double dot, ih, indicates restriction to active orbitals (this choice
of notation is meant to remind readers familiar with the diagrammatic notation of the
"circled arrow" used to denote inactive lines, and the double arrow for active lines).

(HNeff)k = f r + fme mn \e /r + mn\ ef (4.22)

0= fj i+f- mnke)rm + mnlef)r +{(HN,eff)iH (4.23)
fkrm + k 4h+ fm' m 2 4kmn k

0 = (iaIki) + fea (-1)P(ki)fmirT (rhaIki)rm

+ 2(-l),P(ki)(mn|Lki)rC 7 (ma rke)-e i 2(manfikie)re.

+ fmer + -(-1)PP(ki)(mnf]ke)ri + I(mnalef)rkmi + 4 (HN,eff)i

0 = ilabI\kij) (-)PP(k/ij)(mb ij))rka + (-)'P(ki/j)(ablej})T,

+ (-)PP(ab)(man[Ikij)rTn (-)PP(ki/j)(mibi)kie)rEn (rhabIkij)r}

+ (-) P(ab)fberkij (-)PP(ki/j)fm,7jTki + (-) P(ki/j)(mnIrki)r,

(-)'P(k/ij,ab)(ma\lke)rL +ki(HN,eff)l

4.6 Approximations to the Full FSCCSDT Model

In the single reference CC world, the full CCSDT model has seen limited use for

a number of reasons

* Complexity of the method the triples equations have many more terms than single

and double excitations, thus requiring substantially more effort to code them into a

computer program.
* Primary and secondary memory (i.e. physical memory and disk space) requirements

for storage and manipulation of the cluster amplitudes, which are significantly more


numerous than the double excitation amplitudes a factor of n(N-n) in the language

of section 3.6.

Time constraints, since many more amplitudes must be determined, the time required

increases in a similar fashion.

As a result, a number of approximations to the full CCSDT method have been proposed

and have proven quite effective at incorporating important effects of triple excitations

while at the same time minimizing the above concerns [33,47-52]. These concerns

might also be raised with respect to the full FSCCSDT method, so in the same spirit,

we propose several approximations of FSCCSDT as a compromise between the cost and

performance of the full model.

There are two different approaches to approximating FSCC triples which may be

considered. One, in analogy to the approximate triples treatments of single reference CC,

involves truncations and perturbative approximations to the triples equations. The second

method is to restrict the triples amplitudes which are evaluated to something less than

the full set. This approach, which has not been employed in CC theory, is suggested by

the earlier observation that spectator triples must be included to, in some sense, duplicate

the expansion space of a single reference CCSD calculation. The two approaches are not

exactly orthogonal, but can be considered, by and large, independently.

In single reference CC theory, there are no self-evident subsets of the triples

amplitudes, so this approach has not been tried (or at least not reported in the literature).

In FSCC, however, there are two obvious subsets which might be considered. The

smallest set would be the spectator amplitudes, as suggested by the discussion in the


previous section. A larger set would be to include all active-active amplitudes. This set,

which includes the spectators, allows excitations from any active orbital into any other

(spectators involve an active orbital scattering into itself, and the full set of triples allows

scattering from inactive orbitals to occupied orbitals as well as active to active).

As we saw in the previous section, it is necessary to introduce spectator triples in the

Fock space approach in order to match the expansion space of a CCSD treatment of each

of the model functions in turn. On the other hand, since FSCC describes all states of

interest with a single set of amplitudes (with diagonalization of the effective Hamiltonian

distinguishing among individual states), it can be argued that it is impossible for it to

be equivalent to a separate CC calculation for each model function, and the differences

in the expansion space between the two approaches should not be a concern. Indeed, a

number of calculations in the literature [53,54,43,55, for example] support the idea that

FSCCSD is well-balanced and the results obtained are numerically on a par with other

methods. The CCSD method has been widely applied, and its behavior is reasonably well

understood certainly there is far more experience with CCSD than with any FSCC

method and this, more than any special quality of CCSD underlies the desire for a

FSCC method that can mimic it. Only a critical examination of the proposed spectator

and active-active approximations for a wide range of systems can determine if they will

improve the results over FSCCSD or destroy the balance of the FSCC method.

The second approach to reducing the cost of a FSCCSDT calculation is to approximate

or ignore the most expensive terms to evaluate. In this sense, the simplest possible

approximation would be to evaluate the triple excitation contribution perturbatively from


converged FSCCSD amplitudes. The single reference analog to this approach is the

CCSD+T(CCSD) method, in which the triple excitation contribution is correct through

fourth order the first order at which it appears.

For the FSCC perturbation expansion, we follow the choice of Pal, Rittby, and Bartlett

[56] which was also implicit in the work of Haque and Kaldor [57], of counting orders

in terms of H, the reference state Hamiltonian. Alternatives, such as using H would lead

to inconsistencies between the (0,0) sector definitions and those of higher sectors. For

example, f and W are first order in H, while the higher components, III,... are at least

second order. Counting orders in H would place all of these terms on the same level,

making a perturbation expansion of Eq. 4.21 useless.

Although in both CC and FSCC the triple excitation amplitudes first appear in second

order of perturbation theory, in the CC formalism excitations beyond singles and doubles

contribute only indirectly to the energy via their contributions to the singles and doubles

amplitudes, and thus triples first appear in the energy in fourth order. In FSCC, all clusters

contribute directly to HNeff and thus the energy, so the triples first manifest themselves

in third order. Inclusion of the third order triples effects based on converged FSCCSD

amplitudes was first proposed by Haque and Kaldor [57], and called CCSD+T. Later, Pal,

Rittby, and Bartlett [56] discuss the third order model using the name MRCCSD+T(3)

and propose several other approximations. The MRCCSD+T*(3) model attempts to

include some higher order terms by replacing the f and W required in T(3) with their

counterparts from the effective Hamiltonian, f and W. The nature of this method, which

will be referred to as FSCCSD+T(3) in this work, can be seen from the spin orbital


FSCCSDT equations given above. The terms which contribute to the T(0'1) amplitudes

in second order are those on the first row of Eq. 4.25, with appropriate restrictions on

the HN elements: the WN elements of the second and third terms are simply the bare

two-electron integrals in first order, and only two terms of (lab lkij) appear at this level

(the first two terms of (iabfljkl) in Table 3.3). The second order triples amplitudes in

turn contribute to HN,eff (and thus the energy) in third order thought the final term of

HN,eff in Eq. 4.22.

The extension of this approach to higher orders of perturbation theory is obvious. A

fourth-order method, FSCCSD+T(4), would be obtained by the following steps:

1. Second-order T o1) amplitudes, as evaluated above, give third-order contributions

when inserted in the T,(o') and T(O'1) amplitude equations. This correction to the

T"'1) and To'1) amplitudes can be directly contracted into HN,eff, giving a fourth

order correction to the energy.

2. The third-order T(o1) must be evaluated. This requires the fourth, fifth and sixth

terms from Eq. 4.25, and additional terms from (labllkij) (terms 3-10 of (iab\jjkl)

in Table 3.3). Note that even if the (0,0) reference calculation does not include triple

excitations, approximate T 0'0) amplitudes can be formed at second order in pertur-

bation theory, and therefore must be included in the general case. Computational

experiments will be required to determine the effect of approximate T(0'0) in N-y.

Of course such a perturbation expansion could go on for a very long time, but

in practice, better results are likely to be obtained (based on SRCC experience) by

iterating an approximate triples equation as part of the FSCCSD calculation. The

simplest approximation for which this makes any sense is T(4), described in the preceding

paragraph it is the first perturbative method for which T,("') T(o'1) and T3(0') -

T0",) contributions occur. The procedure is simple: In each iteration of the FSCC

calculation, after the single and double excitation amplitudes have been evaluated, the

T(4) procedure is used to calculate the contribution of triples excitations to HN,eff and

the singles and doubles amplitudes. Instead of this being the end of the calculation,

the new ( 1) and T'1) amplitudes, and HN,eff are fed back into the next cycle of

iteration. Note that there are no T(o'1) T3(o') contributions in this method, so the

triple excitation amplitudes do not need to be kept from one iteration to the next, and

also avoids the most expensive part of the T(''1) equation to evaluate. This approach is the

same as that taken by Urban et al. [48] in approximating single reference CCSDT. In the

hopes of maintaining some consistency in nomenclature between SRCC and Fock-space

approaches, this method should be called FSCCSDT-1.

In the single reference case, the idea has been extended in a series of additional

iterative approximations, CCSDT-2 [48,50], CCSDT-3 [48,50], and CCSDT-4 [51],

which include more and more of the terms due to triples, but continuing to exclude

T3(') To'O), thus avoiding the most expensive terms both computationally and for

storage space. In the FSCC case, the T(o'1) equation has a simpler structure, and only

T3o'1) T(o'1) terms remain after the FSCCSDT-1 approximation.

Over the years, experience with various SRCC approximate triples methods has

led to refinements in several of them. One of the main observations is that in-

cluding selected higher order terms in a perturbative or iterative method can pro-

duce better results. For example, Watts, Gauss, and Bartlett [58] have recently intro-

duced the CCSD[T] method, which is a modification of CCSD+T(CCSD), the sim-

plest SRCC perturbative triples method. They differ in the inclusion of a fifth-

order term and in the requirement that non-Hartree-Fock (i.e. QRHF or ROHF) ref-

erences be treated properly. This requires that terms involving off-diagonal f ele-

ments be considered since they are zero for a Hartree-Fock reference, but nonzero

in the non-HF case. Transforming the reference functions to semicanonical orbitals

[59] has the effect of diagonalizing the occupied-occupied and virtual-virtual blocks

of the Fock matrix, so that terms involving fij or fab can be avoided. This leaves

occupied-virtual Fock matrix elements, fai to be accounted for, which involves shift-

ing these terms down one order of perturbation theory, from first to zeroth or-


In the Fock-space coupled cluster framework, however, the situation is somewhat

different. The transformation to semicanonical orbitals diagonalizes blocks of f, not fN.

Even with such a transformation, the leading terms in fij and fab are first order in the non-

HF case. This would result in FT~'1) T(o' 1) contributions, one of the most expensive

terms in FSCCSDT, beginning with the third-order amplitudes (FSCCSD+T(4)). Thus

as far as approximations to FSCCSDT go, such a requirement would be reasonable only

for the most basic third-order method, FSCCSD+T(3). This would add the FT(0'0) term

of (labllkij) (see Table 3.3) and the FT(o'O) terms of (mbhij) and (abllej) (see Table

3.2), and would require evaluation of the FT(o') T(o1 ")-H H term of Eq. 4.24

since fme is now zeroth order.


There is one further matter regarding approximation to the FSCCSDT equations.

Mattie [60] has proven that for the (0,1) sector of Fock-space, the values of the roots

obtained by diagonalizing HN,eff are independent of the active space used for the

calculation. In other words, if two orbitals, a and b are taken as active, the resulting

ionization potentials would be identical to those obtained from separate calculations with

either a or b alone active. This is a very useful result because it means that we do not

have to worry about choosing the "right" active space for a given calculation in order

to get good results. However the proof of this invariance rests on the FSCC amplitudes

and HN,eff satisfying the Fock-space Bloch equation 4.17. None of the approximate

FSCCSDT methods described above satisfy Eq. 4.17, so the invariance is lost. It remains

to be seen how large this effect will be, or if there is an alternative route to the invariance

condition that will work with approximate solutions.

As a start towards understanding the FSCC triple excitation correction as well as we

now understand the single reference triples correction, the remainder of this work will

focus on the third-order triples approximations, FSCCSD+T(3) and FSCCSD+T*(3). The

implementation of these methods for both closed and open-shell reference functions is

described in some detail in the Appendix, including the final, spin-integrated equations

from which the computer code was written. In the following chapters, these methods are

applied to a variety of molecules in order to assess their performance.


5.1 Introduction

Before applying a new theoretical method to otherwise unknown molecular systems, it

is important to get a feel for its behavior strengths, weaknesses, and overall quality -

by comparing the results obtained with new method to other theoretical and experimental

data for similar types of problems. This "calibration" of a method is a continuous and

dynamic process each time the method is used for a new problem, more information

about its performance is gained, thus modifying the assessment of its accuracy and


Choosing the appropriate standard of comparison for a new method is in most cases

not a simple matter. Experimental results are an obvious choice, but there some potential

problems. Experimental corresponding to readily computed results are not always

available or may not be known accurately. For example, FSCC most readily provides

vertical ionization potentials, but only adiabatic may be available from experiment. For

potential energy surfaces geometries, relative energies, and the like there is usually

very limited experimental data. Structures may be known for minima on the PES, but

experiments which can probe transition states are quite new, and have not yet produced

a sizeable body of data on different systems. The vibrational spectra of a given species


may be only partially known, and in many experiments, assignment of spectral lines to

particular species can be a very hard problem. It is also possible for a calculation to obtain

the right (experimental) answer but for the wrong reasons a fortuitous cancellation of

errors. If such problems are given due attention, experimental results can be a very

useful measure the performance of a theoretical method, but it is worth considering other

standards as well.

Comparison of a new method with experimental results often introduces more vari-

ables than one would like. For example, the choice of basis set may not be adequate to

obtain high-quality theoretical results. Full CI provides the exact answer for a given basis

set and Hamiltonian, and therefore offers an excellent standard of comparison for new

methods. Unfortunately, FCI calculations are very expensive, so the results available for

comparison are somewhat limited. Another option is to compare with methods that, from

a well-defined theoretical standpoint, are better than the method in question. An example

of this would be comparing methods which approximates triple excitation effects to the

full triples model such as CCSD(T) vs. CCSDT, or FSCCSD+T(3) vs. FSCCSDT. In

this work, however, as is often the case, the approximate method has been implemented

before the full method because of the additional complexity of the full model, forcing

such comparisons to wait until the more complete method is implemented. A final al-

ternative is comparison of the new method with one for which the behavior is relatively

well understood, but which cannot be clearly designated as a "better" method. Although

these comparisons are the least definite, they are, of necessity, the most common in the

calibration of any new method.


This chapter is a beginning at the calibration of the FSCCSD+T(3) and

FSCCSD+T*(3), models. The molecules chosen for these calculations are drawn largely

from applications of the FSCCSD method already in the literature, and in this work

extended to include triple excitation contributions. The intent is to get an overview of

the performance of the new methods for a variety of systems, most of which have also

been the subject of single reference CC calculations, thus providing a useful point of


5.2 Effect of Changes in the Active Space

The triples approximations, as was noted in the previous chapter, do not have

the same invariance with respect to the choice of active space as occurs for com-

plete methods. In order to simplify later comparisons, we first examine the effect of

the change in active space on the third-order triples models. The "error bars" ob-

tained in this section may then be kept in mind later on, where a single choice of

active space will be used in comparisons of the new FSCC methods with other re-


Because HN,eff is block diagonal due to molecular symmetry, the active space within

a particular symmetry block can be changed without affecting the others. If, in a series

of calculations, the active space for a particular irrep remains unchanged, the IPs of that

symmetry will not change. This is the case regardless of changes to the active space

in other irreps. Conversely, changing the active space within a given irrep will result

in slightly different IPs.


In order to assess the magnitude of this effect, calculations have been carried out on

a number of molecules using different active spaces. The molecules, listed in Table 5.1,

include a number for which the IPs themselves are studied below, and several additional

systems with fairly low symmetry actual or computational symmetry. The idea behind

these additional molecules was to see if having more occupied orbitals in each irrep (and

therefore a larger set of possible active spaces) would result in larger variations in the

calculated IPs. FSCCSD+triples calculations were carried out on each molecule using a

DZP basis at the DZP/SCF optimized geometry for a set of active spaces consisting of

1. the largest active space for which the underlying FSCCSD could be converged

(referred to as the "reference" active space),

2. single orbital active spaces for each of the orbitals in the reference active space, and

3. a series of active spaces obtained by starting with the lowest IP and successively

adding the next highest IP.

This choice provides a slightly larger set of data on the variation of lower IPs, but those are

usually the roots of greatest interest anyway. Analysis of these calculations involved first

reducing the data to obtain a set of unique active spaces within each symmetry block,

then computing for each IP, the mean, range, and mean difference from the available

data. The mean was chosen as the reference point because there is no theoretical basis

for preferring one active space over another as the "correct" one, except perhaps for

having all occupied orbitals active, which proved impossible to converge in most cases.

Table 5.1 presents the largest range and mean difference observed among the various IPs

Table 5.1. Effect of changes in active space in the T(3) and T*(3) models. for a variety
of molecules. Ranges and mean differences of IPs are given in mH.

Point Active FSCCSD+T(3) FSCCSD+T*(3)
Point Active
Molecule nC) Mean Mean
group) orbitalsb) Ranged) Ranged)
Diff.e) Diff.e)
C3H60 Cs 3 3 0.023 0.010 0.023 0.010
CHg) C2v 3 3 0.078 0.035 0.065 0.029
CH2NHg) Cs 4 4 0.744 0.220 0.545 0.165
CH2PHe) Cs 4 2/4 0.103 0.052 0.074 0.028
CH3CIf) C (C3v) 4 4 0.008 0.004 0.006 0.003
CH3CNf Cs (C3v) 4 4 0.354 0.133 0.284 0.107
CH3NCf) Cs (C3v) 4 4 0.044 0.022 0.044 0.022
CH3NH20 Cs 4 4 0.226 0.083 0.175 0.066
H2COg) C2v 2 2 1.118 0.559 0.931 0.466
HOC10 C, 3 2 1.844 0.921 1.581 0.790
HOFO Cs 3 2 3.467 1.734 2.757 1.379
N2g) D2h (Dooh) 3 3 1.297 0.576 0.861 0.383
a) Computational point group (full point group).
b) Largest active space used in the irrep containing the root showing the largest variation.
c) Number of different active spaces sampled for the root showing the largest variation.
Same for T(3)/T*(3) unless two values are given.
d) Largest range (difference between largest and smallest values) for a single root. (in
e) Mean difference (I =1 IP, IP|) for the root showing the largest variation (in mH).
0 Calculated with a DZP basis at the DZP/SCF optimized geometry.
g) Calculated using the geometry and basis used in the ionization potential calculations
in the next section.

for each molecule, along with an indication of the number of distinct active spaces in

which that IP was computed.


The largest variation is observed for HOF, with HOC1, N2, and H2CO close behind.

For a number of other species, however, the variation is quite small. Interestingly, the

largest variations are seen among molecules with higher symmetry and smaller active

spaces rather than large ones. This obviously suggests that the number of occupied or-

bitals is not the most important determinant of the variation, however examination of the

FSCCSD and reference CCSD results provided no clues to distinguish molecules like

HOF or H2CO from ones with a small variation, such as CH3Cl.

Comparison of the T(3) and T*(3) results shows that the range and mean difference

for T*(3), which includes some higher-order terms, are at least as good as, and generally

better than, those for T(3). Although there are a few cases where going to the T*(3)

approximation changes which of the IPs has the largest variation, generally it is the same

ionization for both.

In order to verify that the variation is not due to an inadequacy of the DZP basis sets

used for most results in Table 5.1, a second set of calculations was carried out on HOF

using Dunning's correlation consistent PVQZ basis [61], which is 5s4p3d2flg on O and F,

and 4s3p2dlfon H for a total of 125 basis functions (using spherical harmonic functions),

compared to 4s2pld and 2slp for the DZP basis, totalling 37 basis functions (using

cartesian functions). Table 5.2 presents the mean ionization potential and its variation

with active space for the lowest five IPs of HOF calculated using both basis sets. There are

changes in the range of variation from one basis to the other some roots have a larger

range with the PVQZ basis, others have a smaller range. In general, however, the they are

roughly the same, suggesting that the variation is largely independent of the basis set.

Table 5.2. Effect of changes in the active space on the calculated ionization potentials
(in mH) of HOF for two basis sets.

DZP basis PVQZ basis

Sym. n Mean Mean Mean Mean
Range Range
IP Diff. IP Diff.
a' 3 545.7 1.158 0.608 559.0 1.176 0.480
a' 3 672.8 3.186 1.064 686.0 3.529 1.281
a' 2 775.7 3.467 1.734 765.3 3.714 1.856
a" 2 466.5 0.171 0.085 482.0 0.118 0.059
a" 2 657.6 0.138 0.069 674.2 0.093 0.047
a' 3 538.0 1.391 0.547 550.5 1.076 0.431
a' 3 665.9 2.465 0.853 678.1 2.872 1.033
a' 2 749.8 2.757 1.379 758.4 3.080 1.540
a" 2 459.7 0.073 0.037 474.3 0.033 0.017
a" 2 646.9 0.044 0.022 662.5 0.012 0.006

Given the data above, a reasonable estimate of the error bars due to changes in the

active space is 2 mH for both methods. (It must be pointed out that there are certainly

to be systems with larger variations, but more experience will be required to understand

the origin of the large variations and therefore be able to identify the worst offenders.)

These correspond to about 0.05 eV or 1.3 kcal mol-1. Although this may seem like

a fairly large uncertainty, particularly for comparing energies at different points on a

potential energy surface, there are several important factors worth noting. First, this

is a conservative estimate for entirely unknown molecules. Many of the molecules in

Table 5.1 are more than an order of magnitude better. Also, the table presents only


the largest variation many others were significantly smaller. Secondly, it must be

remembered that for molecules with symmetry, many of the IPs exhibit no variation

because there is only one possible choice of active space in that symmetry. Finally,

it is a fairly simple matter, and not computationally expensive, to obtain a feel for the

variation of any molecule of interest. Only one (0,0) sector calculation is required to

try out several FSCC active spaces, and the FSCC calculations themselves are very

inexpensive compared to the (0,0) reference.

As a result, while the uncertainties associated changing the active space in T(3) and

T*(3) calculations for the general case might be considered large for some applications

much tighter limits may often be obtained with minimal effort. For the remainder of this

chapter, the active spaces used will correspond directly to those results which are reported.

Uncertainties due to the active space will be introduced into discussions as appropriate.

5.3 Ionization Potentials

Although it can be used in many other ways, the "natural" result of a (0,1) sector

FSCC calculation is the energy of the vertical ionization process extracting an electron

from a molecule without allowing the geometry to relax to that of the resulting ion. Cal-

culation of ionization potentials is, thus far, the predominant application of the FSCCSD

method and such calculations will provide the bulk of the data for characterization of the

approximate triples methods. These results will be compared with SRCC calculations,

experimental results, and a benchmark full CI calculation for the methylene molecule by

Bauschlicher and Taylor [62].

Table 5.3. Vertical ionization potentials of 'A, and 3B1 CH2 (in eV) calculated with a
DZP basis.

Method A 3B
12A, 22A1 2B2 2A1 2B1
FSCCSD 10.21 22.54 14.92 10.14 11.30
FSCCSD+T(3) 10.32 22.74 15.00 10.26 11.43
FSCCSD+T*(3) 10.28 22.68 14.96 10.23 11.39

Full CI 10.26 22.14 14.85 10.16 11.31
QRHF-CCSD 10.20 22.25 14.84 10.14 11.30
QRHF-CCSD[T] 10.25 22.28 14.85 10.16 11.31
QRHF-CCSDT-1 10.25 22.28 14.85
QRHF-CCSDT 10.25 22.15 14.85
Note: Geometry, basis set, and full CI results taken from Bauschlicher and Taylor
[62]. QRHF-CC results for 'A1 ionizations from Watts, Gauss, and Bartlett [58]. 3B
calculations used a UHF-CCSD reference function.

Table 5.3 presents FSCC calculations on ionization potentials of methylene in com-

parison with both the full CI data from Bauschlicher and Taylor, and several QRHF-CC

calculations. At the FSCCSD level, four of the five calculated IPs are within 0.1 eV

of the full CI number, with the 'A,--+22A result 0.4 eV off. The effect of the T(3)

correction to FSCCSD is to move all of the results away from full CI, but four of the

roots are still within 0.1 eV, while the A1i-+22A, number has changed to 0.6 eV away

from full CI. Introducing some higher order effects via the T*(3) model shifts all of the

results back toward full CI, but not as close as the FSCCSD result. These differences

from full CI cannot result from uncertainties due to the choice of active space, which

from Table 5.1 are about 0.001 eV (0.04 mH) for this particular molecule. By contrast,


the single reference QRHF-CC calculations start out with all five ionizations within 0.11

eV of full CI, and the inclusion of triple excitations at various levels improves the results

to match full CI within 0.01 eV. This behavior suggests that the single reference CCSD

provides a good approximation to the various states in question, and that triple excitations

are relatively unimportant, which would seem to conflict with the FSCC results.

It is important to recall, however, that there is a fundamental difference between

FSCC and SRCC with respect to the nature of their excitation operators. As discussed in

the previous chapter, spectator triple excitations must be added to the FSCCSD method

to achieve an expansion space that is in some sense equivalent to that of the SRCC. But

even then, the FSCC amplitudes are restricted by the necessity to represent several ion-

ized states simultaneously. Thus, we should not expect the FSCCSD to produce results

identical to a single reference CCSD, nor should we expect triple excitation corrections

to behave in the same way. The differences between FSCC and SRCC shown in Table

5.3 and others, below, should be interpreted as a manifestation of these differences.

These differences between single reference and Fock-space methods can also be seen

in the other calculations presented here (Tables 5.4-5.9). The sign of the FSCC triples

correction is frequently opposite that of the SRCC triples corrections, and their magni-

tudes do not appear to be related a large SRCC triples contribution does not necessarily

imply a large FSCC triples correction, and vice versa. The size of the SRCC and FSCC

triples contributions are rather different as well. The largest SRCC effect is about 0.2

eV, while the FSCCSD+T(3) correction is as large as 0.8 eV, and 0.5 eV for T*(3). All

this is further support for the idea that the nature of the excitation operator hierarchy is

Table 5.4. Vertical ionization potentials of CH2NH (in eV) calculated at the experimental
geometry with a DZP basis.








Experiment [63]
Experiment [64]
Experiment [65]










Notes: Geometries, basis sets, and QRHF-CC results taken from [66].

Table 5.5. Vertical ionization potentials of CH2PH (in eV) calculated at the experimental
geometry with a DZP basis.

la" 5a' 4a 3a
FSCCSD 10.08 10.28 13.17 15.14
FSCCSD+T(3) 10.01 10.33 13.25 15.41
FSCCSD+T*(3) 9.98 10.28 13.25 15.27

QRHF-CCSD 9.90 10.22 13.15 15.15
QRHF-CCSDT-1 10.05 10.24 13.15

Experiment [67] 10.30 10.70 13.20 15.00
Notes: Geometries, basis sets, and QRHF-CC results taken from [66].

different in single reference and Fock-space coupled cluster methods.

Table 5.6. Vertical ionization potentials for 3EEg 02 (in eV) at the experimental geometry.

Method Basis 02 Final State
Method Basis
2rIg 4E- 4iu 4Eu-
FSCCSD DZP 11.76 17.75 16.40 24.16
FSCCSD+T(3) DZP 12.13 17.99 16.30 24.98
FSCCSD+T*(3) DZP 12.08 17.87 16.29 24.66
UHF-CCSD DZP 12.13 17.80 16.15 24.64
UHF-CCSD(T) DZP 11.95 17.77 16.26 24.40

FSCCSD PVQZ++ 12.37 18.34 16.94 24.74
FSCCSD+T(3) PVQZ++ 12.82 18.48 16.74 25.43
FSCCSD+T*(3) PVQZ++ 12.57 18.34 16.71 25.11
UHF-CCSD(T) PVQZ++ 12.39 18.26 16.76 24.82

Expt. (Adiabatic) 12.07 18.17 16.11 24.58
Expt. (Vertical) 12.35 18.33 16.85 24.66

Notes: Geometry taken from Herzberg [68]. Basis sets, SRCC, and experimental result
taken from [69].

Examining the various tables, we see that the SRCC triples correction usually

improves agreement with experiment, while the FSCC triples corrections appear to shift

the calculated result away from the experimental one as often as towards it. It is also

worth noting that the T*(3) correction is generally smaller than T(3) and gives better

agreement with experimental results. Overall, FSCCSD and FSCCSD+T*(3) seem to

do about the same in comparison with experiment, with an average absolute difference

from the experimental results of roughly 0.2 eV. The single reference CCSD results

have an average absolute difference of closer to 0.3 eV, which is improved to 0.2 eV

Table 5.7. Vertical ionization potentials of H2CO (in eV) evaluated at the experimental
geometry with a [5s3pld/2slp] basis.
Orbital of ionized electron
2b2 lbl 5al
FSCCSD 10.51 14.36 15.75
FSCCSD+T(3) 10.96 14.49 16.30
FSCCSD+T*(3) 10.77 14.40 16.06

QRHF-CCSD 10.59 14.19 15.77
QRHF-CCSD[T] 10.63 14.37 15.83
QRHF-CCSDT-lb 10.63 14.39 15.83
QRHF-CCSDT-3 10.58 14.31 15.76
Experiment [70] 10.88 14.38 16.00
Notes: Geometry taken from Herzberg [68]. The basis set consisted of Dunning's
(9s5p)--[5s3p] set on C and 0, and his (4s)--[2s] set with exponents scaled by 1.44
on H [71], with polarization functions exponents of 0.654 (C), 1.211 (0), and 0.7 (H)

by the addition of triple excitation effects. Thus, it would appear that FSCCSD and

FSCCSD+T*(3) perform about as well as single reference CCSD(T) in comparison with

experiment, but the single reference calculations (so far) have an edge in the predictability.

In order to give proper weight to the performance of these methods relative to

experiment, we must consider the potential problems inherent in experimental data, as

mentioned in the Introduction. Most of these problems are displayed by one or more of

the experimental results in this chapter.

The natural result of a FSCC calculation is a vertical ionization potential, while

experiment may provide vertical, adiabatic, or both. Adiabatic IPs are not directly

comparable with calculated vertical ionization potentials, since they depend on the change

Table 5.8. Vertical ionization potentials of N2

(in eV) at R=2.0693 au with a [5s4pld]

Orbital of ionized electron
3aCg l 2au
FSCCSD 15.44 17.14 18.64
FSCCSD+T(3) 15.58 16.86 19.04
FSCCSD+T*(3) 15.44 16.88 18.81

QRHF-CCSD 15.39 16.69 18.71
QRHF-CCSDT-1 15.36 16.87 18.57
Experiment 15.60 16.98 18.78

Notes: Geometry, basis, and QRHF-CC and experimental results taken from Rittby and
Bartlett [73].

Table 5.9. Vertical ionization potentials of ketene and diazomethane (in eV) calculated
at experimental geometries calculated with a DZP basis.
Ketene Diazomethane
2bl 2b2 2bl 2b2
FSCCSD 9.40 14.25 8.60 13.79
FSCCSD+T(3) 9.39 14.49 8.52 13.85
FSCCSD+T*(3) 9.35 14.38 8.47 13.77

QRHF-CCSD 9.29 14.31 8.43 13.78
QRHF-CCSD(T) 9.38 14.24 8.54 13.73

Experiment 9.64 (a) 13.84 (a) 9.00 (a) 14.13 (v)

Notes: Geometries, basis sets, and experimental results taken from Rittby, Pal, and
Bartlett [55]. Experimental IPs are adiabatic or vertical, as indicated in the table.

in geometry from the neutral to the ion. If the ion and neutral geometries are very similar,

the vertical and adiabatic IPs will also be nearly the same, but in general, the adiabatic


IP is only a lower bound for the vertical one. In cases where the experimental value

is adiabatic or unspecified, a larger difference from the computed vertical IP must be


A second problem with experimental results is that, as with theoretical methods,

different experimental techniques can give somewhat different results. This is illustrated

in Table 5.4, which lists three sets of experimental data for methyleneamine, differing

among themselves by as much as 0.2 eV.

A third complication of using experimental data is that details of the calculation

other than the correlation method become important, particularly the basis set. Most

of the calculations presented here use a DZP or slightly better quality basis set. While

very useful for comparing one theoretical method to another, these basis sets may not be

adequate, however, to obtain good results compared to experiment. Table 5.6 shows the

ionization potentials of 02 as calculated using both DZP (32 functions) and an augmented

PVQZ basis (126 functions) in order to demonstrate the effect of basis set on both

the FSCC and SRCC IPs. It is clear that the PVQZ++ basis provides a much better

comparison with the experimental results than the DZP basis, for both FSCC and SRCC

calculations. The change in FSCC results due to increasing the basis set is about 0.5-0.6

eV, and for the SRCC 0.4-0.5 eV. Although there is no experimental comparison, a

similar basis set effect can also be seen in the HOF results in Table 5.2, where the

change from a DZP to a PVQZ basis changes the IPs by 0.2-0.4 eV. It would seem,

then, that correlated calculations with DZP-quality basis sets should not be expected to

reproduce experimental results too well.


Given that the basis sets used here are generally too small, coupled with the other

problems of the experimental data themselves, it would seem that an error of 0.2 eV for

these calculations is not unreasonable, and is not necessarily an indication of intrinsic

problems with the method. Indeed, the fact that the SRCC results show a similar margin

of error with respect to experiment is a good sign for FSCC. It would seem that the

biggest difference in this regard between the single reference and Fock-space results is

the relative lack of consistency of the FSCC triples corrections, which can probably be

better understood with more calculations.

5.4 Potential Energy Surfaces and Related Properties

Another use of FSCC methods is for the study of potential energy surfaces, especially

with the ability to look at several surfaces in the same calculation. The basic result of

the FSCC calculation is still an energy difference, but when added to the absolute energy

of the (0,0) sector reference, it provides an absolute energy for the state of interest.

In the (0,1) sector case, the state of interest has one electron removed relative to the

reference (0,0) state. While this may seem a round about way of calculating the PES

of a molecule, it has practical advantages in many cases. For calculations on molecules

with low-lying electronic states or symmetry breaking phenomena, direct calculation of

the PES by traditional single reference methods frequently suffers from convergence

and other problems due to the proximity of other states to the one of interest. On the

other hand, the same molecule with another electron is often well-behaved at the same

geometries. This the same as the ionization process, but from the opposite point of


view in this case, we want to talk about the final state, while in IPs, it is traditional

to look at it from the point of view of the initial state.

This application of FSCC methods is a very recent development, and there are only

a few such applications in the literature. Moreover, because they deal with symmetry

breaking problems, there is little conclusive experimental evidence to allow concrete

conclusions to be drawn about the performance theoretical methods. Comparisons are

possible with single reference calculations, though as we shall see, the comparison

of high-level SRCC calculations with Fock-space results does not always lead to an

unambiguous determination of which approach is more reliable. Nevertheless, FSCC is

very well suited to such applications and will undoubtedly see increasing use in the study

of potential energy surfaces (beginning with the next chapter of this work), so it is worth

trying to get a feel for its performance.

Engelbrecht and Liu, in a 1983 paper [74], used multiconfiguration self-consistent

field (MCSCF) plus CI calculations to characterize the 3A2 state of CO2. This state,

which is bent (C2v symmetry), has the SCF occupation

(core)3a12b64a13b25a1 b la24b26al

and interacts strongly with the configuration

(core)3a 2b 4a 3b 5a lb 1a 4b 6a

because both the a, and bl orbitals transform as a" in C, symmetry. As a result, SCF

calculations can obtain a lower energy by going to a Cs structure.


Asymmetric Stretch Potential of 3A2 CO2
12 1 ---
10 10
8 FSCCSD+T*(3)
*-' 6

-2 1
-0.10 -0.05 0.00 0.05 0.10
AR (A)

Figure 5.1. Potential curves for the asymmetric stretch mode (AR = !(R1 R2)) of
3A2 CO2 at various levels of theory. The C2v structure is the DZP/FOCIB geometry of
Engelbrecht and Liu [74].

With suitable methods, such interactions can be treated properly, and reason-

able potential energy surfaces obtained, but this requires the ability to describe

inherently multireference states. Most such methods are not many-body in na-

ture, and of course we prefer a size extensive approach if possible (Chapter


It has been observed that QRHF-based coupled cluster methods are capable of han-

dling many multireference problems when sufficient triple excitation effects are incor-

porated [34]. In Figure 5.1, we see two QRHF-CC potential curves for the asymmet-

ric stretch of CO2, which were calculated from closed-shell CO22- SCF orbitals. The

QRHF-CCSD curve has a double minimum shape, while the QRHF-CCSD(T) method


gives a stable C2v structure. This is consistent with the MCSCF+CI calculations of

Engelbrecht and Liu and their conclusion that the equilibrium structure of this state is


The results of three sets of Fock-space coupled cluster calculations using 2A1 CO2-

as the (0,0) reference are also displayed in Figure 5.1. FSCCSD gets the proper shape

of the PES even without triple excitation effects, and the inclusion of triples makes

small adjustments to the shape of the surface. (Note that, as in the IP calculations, the

FSCCSD+T*(3) results lie between FSCCSD and FSCCSD+T(3).) This behavior is what

we might expect from a method that is designed in a multireference framework.

Another symmetry-breaking problem to which a number of single and multiref-

erence approaches have been applied is the nitrate radical, N03. Different exper-

iments have been interpreted in terms of either a totally symmetry D3h equilib-

rium geometry, or a lower symmetry C2v structure with one N-O bond length dif-

ferent from the other two. However, as Kawaguchi and coworkers [75] comment,

no one has advanced a model which explains all of the available spectroscopic


In the most recent theoretical examination of the problem, Stanton, Gauss, and Bartlett

[76] located three stationary points on the NO3 surface using the QRHF-CCSD method.

The D3h structure was determined to be a local minimum, and two C2v structures were

found. They are characterized by the relationship between the unique N-O bond distance

and the other two: one long/two short (1L2S) or one two long/one short (2L1S). The 1L2S

structure was also found to be a local minimum, and the 2L1S a transition state 192 cm-'

Table 5.10. Calculated relative energies (in kcal mol1') of C2v and D3h NO3 structures.

DZP basis PVTZ basis
1L2S-D3h 2L1S-D3b 1L2S-D3h 2LIS-D3h
FSCCSD 2.50 1.05 2.94 1.05
FSCCSD+T(3) 3.03 1.64 3.72 1.73
FSCCSD+T*(3) 2.74 1.35 3.36 1.41

QRHF-CCSD -2.59 -2.04 -2.96
QRHF-CCSD(T) 0.34 0.57
Notes: Geometries, optimized at the DZP/QRHF-CCSD level, taken from [76]. The C2v
geometries are characterized by one of the NO bonds being longer (1L2S) or shorter
(2L1S) than the other two. QRHF-CC relative energies taken from [77].

higher in energy than the 1L2S state; both C2v structures being lower in energy than the

totally symmetry structure. Their conclusion is that the molecule has a C2v equilibrium

geometry, but probably has an extremely low barrier to pseudorotation, which would give

spectra with characteristic of a D3h system.

In a later paper focusing on different orbital choices for single reference treatments of

symmetry breaking problems, Stanton et al. [77] calculated relative energies of the D3h

and C2v minima using more sophisticated methods than in their first paper. Their QRHF-

CC results are shown in Table 5.10 along with FSCC results using the same basis sets

and geometries. Using both DZP and PVTZ basis sets, QRHF-CCSD calculations favor

the C2v minimum by about 3 kcal mol-1. As soon as triples corrections are introduced,

either using the noniterative CCSD(T) approximation or the iterative CCSDT-3 model,

the energies shift in favor of the totally symmetric structure, but by only 0.5 kcal mol'1.


FSCC, on the other hand, predicts the D3h structure to be lower for all three methods

shown, but by roughly 3 kcal mol1.

Given previous experience with QRHF-CC methods applied to problems of this sort,

we should have greater confidence in the QRHF-CCSD(T) and T-3 results than in results

without triple excitations, but with the available data, it is impossible to tell whether

the single reference numbers including triple excitations or the FSCC results are more

reliable. As Stanton et al. [77] observe, "An extremely sophisticated (and expensive)

calculation would be required to firmly establish the nature of the D3h stationary point."

C3+ is another molecule for which both experimental and theoretical studies have

produced differing opinions on the nature of the ground state. Coulomb explosion data can

be interpreted in terms of either a floppy linear molecule or a bent structure, depending on

the vibrational temperature of the sample [78,79]. On the theoretical side, the consensus

is that the ground state is a bent 2B2 with the linear 2 u' higher in energy [80]. The

energy difference between these two states is still not a settled issue. This is important

because if the barrier is small enough, the molecule may be quasilinear although the

true ground state is bent, a large amplitude bending motion could bring the molecule

through the linear structure.

Tables 5.11 and 5.12 present the results of FSCC optimizations, vibrational frequency

calculations, and relative energies in comparison with single reference CC results from

Watts, Stanton, Gauss, and Bartlett [80] which include full CCSDT calculations. All three

Fock-space methods give essentially the same geometry for both structures, comparing

most closely with the single reference CCSD results. For the vibrational frequencies of the

Table 5.11. Optimized geometries and

vibrational frequencies of 2B2 C3+, using a DZP

Metd R () 0 () wl (a,) w2 (al) w3 (b2)
Method R (A) (0) c ( (m'
(cm- ) (cm- ) (cm-1)
FSCCSD 1.328 71.6 1685 583 1269
FSCCSD+T(3) 1.336 70.2 1643 614 1221
FSCCSD+T*(3) 1.335 70.5 1651 603 1225

CCSD 1.337 68.3 1668 687 1215
CCSD(T) 1.350 68.3 1601 515 1074
CCSDT-1 1.350 69.0 1558 569 1078
CCSDT-2 1.347 69.0 1612 595 1213
CCSDT 1.350 68.0 1601 638 1021
Notes: Basis set and SRCC results taken from Watts, et al. [80].

Table 5.12. Optimized geometries, and vibrational frequencies of 2Eu+ C3+, and energy
relative to 2B2 C3+, using a DZP basis.

E(2 + -
(Method R) w(ag) w(ru) E(2 -
Method R (A) (cm-1 (cm-1) (cm-1) 2B2)
(cm4) (cm') (cm-')
(kcal mol')
FSCCSD 1.312 88 1227 492i 7.1
FSCCSD+T(3) 1.315 71i 1203 1911 10.9
FSCCSD+T*(3) 1.316 532 1205 1462 8.6

CCSD 1.314 2500i 13.2
CCSD(T) 1.327 1431i 2.9
CCSDT-1 1.339 1135 -1.7
CCSDT-2 1.329 1595 2.2
CCSDT 1.327 451i 4.9
Notes: Basis set and SRCC results taken from Watts, et al. [80].
2B2 state, the FSCC results are also basically the same, and correspond most closely with


the CCSDT-2 approximate iterative triples model among the single reference calculations.

On the basis of several calculations with a larger basis set, Watts and coworkers suggest

that the true harmonic frequencies should lie in the ranges 1590-1629, 592-638, and

1094-1173 cm-1, and compared to these, the FSCC wl and w3 are slightly high.

The asymmetric stretching frequency (au) of the 2Eu+ state (Table 5.12) is a more

complex problem. An imaginary frequency in this mode indicates that the Doh structure

is a transition state, and a lower energy can be obtained by changing the bond lengths

asymmetrically. For both single reference and Fock-space methods, the calculated

frequency varies widely. Also, the bending mode (Tru) changes dramatically with the

choice of FSCC method. Such large changes due to different correlation methods are

usually an indication that important aspects of the problem are not being taken into

account, and a more detailed examination, using both single reference and Fock-space

methods would be in order.

Table 5.12 also presents the energy differences between the linear and bent structures.

Note that the single reference CCSDT-1 method places the linear structure lower in

energy than the bent one. Watts et al. attributed this to an apparent tendency for the

iterative nature of the CCSDT-1 method to overemphasize problems in the closely-related

CCSD(T) model in difficult cases, thus producing poorer results than the perturbative

CCSD(T) method.

Their CCSDT AE of 4.9 kcal mol-1 is in reasonable accord with CISDTQ calculations

of Grev, Alberts, and Schaefer [81,82] indicating a 4.1 kcal mol-' difference. Taylor et

al. [83] used multireference CI methods to obtain a separation of 1.68 kcal mol-' with


a DZP basis set, and claim that this should be very close to the full CI limit. These

authors also show that multireference effects are important for both the bent and linear

structures. Despite the multireference character of the FSCC approach, all three methods

presented here give a significantly larger energy gap than Taylor's calculations. The

problem with the relative energies of the two states may be due, at least in part, to the

apparent problems in properly describing the 2Eu+ state. It would seem, on this minimal

evidence, that where single reference methods including triple excitations have problems,

Fock-space methods do as well.

5.5 Conclusions

A variety of results has been presented using the Fock-space coupled cluster methods

FSCCSD, FSCCSD+T(3), and FSCCSD+T*(3), for both ionization potential problems at

fixed geometries, and in cases where the potential surfaces must be explored.

Overall, the performance of the FSCC methods is good. FSCCSD results usually

compare well with single reference CCSD(T) calculations. The FSCC triples corrections

seem initially to be less predictable, than the triples correction in single reference methods.

The Fock-space triple excitations have a different nature than their SRCC counterparts,

and it will require more experience with these methods to understand them and discern

the system they follow.

Including some higher-order contributions in the T*(3) method seems to have a

"dampening" effect compared to the basic T(3) model, which causes the T*(3) results to

compare better with experimental results or other calculations. This is probably similar


to the single reference coupled cluster case, where the simplest perturbative inclusion

of triples, CCSD+T(CCSD) was found to be significantly improved by the inclusion

of certain higher order terms, resulting in the CCSD(T) model. For the present, the

FSCCSD+T*(3) method seems to be preferable to the purely third-order FSCCSD+T(3).

In the long run, it is probably worthwhile to investigate more systematic ways of

incorporating higher-order terms in the FSCC triples contribution, as discussed in the

previous chapter.

Numerically, effects of changing the active space in the FSCCSD+triples methods

have been demonstrated to introduce uncertainties as large as 2 mH in calculated

energies (IPs or state energies), however most of the molecules examined have much

lower uncertainties. On the whole, both FSCCSD and FSCCSD+T*(3) seem to do about

as well compared to experiment as SRCC methods including perturbative triples, and

better than CCSD. For very accurate comparison with experiment for both FSCC and

single reference methods, however, a large basis is required.

For energy differences on potential energy surfaces, where the requirements are tighter

than on IPs, Fock-space results do not seem to compare as well with single reference

methods as for IPs. Although the signs of relative energies evaluated using FSCC seem

to match higher-quality single reference calculations, the Fock-space energy gaps are

noticeably larger than in the SRCC case. The cases examined are hard problems, both

experimentally and theoretically, and it is not clear that either approach can be claimed

to be better. Consequently, great care must be taken with such results, and more study

would be most useful.


The FSCC methods all seem to do quite well in determining local properties of

potential surfaces geometries, vibrational frequencies, etc. In most of these cases,

there does not appear to be a substantial difference among the three FSCC methods

examined, except for linear C3*, for which both single reference methods including triple

excitations and Fock-space methods seem to have problems.

In addition to the quality of results of a new method, one should also consider its ease

of use and cost in comparison with competitive methods. By this measure, FSCC methods

have a decided advantage over the single reference approaches. An FSCC calculation,

even including triple excitation effects, amounts to a small fraction of the cost of the

(0,0) sector CCSD reference calculation. In contrast, each of the IP calculations required

separate SRCC calculations for the neutral, and for each ion state, and at least one of

those on an open-shell molecule, for which the SRCC is about four times more expensive

than a closed shell calculation. In most cases, the reference CCSD calculation used in

the Fock-space approach is a closed shell, and all of the IPs can be computed in a single

calculation. For potential energy surfaces, the ability to explore several surfaces at once

is also a distinct advantage over single reference methods.


6.1 Background

Although there is relatively little experimental data on C3+, discussed in the

preceding chapter, there is even less known about Cs+. The experiments which

have observed Cs+ have done so with a mass spectrometer in reactivity [84], frag-

mentation [85-87], or mobility [88] studies. The essential conclusion of these pa-

pers (with respect to C5+) is that there is a single isomer (multiple isomers show

up starting with C7) and that its reactivity and mobility are similar to those of

its neighbors, C4+ and C6+, and the whole set of clusters from C2z to C6 are

commonly labeled as "linear" molecules [84,88]. (C3+, which theoretical results

agree is actually bent, fits the trends for "linear" clusters in the experimental data

as well as C5+ does, suggesting that the "linear" label cannot be interpreted rigor-


Recently, matrix-isolated infrared spectroscopy experiments, using isotopically

mixed graphite, produced an IR band that could be attributed to a planar pentago-

nal (Dsh symmetry) C5+ molecule [89]. The neutral Cs cluster has been well es-

tablished, both theoretically and experimentally, to have a linear ground state [90],

with cyclic structures much higher in energy. In this situation, it would seem un-

likely that the cation would be so different from the neutral that a cyclic cation


would be observable, but the existence of isoenergetic linear and cyclic structures

of neutral C4, which is now generally accepted [90,91], was once also thought un-


A series of SCF and MBPT(2) calculations [92] established that D5h C5+ could not

be responsible for the observed band, and a reanalysis of the experimental data led to

assignment of the band to a C3*H20 adduct, in agreement with a later assignment by

Ortman et al. [93]. The low-level calculations pointed toward a linear ground state

for Cs5, but there were major problems with spin contamination of the UHF-based

calculations, stability of the SCF solution, symmetry breaking, and failure to converge

calculations on some states. In this chapter we attempt to understand the source of the

troubles with lower-level calculations and to establish the true ground state structure of

Cs5 using high-quality many-body methods.

6.2 Computational Methods

The majority of the calculations in this work used a DZP basis set comprised of

Dunning's (9s5p)/[4s2p] double-zeta set [71] and a polarization function (a--0.654) as

optimized by Redmon, Purvis, and Bartlett [72]. All six cartesian d components were

retained, giving a total of 80 contracted basis functions. At selected points on the

various potential energy surfaces, a larger PVTZ basis [61] was used to compute the

energy. This is a (10s5p2dlf)/[4s3p2dlJ] generally-contracted basis, which was used with

spherical harmonic polarization components (5d and 7J), for a total of 150 contracted

basis functions.


Both single reference QRHF-based and Fock-space coupled cluster methods have

been used and in both types of calculations, triple excitation effects have been con-

sidered, mostly through perturbative corrections: QRHF-CCSD(T), FSCCSD+T(3), and

FSCCSD+T*(3). Both QRHF-CC and FSCC calculations used the closed shell neutral Cs

molecule as a reference. The reasons for using both SRCC and Fock-space approaches

are both historical and practical. When work on this project was first begun, no effi-

cient FSCC program was available and QRHF-based approaches were the best choice to

deal with the problems observed in the UHF calculations mentioned earlier. When the

newly-completed ACES II FSCCSD code was applied to this molecule, the nature of the

problems became apparent, and it was clear that understanding the Cs5 molecule was

a problem well-suited to this method because of its closely-spaced, strongly interacting

electronic states.

In some calculations, the carbon Is orbitals and the corresponding virtual orbitals (five

occupied and five virtual for the DZP basis with cartesian angular momentum functions,

five occupied and four virtual for the spherical harmonic PVTZ basis) were not correlated

in order to reduce the computational cost. The effect of dropping these orbitals on the

relative energies is negligible.

Geometry optimizations have been carried out using QRHF-CCSD and all three

Fock-space methods. Analytic gradient techniques were used with the QRHF-CCSD

calculations, but gradients have not been formulated for FSCC yet, so optimizations

were performed via numerical differentiation of FSCC energies. Similarly, vibrational

frequencies were evaluated from the numerical second derivatives of the FSCC energy,