vao

a
. I!
.1...

1%

F thl‘h—ilcvl .5Hh- rdl.fil .

.
9...

«mm.

262%.
‘.'l

a

.0
13;.

t

.u .sdl 2 7

A 53.9.”. .3271
:5.‘ My}...
3.7.3.

1.5:...
93.3.1..-
urns; »
\1! I .
{25.15 . ..
2... .v

.733 . ‘ ;

16:25an7. 2 «Sin

| » . . 35.32.. :-
1<Vl".uh‘
i 5.53.... ,
V :3... $1....

"K
9%!
on.

a... .
130.....r
I—lllnl'.

3:53.233 . . .
1. bu. mm. :
nib

. v.8: .a1‘
15“. \f:

.'.l .9\
In): I;

2:}

‘ 1%... 4 r

011.3'uv

sh
5'51:

L vv bl.-. ~v
IvI-vlbi.

‘3‘...
.11...‘ .
it"):
v! ....

11.4.? n»?

327

H F»...

i1!«v

atlv. ~u\

x
3

s:
‘0 o.."v
.‘ihlf...

2...... as.
. 2

~ \txnt
Viz-5....

:5

{u 7.. .

.2‘. .pxt‘

6:: I.. v.
..I\:....v0
A v.‘

‘ villi...
{11A . . VI. i».

0

'Iv -\l:
I li‘cl‘

%\ -
kl (.an
.1. . .r:
I H‘Kivtf.
3|...“ l!
.1! 31”- ¢ 8

|l'

 

THESIS

JIlllllllllllllllllillllllll llllllllll L

3 1293 01565 9299

 

LIBRARY ‘
Mlchigan State 1
Unlversity

 

 

 

This is to certify that the
dissertation entitled

Advances and Studies

in Computational Chemistry

presented by 2

David Carl Young

has been accepted towards fulfillment
of the requirements for

Ph . D degree in Chemistry

 

 

7%

'Major professor

Date 17 March 1997

 

MS U Lt an Affirmative Action/[q ual Opportunity Institution 0-12771

 

PLACE IN RETURN doxro roman this checkout from your record.
TO AVOID FINES return on or More data duo.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU loAn Affinnatlvo ActionlEcpol Opportunity Institution
mm:

 

ADVANCES AND STUDIES
IN COMPUTATIONAL CHEMISTRY

By

David Carl Yotmg

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Chemistry

1997

ABSTRACT

ADVANCES AND STUDIES
IN COMPUTATIONAL CHEMISTRY

By

David Carl Yormg

Chapter I is an ab initio study of the electronic structure of the mono
hydroxide and oxide cations of Ti, V and Cr. The bonding structure and
energetics of these systems are characterized at the MCSCF and
MCSCF+1+2 levels of theory. Limitations and advantages of other levels of
theory are also discussed. Analysis of bonding structure is presented in the
form of Mulliken population analysis, contour diagrams and qualitative
descriptions of the natural orbitals. Low energy excited states and
dissociation thermochemistry are examined at the MC SCF+1 +2 level of
theory.

Chapter H attempts to answer the pedagogical question of whether

there is a physical reality behind the model of electron lone pairs. The

limitations of orbital descriptions are discussed. Total electron density,
molecular electrostatic potential, electron density difierence and the Laplacian
of electron density are analyzed with two and three dimensional scientific
data visualization techniques. The molecular electrostatic potential is found
to be the most consistent, physically meaningful description of . lone pair
electrons.

Chapter III is an evaluation of the potential for variational quantum
Monte Carlo (VQMC) techniques to be used in describing transition metal
systems. Tests of the efl‘ect of changing various aspects of the algorithm are
presented. A measure of calculation efficiency, called computational work
(CW) is proposed. Several new algorithms for the efi'rcient integration of
heavy element wave frmctions are proposed. A comparison of the efficiencies
of a number of all electron and core potential methods is made.

Chapter IV is a discussion of the existing theorems describing the
behavior of the total electron density. All existing theorems which describe
the electron density for arbitrary systems are presented. General descriptions
of electron density and exceptions are given for systems not covered by
formal theorems. A discussion of other molecular properties, which are

sometimes confused with the electron density, is presented.

To my wife, Natalie

ACKNOWLEDGMENTS

First and foremost I would like to thank my wife, Natalie. I would not
have been able to pursue a PhD. without her support both emotionally and
financially. I hope to make this up to her by giving her the chance to write
the book that she has been hoping to for so many years.

My children Gregory and Ariel can not remember a time that I was not
in school. They have learned the meaning of patience while waiting for their
turn on the computer.

Dr. James F. Harrison, my advisor, has been most patient as I delved
into one area of computational chemistry after another. I would like to thank
him for letting me run with ideas even when he was skeptical of them.

Ed Tisko is my fellow mathematician amongst the sea of chemists. I
have benefited greatly from our many long discussions.

Jesse Edwards is our resident expert on core potentials. I would like to

thank him for making that portion of my work so easy.

As I decide upon the teaching techniques that I like best, I find that
most are those that I learned from Dr. Schwendeman. I would like to thank
him for setting such a professional example. I must admit that I have
occasionally been at odds with other faculty members when I expected them
to live up to the standards that Dr. Schwendeman sets for himself.

The funding for these studies was provided by an NSF equipment grant
CHE 9321436, my receipt of the Dow Graduate Fellowship and an
educational support grant from the Camille & Henry Dreyfus Foundation.
Finally, I would like to thank Nichols Research Corporation which was
willing to hire me with the imderstanding that they would accommodate my

finishing my Ph.D. in abstentia.

TABLE OF CONTENTS

LIST OF TABLES ....................................................................................... xi

LIST OF FIGURES .................................................................................... xiii

LIST OF SYMBOLS AND ACRONYMS ............................................... xviii

INTRODUCTION ....................................................................................... 1
CHAPTER I

METAL OXIDE AND HYDROXIDE CATIONS ....................................... 3

A. Introduction ................................................................................... 3

B. Qualitative description of bonding .................................................. 5

C. Thermochemistry of dissociation .................................................. 13

D. Quantitative results ...................................................................... 16

E. Wave fimction analysis ................................................................ 21

F. Low lying states ........................................................................... 30

G. Wave fimction tests ..................................................................... 36

1. SCF calculations ................................................................. 37

2. MCSCF calculations ........................................................... 4O

3. GVB calculations ............................................................... 40

4. MCSCF+1+2 calculations .................................................. 4]

5. SCF+1+2 calculations ........................................................ 41

H. F inal wave flmction description ................................................... 42

1. Conclusions .................................................................................. 44

List of References ............................................................................. 46
CHAPTER II

LONE PAIR ELECTRONS ON THE WATER MOLECULE ..................... 48

A. Introduction ................................................................................. 48

B. Visualization Techniques ............................................................. 50

C. Orbital Representations ................................................................ 54

D. Total Electron Density ................................................................. 60

E. Molecular Electrostatic Potential .................................................. 63

F. Electron Density Difference ......................................................... 71

G. Laplacian of Electron Density ...................................................... 73

H. Computational Details ................................................................. 75

1. Conclusions .................................................................................. 79

List of References ............................................................................. 81
CHAPTER III

VARIATIONAL QUANTUM MONTE CARLO METHODS .................... 83

A Introduction ................................................................................. 83

B. Random numbers ......................................................................... 89

1. Importance ......................................................................... 89

2. Generation .......................................................................... 91

3. Testing ............................................................................... 93

C. Theory ......................................................................................... 97
D. Computer implementation .......................................................... 109
E. Results and discussion ............................................................... 126
1. Hydrogen atom tests ......................................................... 127

2. Helium atom tests ............................................................. 134

a. Literature survey ..................................................... 134

b. Correlation fimction studies .................................... 138

c. Other tests .............................................................. 155

3. Scandium atom tests ......................................................... 160

a. Simple Metropolis sampling ................................... 165

b. Split 1 technique ..................................................... 165

c. Variable r technique ............................................... 166

d. Split 1 b technique .................................................. 167

e. Split 1, with importance sampling and rejection ...... 169

f. Split 1.’ b, correlating 3S 3p 3d 43 .............................. 171

g. Split ‘t b, correlating 2s 2p 3s 3p 3d 45 .................... 172

h. [Ne] core ECP ........................................................ 173

i. [Ne] core ECP with Jastro fimction ......................... 174

j. [Ar] core ECP ......................................................... 174

k. [Ar] core ECP with Jastro fimction ......................... 17 5

1. Accurate VQMC calculations .................................. 176

F. Conclusions ................................................................................ 179

List of References ........................................................................... 181

CHAPTER IV

. PROPERTIES OF TOTAL ELECTRON DENSITY ................................. 185
A. Introduction ............................................................................... 185

B. Misconceptions about electron density ....................................... 185

C. Description ................................................................................ 187

D. Theorems .................................................................................. 189

E. Cases not covered by Theorems ................................................. 192

F. Other molecular properties ......................................................... 194

G. Conclusions ............................................................................... 199

List of References ........................................................................... 201

1.1

1.2

1.3

1.4

1.5

1.6

1.7

2.1

2.2

3.1

3.2

3.3

3.4

3.5

LIST OF TABLES

Metal - Oxygen bond lengths and dissociation energies ..................... 18
Oxygen - Hydrogen bond lengths and dissociation energies ............... 19
MCSCF Net Atomic Charges ............................................................ 21
Metal valence orbital populations (MCSCF) ...................................... 22
Oxygen and hydrogen valence orbital populations (MCSCF) ............. 25
Total electronic energies (in hartrees) ................................................ 35
Total SCF electronic energies (in hartrees) ........................................ 39
Oxygen atom basis set ....................................................................... 7 7
Hydrogen atom basis set .................................................................... 7 8
Time complexity of computational chemistry techniques ................. 110
Time complexity of VQMC steps .................................................... 116
Efl'ect of changing initial box size .................................................... 130'
Effect of changing box size .............................................................. 132
Dependence upon number of iterations ............................................ 133

3.6

3.7

3.8

3.9

CPU time test for optimized algorithms ........................................... 158

Sample frequency tests .................................................................... 169
Energies using ECP ......................................................................... l7 7
Computational work and energy ...................................................... 17 8

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

1.10

1.11

1.12

2.1

LIST OF FIGURES

TiOIF MCSCF+1+2 energy level diagram ........................................ 14
VOH+ MCSCF+1 +2 energy level diagram ......................................... 15
CrOH” MCSCF +1 +2 energy level diagram ........................................ 17
TiO+ MCSCF sigma bond contour plot .............................................. 23
TiO+ MCSCF pi bond contom plot .................................................... 24
TiOH" MCSCF sigma bond contour plot ........................................... 26
TiOH" MC SCF sigma OH bond contour plot ..................................... 27
CrO” MCSCF pi bond contour plot ................................................... 28
CrO+ MCSCF pi" bond contour plot .................................................. 29
CrOI-F MCSCF sigma bond contour plot ........................................... 31
CrOH+ MCSCF sigma OH bond contour plot .................................... 32
CrOH" MCSCF 0 LP contour plot .................................................... 33
Contour plot of the electron density in the plane of the HZO

molecule. The innermost contour corresponds to 0.142 e/A3

and subsequent contours decrease imiformly by 0.00570 e/A3 ........... 51

2.2

2.3

2.4

2.5

2.6

2.7

2.8

Contour plot of the electron density in the plane perpendicular to

the molecular plane and bisecting the HOH angle. The

innermost contour corresponds to 0.142 e/A3 and subsequent

contours decrease uniformly by 0.005 70 e/A3 .................................... 52

Mesh plot of electron density in the plane of the water molecule,
oriented with the hydrogens towards the viewer. The peak on

the oxygen atom has been tnmcated in order to make the

hydrogens visible ............................................................................... 53

Colorization of the total electron density in the plane of the

water molecule .................................................................................. 55

Colorization of the total electron density at right angles to the

plane of the water molecule ............................................................... 56

Colorization of the total electron density both in the plane of and

perpendicular to the water molecule .................................................. 57

Three dimensional isosurfaces of the water total electron density ....... 58

The 0.0327 e/A3 isosurface with an inscribed sphere (green)

centered at oxygen. The radius of the inscribed sphere is just

large enough to break through the isosurface ..................................... 62

2.9

2.10

2.11

2.12

2.13

2.14

2.15

2.16

2.17

Contours of the electrostatic potential in the plane of the I120

molecule .......................................................................................

Contours of the electrostatic potential in the plane perpendicular

to the plane of the HZO molecule ..................................................

Colorization of the electrostatic potential in the plane of the

water molecule .............................................................................

Colorization of the electrostatic potential in the plane

pependicular to the plane of the water molecule ............................

Colorization of the electrostatic potential in the plane of and in

the plane perpendicular to the plane of the water molecule ...........

Three dimensional nature of regions of low electrostatic potential
The molecular electrostatic potential (MEP) in water. The blue

surface corresponds to a MEP of +27.2 eV while the yellow

sm'face corresponds to a MEP of -2.041 eV .................................

The 0.27 e/A3 isosm'face of the density difference for water

shows regions of increased electron density relative to the

separated atoms in the plane of the molecule behind the oxygen

atom (top lobe) and in the bonding regions (bottom two lobes) .....

The -3.973 isosurface of V2 p .......................................................

..... 64

..... 65

..... 66

..... 67

..... 68

..... 70

..... 72

..... 74

..... 76

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

4.1

4.2

4.3

Correlation fimction f(ry .................................................................. 104

 

Step size bias .................................................................................. 122
Plot of the function flr12)= exp( l/r‘2 ) ............................................ 145
Plot of the fimction f(r]2) = exp( - r12 ) ............................................. 146
Plot of the flmction f(r,2) = exp( -1/r,2 ) ........................................... 146
Plot of the fimction firm) = exp( r12 ) ............................................... 148
Plot of the function f(r,,) = exp(—02 r5”) .......................................... 148
Plot of the function f(r,,) = exp(—0.075132” + 0.3 64499r,‘;-”) .................... 150
Plot of the fimction f(r,,) = exp( 0'5 r” J ......................................... 150
1+rsq,
Electron density for the ground state of the H atom .......................... 188

Electron density extending away from the nucleus of a hydrogen
atom in the excited state with the electron in a 2p orbital ................. 190
Total electron density in the plane of the C2114 molecule plotted
as the height above the plane. The density was obtained fi'om a

CISD calculation with a 6-311G”++ basis set ................................ 193

4.4

4.5

Contour diagram of the electron density of the CO molecule.
The carbon atom is on the left. Contours are plotted for density
values of 0.1, 0.2, 0.3, 0.4 0.5, 1.0 and 5.0 electrons/bohr3. The
density was obtained from a configuration interaction calculation
with a 6-311G*+ basis set. The computed bond length was
1.125 Angstroms ............................................................................. 195
Plot of the electron density along the internuclear line of the Li;
molecule illustrating the magnitude of the maximum in the
middle of the bond. The density was obtained fi'om a CISD
calculation with a 6-311G*+ basis set and a bond length of

2.692 A ........................................................................................... 196

p. p(r)

~68 '6

is:

AREP

CASSCF

LIST OF SYMBOLS AND ACRONYMS

laplacian operator

ionization potential

molecular electrostatic potential

wave length

frequency

electron density

parameter in variable step size algorithms
one electron orbital

electric potential

wave fimction

antisymmetrizing operator

electric field

average relativistic efl'ective core potential

complete active space self consistent field

CC
CI
CPU
CSF

CW

DHS

DMC

(E)

ECP
EO
eI(x)

f(r,-,)

coupled cluster
configuration interaction
central processing unit
configtuation state fimction
computational work

Slater matrix

determinant of alpha spin fimctions
determinant of beta spin fimctions

element ij of the Slater matrix

direct inversion of the iterative subspace
difiusion Monte Carlo

exact energy

analytically integrated expectation value
expectation value fiom N iteration VQMC

effective core potential
equivalent orbitals
local energy

correlation fimction

probability distribution

GFMC

GTO

GVB

t,

MCSCF

Green's function Monte Carlo
Gaussian type orbital
generalized valence bond
Hamiltonian operator

Hartree - Fock

value retruned by a linear congruential random number
generator on iteration j

Multi-Configuration Self Consistent Field

MCSCF+1+2 Single and double excitation CI from the reference space

MEP
MRCI
N

N

N(a)

W>

'm.

P):

of an MCSCF calculation
molecular electrostatic potential
Multi-Reference Configuration Interaction
number of electrons
number of iterations
normalization constant
projection operator

distance from center of mass to the nucleus farthest from the
center of mass

momentum

weight of determinant k in a many determinant wave function

QMC

if

m

SCF
SMO

STO

VQMC

VSEPR

quantum Monte Carlo

charge

distance between two electrons

distance between an electron and a nucleus
distance between two nuclei

self consistent field

spectroscopic molecular orbital

Slater type orbital

CPU time

expectation value of kinetic energy
expectation value of potential energy

variational quantum Monte Carlo
valence shell electron pair repulsion
variance

charge on nucleus n

INTRODUCTION

This thesis is arranged as a series of studies each of which is a separate
chapter. Each of these can be taken as an independent individual study.
However, as a collection these studies present representative samples of four
different aspects of quantum chemistry.

Chapter I presents a study of the electronic structure of the mono
hydroxide and oxide cations of Ti, V and Cr. This is a fairly traditional
quantum chemistry study although a difficult case. Open shell, coordinatively
unsaturated transition metal compounds present significant technical
difficulties in applying quantum chemistry techniques. The automated
algorithms found in popular software packages such as Spartan and Gaussian
do not correctly handle these systems. Simpler levels of theory, such as I-IF
theory, do not give even qualitatively reasonable results for these compornrds.
Thus more complex calculations must be employed.

Chapter II focuses on the nature of electron lone pairs, and employs
scientific visualization techniques. These techniques are a more recent form

of data analysis utilizing object oriented graphics software. This study

1

2

provides a good overview of the capability of visualization techniques as well
as attempting to answer a question of scientific and pedagogical importance.

Chapter III is an evaluation of the potential for variational quantum
Monte Carlo (VQMC) techniques for replacing high order Moller-Plesset
perturbation theory and other many determinant techniques for high accuracy
descriptions of transition metal systems. The theory behind VQMC and
explicitly correlated wave functions as well as issues related to the efficient
computer implementation of these calculations is discussed. Several new
algorithms for the eflicient integration of heavy element wave frmctions are
proposed.

Chapter IV is a discussion of the existing theorems describing the
behavior of fire total electron density. The most general and powerful form of
any mathematical treatment is its theorems and their proofs. This provides an
excellent counter point to the visualization of electron density for specific
molecules illustrated in chapter 11. These theorems represent the most
rigorous description of quantum chemistry. This is the opposite extreme from
simplified models such as VSEPR and the artist’s renditions often found in

text books.

Chapter I
TRANSITION METAL OXIDE AND HYDROXIDE
CATIONS

A. Introduction

The frmdamental understanding of chemistry is built on the systematic
characterization of categories of chemical compounds. Ab initio techniques
provide information about compormds which is extremely difficult to obtain
experimentally. This ab initio study examines the electronic structure of a
number of compounds which are not well described by traditional models
such as Lewis dot structures and VSEPR theory [1]. This characterization is
important not only for understanding these gas phase molecules, but for its
eventual extrapolation to coordination chemistry, materials science and
bioinorganics.

A number of early transition metal oxide and hydroxide cations are
characterized at the Multi-Configurational Self Consistent Field (MCSCF)
and Multi-Reference Configuration Interaction (called MRCI or specifically
MCSCF+1+2 in this study) levels of theory [2]. The qualitative description is

primarily obtained from the MCSCF calculations. More accurate quantitative

results are obtained by doing a configuration interaction calculation from this
reference space (MCSCF+1+2).

Prior to this study, Tilson and Harrison examined the ScO+ and ScOH+
ions at the same levels of theory that are used here [3]. They found that the
ScO" ion has a 12+ grormd state. It is characterized as having a triple bond
formally consisting of two electron pair bonds and one dative bond. The
ScOIF ion has a 2A grormd state in which the Sc - 0 bond is formally a triple
bond consisting of one electron pair bond and two dative bonds.

The CrO+ ion has been examined in a number of prior theoretical
studies [4-6]. New CrO+ calculations were done for this study in order to
have all results at a consistent level of theory. Theoretical studies of the
neutral oxide molecules have been published [7,8], but the studies in
references 3—6 are the only previous cation calculations.

The bond energies of these compormds have been measured
experimentally using mass spectrometric techniques [9-11]. However, there
has been no experimental geometry determination. Other than these
experiments and the theoretical results mentioned in the previous paragraph,

no previous studies have been made of these molecules.

The linear molecules are all arranged with the atoms on the Z axis,
while the bent molecules are in the XZ plane. For the linear molecules, the d
orbitals are labeled according to the symmetry of the orbital (e.g., (1,,z would
be labeled xx). A number of electronic configurations were computed for
each molecule in order to determine the lowest energy electronic structure.
Section G describes the results of tests of a number of wave functions.
Section H gives a detailed description of the basis set and wave frmction used
to obtain the final results.

B. Qualitative description of bonding

Considering how the hydroxide radical might bond to a metal leads to

three possibilities. First is the bent structure prevalent in organic chemistry

which is pictorially represented as

 

A second possrbility is the formation of 1: bonds. Forming a dative 1: bond
with the lone pair electrons on oxygen gives rise to the linear double bonded

structm'e

 

 

 

The third possibility is sp hybridization of the oxygen atom to give the
maximum possrhle bonding interaction, a triple bond with one strong shared

bond and two weaker dative bonds pictorially represented as

 

 

 

The natural orbitals of the MC SCF wave function are used to obtain a
description of the bonding within the hydroxide cations. Note that there is

not a rigorous way to assign which bonds are dative and which are electron

pair bonds. For the linear molecules, the 1t bonds must be equivalent by
symmetry and are kept so in the calculation.

The titanium atom has a ground state electron configuration of 4sz3d2,
while the lowest energy state of Ti” is 4s'3d2. The Ti’ ion has an excited
state with a 3:13 configuration 33 mH higher in energy [12]. When titanium is
bonded to other atoms, a detailed analysis seldom leads to the assignment of a
[significant amormt of 45 character in either the bonds or lone pair electrons of
the grormd state [13,14].

The TiOH+ ground state was formd to be a doubly degenerate 3A state
with a linear geometry. The pictorial representation and Lewis structure for

thisstateare

 

 

 

 

The oxygen atom exhibits sp hybridization as discussed above. Considering
the formation of the Ti-O bond fiom Ti+ and an OH radical would lead to the
characterization of the bond as a triple bond with one electron pair bond and
two dative bonds with both electron pairs coming from the oxygen. We make
the bond assignment based on Ti” since the dissociation to Ti+ is much lower
in energy than dissociation to Ti or Ti”. This triple bond characterization
gives the maximum possible interaction between the Ti and O atoms
assuming energetically accessible orbital occupations. A very low lying 32‘
state with the 8. and 8. orbitals singly occupied was also found.

A doubly degenerate 2A ground state is predicted for the TiO+

molecule. This state is represented as

 

 

 

A low lying 2): state with the do orbital singly occupied was also found.
Without the presence of the hydrogen atom, there is no energetic advantage to
hybridizing the O atom as evidenced by very little mixing of the 0 2s orbital
with it’s p orbitals. Based on having three valence electrons around a Ti” ion
and four p electrons around an O atom, the bond is characterized as a triple
bond with two electron pair bonds and one dative bond from the oxygen.
This bonding again gives the maximum accessrble interaction between the
atoms.

In forming a VOH+ structure analogous to the TiOH+ structure, the
additional electron could be put in the empty d orbital, paired with one of the
unpaired electrons or put in the 4s orbital. The lowest energy state formd was
a 42’ state obtained by filling the last available d orbital. The pictorial

representation and Lewis structure for the ground state are

 

 

 

10

8,. +
00V :ZO—H

5+ c
The ScOI-F molecule studied previously follows the same metal - oxygen
bonding structure as seen here in TiOH+ and VOH+ [3].
A VO” compound with a metal oxygen bonding structure similar to that

of TiO+ will have two unpaired electrons, similar to TiOH“. However, the

computed grormd state is a 32' state, represented as

 

 

 

 

 

Forming a CrOH+ compound with bonding similar to that of VOH+
requires either occupying the 4s orbital on Cr to obtain a pentet state or
pairing electrons on the metal center to get a lower spin triplet state. The

lowest energy structure of the latter type is a 32‘ state with the do orbital

ll

doubly occupied. This is higher in energy than the separated Cr+ (6S) ion and
OH radical (see section G).

Two stable states of CrOH+ were found which do not follow the same
bonding pattern that was observed in the previous hydroxide structures. Note
that the ground state electronic configuration of the Cr atom is 4s‘3d5 rather
than 4s23d4 indicating an additional stabilization due to having a half filled 3d
shell. Likewise, the ground sate of Cr+ has a 3d5 configruation. The
preservation of the halffilled 3d shell around the Cr atom leads to a linear 5’1']
structure with only one 1: bond and a bent 5A; structure with a single 0' bond.
These two states were found to be very close in energy with the 5A. structure
having a slightly lower energy. The pictorial representation and Lewis

structure for the bent grormd state are

cl.
é
coco

 

12

dyz. +

o,, . 7

Cr ---- O
(1,2,,2 o
dxz - 137°\'\
H

The CrO+ structure also deviates from the trend of the previous oxides
in order to maintain a half filled 3d shell around the Cr atom. Two structures

were seen, a 42‘. structure and a 411 structure slightly lower in energy,

represented as

 

5- o
75x OCT: 6—“...
5+.

This molecule has been studied a number of times previously [4-6]. Our

results agree with the most accurate of the previous studies.

13

C. Thermochemist_ry of dissociation
If the oxides could be formed from the corresponding hydroxides by

the tmcoupling of the electrons in the OH bond of the ground state hydroxide,
the result would be a 4A state for TiO+, a 52 state for VO+ and a 6IT state for
GO“. Each of these is a higher spin state than is predicted for the ground
state electronic structure of these molecules.

Energy diagrams presented in F igrues 1.1 and 1.2 depict the hydroxide
at the bottom and the separated atoms at the top. Two routes are presented to
reach the separated products either by first dissociating the metal - oxygen
bond or first dissociating the oxygen - hydrogen bond. There are two levels
associated with the MO+ + H system. The higher energy structure is the
high spin state of the oxide as discussed in the previous paragraph. The
lower energy structure is the ground state electronic configrnation. Both
states reflect optimized bond lengths.

Both TiOIF and VOH+ show the metal - oxygen multiple bonding
structure which was first observed in ScOI-F [3]. The metal oxygen bond in
the grormd state of TiOI-I+ and VOIF is weaker than in the oxide as in the Sc

analogue. This trend is consistent with the characterization of the oxides

I4

 

+IO_._. >

 

IE madn—

I + <~ +O_._./ #

 

I + <v +02. a IO + +_.__..

SE madm—

 

<
I+O++_._.

canoe .32 none $1.582 be: - 2 occur.—

Ev 33am was

r 8.3a-

.. flung-

r 2.2a-

r medmm-

r 8.3a.

a 2.3a-

 

r andmor

15

 

 

 

 

+IO> >
:6 MN:
I + Na +O>/ IO + +>

I .7 Wm +0)

IE 3.3—

 

I+O++>

sense .32 noon mirage: +30> - 3 occur.—

AE 35cm .38.

l and—c7

r and—2-

r end—2..

rad—2-

rem—S-

r 3.32.

Iced—2.

 

r 3.22.

16

having two strong shared bonds and one weaker dative bond while the
hydroxides have one shared and two dative bonds.

Due to the differences in bonding structures between the Cr
compormds and those of Ti and V, the Cr energetics are not expected to
follow exactly the same pattern. However, a similar energy level diagram can
be obtained for the chromium compounds as shown in Figure 1.3 . The (TI
state of CrO+ maintains the half filled shell environment around Cr, but does
not give the maximum possible interaction between the two atoms as formd in
the ground state.

D. Entitative results

The calculated bond lengths and dissociation energies are shown in
tables 1.1 and 1.2 . Past studies have shown that the MCSCF predictions can
be expected to follow the correct trends in bond energy and give good
geometric results [3,14,15]. The MCSCF+1+2 results are expected to give

more accurate bond energies.

17

 

 

I .7 EV +90

 

sense .32 Rocco $1,582 E06 - 3 2:5

ch >326 .83.

r 2..»— _ T

r «6.3 _ T

r 8.3 _ T

r 2.2:.

I GO»— _ T

 

r 3.2 _ T

18

Table 1.1 - Metal - Oxygen bond lengths and dissociation energies

molecule

Scorr+
rion”
VOH
CrOH
ScO
TiO
vo

CrO

r(M-0> ,
MCSCF

3.5093
3.428
3.372
3.403
3.0953
3.044
3.011

3.101

“M02.
MRCI

3.5053
3.378
3.335
3.386
3.1203
3.076
3.048

3.124

162.43
140.8
118.5
45.0
214.13
144.0
95.8

45.3

De(M'O)
MRCI"

172.33
155.4
142.4
63.3
232.93
172.7
145.3

78.65

Values are for the predicted grotmd state configurations

ScO+
'I‘iO+
VO+

CrO+

* MCSCF calculations have the 7: bonds correlated

1 2+
273
3 z-

4n

ScOH"

TiOH+

CrOI-I+

lengths in bohr, energies in mH, references superscripted

De(M-O)
MCSCF'

D0(M'O)
Experiment"

190°, 140‘0
1779, 18010
1659, 171 ‘0
119“), 116“
2639, 254“
2559, 257'1
2209, 209‘1

136“

" MRCI calculations are single and double excitations from the MCSCF

"‘ experimental results from mass spectrometry

19

Table 1.2 - Oxygen - Hydrogen bond lengths and dissociation energies
lengths in bohr, energies in mH, references superscripted

«0m , Dre-r0, Dram ..
molecule MCSCF MC SCF Experiment

 

Soon+ 1.8483 91.63 4810
TiOH+ 1.798 114.3 85‘0
VOH+ 1.798 140.1 11310
CrOH” 1.793 116.6 14410
OH 1.800 117.4 1629

Values are for the predicted ground state configruations

ScOH+ 2A
TiOH+ 3A
VOH+ 42'
@011+ 5A.
OH 211

* MCSCF calculations have the metal - oxygen n bonds correlated, giving an
SCF description of the OH bond
" experimental results from mass spectrometry

20

Both the MCSCF and MCSCF+1 +2 calculations reflect a shorter metal
- oxygen bond length in the grormd state of the oxides than in the
corresponding hydroxide. This shorter bond is accompanied by an increase in
the bond dissociation energy.

Crystal structures for compounds containing titanium - oxygen bonds
yield bond lengths in the range 3.06 - 4.17 bohrs [16,17]. Since the oxide
falls at the low end of this range it is reasonable to assign a triple bond to it.
The hydroxide falls in the middle of the experimental range consistent with
it’s weaker titanirnn - oxygen bond.

Experimental vanadium - oxygen bond lengths fall in the range 2.93 -
4.29 bohrs [16,17]. Having the oxide bond length at the low end of the range
and hydroxide bond length in the middle of the range supports having an
electronic structure similar to that of the titanium analogues.

Chromium - oxygen bond lengths fall in the range 2.86 - 4.03 bohrs
[16,17]. The shorter oxide bond length is consistent with the double bond
assignment in this molecule. Likewise the hydroxide bond length is

consistent with the assignment of a single 0 bond.

21

E. Wave function analysis

A Mulliken population analysis and contour plots were produced from
the natural orbitals of the MCSCF calculations. The Mulliken population
analysis scheme is expected to be reasonable for these calculations since the
basis set did not contain diffuse functions. The net charges show a pattern of
a half electron polarization of the bonds (Table 1.3). This half electron
polarization has been seen in studies of many other transition metal - main

group diatomic and triatomic molecules [3,14,15].

Table 1.3 - MCSCF Net Atomic Charges

Molecule |_S_t_a_t_§ Meg Qtyggn Hydrogen

ScOH+ 2A +1.41 0.77 +0.36 ref3
riOH+ 3A +1.52 -0.82 +0.30

VOI-F 4:- +1.49 -079 +0.29

CrOH+ 5Al +1.45 -0.76 +0.31

Sco+ 12* +1.28 -0.28 ref3
TiO” 2A +1.50 050

v0" 32' +1.44 -0.44

Cro+ 4n +1.43 .043

22

For the diatomics, the a bond is primarily composed of the 2pZ orbital
on oxygen. A contour plot of the TiO+ a bond natural orbital is shown in
Figure 1.4 . Ti has is on the left in all contour plots. The other oxides give a
similar picture. The 1: bonds for the oxides are primarily composed of the p1t
orbital on oxygen bonding to the d, orbital on the metal (Figure 1.5). All of
the bonding orbitals are polarized towards the oxygen with the a bond more

polarized than the 7: bond as shown by the population analysis in Tables 1.4

 

and1.5.

Table 1.4 - Metal valence orbital populations (MCSCF)
State 48 4p, 4p" 4p], 3do 3dr:x 3d7ry 3d8. 3d5+
Scori+ (2A)ref3 0.03 0.10 0.03 0.03 0.10 0.12 0.12 1.0
TiOH+ (3A) 0.13 0.01 0.04 0.04 0.96 0.15 0.15 1.0
VOH+ (42') 0.08 0.02 0.04 0.04 1.01 0.16 0.16 1.0 1.0
CrOH+ (5A1) 0.10 0.05 0.04 0.05 0.31 0.95 1.04 1.0 1.0
Sco+ ('z*)ref3 0.02 0.19 0.05 0.05 0.55 0.43 0.43 1
"HO" (2A) 0.01 0.00 0.03 0.03 0.34 0.54 0.54 1.0
v0+ (32') 0.01 0.00 0.03 0.03 0.30 0.60 0.60 1.0 1.0
Crot (411) 0.02 0.00 0.02 0.05 0.65 0.78 1.05 1.0 1.0

23

Figure 1.4 - TiO+ MCSCF sig bond

 

I.

 

I 1.

II III-

a cii--|009f

I .m/ or I
/. II. I I!

I

   

/ t\ \ \ \ \
I a l' s\§ 0‘: u\ \\ \
.l 1‘ 0‘ .\ \t \

 

 

5.00

 

50

6.

-3.50

2

24

Figure 1.5 - TiO+ MCSCF pi bond

 

 

.".

.0'
----‘ ' '

 

 

5.00

 

 

-5.00

6.50

-3.50

Z

25

Table 1.5 - Oxygen and hydrogen valence orbital populations (MCSCF)

 

State 28 2pc, 2px 2py Hydrogen

ScOH’ (2A) 1.74 1.37 1.83 1.83 0.64 ref3
TiOH+(3A) 1.81 1.41 1.80 1.80 0.70

VOH' (42') 1.81 1.39 1.79 1.79 0.71

CrOH+ (5A1) 1.77 1.43 1.66 1.90 0.69

Sco+ (1):“) 1.80 1.44 1.52 1.52 ref3
r10+ (2A) 1.97 1.68 1.42 1.42

vo+ (32-) 1.98 1.71 1.36 1.36

Cro+ (Ti) 1.99 1.34 1.19 1.89

The two linear hydroxides, TiOIF and VOH+ have 1: bonds with

electron distributions similar to those of the oxides. The 0 bonds in these two

compormds show a dative M-0 0 bond and a covalent OH 0 bond such as

the ones for TiOH+ in Figures 1.6 and 1.7 . In these compormds, the 7: bonds

are more polarized than the 0 bonds. The node in the a bond of Figure 1.6

can reasonably be considered a polarization near the nucleus.

26

Figure 1.6 - TiOH+ MCSCF sigma bond

6.00 l l L

 

 

 

 

 

~6.00 6.00

zaxis

27

Figure 1.7 - TiOH+ MCSCF sigma OH bond

 

anlS

 

 

-6.00
-6.00 6.00

 

 

zaxrs

28

Figure 1.8 - CrO+ MCSCF pi bond

 

       

.°‘~

 

 

 

 

6.50

-3.50

IS

Z

29

Figure 1.9 - CrO+ MCSCF pi* bond

 

 

 

’0' [I I
O
0.0. t/olr . .

 

   

’
. "-.-.--/

\

6.50

 

 

5.00

-3.50

 

-5.00

30

The CrO+ 7t system in the y direction contains 3 electrons, two in the
bonding orbital, which is polarized even further towards the oxygen (Figure
1.8) and one in the 7r“ which is polarized towards the Cr atom (Figure 1.9).

The bent CrOI-F 5A] state has three natural orbitals in the valence a
system, which describe a o(Cr-O), a 0(O-H) and a lone pair on oxygen.
(Figures 1.10, 1.11 and 1.12 respectively.) All three figures are in the plane
of the molecule with oxygen at the origin, Ti at z=-3.403 bohr and H at
z=1.311 bohr, y=l.223 bohr. Figures 1.4 through 1.12 all depict the positive
contours with the solid lines and negative contoms with dashed hues for
contour values +0.002, +0.005, +0.01, +0.02, +0.04, +0.08, +0.16, +0.32,
+0.64, +1.28.

F. Low lm’ states

Low lying excited state wave frmctions were obtained through the use
of symmetry constraints and appropriate initial guesses. The energies of low
energy excited states are summarized in Table 1.6 .

The two TiOH+ states are only different by the occupation of the
impaired electrons. We were tmable to get the MCSCF calculation to
converge to the 32 state. This suggests that the 32 state is not the predicted

ground state at the MCSCF level of theory.

31

Figure 1.10 - CrOH+ MCSCF sigma bond

 

 

/.

 

 

6.00

 

-6.00

6.00

-6.00

32

Figure 1.11 - CrOH+ MCSCF sigma OH bond

 

 

 

6.00

 

-6.00

 

-6.00

33

Figure 1.12 - CrOH+ MCSCF O LP

 

 

 

 

6.00

 

 

6.00

34

The MCSCF+1+2 energies for CrOH+ are too close in energy to be
certain whether the ground state is the bent 5A] state or the linear 52 state.
The chromium oxide cation also has a very low energy excited state [4-6].

The CrO+ ground state has been a matter of debate in the literature [4-
6]. These results agree with the most accruate methods applied thus far, but

are not significantly better than the previous works.

35

Table 1.6 - Total electronic energies (in hartrees)

 

molecule state energy level titheorv
norr: :A -923.74605 MCSCF
non 3A -923.79362 MCSCF+1+2
non+ 2‘ -923.75673 MCSCF+1+2
VOIF 42- 4018.20456 MCSCF
VOH“ 4): -1018.26626 MCSCF+1+2
CrOH: :Al -1118.58523 MCSCF
CrOH+ 5Al -1118.67616 MCSCF+1+2
CrOH+ 511 411857381 MCSCF
CrOH 11 -1118.67476 MCSCF+1+2
no“ 2A -92313249 MCSCF
no“ 2A -923.17622 MCSCF+1+2
no“ ‘11 92302322 MCSCF
r10“ “A -923.14174 MCSCF+1+2
v0“ 32: -1017.56520 MCSCF
vo“ 3): -1017.62547 MCSCF+1+2
v0“ 52 -1017.51691 MCSCF
vo“ 5: -1017.60985 MCSCF+1+2
00: :11 411796931 MCSCF
Cro+ 4r1 -1118.04068 MCSCF+1+2
00+ 42 -1117.94804 MCSCF
CrO 2 -1118.02701 MCSCF+1+2
cm“ 611 4117.93228 MCSCF
etc“ 611 -1118.06084 MCSCF+1+2

36

G. Wave function tests

All of these calculations have been done with an extension of
Dunning’s reoptimization of the Wachter basis set [19]. The precise
description of this set is given in the next section. This set was chosen
because it has proven to be reliable in many previous studies [3,14,15]. This
also has the advantage of keeping the results of this study comparable to
these previous studies. The more difficult question was how to include
correlation in the wave frmction.

It has been well documented that correlation is crucial to obtaining
even a qualitatively correct description of similar chemical systems
[3,4,6,14,15,18]. This of course immediately opens the question of how to
construct a good wave frnrction.

Excellent results could be obtained with heavily correlated methods
such as firll CI’s or complete active space self consistent field (CASSCF)
calculations. Unfortrmately these methods are too computationally intensive
to be done on the crurent generation of work station class computers. In HF
calculations, the difficulty is in computing the integrals which scale as N4
where N is the number of basis frmctions and calculations with the same basis

set are being compared. In many-determinant calculations the number of

37

integrals has not changed but becomes a relatively minor problem compared
to the computation of the configmation list (list of individual determinants)
which can scale as poorly as N8.

In the course of this study a number of types of wave frmctions were
used. Some of these are useful only as intermediate steps and thus omitted
from the previous discussion. For the sake of completeness all of these levels
of theory are discussed here.

1. SCF calculations

SCF or HF calculations (the terms are synonymous in this day) are
known to be ofien qualitatively incorrect to the point of predicting the wrong
ground state and ordering of states for metal - main group systems. However,
SCF calculations are still a useful first step in the process of arriving at the
final results. SCF wave frmctions are often used as an initial guess for more
complex calculations. Also, although the ordering of nearby states may be
incorrect, states which show a very large difference in energy at the SCF level
can be expected to maintain a 1arge spacing with more accurate wave
functions.

Table 1.7 lists the energies of various electronic states at the SCF level

of theory. Bonding orbitals are designated by Greek letters (1! or o) and

38

orbitals centered only on the metal are designated by English letters (8 or (1).
These initial calculations did not include a geometry optimization although an
effort was made to use as reasonable a geometry as possible.

All of the metal hydroxides listed in Table 1.7 were calculated as linear
molecules except for the 5A1 CrOH+ calculation, which had a 105 degree
bond angle. The TiOH+ calculations use a Ti-O bond length of 3. 12 bohr and
an O-H bond length of 1.78 bohr. The VOHI+ calculations use a V-0 bond
length of 3.43 bohr and an OH length of 1.80 bohr. The CrOI-I+ calculations
use an O-H length of 1.79 bohr and a Cr-O bond length of 3.37 bohr, except
for the 5A1 which used a 3.369 bohr length and the 5IT which used a 3.323
bohr bond length. All of the hydroxide calculations showed Ti-O and O-H 0
bonds (included in the core) except for the 5 IT CrOH+ calculation which had
one 1! bond with no a bond.

The TiO+ calculations were done with a 2.93 bohr bond length. The
VO+ calculations used a 3.372 bohr bond length. The CrO+ calculations were

done with a 3.0 bohr bond length.

39

Table 1.7 - Total SCF electronic energies (in hartrees)

molecule
TiOH+

non"
TiOH*
TiOI-F

VOIF

VOH”
VOH*

Crorl“ (bent)
CrOH+
CrOH”

CrOH+
CrOI-F

TiO+
TiO+
TiO+

V0+

VO+

CrO+
Cr0+
CrO+

m
3A-
32-
'2
'2

m

2
(core) 7:,‘

2
x

(core) 7:

2
(core) 7:,

 

(core) ”2

I

2
(care) It,

2
(care) It,

 

2
(care) It,

2
(care) It,

2
(care) It,

 

2
(core) 7r,r

!(core) 71':

2
(core) 71',r

2
(care) It,

 

2
(care) It,

2
(care) It,

2
(care) It,

 

2
(care) It:

2
(core) 7:,

2
(care) It,

 

2
(care) 7:,

y ‘9' z
a; d;,_y. dj.)
7:; d3)
71; df=)

7:2 d' d',)

n2 d', , d‘, 41,)

y x—y :

.212. . 1,)

y x-y

It: 43' 61:, ally)

(core) 02 d; d:,_y, d; (11),)
d; 0112-..: “C: ally)
a; d}, d;,_y. egg)

2 I l I I
7:, 4s dirty, d .2 d”)

7!: d3, 43' (1;)

”i 11;)
7’: d1.)
1148‘)

7:2 d',_y, lily)

y x

7:2 d‘, , d:,>

y x-y

2 2
:5 dirty)

7:2 d'2 , d'2 d;

y x—y z
7:2 d‘ d', 4,)

I”:

2 2 l
a, dx,_y, or”)

energy
-923.67796

-923.67319

-923 .633 l 6
-923.57630

-1018.l4967

-1018.03075
-1018.02580

-1118.55670
-lll8.55520
-lll8.44184

-lll8.41322
-1118.26561

-923.04632
-922.97239
-922.96448

-lOl7.39671
-1017.36484
-1017.34811

-1017.79055
-lOl7.73586
~1017.72798

2. MCSCF calculations

The titanium and vanadium MCSCF calculations had both metal -
oxygen 7: bonding orbitals correlated. Correlating the more polarizable 7r
bonds while keeping an SCF description of the 0' bond was found to be
optimal for the scandiurn analogues by Tilson and Harrison [3]. Test
calculations with the sigma bond correlated as well, verified the validity of
this assmnption for these molecules. The MCSCF wave frmctions are
descrrbed in more detail in section H.

3. GVB calculations

Generalized valence bond calculations (GVB) use only two
determinants per pair of correlated electrons. These calculations were
explored as a reference space for a MRCI calculation in order to make the
larger calculations easier to do. Test calculations for the chromium
compounds showed a loss of 85 ml-I compared to the MCSCF calculations.
Since it was possible to obtain suficient disk space to use the MCSCF
reference space, the GVB calculations were not pursued due to their

questionable accrnacy.

41

4. MCSCF+1+2 calculations

Based on previous work, multi-reference configuration interaction
calculations are expected to be the most quantitatively acclnate description of
the geometry and energetics of these systems. A valence CI was done from
the MC SCF reference space. For the titanium and vanadium compounds, the
1t bonds and 1mpaired electrons were included in the active space. For the
chromium compormds the metal - oxygen sigma bond was also included in the
active space. This wave frmction gave the optimal acclnacy with a
reasonable amount of computational resources. These calculations are
discussed further in section H.

S. SCF+1+2 calculations

Inorganic chemists often think of metal main group bonds as being
entirely ionic coordinative bonds. If this were indeed an adequate
description, configuration interaction calculations from an SCF reference
space would be expected to perform well. This is due to the fact that the SCF
wave function usually describes ionic systems well.

For the grormd state of TiOH+, the SCF energy is 68.1 mH higher than
the MCSCF energy. An SCF+1+2 calculation with an optimized geometry

gave a total energy of -923.78739 H, only 6.2 mH higher than the MRCI.

42

This result indicates that this system is very sensitive to electron correlation
and that SCF+1+2 calculations give a reasonably good description of the
molecule. The MCSCF+1+2 calculations are used for the rest of this study
due to their increased accruacy, ability to dissociate to the adiabatic products
and quality of population analysis from the MCSCF orbitals.

H. Final wave function dggription

The metal atom basis contraction was (14sllp6d) -> (5s4p3d). The
basis set used was Dunning’s [l9] reoptimization of the Wachter [20] basis
set with the 6th d exponent recommended by Hay [21].

The oxygen basis was contracted as (11s7p1d) —-) (433p1d). The
Duijneveldt [22] 1157p was used with an additional diffuse d fimction with an
exponent of 0.85 .

The hydrogen contraction was (481p) —) (Zslp) using Huzinaga’s [23]
45 basis with Drmnings [24] 2s contraction coefficients. An additional p
exponent of 1.0 was added as suggested by Harrison [25].

Calculations at the MCSCF level correlated some of the bonding
orbitals against their antlbonding orbitals. Configuration interaction

calculations were performed using the MCSCF active space as the reference

43

space. This MC SCF+1 +2 method has been shown to give fairly accurate
bond lengths and bond energies with a minimmn amount of computation.
In the titanium compound MCSCF calculations, pairs of orbitals were

correlated as follows.

no" 2A’ ] (core)(7rr,7r32)2 ("ya”b2 51—)

3

TiOH+ A' I (care) (unit; )2 (zy,7r:,)2 61. 01)

The 0'1 refers the do (dzz) orbital. By keeping the 0' bonds in the core,
we are using an SCF description of these bonds even at the CI level.

The description of the vo+ 32' and von+ ‘2‘ differ from their titanium
analogs only in the presence of the additional unpaired electron in the active
space.

The difi‘erence in bonding in the chromium compounds required a
consistent correlation scheme to be constructed as

4

CrO+ II

 

(core) (030*)2 (Irx,7r;)2p§ 61.61 d;>

5

CrOH+ A |(core) (030")2 61 6: d; 61;.)

l

The TiOH+ MCSCF wave function contained 25 configuration state

functions (CSF). The no" MCSCF contained 17 CSF’s. The VOH” and

44

vo+ MCSCF’s contained 34 csr’s and 25 CSF’s respectively. The CrOH+
and CrO+ MCSCF wave fimctions contained 7 and 34 CSF ’s respectively.

The non" MCSCF+1+2 wave function contained 36,686
configuration state functions. The no+ MCSCF+1+2 contained 10,748
CSF’s. The VOH+ and vo“ Cl’s contained 64,734 CSF’s and 27,927 CSF’s
respectively. The CrOH+ and CrO+ MCSCF+1+2 wave functions contained
61,574 and 59,034 CSF’s respectively.

All calculations were carried out with the Columbus software suite
[26].

1. Conclusions

A number of conclusions about these molecules can be drawn fiom this
study.

The early transition metal mono oxide cations exihibit nearly the
shortest bonds known for oxygen with these metals. The oxide bonds are
very strong bonds which can reasonably be assigned to be triple bonds.

The mono hydroxide cations of Sc, Ti and V have linear grormd states.
The metal - oxygen bonds in these molecules are short compared to the range
of experimentally determined bond lengths. These bonds are characterized as

being triple bonds consisting of one electron pair bond and two dative bonds

45

in order to give the maximum possible interaction between the metal and the
OH ligand.

The most stable electronic states of chromium compormds are those
states which preserve the half filled 3d shell environment around the Cr atom.
This leads to electronic and geometric structures considerably different from
those of the corresponding vanadium compounds.

Contour plots and population analysis indicate very polar bonds
characterized typically by a 0.5 electron polarization. This characterization
makes these systems nearly perfectly centered between completely covalent

and completely ionic behavior.

LIST OF REFERENCES

10.

ll.

12.

13.

LIST OF REFERENCES

D. D. Ebbing “General Chemistry” Houghton Mifllin (1987)
J. F. Harrison, D. C. Yormg submitted to J. Phys. Chem.

J. L. Tilson, J. F. HarrisonJ. Phys. Chem. 9_5, 5097 (1991)

J. F. Harrison J. Phys. Chem. 20, 3313 (1986)

P. G. Jasien, W. J. Stevens Chem. Phys. Lett. 151, 72 (1988)

Y. Takahara, K. Yamaguchi, T. Fueno Chem. Phys. Lett. 1_5_8_, 95
(1989)

C. W. Bauschlicher, Jr., P. Maitre Theor. Chim. Acta. .9_0, 189 (1995)

R. Bergstrbm, S. Limell, L. Eriksson Int. J. Quantum Chem. 52, 427
(1996)

D. E. Clemmer, N. Aristov, P. B. Armentrout J. Phys. Chem. 21, 544
(1993)

T. F. Magnera, D. E. David, J. Michl J. Am. Chem. Soc. 111, 4100
(1989)

H. Kang, J. L. Beaucltamp J. Am. Chem. Soc. 198, 5663 & 7502
(1986)

C. E. Moore “Atomic Energy Levels ” National Bureau of Standards
(1949)

F. A Cotton ‘O‘Chemical Applications of Group Theory ” John Wiley &
Sons (1990)

14.

15.

l6.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

47

J. F. Harrison, K. L. Kunze “Organometallic Ion Chemistry” Kulwer,
89 (1996)

D. B. Lawson, J. F. Harrison J. Phys. Chem. 199, 6081 (1996)

F. A. Cotton, G. Wilkinson “Advanced Inorganic Chemistry ” John
Wiley & Sons (1988)

A. G. Orpen, L. Brammer, F. H. Allen, 0. Kennard, D. G. Watson, R.
Taylor J. Chem. Soc. Dalton Trans. S1 (1989)

B. H. Botch, T. H. Dunning, J. F. Harrison J. Chem. Phys. 15, 3466
(1981)

T. H. Dlmning private communication.
A. J. H. Wachter J. Chem Phys. _5_2, 1033 (1970)
P. J. HayJ. Chem. Phys. 66, 4377 (1977)

F. B. Duijneveldt IBM “Technical Research Report No. Rl-945 ” IBM
Research Laboratory, San Jose, CA (1971)

S. Huzinaga “Approximate Atomic Functions 11.” Research Report;
Division of Theoretical Chemistry, Department of Chemistry,
Univeristy of Alberta (1971)

T. H. Drmning J. Chem. Phys. 53, 2823 (1970)

M. Rivera, J. F. Harrison, A. Alvarado-Swaisgood J. Phys. Chem. 24,
6969 (1990)

R. Shepard, I. Shavitt, R. M. Pitzer, D. C. Comeau, M. Pepper, H.
Lischka, P. G. Szalay, R. Ahlrichs, F. B. Brown, J. G. Zhao. Intemat.
J. Quantum Chem. S2; 149 (1988)

Chapter II

LONE PAIR ELECTRONS ON THE WATER MOLECULE

A. Introduction

Lone pairs in water and related systems are frequently described as
being regions of enhanced electron density [1]. Indeed, the familiar "rabbit
ears" representation, also called equivalent orbitals (E0), epitomizes this idea
and allows the systemization and rationalization of many chemical properties

[2].

 

E0 SMO

The lone pair electrons on the water molecule
can be modeled as spectroscopic molecular
orbitals (SMO's) or equivalent orbitals (EO's).

 

 

 

48

49

Alternatively, the lone pairs can be viewed as two different orbitals
called spectroscopic molecular orbitals (SMO). This is the view given by the
symmetry adapted orbitals in most computational chemistry programs. This
model has in the past been cited for the existence of two peaks in the
photoelectron spectrrnn, but it is now recognized that both models predict two
peaks if examined closely [3].

In spite of the enormous appeal and success of these representations,
little is known about the physical nature of the lone pairs. Most discussions
invoke the molecular orbital model and focus on the nature of individual
orbitals [4,5]. Since these orbitals are not unique, these discussions are
somewhat subjective [6].

In this chapter, the physical manifestation of the lone pairs in an
isolated, gas phase 1-120 molecule is explored in terms of several aspects of
the wave function. First the limitations of orbital based models are discussed.

Then two invariants of the 1-120 wave frmction with well understood physical
interpretations are examined: the electron density, p(?) and the molecular
electrostatic potential (MEP) derived from this density. Finally, two
properties obtained by mathematical manipulations of the wave frmction are

discussed: density difl'erence plots and the Laplacian of the electron density.

50

The electron distrrbution to the "rear" of the oxygen (in the "lone pair"
region) is examined to see if there is any observable structure to support the
idea of the lone pairs as regions of enhanced electron density.

B. Visualization Technig us

This project utilizes a nmnber of visualization techniques. These
techniques are computer graphics extensions of traditional graphs. All of
these images are views of actual data, not artist’s renditions as are sometimes
encountered.

A typical graph of f(x) versus x can be viewed as having one
coordinate dimension, x (where x is being thought of as a physical distance),
and one data value, f(x), at each point in space (corresponding to each value
of x). This traditional plot of Ex) would be considered a one dimensional
graph in the language surrormding visualization techniques.

In contrast a plot of land elevation would be thought of as having two
coordinate dimensions (latitude and longitude) and one data value (the
elevation) at each point and is thus called a two dimensional data set. The
traditional way to plot two dimensional data sets is as a contour plot as shown
in Figure 2.1 and 2.2 . A more intuitively quantitative plot, called a mesh

plot, is obtained by using the data value as a displacement in the third

51
Figure 2.1 - Contour plot of the electron density in the plane of the HZO

molecule. The innermost contour corresponds to 0.142 e/A3 and subsequent

contorn's decrease 1miforrnly by 0.00570 e/A3.

 

 

 

 

 

52
Figure 2.2 - Contour plot of the electron density in the plane perpendicular to
the molecular plane and bisecting the HOH angle with the hydrogen atoms to
the left. The innermost contour corresponds to 0.142 e/A3 and subsequent

contours decrease uniformly by 0.005 70 e/A3.

 

 

 

 

 

53

Figure 2.3 - Mesh plot of electron density in the plane of the water molecule,
oriented with the hydrogens towards the viewer. The peak on the oxygen
atom has been truncated in order to make the hydrogens visible.

 

1“]
l‘~ [\IV
l\/\/

xI‘K‘I

5.1"

 

 

 

 

 

 

 

 

 

 

 

 

 

 

54

dimension. Figure 2.3 is a mesh plot of the electron density for the water
molecule in the plane of the molecule. A third visualization technique for two
dimensional data is to assign colors to data values in order to give a
colorization plot as shown in Figures 2.4 - 2.6, which represent high densities
as red and low densities as blue.

In the case of electron density, there is a density associated with every
point in three dimensional space. The three dimensional analogue to a
contour of electron density would be a che, called an isosru'face, as shown
in Figrne 2.7 .
C. Orbital Renrggngtions

Wave frmctions for many electron systems are often constructed from
one electron frmctions, called orbitals. For example, the Hartree-Fock wave
filnction is a single determinant of one electron orbitals. Orbital based
descriptions have proven to be useful for the derivation of many simplified
models such as VSEPR theory [2] and the Woodward-Homnan e rules [7].
However, orbitals are not unique.

Most ab initio and semiempirical programs use symmetry adapted
orbitals. These are the linear combinations of basis frmctions necessary to

give orbitals which transform according to the irreducible representations of

55

Figure 2.4 - Colorization of total electron density in the plane of the water
molecule

Water total electron density

 

o ‘ I.-
.. .' t 2
. ' ‘ I
’_ ...., ‘ e '

 

56

Figure 2.5 - Colorization of the total electron density at right angles to the
plane of the water molecule

Water total electron density

3

 

57

Figure 2.6 - Colorization of the total electron density both in the plane of and
perpendicular to the water molecule

Water total electron density

\

 

58

mvfiood

.ommmumwm.
um...“ .. .1 .u .

 

s84

.. ... . t... .

1. —. ...- .1... . .

a .15»... . . . . . .3... WA...
1&— v . . .. . C 3.3+. .
.c. . N c .. .. A.

. . . .

.

  

mgod

a. . .
n...a

...
|.. V r
a .. ,.
c.. .

r.
.7.

  

3.0.0

 

AncchammémcobuoE
mam—Bu 5.5020 :33 .833

vvvod

 

3 2.9m

59

the molecular point group. These orbitals reduce the cost of computation by
giving a block diagonal Hamiltonian matrix with one block for each
irreducrble representation [6].

Another set of useful orbitals are localized orbitals. Localized orbitals
may be constructed by finding the set of orbitals with the minimum
Coulombic repulsion between different orbitals. This results in having one
orbital for each pair of bonding, lone pair and core electrons. This intuitive
description is often invoked when discussing molecular structure [2].

There are two lone pair orbitals in the water molecule 1mder either of
these systems. The symmetry adapted orbital system gives an A. and a B1
orbital (SMO’s). The localized orbital system gives an identical orbital for
each of the lone pairs (EO’s). Two equivalent orbitals can be obtained
merely by taking positive and negative linear combinations of the symmetry
adapted orbitals.

Past discussions of lone pair electrons have described the lone pairs in
terms of these orbital descriptions [4,5]. Since orbitals are truly only
convenient ways of partitioning a wave fimction into subfirnctions, arguments
over which orbitals are correct are tantamormt to arguing over whether 4 is

2+2 or 3+1.

D. To lElectron Densi

The total electron density is an observable obtained from the wave
flmction [8]. The electron density is the number of electron charges per unit
volume (e/A3 in this study). It is experimentally observed via X-ray
crystallography, scanning microscopy techniques and electron scattering [9-
12]. To date, easily computed electron densities are more accurate than any
experimentally obtainable densities.

Many texts have traditionally described the lone pair electrons as
regions of high electron density. The electron density arising from a single
orbital may indeed match this description. However, the experimentally
observable total electron density is c0mputed as the sum of all of the orbital
densities rather than a single orbital density.

Figures 2.1 and 2.2 show two views of the electron density in the 1120

molecule. The innermost contour corresponds to 0.142 e/A3 and subsequent

contorus decrease uniformly by 0.005 70 e/A3. Clearly, the electron density is
rather featureless, especially in the lone pair region.
Figure 2.7 shows a sequence of three dimensional isosurfaces of the

density beginning at 0.148 e/A3 and decreasing to 0.00148 e/A3. We see that

contour levels 0.148 and 0.0739 e/A3 represent density essentially localized

61

on oxygen with the hydrogens beginning to contribute at 0.0592 e/A3. By the

time we get to the contour level 0.00148 e/A3 we have a rather structureless
electron distribution. The size of this last che is, at most, 3.23 A across.

Although the density is featureless, it is not spherical in the lone pair

region. Figure 2.8 shows the 0.0327 e/A3 isosurface with an inscribed sphere
centered about the O atom. The inscribed sphere is colored green and has
been increased in size 1mtil it just breaks through the isosurface. If the
isosru'face was equidistant from the oxygen, the inscribed sphere would break
the srn'face at all points simultaneously. Consequently the white areas are
firrther from the nucleus than the green areas. If we increase the radius of the
inscribed green sphere very slightly we see the green sphere completely
outside of the isosurface. The electron distribution is very nearly spherical
and we See no pronomrced variation in the electron density that reminds us of
the lone pairs. Note that the ragged line where the sphere meets the density
isosurface is due to the visualization software approximating a sphere by a
high order polygon. However, the slight asymmetry in the electron
distribution, as shown in Figure 2.8, is indeed consistent with the notion of

lone pairs oriented approximately tetrahedrally with respect to the OH bonds.

62

Figure 2.8 - The 0.0327 electron per cubic Angstrom isosurface with an
inscribed sphere (green) centered at oxygen. The radius of the inscribed
sphere is just large enough to break through the isosurface.

 

63

There is now a large body of computational, experimental and
axiomatic evidence to contradict the idea of lone pair electrons being regions
of high electron density. See chapter IV for fruther discussion of the total
electron density.

E. Molgular Electrostatic Potential
The molecular electrostatic potential (MEP), 41, is an observable,

intimately related to the electron density p(?) [13]. The MEP at any point in

space, R , arormd the HZO molecule is given by

 

 

MFMV w“ Zr
R - - _
¢()= I]r-1‘é|+ H le_R|

where, p(?) is the electron number density at r, RJ- is the vector to and Zj is
the charge on the jth nucleus. A point particle with an infinitesirnally small
charge, q, will have a potential energy q ¢(R) at the point R. There are
regions arormd the water molecule where AR) is positive and regions where

it is negative. Positive point charges will be attracted to the regions of
negative potential and negative charges to the regions of positive potential.

AR) was evaluated using the same electron density function as is

displayed in Figures 2.1 - 2.8 . Figures 2.9 - 2.13 show contours and

Figure 2.9 - Contours of the electrostatic potential in the plane of the HZO

molecule.

 

 

 

 

 

 

 

65

Figure 2.10 - Contours of the electrostatic potential in the plane

perpendicular to the plane of the 1-le molecule.

 

 

 

 

 

66

Figure 2.11 - Colorization of the electrostatic potential in the plane of the
water molecule.

Water electrostatic potential

 

 

. a... s. ;,r i. -

1:3 “xiv“

67

Figure 2.12 - Colorization of the electrostatic potential in the plane
perpendicular to the plane of the water molecule.

Water electrostatic potential

    

68

Figure 2.13 - Colorization of the electrostatic potential in the plane of and in
the plane perpendicular to the plane of the water molecule.

Water electrostatic potential

 

_w,.:¢l¢nf "Vi-a3 . .

4’."- A 1‘2! “the"
A i 4,... r ~
~. 1 1": ‘ ‘ ~'A'

9.

‘3‘” 7 , jam“?

.1!

       

69

colorizations of d) in the molecular plane and in the plane bisecting the HOH
angle. In Figures 2.11-2.13, the black region around the nuclei is a region of
positive electrostatic potential, which becomes red as large positive values
are encountered near each nucleus. The colored areas are regions of negative
electrostatic potential with blue being near zero and the largest magnitude
negative regions represented by a black patch inside of a red region.

There are several interesting observations to be made. First, we see
that the negative contours are clustered to the rear of the water molecule
away from the positive hydrogens. Secondly, there is considerable structure
intheMEPandwhilethereareminimainbothplanesthe globalmim'maare
in the plane bisecting the HOH angle. These two regions, symmetrically
located above and below the plane containing the nuclei, are 1.26 A from the
oxygen nucleus and subtend an angle of 98° with the oxygen.

Figure 2.14 shows the three dimensional nature of these regions of low
electrostatic potential. The contour values are given in energy units (eV),
which corresponds to the energy a tmit positive charge would have at these
points. The -1.90 eV contour encloses a volume which extends continuously

from one side of 1120 to the other. At -1.97 eV we see the regions beginning

to separate and at ~2.04 eV the minima in the MEP are symmetrically

70

Figure 2.14 - Three dimensional nature of regions of low electrostatic
potential.

Water electrostatic potential

 

—1.973 eV

   

—2.041 eV

71

arranged above and below the molecular plane. The absolute minimum in

¢(R) is approximately -2.27 eV at which value the enclosed volumes are

zero.

Figure 2.15 shows a positive value of the MEP in blue (+272 eV) and
a negative region in yellow (-2.041 eV). These minima define directions in
space not too dissimilar to those associated with the approximately
tetrahedrally located lone pairs often associated with water. This orientation
is also consistent with known hydrogen bonding in gas phase clusters and
liquid water as well as the crystal structure of ice.
F. Electron Densm' Diflerence

A density difference plot is a plot of the electron density in the
molecule minus the electron densities of the isolated atoms. Thus positive
regions represent a gain in electron density as might be expected for chemical
bonding regions and lone pair regions. Likewise negative regions represent a
net loss of electrons.

In order to obtain a density difference plot, a judgement must be made
as to how the individual atoms are oriented. For example, the two unpaired
electrons on an oxygen atom could be pointed towards the hydrogens of the

water molecule, or the atomic electron densities could be spherically

72

Figure 2.15 - The molecular electrostatic potential (MEP) in water. The blue
surface corresponds to a MEP of +27 .2 eV while the yellow surface
corresponds to a MEP of -2.04l eV.

 

73

averaged. The choice of orientation will afi‘ect the behavior seen in the
density diflerence making this method somewhat subjective.

Analysis of the density difference indicates an enhanced density
relative to the separated atoms in the bonding regions and a single lone pair
like region in the plane of the molecule as shown in Figure 2.16 . Although
this is consistent with the SMO description of lone pair electrons, it is a
mathematically derived quantity rather than a physical property of a single
isolated molecule.

G. La cian of E n

The Laplacian of the electron density, \72 p , shows the behavior of the

electron density relative to an exponential decay. A perfect exponential
decay, fix) = exp( -x ), will have all of it’s derivatives continuous and
monotonic giving no peaks in plots of derivatives of any order. If the electron
density is a nearly exponential decay with an additional shoulder, that
shoulder will appear as a peak in a plot of the second derivative. Thus there
will be a negative minimum in the Laplacian of the electron density where
there is the maximrmr positive deviation from exponential decay. Bader and
coworkers have utilized this property to interpret electron lone pairs as

maxima in plots of -V2p [14]. Note that this is a maximmn relative to an

74

Figure 2.16 - The 0.27 e/A3 isosurface of the density difference for water
shows regions of increased electron density relative to the separated atoms in
the plane of the molecule behind the oxygen atom (top lobe) and in the
bonding regions (bottom two lobes).

 

75

exponential decay, but is still a region of monotonically decaying electron
density.

Figm'e 2.17 shows that the lone pair regions in plots of the Laplacian of
electron density have a shape that is similar to the regions of negative
electrostatic potential. Note however, that these regions are 0.33 A from the
oxygen nucleus, significantly closer than the 1.26 A distance to the minimum
in the electrostatic potential (Figures 2.14 & 2.15). Like the density
difl'erence plot, the disadvantage of this technique is that it is mathematically
derived rather than a direct physical interpretation.

H. Computational Details

The electron density and the electrostatic potential have been
calculated from high quality ab initio wave fimctions at the self consistent
field (SCF) and configuration interaction (CI) levels [6]. The oxygen atom
was described by a Duijenveldt lls 7p basis set [15] contracted to 65 5p
along with three tmcontracted d primitives [16] with exponents of 2.25, 0.75
and 0.25. The 5 basis had the 7th, 8th, 10th and 11th primitives uncontracted.
The p basis had the 2nd, 3rd, 6th and 7th primitives rmcontracted. The
hydrogen atom is described by a Polarized Valence Double Zeta (PVDZ)
basis set [17] with the last two primitives tmcontracted and a single p fimction

with an exponent ofl.0 . The basis sets used are listed in Tables 2.1 and 2.2.

76

Figure 2.17 - The -3.973 isosurface of V2 p

 

77

Table 2.1 - Oxygen atom basis set

a contraction p contraction
W Eur-sisal firms-14 Ema
31195.6 ’ 0.000210 114.863 0.002266
4669.38 0.001628 26.8767 0.017192
1062.62 0.008450 8.32077 0.07534]
301.426 0.034191 2.97237 0.212775
98.5153 0.110311 0.423600 0.398523
35.4609 0.269493 0.150740 0.184703
13.6179 0.423545
5.38618 0.283038 p contraction
1.53873 0.027484 Exponent Qgircicnt
26.8767 1.000
a contraction
mm 2mm 9 contraction
31 195.6 0.000048 £529.19!!! goefi_c_ient
4669.38 0.000369 8.32077 1.000
1062.62 0.001936
301.426 0.007860 p contraction
98.5153 0.026719 13m ELM-1912!!
35.4609 0.070298 0.423600 1.000
13.6179 0.146011
5.38618 0.147898 p contraction
1.53873 -0.238701 m m
0.150740 1.000
a contraction
m We!
13.6179 1.000
a contraction d contraction
Man! W W
5.38618 1.000 2.25 1.0
s contraction d contraction
Liane-rm Meet W germ
0.605500 1.000 0.75 1.0
a contraction d contraction
cm mm mm
0.220540 1.000 0.25 1.0

78

Table 2.2 - Hydrogen atom basis set

s contraction
Elam we

13.01 0.019685
1.962 0.137977
0.4446 0.478148
0.1220 0.501240

8 contraction
Exmnent Coeficient
0.4446 1.000
s contraction
Exmnent oeficient
0.1220 1.000
p contraction

W Coeficient
1.000 1.000

The SCF total energy of -76.06160H and the CI total energy of
-76.34446H compare favorably with previous calculations [16] which predict,
at best, -76.06429 H for the SCF total energy and 46.35328 H for the CI
total energy. The optimized geometries used were a bond length of 0.955 A
and a bond angle of 105.1° for the CI wave fimction and a length of 0.942 A
and angle of 106.2° for the SCF wave fimction. The experimental geometry
is a bond length of 0.9587 A and an angle of 103.9° [18]. The reported

results were obtained with the CI wave fimction but we find essentially the

79

same results using the SCF function. The calculations were performed using
the Colrnnbus [19] and Gaussian 90 [20] programs and the subsequent
visualizations used the Advanced Visualization System (AVS) [21].
1. Conclusions

Although useful, orbital descriptions are inconclusive due to being non
tmique. The electron density, which has been cited as the physical
manifestation of lone pair electrons, does not have a maximum in the region
to the rear of the water molecule. Density difference plots and the Laplacian
of electron density show structm'e consistent with lone pair electrons, but do
not have a well defined physical interpretations in themselves. The
molecular electrostatic potential bears both a structure and physical
interpretation consistent with the definition of lone pair electrons. Subsequent
studies have formd the MEP to be consistent with the lone pair model in a
variety of other molecules [22].

These results suggest that it might be more physically accmate to

consider the lone pairs in H20 as a derived property of p(?) e.g., the MEP,
KR) , rather then as the electron density itself. While the concept of a lone

pair is extraordinarily valuable in interpreting, predicting and systematizing

many aspects of chemistry, its characterization, as regions of space with

80

significantly enhanced electron density is a considerable oversimplification.
Indeed, as we have shown, lone pairs (at least in H20) correspond to regions
of minimum electrostatic potential rather than regions of high electron

density.

LIST OF REFERENCES

10.

ll.

12.

13.

LIST OF REFERENCES

D. Eisenberg, W. Kauzmann “The Structure and Properties of Water”
Oxford University Press (1969)

See for example any text that discusses the VSEPR theory of molecular
geometry. A representative one being, P. W. Atkins, J. A. Beran
“General Chemistry” Scientific American Books (1992)

J. Sirnons J. Chem. Educ. 62, 522 (1992)

R. B. MartinJ. Chem. Educ. 66, 668 (1988)

M. J. Laing J. Chem. Educ. 64;, 124 (1987)

I. N. Levine “Quantum Chemistry” Prentice Hall, (1991)

J. March “Advanced Organic Chemist ” John Wiley & Sons (1992)

V. H. Smith in “Electron Distributions and the Chemical Bond” (Eds:
P. Coppens and M. B. Hall), Plentnn, 3 (1981)

S. Shrbata, F. Hirota, N. Kakuta, T. Muramatsu Int. J. Quant. Chem.
16, 281 (1980)

D. A. Skoog “Principles of Instrumental Analysis” Satmders College
Publishing (1985)

C. M. Lieber Chem. Eng. News M 28 (1994)
S. N. Magonov, M.-H. Whangbo Adv. Mat. 6, 355 (1994)
P. Politzer, J. S. Mm'ray in “Transactions of the American

Crystallographic Association” (Ed: R. Blessing) 26, 23 (1990)

81

14.

15.

16.

17.

18.

19.

20.

21.

22.

82

R. F. W. Bader “Atoms In Molecules; A Quantum Theory” Oxford
(1994)

F. B. van Duijneveldt “IBM Research Report RI 945” (1971)

R. R. Lucchese, M. P. Conrad, H. F. Schaefer, 111 J. Chem. Phys. 68,
5292 (1978)

T. H. Dtmning, Jr., P. J. Hay in “Methods in Molecular Electronic
Structure Theory” (Ed: H. F. Schaefer, III) Plenlnn, (197 7)

R. L. Cook, F. C. De Lucia, P. Helminger J. Mol. Spectrosc. _56, 62
(1974)

R. Shepard, I. Shavitt, R. M. Pitzer, D. C. Comeau, M. Pepper, H.
Lischka, P. G. Szalay, R. Ahlrichs, F. B. Brown, J. G. Zhao. Internat.
J. Quantum Chem. _S_2_2_, 149 (1988)

J. J. Frisch, M. Head-Gordon, G. W. Trucks, J. B. Foresman, H. B.
Schlegel, K. Raghavaohari, M. Robb, J. S. Binkley, C. Gonzales, D. J.
Defrees, D. J. Fox, R. A. Whiteside, R. Seeger, C. F. Melius, J. Baker,

R. L. Martin, L. R. Kahn, J. J. P. Stewart, S. Topiol, J. A. Pople
Gaussian 90, Revision J, Gaussian, Pittsburgh, (1990)

A VS Advanced Visualization Systems, Waltham, (1992)
L. Appleby “The Study of Lone Pair Electrons Using Electron Density

and Molecular Electrostatic Potential” (thesis) Michigan State
University (1996)

Chapter III

VARIATIONAL QUANTUM MONTE CARLO METHODS

A. Introduction

The wave mechanics formulation of quantum mechanics was originated
in 1926 by Erwin Schrodinger [1]. This theory gives a mathematical
description of the behavior of electrons in a molecule. Quantum chemical
computations typically accormt for the kinetic energy of the electrons, the
Coulombic attraction of electrons to the nuclei, the repulsion between the
electrons and the repulsion between the nuclei. The equation describing this

system is called the Schrodinger equation and is written, within the Bom-

     

Oppenheimer approximation, as
[T‘P=E‘-P (1)
elect _V2+ elect nucl_elect1
H=Z— "+2 2,1 +2 Z—"— (2)
i<j r'tj n<m rnm

Equation (1) is an eigenfimction equation in which E is the total energy of the

system and ‘1’ is the wave fimction. The wave fimction is a flmction of the

83

spatial coordinates of the nuclei and the space and spin coordinates of the
electrons. The wave function is a probability amplitude, meaning that the

probability/(multiple volume) of the particles being formd at a particular set of
coordinates is given by |‘I’|2 evaluated at those coordinates. Once ‘P is

known, all of the properties of the molecule can be obtained relatively easily.

Equation (2) defines the Hamiltonian operator, it. The first
summation, containing the Laplacian operator, is the kinetic energy operator
of the electrons. The second summation is the operator for the Coulombic
attraction of the electrons to the nuclei. The third summation is the electron -
electron repulsion operator. The final smnmation is the nuclear - nuclear
repulsion. This last term is not necessary for the computation of the wave
frmction, and is ofien omitted for the calculation, but is always added at the
end of the calculation.

The Schrbdinger equation has been solved exactly for the hydrogen
atom [1,2] and a mnnber of other one electron systems [3-5]. No exact
(analytical) solution has ever been found for any atom or molecule containing
two or more electrons. Since mathematics has not yet ofi‘ered any technique
amenable to solving these equations, theoretical chemists have applied a

number of methods for approximating the solution.

85

The accuracy of approximate calculations can be a matter of great
concern. The total energy which is computed is the energy required to
remove all of the electrons fi'om the molecule and separate the nuclei. For a
very small transition metal complex this energy may be on the order of
several thousand hartrees whereas the energy of the chemical bonds would be
around 0.1 hartrees and Van der Waals energies are on the order of 0.001
hartrees. Thus an error in the energy of 1% may be completely unacceptable
in describing chemical behavior. 2

Quantum chemisny approximation methods can be roughly classed into
forn categories; pertin'bation theory, variational methods, approximate
Hamiltonians and density functional theory [3,6]. Only the variational
method will be discussed here.

The calculus of variations is the branch of mathematics concerning
extreme mathematical cases such as largest, smallest, most, least, etc. This is
applied to quantmn chemistry by setting up an approximation scheme in
which the exact energy is a lower bound for the entire space of approximate
energies. Specifically, a fimctional form for the wave fimction is assmned,

then the expectation value of the energy, (E), is computed where

86

 

, ‘ \rt‘fiwt
<E>=<H>=<T$TJ>=W ‘3)

and the integration is over all space and spin coordinates. The expectation
value of the energy will always be greater than or equal to the exact energy as
long as the boundary conditions are satisfied. Thus the calculation with the
lower energy is closer to the exact solution. Typically the wave fimction
parameters are optimized by determining the values that give the lowest
energy.

A number of fimctional forms for approximate wave fimctions have
been proposed. One common technique is to use one electron fimctions
which have no explicit dependence upon the electron - electron distance. A
wave function formed fiom a single determinant of one electron fimctions
takes into account the average value of the electron - electron repulsion only.
This is the Self Consistent Field (SCF) method [3,7].

Hartree-Fock calculations can now be done routinely for small organic
molecules. The difference between the HF energy and the exact energy is on
the same order of magnitude as the energy of a chemical bond. HF
calculations can give reliable results where there is a systematic canceling of

errors (i.e. difl'erence in energy between two conformers). However, even a

87

very good HF calculation is not always accurate enough to give the desired
information.

A HF calculation can be used as the starting point for a more complex
calculation such as a Moller-Plesset perturbation theory method and other
multiple determinant methods such as a configuration interaction calculation
[3]. Unfortrmately these more complex calculations require far more
powerful computational resources than do HF calculations and may still not
be as accurate as desired.

Another fimctional form for the wave fimction was proposed by
Hylleraas [8] for two electron atoms. He constructed a wave fimction from
one electron ftmctions multiplied by a correlation ftmction. The correlation
fimction is a fimction of the electron - electron distance. In this scheme the
electron - electron repulsion is being accormted for more accurately than in
HF theory.

In the early 1960’s both styles of variational wave fimction were being
examined by many workers [9]. It was shown that the Hylleraas style
calculations were more accurate and converged faster than the HF method
and its extensions [10]. Unfortunately the Hylleraas wave fimctions resulted

in integrals too complex to evaluate analytically for any systems other than

atoms and homonuclear diatomic molecules [1 1]. Due to the generalizabflity
of the one electron calculations, explicitly correlated wave functions were all
but completely abandoned.

The reason for a renewed interest in explicitly correlated calculations is
the development of a numerical Monte Carlo technique for evaluating the
necessary integrals. The Monte Carlo method was first proposed by

Metropolis, Rosenbluth, Rosenbluth, Teller & Teller (casually referred to
asM(RT)2 due to the initials of the authors) in 1953 [12].

The early application of the Monte Carlo technique to molecular
quantum mechanics calculations was pursued by Conroy in the 1960’s [13-
15]. These early calculations could be done on only the smallest systems due
to the ineficiency of the computer algorithms and capability of the computers
available at that time.

Quantum Monte Carlo calculations really became feasrble for
molecular systems with the creation of a greatly improved computer algorithm
by Ceperley et. al. in 1977 [16].

All numerical integrals replace analytic integration with a smnmation of
finite elements. If the locations of these finite elements are chosen on a

lmiform grid in each dimension, the number of points to be sampled would be

89

the mnnber of points in one dimension raised to a power, which is the amber
of dimensions. In the case of quantlnn chemical calculations, the number of
dimensions is three times the number of electrons. This makes for a
prohibitively large nmnbcr of points.

In a Monte Carlo integration, the locations of the finite elements to sum
are chosen randomly to give an approximate integral. The error incrnred in
this statistical integration can be computed as a standard deviation between
successive “identical” integrations. As the number of points chosen becomes
very large, the answer converges to the exact integral (standard deviation
goes to zero).

B. Rand6m numbers

1. Importance

Random numbers are treated mathematically by the fields of statistics
and probability. These fields give mathematical descriptions of how random
systems will behave based on mathematical descriptions of what “random
nlnnbers” are. Results of quantum Monte Carlo calculations will not be
reliable lmless there is a method available to obtain random mnnbers which

behave in the ways which statistics assumes them to behave.

The best way to describe what “good” random nmnbers should be like

is to describe behavior that should and should not be exhibited. The behavior

that should be seen is, as follows:

A sequence of random numbers should give every possible number in its
range of values and on average each possrble nmnber should appear the
same number of times.

A mnnber should appear in a random sequence, on average, as often as
there are possible values retln'ned (e. g. the mnnber 7 in a sequence of
random nlnnbers in the range 0 - 32000 should appear with the frequency
of about once in every 32000 entries).

The space of random ntnnbers should be spanned evenly (e.g. in a
sequence of random numbers in the range 0 - 32000, a number in the

range 0 - 4000 should on average appear every eighth iteration).

The behavior that should NOT be seen is, as follows:

Random numbers should not repeat a set sequence.

Random numbers should not show an oscillation where every other
number is high or low.

Random mnnbers should not give long strings of high or low values.

Serial subsequences of numbers should not repeat.

91

0 Random numbers should not give the same sequence of numbers every
time they are used.

The last point is usually handled by starting from some initial nmnber,
called a seed, which is taken from the system clock. This ensures that no two
sequences start at the same number because it is never the exact same date
and time twice. The generation and testing of random nlnnbers will be
discussed in the following sections.

2. Generation

For the generation of random numbers, the seminal text is Knuth’s
work [17] which is old but still a classic work of computer science and
computational mathematics. This is a series of texts titled “The Art of
Computer Programming”. Random number generation and testing are
addressed in chapter 3 of volume 2, titled “Seminumerical Algorithms”.

One of the earliest computer systems obtained random mnnbers by
using an analog to digital converter to sample the noise on the incoming
power lines. While this is a rmique and inventive solution, it is not really the
best. Electric utilities use power generating equipment with a slight
fi'equency droop, which means that the frequency of AC electricity is slightly

lower during peak use times than dming light usage times. Thus even the

92

noise may be skewed to deliver statistically higher or lower values depending
upon the time of day, day of week, outdoor temperature, etc.

Most random number generators are pseudorandom number generators.
These are algorithms which generate numbers by nonlinear mathematical
equations in which the results of the previous iteration are used as the inputs
for the next iteration. These are not perfectly random because the same
sequence results a second time if the same initial starting point is used.
However, as long as a mechanism is available to prevent starting fi'om the
same point, these generators can generate sequences of numbers which obey
all of the conditions which may be applied to defining random numbers.

The most common type of random mnnber generators are called “linear
congruential generators”. These typically generate mnnbers between 0 and
m-l by using the results of one iteration, 1,, to calculate the next result, I,» 1,
according to the formula

1,-+1=(a1,°+c)modm (4)
where a and c are constants. These random number generators can give
extremely poor or extremely good results depending upon the choice of

values for the constants [17,18]. One fact that is always true is that a linear

93

congruential generator must eventually repeat its sequence of random
numbers.

For some applications, no additional acclnacy can be gained once the
random number sequence begins to repeat. An example of this is an
integration in which the random number is used as a location at which to
evaluate the flmction. Once the sequence begins to repeat, the space has been
spanned as thoroughly as possible with the algorithm.

In other applications, the random number represents a relative
displacement fi'om the current location of a particle. In this formulation, the
results are not being repeated unless the random numbers begin their repeat
cycle when the particles are exactly at their original locations, as is the case
for variational quantum Monte Carlo calculations.

A second popular random mnnber generation scheme is based on a
Fibinacci sequence. The most well known Fibinacci generator is the ran3
generator given in “Nlnnerical Recipes” [19]. It uses a Fibinacci series which
is lagged by 55 terms. This generator returns floating point values in the
range 0 to 1 exclusive, (0,1) . This gives excellent results for many types of
calculations without repeating its sequence. The disadvantage of using a

Fibinacci generator is that the results can be extremely poor if the process

being simulated can be mathematically mapped onto a F ibinacci sequence
[20].

3. Testing

There are two classes of tests for random number generators. The first
is a mathematical analysis of the equations used to generate the mnnbers.
This is discussed by Knuth [17] but will not be pursued in this discussion.
The second class of tests is an analysis of the sequence of numbers produced
by the random number generator.

All random number analysis methods are based upon testing for
confidence intervals. For example, the results of a set of identical trials
should be within one standard deviation of the analytic answer in 68% of the
trials, within two standard deviations in 95% of the trials and within three
standard deviations in 99.73% of the trials. An alternative confidence interval
test utilizes the chi squared intervals. The standard deviation and chi squared
tests yield mathematically equivalent information. The choice of a confidence
interval test is usually made based on mathematical expediency. The
standard deviation method was chosen for this work because it is the most

commonly used in chemical applications.

95

The random number generators tested were the ran3 generator and
ones built into the GNU C++ and Silicon Graphics C++ compilers.
Traditionally, the random number generators provided with compilers were
very poor linear congruential generators. These were created in the
anticipation of needing random elements in video games and such but never
intended to give statistically valid results. All three random number
generators gave excellent results. The ran3 generator was finally chosen for
the sake of future portability of the software.

Therandomnmnberteststhaterdstareeachaway ofdirectly testing
foroneofthe properties ofrandomnmnbersetsthatarelistedinthe section
on the importance of random numbers. A comprehensive listing of well
rmderstood tests is presented by Knuth [17].

Some of these test are merely alternative ways of verifying the same
property. When a well known random number generator, such as ran3, is
used there have already been documented tests of desired properties. It is
however still the case that researchers concerned with statistical results

should have an understanding of random number testing and verify the tests

appropriate to their application.

Eleven of the tests presented by Knuth were examined. Those that
were not redlmdant were performed in depth. Additional spot checks of
redlmdant properties were done.

The single most important test is to verify the results for a known case
with the algorithm and code that will be used for experimental cases. This is
illustrated by the ran3 generator which gives excellent results for many
applications, but gives very poor results when applied to Ising lattice
calculations [20]. The reliability of the ran3 generator was verified
extensively for use with variational quantum Monte Carlo calculations using
both Gaussian type (GTO) basis sets and Slater type (STO) basis sets. This
was done by nmning calculations which can be done without the use of

random number generators to verify the results.

97

C. Theory
Quantum Monte Carlo (QMC) methods are broken into three

categories; variational quantum Monte Carlo (V QMC), difi'usion quantlml
Monte Carlo (DMC) and Green’s fimction quantlmr Monte Carlo (GFMC).
Only the VQMC method uses an analytic wave function and obeys the
variation principle [21].

DMC calculations use a numerically defined wave fimction rather than
integrating an analytic flmction. This wave function is approximate because
exact bolmdary conditions, in the form of nodal surfaces, are not known. The
energy expression is expressed as a diffusion equation in imaginary time.

The Green’s frmction Monte Carlo technique uses an approximate
Green’s flmction integration as well as a numerical wave fimction.

For the pmposes of this study, there are strong advantages to using
approximate analytic wave functions rather than numerical wave frmctions.
Specifically, it is possrhle to examine the distribution of electrons within the
wave fimction by examining the orbital coefficients. This gives a means to
analyze bonding, hybridization and population analysis. Also the variational

principle gives a measrn‘e of the suitability of a wave flmction for describing a

molecular system. The rest of this chapter will focus on the variational
quantum Monte Carlo technique only.

Expectation values are computed as an integral over the wave flmction
(equation 3). An alternative form of this integral is to use the probability
distribution (the square of the wave flmction) instead of using the wave
fimction to perform the integration. Thus the expectation value of energy,

(E), would be computed as follows

<E>= If<x>e1<xldx (5)

f (x) = I‘I’I2 is the probability distribution (6)
el(x) = £3: is called the local energy (7)

Thus f (1:) evaluated for a specific location of the electron, x, gives the
probability of finding the electron at that location. Likewise el(x) gives the
contribution to the total energy of the system for the electron at a specific
location. This method is entirely generalizable to a system of N electrons
giving a calculation in 3N dimensions.

The advance that was presented by Menopolis, Rosenbluth,

Rosenbluth, Teller and Teller [12] was a numerical form of the above

integration. In this form, the continuous integral is replaced by a discrete

smnmationoftheform
1 N
E =_ I 8

where el,, is the expression from equation (7 ) evaluated for the electron
positions on iteration n. The choice of particle positions is made in
accordance with the probability distribution, f (x) , giving a numerical
integration. The Metropolis algorithm, which ensm'es that the probability
distribution is sampled, will be presented in the following section on the
computer implementation of these techniques.

The advantage of this procedure is that the approximate value
converges asymptotically to the analytic integral as the number of iterations,
N, becomes very large

limN-m EN = (E) (9)
Thus the desired computation time could be chosen and the answer would be
obtained along with a standard deviation to indicate how accmately the
approximation was converged. Alternatively, a desired accuracy could be
specified and the only barrier to obtaining that accuracy would be the

availability of the necessary computational time. Note that this discussion is

100

addressing only the acclnacy of integration, not the accuracy of a specific
wave fimction for describing a specific chemical system.

The question then becomes one of “What form for the wave function
will give an accurate result?”

This method could be used to do 3 HF calculation with Gaussian type
orbitals (GTO) of the form exp(-C; r2 ). Indeed this is done to verify the
correct fimctioning of the computer program. However, this is not the most
efficient way to do such calculations since HF GTO wave fimctions can be
integrated analytically.

An improvement over the acclnacy of GTO orbitals can be gained by
the use of Slater type orbitals (STO) of the form exp(-Q r ). The exact
solution to the hydrogen atom is a STO flmction and STO primitives are
known to give a more accmate representation of molecular charge densities.
Also it takes fewer STO flmctions than GTO frmctions to achieve the same
accruacy of representation of the charge density. The disadvantage of STO
fimctions relative to GTO flmctions is the necessity of doing numerical
integrations. Since QMC is, by definition, a numerical calculation, there is no
disadvantage to using STO fimctions for QMC calculations. Indeed, the need

for fewer flmctions increases the speed of the computation making STO

10]

functions the better choice for QMC calculations. Almost all VQMC
calculations are done with STO primitive fimctions.

The use of STD primitives alone does not result in a sufiicient increase
in acclnacy to choose QMC calculations over HF GTO calculations. There
are numerical techniques for doing STO HF integrals more quickly than is
possible with the QMC technique [22-25]. Also, the use of STD primitives
does not correct for the inacctnacies in the HF method itself.

’ The real advantage of QMC techniques is that the electron correlation,
which is lacking in HF calculations, can be put into the calculation at little
additional computational expense. This is usually done by using a wave
flmction of the form

i=1)“ D” 11%) (10)

i< j
Where D“ and D’9 are determinants of alpha and beta electrons respectively.
The use of this two determinant fimction gives the mathematical convenience
of canceling out the spin terms in the local energy flmction. This does not
introduce any approximation into the energy computation as long as the
Hamiltonian does not contain any terms dependent upon the electron spin

states (i.e. spin-orbit coupling).

102

The factor f(r,-,-) in the above equation is called a correlation fimction.
This fimction gives the electrons a lower probability of being very close
together than far apart thus accormting for the electron correlation.
Correlation fimctions are always a fimction of the electron - electron distance.
Correlation fimctions, which are a fimction of electron - nucleus distances
also, have been used for some calculations (see the “literature survey”
section). There is no proof that the correlation flmction must be a flmction of
the electron - nucleus distance. Indeed, correlation fimctions which depend
upon the electron - electron distance only are readily generalizable to any
arbitrary molecular system and have seen the most wide spread use.

There are five known conditions which a correlation function must
obey. Calculations can be done with correlation fimctions which do not obey
these conditions as long as the bormdary conditions on the total wave flmction
are not violated. However, only correlation functions which obey the
following conditions have been found to be generalizable to all molecular
systems and give consistently reasonable results [21].

1. The value of the correlation flmction at zero electron - electron distance

should be finite [26]. Thus «0) = c.

103
2. The value of the correlation flmction at infinite electron - electron distance

should be finite [26]. Thus {(00) = c...

3. The slope of the correlation fimction at zero inter-electron distance should

 

be finite [27]. Thus d is finite, where d=[df:l Once a
r..
'1 rv=0

normalization scheme is defined, the value that d should have for an exact
representation can be determined [28].

4. The correlation fimction must be continuous everywhere [29].

5. The correlation fimction should be monotonically increasing from on to c...
Thus given two electron - electron distances, r, and r; such that 0 < r, < r2
< co the values of the correlation functions evaluated at these points will
obey0<co<f(r,)<f(r2)<cw<oo. f

These constraints give the correlation flmction the shape shown in

Figure 3.1 . Frmctions obeying these conditions should yield reasonable

results. However, the accuracy of the calculation will still be very dependent

upon subtle variations in the slope and cln'vature of the correlation flmction as
well as the ratio of co to co
As mentioned above, QMC calculations do a statistical integration with

an accuracy represented as a variance or standard deviation. This variance is

104

Figure 3. 1 - Correlation function fir”)

 

 

L6 1’

1.4 ,_

1.2 4

0.8 T

3..

0.6 4~

0.4 «r

0.2 4-

ad
«N
ha
o."
md
in
MN
Nd
3a..

m.—
w.—
5..
o.—
n.—
v.—
m.—
N.—
fi—

06
ad
ho
v.9
vie
Yo
Md
N6
_6

l05

intrinsically linked to the speed and accruacy of the calculation. Thus it is
useful to examine the behavior of the variance.

There are two variances associated with a variational QMC
calculation. The comparison of several identical calculations gives a variance

in the total energy denoted as var{EN}. There is also a variance in the local
energy denoted var {eln }. In theory these two variances are linked by the

relationship [30]
var{EN} =-]lV-var{el,,} (11)

In practice, even though a calculation may span the entire space of
electron coordinates evenly it may not converge at the same rate in all regions
of this space. Also, the contribution to the total energy is different from one
region of space to another. Thus there is an efi‘ective relationship between

the two variances of
C
var{EN} = N var{el,,} (12)

where C is a constant.
For a given algorithm, the variance in the local energy will converge
fairly quickly and stay at a constant value. The variance in the total energy

will continually decrease as the mnnber of iterations, N, increases. This is a

106

usefirl result because it makes it possible to predict the variance of results for

a given number of iterations.

Both C and var {eln} are affected by the choice of algorithm. This

points to the fact that the calculation speed could be improved by devising
more efficient computer algorithms. Specific algorithmic concerns will be
addressed in the next section.

An implication of the above discussion is that a reliable measrn‘e of

efficiency should be based on var{EN} in order to measrne the overall

efficiency of the algorithm instead of one component of the calculation.
Improved algorithmic techniques typically achieve a decrease in

var{EN} at the expense of an increased CPU time, T. This is only an
overall improvement if it decreases the product which we have proposed
T var{EN} = computational work = CW (13)
Note that this does not give units of work. It is rather something that is best
minimized to obtained the greatest efficiency.
This does not, however, give any indication of the accrnacy of results.
The best way to choose among a number of methods is to list the methods

with acceptable accrnacy then select the one with a minimrnn value of CW.

107

Note that acceptable acclnacy will still be dependent upon the chemical
system and the type of information desired.
The energy computed by a VQMC calculation is a variational upper
bound to the exact energy. Thus it satisfies the condition
E, 2 E (14)
where the equality is true in the case of an exact wave fimction. A capability
of VQMC calculations which is not present in conventional ab initio
calculations is the ability to get a variational lower bound to the exact energy
via the relationship [21]
E, +vn{e1,} 2 E, 2 E 2 E, -m{et,} (15)
Expectation values of electron positions can be computed easily
because of the way in which the Metropolis algorithm samples the probability
distribution. For example, the expectation value of electron - electron

distance, r,-,-, would be computed as
1 N
«kW—gt,- (16)

Likewise, the electron density in a region of space could be obtained
by finding the average number of electrons in that region at any given

iteration of the simulation.

108

The computation of a Fermi contact density is somewhat more difficult
since it requires knowing the electron density at one exact point rather than in
a region of space. We have computed Fermi contact densities by computing
the electron density within two spherical regions arolmd the nucleus then
fitting these densities to an exponential in order to extrapolate to an
infinitesimal point on the nucleus. I

Other properties, such as excited state energies and spin interactions
have been examined using Monte Carlo methods. For a discussion of these

types ofcomputations the reader is directed to reference [21].

109

D. Computer implementation

The usefirlness of a computational chemistry technique depends upon
both its accmacy and eficiency. Calculations which are too inaccurate to
predict experimental results are of little practical use. Extremely accurate
calculations which can only be done on very small molecules using super
computers are also of only limited usefulness. Quantum Monte Carlo
calculations tend to fall in the latter category.

An enormous amormt of work has been done to find ways to do very
acctnate quantum chemistry calculations. Of all of the known methods for
approximating the solution to the Schrodinger equation, the quantum Monte
Carlo methods are the most accurate known to date. Anderson and
coworkers have done QMC calculations with an acclnacy exceeding that of
fill] CI calculations on the same molecular systems [31]. However, these
calculations required access to massively parallel computer facilities.

In light of these previous studies, this study will focus on making
improvements to the efficiency of calculations rather than on improving the
acctnacy. Emciency can be related to three areas; disk usage, memory
requirements and CPU usage. Disk usage is completely negligible for QMC

methods. Memory requirements are also far less than for conventional ab

llo
initio calculations. CPU usage is a major concern for QMC methods and will

be the focus here.

The most important measrn'e of CPU usage is time complexity. Time
complexity is the way in which resource usage scales as the calculation size
increases. For example, Hartree-Fock calculations have a time complexity of
N4 where N is the number of electrons. This means that a calculation with
twice as many electrons will take sixteen (24) times as long to am. Table 3.1
lists time complexities of a number of common computational chemistry

techniques.

Table 3.1 - Time complexity of computational chemistry techniques

 

t hnigue time complexig
Semiempiricals N2
VQMC (Ceperly) N3

Density fimctional theory N
Hartree-Fock NA
Configlnation interaction N3

VQMC (original) Nl

111
The time complexity listed is by definition the scaling obtained for

large values of N. A more in depth analysis may list a time complexity as
“aN’ + or? + cN + d”. The N3 behavior would be dominant for large values
of N, but for smaller calculations the other terms or even d, the over head,
may be the behavior observed. Thus the calculation with the smallest time
complexity may not necessarily be the fastest unless a large number of
electrons are being modeled.

The N3 time complexity predicted for VQMC applies to varying the
number of atoms in a molecule. As heavier atoms are computed, there are
difficulties sampling both the core and valence regions eficiently. This leads
to a time complexity of Z’" for varying the nuclear charge while maintaining a
consistent acclnacy. Valence only calculations for heavy elements have an
effective nuclear charge, Zea, equal to the mnnber of valence electrons
modeled. These calculations are predicted to have a time complexity of 2:;
[21].

Implementation of VQMC as written requires the evaluation of the
local energy expression (equation 7). This requires the evaluation of
determinant expansions consisting of N! terms. Such calculations were done

for very small chemical systems by Conroy and coworkers in the 1960's [13-

112

15] but the CPU time requirements for large molecules are prohibitive even
for today’s computers. An improved algorithm was formulated by Ceperely
and coworkers in 1977 [16].

The Metropolis Monte Carlo method (also called the random walk
method) is a technique for the statistical evaluation of a function by sampling
the probability distribution. In this case it is being used as a numerical
integrator. Monte Carlo methods are often the most eflicient method for
doing numerical integrations with a large number of dimensions. In this case
the number of dimensions is three times the mnnber of electrons. The

Metropolis method, as applied to VQMC, is as follows:

113

H

. Select a trial wave function.
2. Pick a random location for each electron (x).
3. Determine the electron - electron and electron - nucleus distances (r).

4. Calculate the local energy.

e=——- (17)

Ur

. Repeat steps 5.a through 5.e for N iterations.

5 .a. Randomly pick a new location (xm) for one of the electrons.
5 .b. Compute inter particle distances (rm).

5 .c. Compute the acceptance ratio (q).

__ Prom] (1 8)

 

5.d. Get a random number (rand) between 0 and 1.

If q > rand then keep the new positions, r = rm, x = x

new'

5.e. Evaluate and add in the local energy.

19‘?
e-e+—‘p— (19)

6. Divide by the number of iterations.

E. = {7 (20)

114
A variation on this is the use of a sampling frequency. This means that

step 5 .e is only performed on every mth iteration and the result of step 6 is
multiplied by m. This is an attempt to avoid over sampling the region along
the path of the particles Monte Carlo walk relative to other configurations.
This is of the most concern when very small step sizes are used as in DMC
calculations.

Importance sampling techniques replace steps 5.3 c & d with equations
utilizing a knowledge of the fimctional form of the probability distribution.
These techniques will be discussed later in this section.

All of the significant CPU time in the calculation is in the loop in
step 5. Any improvement in the speed of the steps in this loop will be an
improvement in the over all speed of the calculation. Any additional work
that can be done out side ofthis loop will be more than made up for by
improvements in the CPU time.

As mentioned previously, the algorithm on the previous page is not
practical due to a N! time complexity. In 1977 Ceperley et. a1. created a
greatly improved algorithm in terms of the inverse Slater matrix [16]. . The
Slater matrix, D, is the matrix of occupied orbitals evaluated at electron

positions as follows

115
Dij = (Pj(xi) (21)

The inverse Slater matrix is the matrix of cofactors of the Slater matrix

divided by the determinant of the Slater matrix, which is the wave flmction.

.1 _ Cii
n _ 1F (22)

Once the inverse Slater matrix is known, the acceptance ratio and
evaluations of terms in the local energy flmction become a sum down a

column of the inverse Slater matrix according to the following formulas

 

 

 

'73,) an on
V??? = Evil/50.19;! (r) (24)
$552, {$3142 Bit) (25)
31%,): Evian-(0920) (26)

Further discussion of equations 23, 24 and 26 is in appendix B] of reference
[21]. Equation 25 was not found in the literature.
This is implemented by moving one electron at a time, so only one

column of the Slater matrix must be recomputed on any iteration if the move

is accepted.

no
The most time consuming step of this algorithm is taking the inverse of

the Slater matrix when the move is accepted. This step scales as N3, however
each of these Operations are very fast. In tests of up to 21 electrons, the N2
behavior of computing the local energy and gradients was still the dominant
factor in determining the calculation time (see the results and discussion
section).

Table 3 .2 lists the time complexity of selected steps in the VQMC

calculation for the algorithm described here

Table 3.2 - Time complexity of VQMC steps

 

opfltion time complexity
Matrix inversion N3

particle - particle distances N‘2

Energy expression MN2
Random number calls 3
Parameter gradients P

N is the mnnber of electrons

M is the number of occupied spatial orbitals

P is the number of parameters

117

These time complexities are seen for a constant number of iterations in the
Monte Carlo procedlue. The CPU time varies linearly with the number of
Monte Carlo iterations, which is the single biggest factor in determining the
total computation time.

In DMC and GFMC calculations, many virtual particles are used, far
more than the number of electrons. For these calculations a more complex
procedure for updating the inverse Slater matrix increases the efficiency of
the calculation as described in appendix B2 of reference [21]. For VQMC
calculations with a very large number of electrons, this method would be
advantageous.

In general the variation of parameters is done by determining the

derivative of energy with respect to the parameter by [16]
g:(r(x)el(r))-(r(x))(ez(r)) (27)
I’(x)=2%ln(‘l’(x)) (28)

This is usually used in a gradient based optimization method such as a
steepest descent method. Improved performance can be obtained by saving

the gradient from the previous iteration and using the two gradients to

113
compute a cruvatrue for use in a quasi-Newton formulation. Non gradient

based methods, such as a simplex optimization can be used for a very small
number of parameters but are not as efficient as the gradient based methods
[19]. An adaptation of Pulay's direct inversion of the iterative subspace
(DIIS) method has not yet been implemented for QMC calculations [30].

AnimportantfactorinthespeedofaQMCprogramisthe
preprocessing of needed mnnbers. The primitive and contraction coefficient
normalizations can be done before starting the iterative Monte Carlo
procedlne. Likewise, intermediate steps in the computation of gradients and
the local energy expression are saved so that only the column corresponding
to the electron that was moved need be recomputed. These types of
algorithmic changes can increase the computation speed by orders of
magnitude. Specific examples will be given in the next section.

Another valuable algorithmic tool in creating eficient programs is the
use of difi'erent sampling techniques. In the simplest Metropolis sampling
technique the new location for an electron is chosen randomly within a small
cube around the present location of that electron. The acceptance / rejection
test ensrues that this gives a proper sampling of the probability distribution.

This random choice of the electron movement does not take into account any

119
knowledge of where there will tend to be a high probability of finding the

electrons.

There are a number of techniques, called importance sampling
techniques, which do take into account the nature of the probability
distribution when choosing new electron locations. Some of these introduce
additional approximations into the calculations while others are still giving a
numerical integration which must converge the exact integral.

Fokker-Plank importance sampling is an approximate technique which
does not utilize a rejection step [21]. This method is an adaptation of a
statistical mechanics technique. The new location for the electron being

moved is chosen according to the formula

xm = x + step @312: g + step_ size gasdev (29)

 

where gasdev is a random number in a Gaussian distribution with a mean of
zero and a variance of one. Tests of this method, presented in the next
section, showed it to be too approximate to be useful for chemical
computations on even the smallest molecules.

Importance sampling with rejection, also called Metropolis importance
sampling, is an importance sampling technique which does not introduce any

additional approximation into the calculation [21]. This technique uses the

120
Fokker-Plank equation to move the electrons and adds a rejection step. The

rejection step serves to maintain an exact sampling of the probability
distribution for finite step sizes.

When modeling transition metals, an additional problem arises which is
negligible for light elements. There is a significant difl'erence in the behavior
of core versus valence electrons as evidenced by the large difference in the
electron density between these regions. This creates dimculties in choosing a
reasonable step size for the calculation. A step size which is very small
would be reasonable in the inner core region where the electrons are densely
packed. However, a small step size would take an extremely long time to
sample the entire valence region. Likewise, a step size large enough to
efficiently sample the valence would often overshoot the high probability
region in the core thus requiring more iterations to get a reasonable sampling
of the core region. A number of variable step size techniques have been
proposed for application to heavy elements.

Variable step size teclmiques often introduce an additional error into
the calculation referred to as a step size bias. This is a statistical error at the

point where an electron is seeing a change in the step size, which is due to the

121
electron making a move for which there is zero probability of making the

reverse move as illustrated in Figure 3.2 .

The split 1 technique was introduced by Belohorec, Rothstein and
Vrbik in a paper examining CuH vrbrational properties [33]. This technique
uses four diflerent step sizes for difierent electrons as follows

Is 2s 2p 0.016 1

3s 3p 0.04 r
3d 0.12 1.’
4s 1

The assignment of electrons to orbitals was made based on the ordering of the
electrons from the nucleus. Thus the ten electrons closest to the nucleus
would be moved with a 0.016 I step size, etc.

The above technique incorporates a significant statistical error at the
point where two electrons are interchanged between shells and no error other
wise. We have proposed a method dubbed "variable t" utilizing the opposite
philosophy of a very small error at almost every iteration. This method
chooses the step size based on the distance fi'om the nucleus according to

step_size = r( 1.0 - 0.99 exp ( -L r) ) (30)

122

Figure 3.2 - Step size bias

A particle moving (the arrow) from a region with a large step size (the large
square) to a region with a small step size (the small square) has zero
probability of returning to it’s original position. This introduces a statistical

error, called a step size bias.

 

 

 

 

123

L=+ (31)

Where L is chosen to give the single STO that reproduces (r). In practice the
accumulation of very small errors using this technique was larger than the
sum of less frequent larger errors given by the split 1: method.

Based on this information and experience with transition metal
calculations, we proposed a modification of the split 1: method called
"split 1' b". This method varies the step size by one order of magnitude rather
than the two orders of magnitude used in the split 1: method. It also classifies '

the electrons into three groups rather than four according to

ls 25 2p 0.1 1
3s 3p 0.3 r
3d 45 I

As shown in the next section, this proved to be the most efficient method
tested for all electron VQMC calculations on transition metals.

If a core - valence separation is to be utilized to treat heavy atoms, an
obvious step would be to eliminate the core electrons completely by the use

of core potentials.

124
Core potentials use a Gaussian expansion to describe the effective

potential of the nucleus and core region on the valence electrons.
Orthogonality between the core and valence is accormted for by the use of
projection operators.

Hay and Wadt's [34] averaged relativistic efi‘ective core potentials
(AREP) were chosen based on previous experience with HF and correlated
calculations. For Sc - Zn their core potentials are nonrelativistic potentials
fitted to numerical Hartree-Fock wave flmctions. For Y - Hg they report
relativistic potentials fitted to Cowan-Grifin calculations. This distinction is
noted because calculationsusingtheHayandWadtsetareofien seeninthe
literature which refer to the core potentials generically as relativistic
potentials even when first row transition metals are being studied. These are
called averaged potentials because they use only an orbital angular
momentlnn projection [3 = |l)(l| rather than the l and m projectors P = |lm)(lm|
[3 5].

One of the most promising ways to improve the performance of
computer applications is to use a parallel computer or a collection of work
stations to break the calculation into pieces and let difl'erent processors do

part of the computation in parallel. QMC techniques are probably some of

125
the most readily parallelized calculations known. These techniques can be

parallelized in a very course grained way in order to spread the calculation
over separate machines without being hindered by the speed of
communication over a network link.

The most cornse grained way to parallelize a QMC algorithm is at the
level at which the calculation is repeated several times in order to get average
values and their standard deviations. Each of these repeat calculations can be
done on a separate processor. Once each processor has done a complete
calculation, the results are used to determine the actual values (average) and
accuracy of those values (standard deviation). This scheme is useful when a
relatively small number of very powerful processors is available.

In order to utilize a very large number of processors, the iterations in
the Monte Carlo procedure can be broken up and passed to different
processors. The intermediate summations can then be totaled to get the

calculation result.

126
E. Results and Discussion

The reason for lmdertaking a study of the VQMC method is to evaluate
it's potential for high acctnacy calculations on small molecules containing
transition metals. These compounds are extremely sensitive to the effects of
electron correlation and must be described by a reasonably well correlated
method in order to get even qualitatively correct results [36]. These highly
correlated calculations can become extremely computationally intensive
(Table 3 .2). The previous work showing the high acclnacy obtainable by
QMC calculations and the more favorable time complexity indicate that these
methods may provide an efficient way of doing very accurate calculations.

In order to 1mderstand the transition metal results it is important to
1mde the behavior of the underlying VQMC method. The hydrogen
atom will be used to test the VQMC method by testing the effect of changing
various aspects of the simulation. Since the hydrogen atom wave flmction
used can be integrated analytically, this also serves as the single most
important test of the suitability of the ran3 random number generator for
VQMC applications. The helium atom is examined as the simplest system

where a correlation flmction must be used. Finally a mnnber of So atom

127
calculations are presented in order to test the application to heavy atoms. The

extension to molecular systems will not be examined in this work.
1. Hydrogen atom tests

Note that VQMC gives the correct exact energy for the H atom when

A

the exact wave flmction is used. The local energy expression is £31—

(Equation 7 ). When the exact wave flmction is used, the numerator becomes
-0.5 times the wave flmction giving a local energy of -0.5 H, the exact total
energy for the H atom. If many iterations are done, this constant is averaged
with itself which of course gives an average of ~0.5 H. The same result
would be obtained for any molecule if the exact solution were known.

The previous paragraph seems trivial at first, but it illustrates several
points which should be mentioned for comparison to other computational
methods. First, QMC methods are ab initio methods meaning that no
experimental data are used in the formulation of the method. Second, the
method will admit an exact solution if the variational frmction spans the space
of the exact fimction. (Density Flmctional theory has not yet yielded an exact
solution to the hydrogen atom.)

The following is an approximate solution to the hydrogen atom

consisting of a single Gaussian type orbital.

128
‘1! = N(ot) exp( -or r2 ) (32)

N(ot) = ( 2 or / 1r )3“ (33)
The analytic integration gives approximately
(E) = 1.5 or - 1.59576912 or"2 (34)

The local energy for this wave fimction is

A

—’fPE=3ot-2o2r2-1/r (35)

which integrated by VQMC should yield the same result as Equation (34).
Choosing the value of or to be 1.0 should give an expectation value of
(E) = 4.095769122 H. The following were obtained by VQMC
100 tests of 10“ iterations, EN = 0.09 :t 0.03 H
100 tests of 106 iterations, EN = 0.096 a: 0.003 H
10 tests of 109 iterations, EN = 0.09575 :l: 0.00003 H

The number of tests is the number of times that the calculation was repeated.
The energy reported is the average energy from these repeat calculations and

the standard deviation in their values.

129
For a single GTO H atom, a near optimal value of or would be 0.3

giving (E) = 0424039 H. Doing a VQMC calculation of 10 tests of 106
iterations gave E N = -0.4239 1 0.0006 H.

These H atom calculations give a check on several things. First, they
are a check of the correct implementation of VQMC in the computer program
written for this study, called "rwh". Second, they show that the VQMC result

is converging to the analytic result as the number of iterations is increased.

Third, the two tests with 10‘5 iterations illustrate the fact that the standard
deviation for a given size calculation should decrease as the wave fimction is
improved. In fact, calculations have been done where the variation of
parameters was performed by minimizing the deviation rather than minimizing
the energy [21]. Fomth this is a check of the fact that the ran3 random
number generator is appropriate for VQMC applications.

All of the above H atom calculations were done choosing the initial
position for the electron within a cube of side length one bohr centered on the
nucleus. Also, the new location for the electron was chosen randomly within
a box of side length one bohr centered on the present location of the electron.
It would be prudent to test these settings to determine what values will give

the most accurate calculation in the form of the lowest standard deviation.

130
Table 3 .3 gives the results of running 10 tests of 106 iterations for each

of several difl'erent initial box sizes. The settings in all of these calculations

are a = 0.3 and box size for moving the electron of one bohr.

Table 3.3 - Effect of changing the initial box size

 

Lox_s_ize (B) EA (H) Std. dev. (E)

0.01 -0.425 0.001
0.1 -0.424 0.001
1.0 -0.424 0.001
2.0 —0.424 0.002
5.0 -O.423 0.001
10.0 -0.424 0.001
100.0 -0.50 0.05
1000.0 -80 50

Several conclusions may be drawn from Table 3 .3 . First, the
calculation is not extremely sensitive to the small changes in the value of this
setting. Second, in this case a value close to the distance between the nucleus

and the peak in the radial distribution frmction comes out about in the middle

131
of the range of reasonable values. Finally, a very lmreasonable value might

have a very adverse efl'ect on the calculation as with the last calculation in the
table.

The efiect of changing the size of the box within which new electron
positions are chosen was tested. This was done by keeping the same settings
as the above sequence along with an initial box size of one bohr. The results

of these tests are in Table 3.4 .

132
Table 3.4 - Effect of changing the box size

 

hox s_i_ze (B) EN (H) Std. dev. (H)
0.001 -1.3 0.8
0.01 .042 0.08
0.1 .042 0.01
1.0 .0424 0.001
2.0 -0.4238 0.0009
5.0 .0425 0.001
10.0 0425 0.005
100.0 -0.7 0.2
1000.0 -1.1 0.3

Although the calculation will converge with nonoptimal values, it is
worthwhile to use reasonable values.

A comparison of the number of iterations to the calculation accuracy is
given in Table 3.5 These tests were each done with 10 calculations, a = 0.3,

box size and initial box size of 1.0 B.

133

Table 3.5 - Dependence upon number of iterations

 

iterations EN (H) Std. dev. (H)
103 -035 (12
104 -044 0.03
105 -043 0.02
106 -0423 0.003
107 -0.424 0.001
108 -o.4239 0.0003
109 .04240 0.0001
10” -0.42406 0.00006

It is also important to have some indication that a calculation is giving
a reasonable result. One test of a reliable calculation is a virial theorem test.
The virial theorem states that the expectation value of the potential energy
divided by the expectation value of the kinetic energy should equal negan've

two.

(_Vl _
(71 2 (36)

134
The virial theorem check is used in conventional HF calculations to verify that

the wave frmction is a reasonable description of the system. This usage still
holds for VQMC calculations along with a second usage. The virial theorem
test in a VQMC calculation serves as a check that a sufl'lcient nmnber of
iterations have been done to give an adequate statistical sampling of the
probability distribution.

The most extreme example of the need for this test is presented by the
hydrogen atom in which the exact wave frmction is used as discussed at the
beginning of this section. Only a single Monte Cario iteration is necessary to
give the exact energy of -0.5 H. However, all of the other properties of the
atom are still expectation values over the probability distribution. To

illustrate this two tests of 10 iterations were done. The average virial theorem
value was -1.4774780. Doing 10 tests of 104 iterations gave a virial theorem
value of -1.939564 and ten tests of 107 iterations gave -1 .999373 .

2. Helium atom tests

a. Literature survey
A nmnber of Hylleraas style calculations have been done for the He

atom. These are calculations that use an explicitly correlated wave frmction,

similar to VQMC calculations. The distinction between the two is that the

135
Hylleraas calculations used analytic integration of the expectation value

integrals while VQMC calculations use a numerical Monte Carlo integration.
Hylleraas calculations have the advantage of being much faster than
VQMC calculations. The disadvantage is that the flmctional form of the
wave frmction is limited to functions which can be integrated analytically.
This is why Hylleraas calculations are not practical or necessarily even
possible for molecules.
The simplest Hylleraas calculation uses a wave fimction of the form
‘I’=exp(-§r,)exp(-§r2)(l+br,2) (37)
The optimal values of C = 1.849 and b = 0.364 give an energy of 2.892 H [8].
Kinoshita used a correlation function which was a 39 term polynomial
expansion in r1, r; and n; [29]. He computed an energy of -2.9037225 H.
Pekeris used a 1078 term expansion to compute an energy of
-2.903724375 H for the 1S gromld state of the He atom. He also computed an
energy of -1 .47486588483 H for the 3S first electronic excited state of the He
atom [3 7].
The best available He ground state energy seems to be from a 308 term

expansion used by Thakkar and Koga [38]. They computed an energy of

-2.9037243770341144 H.

136
A number of other studies have been done which focused on obtaining

similar accuracy with a smaller expansion form for the correlation function
[39-45]. These studies obtained excellent accuracy with 200 - 300 term
polynomial expansions.

These Hylleraas calculations were all done with correlation functions,
which were polynomial expansions due to the constraint of having to be able
to do the integrals analytically. The choice of correlation flmction for VQMC
calculations is usually based on making the local energy expression easy to
evaluate as well as being an accurate representation of the wave frmction.

Wave flmctions which are most convenient for VQMC calculations are
those of the form given in Equation (10). Furthermore it is most convenient

to evaluate local energy expressions if the correlation frmction, f(rij), has the
form
fro.) = exp( g0,» (38)
The simplest case of this form is illustrated by Tsien and Pack [46]
who used a wave frmction of the form
‘P=exp(-§r,)exp(-§r2)exp(ar12) (39)
They formd the optimal values of the constants to be I; = 1.8581 and

or = 0.2547 to give an energy of -2.8896l8 H.

137
Pritchard used a three term exponential correlation ftmction along with

three STO orbitals that yielded an energy of -2.903148 H [47].
Lester et. al. tried using a 632p GTO expansion for the orbitals with
three s type GTO's for a correlation frmction to get an energy of -2.90039 H

[48].

138
b. Correlation function studies

Correlation flmctions are not as well understood as more well studied
methods such as HF theory. The known boundary conditions have been
stated in the “Theory” section. There are still aspects of these calculations
open to debate and others not yet explored.

Calculations that do not obey all of the bormdary conditions given in
the section on theory were performed. This has been done in order to gain an
tmderstanding of how each condition afl'ects the calculation results.

The calculations that will be presented utilize wave flmctions
consisting of determinants of one electron flmctions times a product of
correlation fimctions as in Equation (10). This immediately opens the
question of what one electron flmctions should be used. There are three
obvious choices.

1. Hydrogenic frmctions.
2. Harlree-Fock frmctions.
3. thctions optimized for use in the VQMC calculations.

The use of hydrogenic frmctions would seem to be indicated by the
flmctional form of the Hamiltonian. The Hamiltonian has kinetic energy and

nucleus - electron attraction terms which depend upon one electron only. The

139

repulsion terms depend upon two electrons. Having a wave function that also
had one and two electron terms that corresponded to these would seem
logical. The one electron flmctions which are eigenfunctions of the
approximate Hamiltonian formed by omitting electron - electron repulsion
terms are the hydrogenic atom flmctions.

The second way to think of VQMC is as a second step correction to
the Hartree-Fock wave frmction as CI and CC calculations are done. This
would give an efficient mechanism for optimizing the one electron functions
in a HF calculation where the integrals can be done more eficiently. Under
this scheme, the electron correlation would be put in twice, first in net efi‘ect
only by the HF wave flmction then explicitly by the correlation function. This
might lead to dificulties due to the correlation flmction having to make up for
the short comings of the HF wave flmction as well as putting in the electron
correlation.

The third option of simply optimizing the one electron fimction for use
in the VQMC wave flmction would admit either of the previous possibilities if
they were the optimal frmctions. However, this is also the route which would

always require the most computation time. Since excessive CPU time is the

140
major draw back of the VQMC method, it would be best to avoid this option

if a more emcient technique can be found.

The three choices of one electron frmctions mentioned here are still a
matter of debate in the literature. There are not yet any rigorous proofs of the
correctness of any of these options. From a pragmatic standpoint any of these
choices can yield accruate energies if the correlation ftmction is sufficiently
- complex to compensate for any inaccmacies in the one electron frmction.
This is a difiicult question to resolve. A marginal correlation function may
yield a better energy than a good correlation frmction if the shortcoming of the
marginal fimctionisofasortthatcanbecompensated forbyashifiinthe one
electron fimction and vice versa.

The rest of this section will be concerned with the shape of the
correlation flmction, firm). As such it is most convenient to use terminology
relative to the correlation function. The equations will continue to be written
in the same notation that has been used thus far. However, the terminology
will reflect a one dimensional system with the origin on one of the electrons
of the helium atom and lengths being distances to the other electron. Thus

f(O) will be referred to as the origin and Koo) will be considered to be infinity.

141
For the sake of the following discussion, the position of the electrons relative

to the nucleus is irrelevant.

The next issue is the value of the correlation flmction at the origin, f(O).
At first glance it might be expected that the value should be zero at this point
in order to give a zero probability of the electrons being at the exact same
point in space (early studies using correlation flmctions assumed, incorrectly,
that this was true [49]). However, there is a theorem stating that the
correlation flmction should have a finite value at this point [26].

There is a theorem stating that the correlation flmction should have a
finite positive slope at the origin thus giving it a cusp [27]. Functions will be
examined which have a zero slope analogous to using GTO functions to
approximate a STO.

Once a normalization scheme has been defined, the exact value of the
slope at the origin can be determined. Some workers have done calculations
which enforce this slope [50]. The studies presented by us will allow the
correlation frmctions to optimize variationally in order to give the best overall
description of the wave frmction. This is the more frequently used technique.

At large distances the correlation function should go to a constant finite

value. This fact was originally presented by Hartree and Ingman [26]. Three

142
later works argued that the correlation function goes to infinity [27,51,52]. It

is now generally recognized that the original work was correct. The
shortcomings inherent in the incorrect proofs will be discussed.

Finally, the correlation flmction should be monotonic (continuously)
increasing fi'om f(0) to fit»). This is usually a tacit assumption in the choice
of correlation frmction. It is possible to obtain a reasonable energy for the He
atom using a correlation frmction that disobeys this condition, as shown by
Roothaan and Weiss [28] and Equation (50).

The following paragraphs summarize the energies which can be
compared to the results presented subsequently for the 1S grormd state of the
helium atom.

The very crude approximation obtained by completely omitting

electron - electron repulsion yields the following equations which are solvable

exactly.
2 2
fi=-_V_I__Y_L__2___7; (4o)
2 2 r1 r2
‘i’=exp(-2rl)exp(-2r2) (41)

E=-4H (42)

143

Using the above wave flmction with the 1/r12 electron - electron
repulsion term in the Hamiltonian, the energy (E) = -2.7495 H is obtained.

Adding more ftmctions with optimized coefiicients and exponents will

eventually result in reaching the HF limit of (E) = -2.86168 H. For this
discussion, Thakkar and Koga's value of (E) = -2.9037243770341144 H will

be considered to be an exact solution.

The most commonly used test of performance is the variation principle
which states that the most accurate energy will be the lowest energy obtained.
In the case of VQMC calculations, both a variational upper bound and a
variational lower bound to the energy can be computed (Equation 15). As the
wave fimction is improved, both of these should converge towards the exact
energy.

There is the possibility of a situation in which the upper bound is closer
for one calculation whereas the lower bound is closer for another calculation.
In this case the calculation with the closer lower bormd probably has a better
overall description of the wave frmction while the other has obtained a better

energy with large positive and negative deviations in various regions of space.

144
A program, called rwhe, was written for performing VQMC

calculations on two electron atoms. The program uses a wave function of the
form given in Equation (10).

In order to compare correlation functions with a range of properties the
one electron frmctions will be held fixed. As a simplest case the following
wave function form will be used

‘P=exp(-2r,)exp(-2r2)f(rlz) (43)
The worst case example would be the correlation ftmction
f( r12 ) = exp( 1/r12 ) (44)
which is plotted in Figure 3.3 . This gives a wave frmction which is not
normalizable due to the correlation frmction going to infinity. Attempting a
VQMC calculation with this frmction resulted in a failed calculation due to a
floating point overflow error.

A function that is finite, but monotonic (continuously) decreasing

(Figm'e 3.4) rather than increasing is
fl r12 ) = eXP( - r12 ) (45)
Doing a VQMC calculation of 10 tests of 106 iterations gave

EN = -0.97 :t 0.03 H. The high energy and large deviation are both the result

 

exp( l/r12 )

145

 

1

Figure 3.3 - Plot of the flmction firm)

r «N

i Nu
i fin
.4. N

r a.—
rr a;
r b.—
r e.—
i n.—
.r v.—
r n.—
i N.—
.r m—
i _

r m6
i ad
r rd
i 0.9
r md
1. ed
r n6
- Nd
r fie

 

1.00E+01 --

1.00E+02 T
OBOE-Kl)

9.0013401 1*
2.00E+01 3r

8.00E+01 ..
7.00E+01 «r
6.00E+01 ~»
5.0013+01 +
4.00901 ..
3.0013+01 «r

h

A
C.’

 

 

exp(-r12)

146

Figure 3.4 - Plot of the ftmction f(rlz)

\

 

 

 

 

0.9 4
0.8 +
0.7 ..
0.6 ~~
0.5 '1’
0.4 it
0.3 it
02 4
0.1 ..

2‘13

 

r m - n
1r a." r QN
e QN r. a."
rr EN 1. EN
r o." \N r 9N
T n.” r1 1 m.”
. e." U . e."
r 0N . MN
.. «a 1m e «a
i —.N e r _.N
r N = fir N
r 04 \m r m.—
r . r r .
.. M.” K ,M.”
. n .
4 0 ~ .0 ma 0 ~
4 2 r m t 2
if V; M r v.—
: @— a1 MA
. e .
r N ~ .m rr N —
r. _.— Cm hr —.—
r — rr _
. 3 M . .3
T as D“. . ad
a he 5 r ho
rfi 0.0 3 r ooc
. 3 e .. 3
L. v.6 W 11 v.9
r M.° o .. M.°
r N.O F r N.O
1T —.° . r —.o
o . at it 4 4. a a 4 o
0 8 7 6 5 4 3 2 J 0
0 0 0 0 o 0 0 0

£23

147

of this being a very poor correlation frmction to the point where it is
detrimental to the calculation rather than an improvement.
The following flmction (Figure 3.5)
t( r12 ) = exp( -l/r12 ) (46)

is monotonic increasing but has both an incorrect value at the origin, f(0) = 0,

and an incorrect slope g! = O. This gave an energy of
r=0

EN = -2.365 :I: 0.003 H.
The opposite test would be the filnction (Figure 3.6)
K ’12 ) = exp( r12 ) (47)
This fimction has a finite value and slope at the origin but goes to infinity at
long distances. This is not a violation of the normalizability boundary
condition because the total wave fimction still goes to zero at large electron
separations. This frmction resulted in the energy E N = -2.435 i 0.006 H.
Tsien and Pack [45] used the similar frmction
f( r12 ) = exp( 0.2547 r12 ) (48)
to obtain an energy of EN = -2.870 :t 0.001 H. Note that this is just slightly
below the HF limit.

Next we will look at a frmction similar in behavior to Equation (46).

 

148

exp( r12 )

Figure 3.6 - Plot of the frmction f(rlz)

.. ad
tr EN
i ON
1.. W.N

a. Na
1. _N

r. a.—
rr a.—
m. E
nr 0.—
rv n.—
n. E
i m.—
.a n;

i ad
1T rd
1. 0d
3. m6
Lt V.°
11 M.°
i nd
4. _d

 

 

25 ~~

20 «>-

15 +

3..

10 ~~

—05
2

1

r

Figure 3.7 - Plot of the frmction f(rn) = exp(-02

. m."
r a."

.. 5.a

r a."
r m."
. 1N
T MN

it N."

r _.N
r N

a m.—
r a.—
r b.—
r o.—
rm.—
. 9.—

rr m.—
.r N._

r m—

r dd
r ad

1r. 5.a

4 9o

4 ad

 

 

8..

 

149
f(r 12) = exp(“ ’1?) (49)

This will have the correct behavior at long distance and incorrectly go to zero
at the origin if a and B are both less than zero (Figure 3.7). This frmction
does, however, have more freedom to find an optimal function shape
compared Equation (46). The optimal values or = -0.2 and [3 = -0.5 were
obtained. This gave an energy of -2.795 1 0.009 H.

Using the correlation flmction from the previous paragraph, holding the
values of 01 and [3 fixed and allowing the orbital exponent to optimize to 1.75
gave an energy of -2.860 :1: 0.003 H which is still slightly higher than the HF
limit.

A flmction with even more freedom to give the correct shape near the
origin would be (Figure 3.8)

f(r,,) = exp(—0.079;; + 0.3 64499r3’”) (50)
This gave an energy of EN = ~2.901 i 0.003 H. This gives 90% of the

correlation energy beyond the HF limit. This is an excellent energy for a two
term expansion. Since the behavior of this frmction far from the origin is

incorrect, it is not adaptable to molecular calculations. It does, however,

0.75
12

+ 03 64499!‘

1.5
12

r

exp(-0.075

150

 

 

 

 

 

.4

Figure 3.8 - Plot of the ftmction f(rn)

1.6 ~-
l.4 ..
12 at
0.8 4r
0.6 4»
0.4 ~~
02 T

A

be

Cr!

 

z- a - a
rr QN rr 0N
4r ad 1r ad
rr «MN 4r 6."
.r 9N \ll/P. 4r QM
.. on n r l. n.
r
t v." 5 .010. 1. an
it MN 0 1H 1.. M.N
i . f K it .
.- u __ t n
L.. Q— \lal. r. d.—
11 w.— rel it a.—
Lr 0.— m 11 °.~
1r W.— T. m at W.—
.. 2 a 2
11 N.— m 4.. N.—
r. —.— Pm hr ”.—
lT — Lil
.. as W. a. as
11 ”.c P hr ”.°
rfi 5.a - hf 5.a
. 9 r .
1.. 06 «m a we
5. m6 e .r md
rr v.6 r tr V6
in Md .m r Md
.4 3 F .T 3
i _d r —d
o n .1 . a a a a o
0 4 2 1 8 6. 4 2 0
.I. l 0 0 0 0

3a

151
serve as an excellent model for the correct shape of the correlation frmction

near the origin.
A correlation flmction often used in molecular computations is the

frmction (Figm'e 3.9)

«or «rm—2:.) 91>

This frmction is called a Jastro frmction, Jastro-Bijl fimction or Pade frmction.
The origin of this frmction is somewhat ambiguous. This is indeed the
flmctional form obtained from a first order Pade approximant. This flmction
does not appear in the papers by Jastro [49] or Bijl [53] which are used as
references to its origin. The earliest paper formd using this frmction is by
Moskowitz and Kalos [54]. Regardless of its origin, this is a simple filnction
that obeys all of the known bormdary conditions on the correlation frmction.
When used for light atoms with HF orbitals, this function consistently gives
50% of the correlation energy. For heavy atoms, problems arise as discussed
in the following section. Typical values of the parameters are a = 0.5 and

B = 1.8 which give the energy EN = -2.7981 $00009 H. The results are

relatively insensitive to scaling both parameters as shown by the fact that the

values or = 4.3322 and B = 4.6717 give the energy EN = -2.8184 :t 0.0009 H.

152
Many correlation filnctions were tried before the correct behavior was

known [28,39,46]. To date, only a handful of fimctions have been tested
which obey all of the known bormdary conditions [21,55]. As such there is

not yet an established systematic way of choosing correlation frmctions of

successively increasing accuracy.

153
The behavior of the correlation frmction at long distances, f(oo), has

been debated in the literature. Hartree and Ingman [26] correctly predicted
that the correlation frmction should go asymptotically to a constant value,
f(oo) = c. They based this conclusion on a combination of physical insight and
mathematical argument. The physical insight used was the fact that the
correlation frmction accounted for the coulombic repulsion between the
electrons; thus, a constant value would give no interaction at long distances.
The mathematical argmnent employed was that a H; wave function would
separate into two H atom wave flmctions at long distances if the constant
value of the correlation frmction at long distances was a non-zero finite value.
The exact value that the wave frmctions goes to would be determined by the
choice of a normalization scheme. Their realization that their method was not
truly rigorous left the way open for other groups to attempt a more complete
analysis.
Baber and Hassé [51] obtained the incorrect result
lim,u_,,f(r,,) = exp(]: r,,) (52)

thus giving f(oo) = oo since they estimated k to be positive for He. They did
not explain how they derived this expression, but did show how they set up

the difi'erential equations. It seems likely that their means for determining

154
which terms in the differential equations were negligible at long distances was

the source of the incorrect result. Green et. al. [52] followed a similar route
in defining a rather complex series fimction for the large separation limit of
firm) which also went to infinity.
Pluvinage [27] used perturbation theory to obtain the incorrect result
lim,u_,,,f(r,2) oc exp(rn) (53)
thus giving Koo) = oo. He had ignored the cross terms of the form 01' f'in
order to make the equations separable for his zero‘11 order approximation. He

was quick to blame his relatively poor energy of (E) = -2.878 H on a poor

zeroth order approximation. He did not deem it necessary to present any
secondary analysis to verify that Equation (53) was qualitatively correct.

The above incorrect results had as their primary verification the fact
that wave flmctions with the specified behavior gave reasonable numerical
results for individual atoms. However, at long distances in an atom the total
wave frmction decays exponentially towards zero. Thus the large separation
behavior of the correlation frmction has very little effect on the total wave
flmction and thus the numerical results.

In a molecule, however, the correlation fimction has a significant

weight in the wave frmction at large distances where it is describing the

155

correlation between electrons centered on two different nuclei. It has been
formd that for molecular calculations, correlation ftmctions obeying the
original conditions set forth by Hartree and Ingman [26] must be used in order
to obtain reasonable results [21].

There is much work yet to be done pertaining to the study of
correlation frmctions. There are not yet llbraries of well tested fimctions
analogous to the binaries of basis frmctions available for HF calculations.
There are also still unanswered questions, as we have pointed out, regarding
the best form of the wave fimction and exact shape of the correlation frmction.

Although correlation flmctions are an important part of VQMC, this
study will move on to exploring the dificulties inherent in applying VQMC
techniques to transition metal atoms. This is being done for pragmatic
reasons. It would be rmwise to spend an enormous amount of time ironing
out difliculties with the He atom just to find some insurmormtable difiiculty
preventing the study of other systems.

c. Other te_st_s

Even though a computer program gets correct results it may not be
doing so as eficiently as possible. Since VQMC calculations can use

enormous amounts of CPU time, it is worthwhile to attempt to optimize the

156
program as much as possible to use a minimal amount of CPU time. This is

done by utilizing functions that are equivalent in result but faster. For
example, computers perform multiplication operations faster than divisions so
it is faster to multiply by 0.5 than to divide by 2. A good set of principles for
designing programs for fast computation speed are
1. Use the algorithm with the best time complexity.
2. Don’t compute the same number twice.
3. Don’t put anything inside of a loop that does not have to be there.
4. Use the fastest operation possible.
The relative speed of operations is as follows
fimctions inside of loops slowest
procedure calls
trigonometric fimctions
conditional statements
division
multiplication
addition or subtraction

assignment fastest

157
The VMCatom program was written first to get the correct answer by

directly transforming equations into C++ som'ce code. It was then rewritten
to optimize performance. As successive improvements were made a time test
was done consisting of a He atom calculation with 10 tests of 10000
iterations using Clementi’s SS STO basis [56] and a Jastro fimction. Table
3 .6 lists the CPU time necessary to do this calculation as successive

improvements to the algorithms were made.

153
Table 3.6 - CPU time test for optimized algorithms

CPU time (minzsec) algorithm changed

 

 

1.28 no compiler optimization

1:28 -01 compiler Optimization

1:11 ~02 compiler optimization

0:57 primitive normalization constants hard coded
0:57 avoid repeat calculation of determinants

0:57 mathematical form of Jastro frmction

0:55 ~03 compiler optimization

0:43 avoid repeat calculation on rejected steps
0:38 recompute only one row of Slater matrix
0:18 normalize all primitives in advance

0:16 save intermediate steps of local energy calc.
0:13 avoid repeat orbital calculations on accepted steps
0:08 implement inverse Slater matrix method

159

Computer science classes no longer put a heavy emphasis on
optimization of algorithms. This is because computers are now fast enough to
make compiler optimization suficient for many applications. Table 3 .6
illustrates the fact that this is not yet the case for quantum chemistry software.

Fokker-Plank importance sampling introduces an additional
approximation to the calculation. In order to test the magnitude of this error,
the He atom test calculation that was used for time testing was performed
with F okker-Plank importance sampling. The energy for this test should be
-2.88 :1: 0.02 H. Fokker-Plank importance sampling with an optimal step size
of 0.2 bohr gave an energy -2.628 :1: 0.008 H. While the standard deviation
did decrease somewhat, the loss in accuracy was too great to be acceptable.

No further use of Fokker-Plank importance sampling was made.

 

160

3. Scandium atom tests

Quantum mechanics calculations on transition metals are particularly
difficult. Transition metal calculations are very sensitive to electron
correlation. A good description of electron correlation can be obtained by the
use of CI calculations, which are extremely time consuming as shown in
Table 3.1 . In order to address many chemical questions it is necessary to
include spin coupling and mass defect terms in the Hamiltonian as well.
These relativistic efi‘ects are usually negligible for first row elements but
become significant for heavy elements.

The choice to explore QMC methods was made based on the favorable
time complexity (Table 3.1). A better time complexity will ensure that a
method is more efiicient for suficiently large systems. For small systems this
may not be true since the overhead in the calculation may be a larger
contributing factor than the over all time complexity.

A second advantage of QMC calculations is that additional terms can
very readily be added to the Hamiltonian. For example, the program
modifications necessary to perform ECP calculations were done in two days
time. Although relativistic terms were not included in this study, the ease of

adding them is an advantage which is expected to be utilized in the future.

161

The original goal of being able to do accurate calculations on
compounds containing transition metals, has translated into the operational
goals of
1. Evaluate the accuracy of VQMC methods.

2. Have the ability to use small test jobs to predict the quality of final results.
3. Evaluate the efficiency of VQMC methods.

The first of these goals has already be well documented in the literature
[21,32]. The second is accomplished with the use of Equation (12). This
leaves the third yet to be addressed by comparing the computational work,
CW, obtained for various methods (Equafion 13).

Although many QMC studies have been done on light elements, the
application to heavy elements presents additional difiiculties [21]. To date
there have been four published QMC calculations on transition metals [33,57-
59].

Christiansen published a paper in which be computed excitation
energies for Sc and Y atoms [57]. These were done as 11 electron
calculations using core potentials created by Hurley et. al. [60] and J 31 ohn et.
al. [61]. VQMC calculations were done using the J astro correlation frmction

(Equation 51).

162

Belohorec, Rothstein and Vrbik published a calculation for the CuH
molecule [33]. The primary focus of their work was the vibrational behavior.
They introduce the split r method, described in the section on computer
implementation, which they used with a damped core calculation. Damped
core calculations use VQMC to describe the core region and diffusion Monte
Carlo to describe the valence. They also utilized data blocking which is
similar in result to the sampling fiequency tested in this study although
somewhat more sophisticated. They did not publish a total energy or valence
energy for their calculation which is highly Imusual for any study regardless
of what properties are being studied.

Mitas did a large study of [Ar] core ECP methods [58]. He used the
same ECP set for HP, VQMC and DMC calculations on Sc, Ti, Mn, Fe, Cu
and Cuz. He used numerical orbitals which were made piecewise smooth
with cube spline fimctions. The correlation frmction used was an expansion
of products of squares of Jastro type fimctions. Mitas computed ionization
energies, electron affinities and the energy of excitation to low lying excited
states. The QMC methods yielded reasonable agreement to experiment. His
QMC calculations were consistently a significant improvement over HF

calculations.

163

In a second study [59] Mitas used [Ne] core potentials for describing
the Fe atom with a range of computational techniques. The methods
examined were HF, VQMC, DMC, CCSD(T), QCISD(T) and CISD (These
last three are variations on the multiple determinant CI method). The energies
computed and reliability of results were the same as in work in the previous
paragraph. The QMC methods yielded more accurate excitation energies
than the large multiple determinant calculations.

The following sections present a mnnber of algorithms applied to the
2D grormd state of the Sc atom. Each technique will be analyzed for both
speed and accuracy.

Two basis sets for the Sc atom are used. The all electron calculations
use Clementi’s lls 6p 5d STO basis set ((E)HF = -759.73555 H) [56]. The
core potential calculations use Hay and Wadt’s eflective core potentials (non-
relativistic for Sc) [34]. The [Ne] core potentials representing the inner 10
electrons and the [Ar] core potentials representing the inner 18 electrons have
been tested.

The quantity E N — var{el,,} is given because it is a variational lower

bormd on the total energy available fi'om VQMC calculations.

164

The acceptance ratio listed is the ratio of the nmnber of iterations on
which the step was accepted to the total number of iterations. An acceptance
ratio very close to one indicates that a larger step size may be more efficient.
An acceptance ratio below 0.5 indicates that a smaller step size may be more
efficient.

A number of small calculations are presented in sections a - k. Each of
these calculations is done with ten tests of one million iterations in order to
compare CPU times. More accurate results are presented in section 1.

Sections a - e present HF calculations (no correlation frmction) with
VQMC integration. These provided a comparison of the speed and accuracy
of various algorithms. HF calculations are done so that the accmacy can be
determined by comparing the results to the known correct value for
Clementi’s basis set.

Sections f and g present correlated all electron calculations. This gives
an illustration of the effect of correlating different shells of electrons.

Sections h - k present core potential calculations with and without
correlation frmctions.

The results of sections a - k are summarized in Table 3.9 in the

“Conclusions” section.

165

a. Simple Metromlis sampling

This is the Metropolis algorithm with the new electron location chosen
within a cube of a set size, referred to as a step size. The following is a
summary of the performance and computational details for this calculation

EN=4w7aosH
EN —var{e1,,} = -12424 H
CW = 184 Hmin

(9) = -

step size = 0.5 B (optimal value)
acceptance ratio = 0.58
CPU time = 4 hr 33 min (SGI with MIPS 4400)

This is the simplest algorithm to implement. Fmthermore it introduces
no additional approximations into the calculation. It may not, however, be
the most efficient computation.

p. Split 1 technique

This technique was original introduced for use in damped core
calculations. It bears testing since this calculation is the only all electron

QMC calculation published thus far. The parameter r was varied to find the

166

value that gave the lowest standard deviation and thus the most efficient

calculation.

EN —var{e1,,} = -14756 H

CW = 2934 Hmin

0%) = -1.86
r = 2.5 (optimal value)
acceptance ratio = 0.90
CPU time = 5 hr 26 min (SGI with MIPS 4400)
The energy is two standard deviations from the correct value -759 H.
The virial theorem check is likewise a bit far from the correct value of -2.
Since the computational work is greater than that for simple Metropolis
sampling, this method does not represent an improvement for VQMC
calculations.
0. Vm'able t technique
This technique uses a step size which is a fimction of the distance from
the nucleus as given by Equations (30) and (31). It was created in an attempt

to circumvent some of the step size bias due to a large discontinuity in step

167

size right at the interface between the core and valence region. The price for
avoiding a large discontinuity in step size was that there was a small step size
bias at practically every iteration rather than only at the steps where core and
valence electrons interchange places.

EN: -738. a 3 H
EN -var{e1,,} = -52534 H

CW = 2835 Hmin

0%) = -134
r = 0.5
acceptance ratio = 0.82
CPUtime = 5 hr 15 min (SGI with MIPS 4400)
The energy of -7 38 H compared to the correct result of -759 H indicates that
this approach introduced far more error than was eliminated. The virial
theorem check of -1.34 gives a second confirmation of this fact.
91. Split 1 b technique
This technique was proposed based on an analysis of the two previous
techniques.

EN: -758.5 a 0.6 H

168

EN —var{el,,} = -l3875 H

CW=104Hmin

0% = -192
r = 2.0 (optimal value)
acceptance ratio = 0.67
CPU time = 4 hr 48 min (SGI with MIPS 4400)
Based on the computational work measm'e of . efficiency, this is an
improvement over simple Metropolis sampling.
This same method was also used to test the efl‘ect of introducing a
th

sample frequency in which the local energy is only summed on every m

iteration. The results of sample frequency tests are given in Table 3 .7

169

Table 3.7 - Sample fiequency tests

 

Sample fimuencv Er_1e_rg1_(H)
1 -758.5 i 0.6
2 -758 :1: l
5 -759 a: 2
10 -759 :1: 1
21 -759 :h 1
42 -759 :l: 1

Less fi'equent sampling did not decrease the standard deviation. Any
improvements to the accuracy of calculation were not visible in the case of
these small test calculations.

e. Split 1, immrtance sampling with rejection

Light element calculations have shown importance sampling techniques
to be a more eficient way of performing numerical integrations. Importance
sampling was used with the split r method because simple Metropolis
sampling and split 1 b were not numerically stable (virial theorem and
energies off by an order of magnitude).

EN: -684. a 5 H

 

170

EN —var{e1,,} = .4951 H

CW = 13050 Hmin

0% = -123
r = 0.5 (optimal value)
acceptance ratio = 0.89
CPUtime= 8hr 42min (SGIwithMIPS4400)
For Sc, this technique does not show the same increase in speed that

has been observed for light elements. The energy and virial theorem checks

are also very- poor.
Tests with lighter atoms verified that this method was implemented

correctly.

171

f. Split 1 b, correlating 3s 3p 3d 4s

Here a Jastro correlation frmction is applied to all pairs of the outer ll
electrons. This increases the accuracy of the calculation by including
correlation. The Jastro frmction parameters, or and B, were varied to give the

lowest total energy.

EN: “760. i 2 H
EN —var{e1,,} = -14297 H

CW = 2156 Hmin

0%) = -197
a= 2.5 (optimal)
fl= 3-0 (Optimal)
r = 2.0
acceptance ratio = 0.69
CPU time= 8hr59min (SGI with MIPS 4400)
The additional CPU time required for computing the correlated energy

expression is evident by the comparison to the 4 hour CPU time in section d.

172

g. Split 1 b, correlating 28 2p 3s 3p 3d 48

A J astro correlation frmction was applied to all pairs of the outer 19
electrons. This gave the following results.

EN: ’758.7 It 0.9 H
EN -var{e1,,} = -13750 H

CW = 787 Hmin

0%) = -1.92

a = 0.04 (optimal)

)5 = 1-84 (Optimal)

r = 2.0

acceptance ratio = 0.67

CPU time= 16hr 11 min (SGIwithMIPS4400)

This not only required more CPU time, but gave a higher energy than

the calculation in section f. This loss of accmacy can most logically be
explained as being due to the inability of a single set of J astro parameters to

represent the correlation in both valence and inner core regions.

173

h. lflelcore ECP
VQMC with simple Metropolis sampling was used to do a HF

calculation describing 11 electrons explicitly. This was done for the sake of
timing comparison.

It is interesting to note that ECP frmctions can be added into the
VQMC energy expression extremely easily. In fact, the increase in speed due
to the use of an ECP basis is greater for VQMC calculations than for
convention HF GTO calculations.

EN = -45.94 i 0.05 H

EN —var{el,,} = -96 H

CW = 0.27 H min

step size = 1.0 B (optimal)

acceptance ratio = 0.67

CPU time = 1 hr 47 min (SGI with MIPS 4400)

The virial theorem test is not reported for core potential calculations
because the virial theorem does not hold for these calculations, QMC or

otherwise.

174

i. INeI core ECP with Jastro flmction

A Jastro function was applied to all pairs of the 11 outer electrons to

obtain

EN = -45.96 a 0.06 H
EN —Var{eln} = '96 H

CW = 0.89 Hmin

a = 0.01 (optimal)
B = 0.5 (optimal)
step size = 1.0 B

acceptanceratio= 0.68
CPUtime= 4hr 6min (SGIwithMIPS4400)

j. |Ar| core ECP
A HF calculation integrating the outer 3 electrons explicitly gave

EN = -153724 :1: 0.02 H
EN -var{e1,,} = -67 H

CW = 0.0084 Hmin

step size = 2.5 (optimal)

acceptance ratio = 0.65

175

CPU time = 21 min (SGI with MIPS 4400)

k. |Ar| core ECP with Jastro ftmction

A Jastro frmction applied to all pairs of the 3 outer electrons gave

EN = -1.55 :l: 0.06 H
EN -var{e1,,} = —64 H

CW = 0.10 Hmin

or = 2.5 (optimal)
B = 3.0 (optimal)
step size = 2.5

acceptance ratio = 0.66
CPU time = 28 min (SGI with MIPS 4400)

176

l. Accurate VQMC Mpg

Due to the necessary CPU time for performing all electron calculations,
even with the most efiicient methods tested, core potential calculations were
used to test the accuracy of the energies obtained by VQMC.

Table 3 .8 lists the valence energies and excitation energies obtained.
The [Ar] core ECP calculations used 10 tests of 107 iterations. The [Ne] core
ECP calculations used 10 tests of 108 iterations for the ground state and 22
tests of 108 iterations for the 4F state. The rest of the calculation details were

thesameasinsectionsh-k.

177

Table 3.8 - Energies using ECP

[Ar] core with J astro valence energies

State Energyflfl
4s2 3d1 21) -1545 :1: 0.004
3d3 4F -1399 :1: 0.008

[Ne] core with Jastro valence energies

Stag Energyfl)
4s2 3d1 21) 45.936 2 0.004
3d3 4F 45.786 :1: 0.008

4s2 3d1 -) 3d3 excitation energy (eV)
HE V MC Exp.
[Ar] core 4.42 3.97 4.19

[Ne] core 4.47 4.08 4.19

178

Table 3.9 - Computational work and energy

CPU time‘

Mm Milan) E_N(Ii) Gauguin)
[Ar] core potential 0.0085 -1.54 i 0.02 0:21
[Ar] core potential with correlation 0.10 -l .55 :1: 0.06 0:28
[Ne] core potential 0.27 -45.94 :t 0.05 1:47
[Ne] core potential with correlation 0.89 -45.96 :1: 0.06 4:06
Split tau b 104 -758.5 d: 0.6 4:48
Simple Metropolis sampling 184 -759.7 d: 0.8 4:33
Split tan b, correlating 19 electrons 787 -758.7 :1: 0.9 16:11
Split tan b, correlating ll electrons 2156 -760. :t 2 8:59
Variable tau 2835 -738. :t 3 5:15
Split tau 2934 -753. i 3 5:26
Split tau, importance sampling 13050 -684. :1: 5 8:42

with rejection

"' SGI Onyx with MIPS 4400 processor

179

F. Conclusions

The computational work for each of the methods discussed is
summarized in Table 3.9 . As expected, the core potential calculations were
the fastest. Although very fast, the [Ar] core potentials are not accurate
enough for many chemical applications.

The “split 1 b” method, introduced here, is the most efficient all-
electron VQMC algorithm tested. Simple Metropolis sampling does
sm'prisingly well for all electron transition metal calculations. The use of
importance sampling, which has seen much use for light elements, was the
least efficient method.

Accurate ECP calculations are presented in Table 3 .8 . As expected,
the correlated [Ne] core frmctions came the closest to the experimental
results. Fmther improvement of the accuracy of these calculations can be
achieved by the use of more complex correlation frmctions and optimization
of orbital coefi'rcients. The [Ne] core potentials may be an optimal balance of
speed and accuracy for many applications.

The VQMC method has given us the ability to correlate 3s, 3p, 3d and
4s electrons for transition metals. The current generation of work station

computers that are readily afi‘ordable by academic institutions typically have

180

several hundred Mb of RAM and several Gb of disk space. For such
machines, VQMC represents the only available way to do such calculations
with this level of accruacy. The most accurate [Ne] core ECP calculations
required in excess of a week worth of CPU time on the fastest SGI work
stations available. As such these calculations are just becoming practical for
small parallel processing machines.

The major downfall of all of the QMC methods is the requirement for
large CPU times. However, this detriment is being overcome by increases in
CPU performance and availability of parallel computing facilities. Indeed,
VQMC calculations have the potential to be the next work horse of the
computational chemistry industry.

At the risk of being wrong, the following listing is an attempt to predict
what advances to the VQMC method will be the most crucial in developing
methods with the accuracy and eficiency already in demand by the chemical
industry. The list is numbered in order of most important developments first.

1. Development of accurate correlation frmctions.
2. Convenient user interfaces.
3. Parallelization.

4. More eficient sampling schemes.

 

LIST OF REFERENCES

 

10.

11.

12.

13.

LIST OF REFERENCES

E. Schrodinger Annln. Phys. :72, 361 and 489; 89, 437; 8_l, 109 (1926)

L. Pauling, E. B. Wilson, Jr. “Introduction to Quantum Mechanics
With Applications to Chemistry ” Dover (1985)

I. N. Levine “Quantum Chemistry ” Prentice Hall (1991)

D. R. Bates, K. Ledsham, A. L. Stewart, Philosophical Transactions
of the Royal Society of London A246, 215 (1953)

W. Bian, C. Deng Theor. Chim. Acta. 22, 135 (1995)

R. G. Parr, W. Yang “Density-Functional Theory of Atoms and
Molecules” Oxford (1989)

D. R. Hartree Proc. Cambridge Phil. Soc. 251, 89, 111, 426 (1928)
E. A. Hylleraas Z. Physik _SA, 347 (1929)

Conference proceedings in Rev. Mod. Phys. volume _3_2, nmnber 2
(1960)

T. L. Gilbert Rev. Mod. Phys. 3_5, 491 (1963)
W. Kolos, C. C. J. Roothaan Rev. Mod. Phys. g, 205, 219, (1960)

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E.
Teller J. Chem. Phys. 2_1, 1087 (1953)

H. ConroyJ. Chem. Phys. 4_1_, 1327, 1331, 1336, 1341 (1964)

181

 

14.
15.

16.

17.
18.

19.
20.
21.

22.
23.
24.
25 .

26.

27.
28.
29.

30.

182

H. Conroy J. Chem. Phys. 41, 912, 930, 5307 (1967)
H. Conroy, B. L. Bnmer J. Chem. Phys. _4_7, 921 (1967)

D. Ceperley, G. V. Chester, M. H. Kalos Phys. Rev. B. 16, 3081
(1977)

D. E. Knuth “Seminumerical Algorithms ” Addison Wesley (1969)
J. Gleick “Chaos ” Viking (1987)

w. H. Press, B. P. Flannery, s. A. Teukolsky, w. T. Vetterling
“Numerical Recipies ” Cambridge (1989)

I. Vattulainen, T. Ala-Nissila, K. Kankaala Phys. Rev. Lett. ‘_7_3_, 2513
(1994)

B. L. Hammond, W. A. Lester, Jr., P. J. Reynolds “Monte Carlo
Methods in Ab Initio Quantum Chemistry” World Scientific (1994)

W. Hierse, P. M. Oppeneer Int. J. Quantum Chem. fl, 1249 (1994)
H. W. Jones Int. J. Quantum Chem. 5_1_, 417 (1994)

O. Treutler, R. Ahlrichs J. Chem. Phys. _lfl, 346 (1995)

I. I. Guscinov, R. F. Yassen Int. J. Quantum Chem. 3, 197 (1995)

D. R. Hartree, A. L. Ingman Mem. Proc. Manchester Lit & Phil. Soc.
_71, 79 (1933)

P. Pluvinage Ann. Phys. 5, 145 (1950)
C. C. J. Roothaan, A. W. Weiss Rev. Mod. Phys. 32, 194 (1960)
T. Kinoshita Phys. Rev. 10_5_, 1490 (1957)

P. PulayJ. Comp. Chem. 3, 556 (1982)

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

183

M. H. Kalos, P. A. Witlock “Monte Carlo Methods Volume 1: Basics ”
John Wiley & Sons (1986)

J. B. Anderson, C. A. Traynor, B. M. Boghosian J. Chem. Phys. 92,
345 (1993)

P. Belohorec, S. M. Rothstein, J. Vrbik J. Chem. Phys. E, 6401
(1993)

P. J. Hay, W. R. Wadt J. Chem. Phys. 82, 270, 299 (1985)
M. Krauss, W. J. Stevens Ann. Rev. Phys. Chem. 3_5_, 357 (1984)

J. F. Harrison, K. L. Krmze “Organometallic Ion Chemistry” 89,
Kluwer (1996)

C. L. Pekeris Phys. Rev. _11_5_, 5 (1959)
A. J. Thakkar, T. Koga Phys. Rev. A _S_Q, 845 (1994)

D. E. Frermd, B. D. Huxtable, J. D. Morgan, III Phys. Rev. A Q, 980
(1984)

J. Midtdal Phys. Math. Univ. Osloensis 32 (1966)

T. Koga, Y. Kasai, A. J. Thakkar Int. J. Quantum Chem. 44, 689
(1993)

H. M. Schwartz Phys. Rev. m3, 110 (1956)
H. M. Schwartz Phys. Rev. m, 483 (1960)
H. M. Schwartz Phys. Rev. m, 1029 (1963)
C. Schwartz Phys. Rev. _12_8, 1146 (1962)

r. P. Tsien, R. r. Pack J. Chem. Phys. 49, 4247 (1968)

 

47.

48.

49.

50.

51.

52.

53.

54.

55.

56.

57.

58.

59.

60.

61.

184

H. O. Pritchard J. Chem. Phys. 49, 1496 (1964)

W. A. Lester, Jr., M. Krauss J. Chem. Phys. fl, 1407 (1964)

R. Jastrow Phys. Rev. 2_8_, 1479 (1955)

S. E. Koonin “Computational Physics ” Benjamin/Cummings (1986)

T. D. H. Baber, H. R. Hasse' Proc. Cambridge Phil. Soc. 33, 253
(1937)

L. C. Green, C. Stephens, E. K. Kolchin, C. C. Chen, P. P. Rush J.
Chem. Phys. fl, 1061 (1959)

A. Bijl Physica 1, 869 (1940)

J. w. Moskowitz, M. H. Kalos Int. J. Quantum Chem. 24, 1107
(1981)

H.-J. Flad, A. Savin J. Chem. Phys. 143, 691 (1995)

E. Clementi, “Ab Initio Computations in Atoms and Molecules”, IBM
Journal of Research and Development 2, 2 (1965)

P. A. Christianst Chem. Phys. _9_5_, 361 (1991)

L. Mitds “Computer Simulation Studies in Condensed-Matter Physics
V” Eds. D. P. Landau, K. K. Mon, H. B. Schr’ittler; Springer-Veflag
(1993)

L. Mitas Phys. Rev. A 42, 4411 (1994)

M. M. Hurley, L. F. Pacios, P. A. Christiansen, R. B. Ross, W. C.
Ermler J. Chem. Phys. 8_4, 6840 (1986)

L. A. LaJohn, P. A. Christiansen, R. B. Ross, T. Atashroo, W. C.
Ermler J. Chem. Phys. 8_7, 2812 (1987)

Chapter IV

PROPERTIES OF TOTAL ELECTRON DENSITY

A. Introduction

In the attempt to determine the physical nature of lone pair electrons
(Chapter 11), one common misconception formd was that these electrons were
visible as a maximum in the total electron density distribution. Orbital
descriptions of bonding, although a useful description of chemistry, have also
led to confusion between orbital densities and the total electron density.
These misconceptions prompted the writing of a paper on electrdn density for
the Journal of Chemical Education [1]. The focus of this work is a
clarification of what theorems describing the electron density exist, what is
typically observed and what exceptions have been formd.
B. Misconceptions about electron density

While students are taught orbital descriptions of bonding and geometry,
they often have misconceptions about the name of the total electron density.

For example, many students and experienced researchers expect to see large

185

186

bumps in the electron density where the lone pair electrons are located which
simply don’t exist.

Some misconceptions are starting to be eliminated through the
availability of programs to graphically display the electron density. However
the most general, and thus most powerful, statement of the form of the
electron density is given by the associated theorems. Theorems are
generalized descriptions of the behavior of all atoms or molecules. As such,
theorems will describe trends which are true for many systems, but will not
predict exact values for specific molecules. Although some scientists may be
well aquainted with these theorems, many misconceptions remain among
chemists at all levels.

The remainder of this chapter will attempt to srunmarize the existing
theorems describing the behavior of the electron density. For cases, where
theorems have not been constructed, the normal behavior and anomalous
cases will be described. We will also give a very brief description of some of

the other properties which are often confused with the electron density.

 

 

187

EMILE

The electron density, p, is the number of electrons expected to be
found per 1mit volume of space, usually expressed in units of electrons/A3 .
This is the number of electrons times the probability of finding an electron.
As such, we do not examine the trajectory of a single electron as we would
the planets in the solar system. In an orbital description this is the smn of the
densities of the occupied orbitals. The electron density is sometimes referred
to as the total electron density in order to difierentiate it from orbital
densities.

For the hydrogen atom, the density is known exactly [2]. For the

ground state of the hydrogen atom

p0) = ire-2' (1)

7:
Thus the density is a maximum on the nucleus and decays exponentially and
monotonically (continuously) to zero at infinite distance as shown in figure
4.1 . The electron density is spread away from the proton in order to obtain
the lowest energy balance between potential energy which is lowest near the
nucleus and kinetic energy which is highest near the nucleus. For excited

states, the decay to infinity is still exponential, but there is not necessarily a

 

188

Figure 4.1 - Electron density for the grolmd state of the H atom

 

 

 

13
14
15
1.64
174»
1a

1.9 1)

 

189

peak on the nucleus, nor is the decay monotonic, as evidenced by the 2p
density shown in Figure 4.2 .
D. Theorems
The wave frmction, and thus p, must decay to zero at long distances
from an atom or molecule in order to satisfy the normalizability condition [2].
An upper bormd for the decay of the electron density to zero at infinite
distance was derived by Hofl'rnann-Ostenhof & Hofinann-Ostenhof [3].
Their proof utilized properties of differential equations, the variational
principle and the Cauchy-Schwartz inequality. For molecules, they found that
p(r) s k r‘2 (r — p)"’“ e”“"”’ for large values of r (2)
a = (28)"2 (3)
where k is a constant, a is the ionization potential of the molecule, Z is the
sum of nuclear charges, r is the distance from the center of mass and p is the
distance fi'om the center of mass to the nucleus farthest fi'om the center of

mass. For atoms this reduces to

p(r) s k r” e’“ (4)

fl=Z/a-1 (5)

 

190

Figure 4.2 - Electron density extending away from the nucleus of a hydrogen
atom in the excited state with the electron in a 2p orbital.

 

 

 

191

This work is based upon the analysis of wave frmction behavior by Ahhiches
[4-6]. A more detailed mathematical wave frmction analysis has been
published by Agrnon [7].

Ta] [8] increased the accmacy of the previous expression for atoms by

deriving

 

fl=Z-n+l_1 (6)

a
where n is the number of electrons. His derivation utilized a density matrix
theory formulation. This is in agreement with the wave frmction analysis of
Katriel and Davidson [9].

For the ground states of atoms and atomic ions, Weinstein, Politzer &
Srebrenik [10] have shown that the density is a maximum at the nucleus and
decays monotonically. This is proven as a statement which they did not put
in equation form. Their analysis shows this to be a direct consequence of
Poisson’s equation. These findings are further supported by the approximate
analysis of Angulo & Dehesa [l 1] based on the Stieltjes-moment-problem
technique.

‘ The preceding results show that the electron density of the grormd state

of all atoms is expected to look qualitatively like the hydrogen atom ls

192

density plotted in Figure 4.1 . For excited states of atoms, the density need
not decay monotonically or be a maximum on the nucleus [10].
Hofinann-Ostenhof & Morgan [12] have analyzed the behavior of one

electron homonuclear diatomic molecules (H; , He;3 , etc.). They have shown

that there can be maxima in the electron density only on the nuclei. They
have also shown that there are no local minima. The density is a minimmn at
the center of the bond along the line connecting the nuclei and a maxirmun
along lines perpendicular to the bond. They derived these results by
analyzing the wave frmction with differential inequality techniques in several
coordinate systems.
E. Cases nqt covered by Theorems

Although no formal proofs have been presented regarding maxima and
monotonicity in molecules, many studies have been done [13-16]. For the
grormd state electron distribution in almost all of the molecules ever studied
the only maxima in the electron density are on the nuclei. This is illustrated
by Figure 4.3 in which the electron density in the plane of the Czl-Ia molecule
is plotted as a displacement in the third dimension. The peaks on the carbon

atoms have been tnmcated in order to make the hydrogen atom peaks visible.

193

Figure 4.3 - Total electron density in the plane of the C211. molecule plotted

as the height above the plane. The density was obtained from a CISD

calculation with a 6-311G*"‘-H- basis set.

 

 

 

 

 

 

 

 

 

 

 

 

J/\/ (7\

 

 

 

 

194

The density around each atom appears nearly spherical when plotted as
shown for the CO molecule in Figure 4.4 . The bonding and lone pair regions
show some elongation. Note that there are no pronounced tear-drop shaped
regions associated with individual bonds or lone pair electrons. Thus the
electron density for a molecule is very similar to the sum of the densities of
the component atoms. The shift in electron density upon forming a molecule
from the separated atoms is typically less than or equal to 10% of the atomic
densities [10].

In a few rare cases, a peak in the electron density is found in the
middle of the bond for the grormd state of a molecule. This maximum in the
bond region is referred to as a nonnuclear attractor. The most pronomrced
and well studied case is the Li; molecule [17]. The magnitude of this peak is
three orders of magnitude smaller than the peaks on the nuclei as shown in
Figure 4.5 . This behavior has been shown to be expected only in very
polarizable systems such as Li and Na clusters.

F. Other molecular progrties
The remaining question is “Where does the chemistry come fi'om?”

Chemistry is energetically a very subtle effect as evidenced by a carbon-

195

Figure 4.4 - contour diagram of the electron density of the CO molecule.
The carbon atom is on the. left. Contours are plotted for density values of
0.1, 0.2, 0.3, 0.4, 0.5, 1.0 and 5.0 electrons per cubic bohr.

The density was obtained from a configuration interaction calculation
with a 6-311G*+ basis set.

The computed bond length was 1.125 Angstroms.

/\
(a
V/

 

 

 

 

 

 

 

196

Figure 4.5 - Plot of the electron density along the internuclear line of the Li;
molecule illustrating the magnitude of the maximrun in the middle of the
bond. The density was obtained from a CISD calculation with a 6-311G“+
basis set and a bond length of 2.692 A

L 1 .
fi' T
‘ _ ‘T...’
____,.__.——-——
“.-
- 537:»

density (olAngstrom"3)
'8 8 8 S 8 3‘ 8 8

 

‘

C G

L

T
//’
/

/

t\
i
i.
(

AAAAAA

-25
.23 I.
-21 i
-1 9 J
-17
-15
-13
-1 1 .
-03 Y
-07 i
.95
-03

...........................

~
0
.......................................................................................................................
a

 

 

a
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

 

197

carbon single bond having an energy of arormd 4.3 eV while 1029.8 eV are
required to remove all of the electrons fi'om a carbon atom.

There are a number of molecular properties, other than the electron
density, which illustrate chemical structure much more graphically. The
following are short descriptions of a number of such properties and their
relationship to the electron density. These are mentioned here only in order
to decrease confusion. It is not the intention of this paper to give a detailed
description of these properties.

Electron density difference plots show the result of subtracting the
atom densities fi'om the molecular densities. This gives an indication of how
the density shifts on bond formation rather than the absolute density [2, 16].

In radical compormds, the spins of all of the electrons do not cancel to
a net zero spin. Spin densities are plots of the net electron spin in radical
compounds [13, l4, 16].

The electrons in an atom are often assigned to K, L and M shells. This
atomic shell structure is seen in the total radial distribution fimction but not in
the electron density [2]. The radial distribution function gives the probability
of finding an electron at a given distance from the nucleus of an atom. The

radial distribution frmction is obtained by integrating over the angular

198

variables in order to obtain an angular averaged frmction, which is a function
only of the radial distance from the nucleus.

Much of our understanding of bonding comes from molecular orbitals
[2, 13, 14]. Orbitals, such as those obtained from Hartree-Fock calculations,
are one electron frmctions, meaning that they show the behavior of one
electron in the presence of all of the other electrons. Orbitals, which are
similar but not identical, can be generated from other calculations such as
configruation interaction and density frmctional theory calculations.
Molecular orbitals can show quite a bit of feature such as sigma and pi
bonding characters or electron lone pairs. The square of an orbital gives the
contribution to the electron density from that orbital. The total electron
density is the sum of the electron densities from all of the orbitals.

The energy of an infinitesimally small charge placed at a point in space
is the charge times the molecular electrostatic potential evaluated at that point
[16]. The electrostatic potential at a point in space is calculatedby summing
the repulsions with all of the nuclei and adding to that the integral of the
attractions to the electron density at all points in space. This is sometimes
used as an indication of hydrophobic vs. hydrophilic regions or equated to the

behavior of lone pair electrons.

199

The gradient of the electron density is a vector quantity. This shows
how electrons would move most rapidly from regions of high electron density
to regions of low electron density [15].

The Laplacian of the electron density, V2 p , shows the behavior of the
electron density relative to an exponential decay. There will be a negative
minimum in the Laplacian of the electron density where there is the maximum
positive deviation from exponential decay. Bader and coworkers have
utilized this property to interpret electron lone pairs as maxima in plots of
—V2 p [15]. Note that this is a maximum relative to an exponential decay, but
is still a region of monotonically decaying electron density.

G. Conclusions

A nmnber of theorems have been constructed to describe the behavior
of the electron density for simple systems. However, existing theorems do
not describe all aspects of the electron density for all systems. Descriptions
of the expected behavior of the electron density in molecules are supported by
a large volume of theoretical and experimental evidence.

The electron density is for the most part featureless. The true nature of
the total electron density has been shown to be devoid of large peaks

denoting bonds and lone pair electrons. 'Featm’es which are sometimes

200

attributed to the electron density are often those of other properties which
focus on the more subtle effects which dictate the chemistry. Most of these
properties are related to the electron density in some way. Much insight can
be gained from these related techniques. However, their behavior should not

be confused with that of the total electron density.

 

LIST OF REFERENCES

10.

ll.

12.

13.

LIST OF REFERENCES

D. C. Young submitted to J. Chem. Educ.
I. N. Levine “Quantum Chemist ” Prentice Hall (1991)

M. Hoflinann-Ostenhof, T. Hofirnmn-Ostenhof Phys. Rev. A ,1_6, 1782
(1977)

R. Ahlrichs Chem. Phys. Lett. ii, 609 (1972)

R. Ahlrichs J. Math. Phys. 1;, 1860 (1973)

R. Ahlrichs Chem. Phys. Lett. _1_8, 521 (1973)

S. Agrnon “Lectures on Exponential Decay of Solutions of Second-
Order Elliptic Equations: Bounds on E igenfimctions of N-Body
Schrodinger Operators” Princeton Univerisity Press (1982)

Y. Tal Phys. Rev. A _1__8_, 1781 (1978)

J. Katriel, E. R. Davidson Proc. Natl. Acad. Sci. USA _71, 4403 (1980)

H. Weinstein, P. Politzer, S. Srebrenik Theoret. Chim. Acta (Berl.) 38,
159 (197 5)

J. C. Angulo, J. S. Dehesa Phys. Rev. A 44, 1516 (1991)

T. Hofi'mann-Ostenhof, J. D. Morgan III J. Chem. Phys. 2;, 843
(1981)

E. R. Davidson “Reduced Density Matrices in Quantum Chemistry”
Academic Press (1976)

201

14.

15.

16.

17.

202

R. G. Parr, W. Yang “Density-Functional Theory of Atoms and
Molecules” Oxford University Press (1989)

R. F. W. Bader “Atoms in Molecules A Quantum Theory” Clarendon
Press (1994)

P. Coppens, M. B. Hall, Eds. “Electron Distributions and the
Chemical Bond” Plenum Press (1982)

G. I. Bersuker, C. Peng, J. E. Boggs J. Phys. Chem. 91, 9323 (1993)

"tannin)“