Molecular multipole moments and their performance in the self-consistent reaction field model


Material Information

Molecular multipole moments and their performance in the self-consistent reaction field model
Physical Description:
x, 111 leaves : ill. ; 29 cm.
Zheng, Xuehe
Publication Date:


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


Thesis (Ph. D.)--University of Florida, 1995.
Includes bibliographical references (leaves 105-109).
Statement of Responsibility:
by Xuehe Zheng.
General Note:
General Note:

Record Information

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

This item is only available as the following downloads:

Full Text







To the Clbhs


I would like to thank Professor Michael C. Zemer for introducing me to the challenges

that lead to this thesis. Professor Zerner has always showed great enthusiasm for the

work presented here, and he let me watch over his shoulder as we chased electrons,

yet still, he is brave enough to give the carte blanche for me to chase on my own. I

also benefited from many other workers that clustered in the Zerner group. I would

especially like to thank Professor Paul Wormer (Nijmegen) for several discussions on

rotational and translational methods used in angular momentum, and pointing to some

important literature; Dr. Jim Rabinowitz (EPA) for some in-depth look at the multiple

moments; Professor Mati Karelson (Tartu University) and Dr. Tamm (QTP, FL) for some

ground work on the SCRF code; Professor Dan Edwards (Idaho and QTP, FL) for much

help with the general structure of MO and CI calculations in the ZINDO package; and

particularly, Professor Feng Jikang (Jilin and QTP, FL) for undergraduate instruction and

early cooperation on some interesting applied work using the ZINDO package. I also

acknowledge the devastating influence of the boys and girls of the legendary clubhouse,

which forever shaped the character and style of this work, and all future works to come.

I thank Professor N. Yngve Ohm for helping to provide the uniquely stimulating

environment of the Quantum Theory Project, and the constant encouragement that enabled

me to regain my sense of purpose for this thesis!


I certainly thank Dr. Marshall G. Cory (QTP, FL) for his help in constructing this


I am most grateful to my family, my mother, brother and two sisters, half a globe

away, for their unfailing love, patience and support. For the same I also thank my

extended family that I met quickly upon my arrival, Dr. E.T. and Mrs. Vam York. I also

thank my wonderful friends: Anne Kassem, for many enjoyable conversations, dinners,

and adventures in exploring natural Florida from its back roads to its shining sea; the

DeWildes and Wei, for always being there; Mrs. Vanderwerf and late Dr. Vanderwerf

for encouragement and friendship; and Eva for her understanding and tolerance of the

many late nights spent on writing this thesis.


ACKNOWLEDGMENTS ...... ... ... ....... ............ iii

LIST OF TABLES ......... ... .................... vii

LIST OF FIGURES ....... .. ............. ........... viii

ABSTRACT ........... ... ............ .... .......... ix


1 INTRODUCTION .................................. 1

OVER SLATER-TYPE ORBITALS ........................ 8
Introduction ........... ... ........ .. .. ........ 8
Preliminary Considerations.............................. 9
M ethod .. .. . .. .. .. 14
Discussion and Conclusions .......................... 33

M OM ENTS ..................................... 39
Definitions of Electric Multipole Moments ................ 39
The Theory and Calculation of Molecular Orbitals ............. 47
Some Calculated Molecular Multipole Moments ............. 57

The Reaction Field ............................. 71
The Interaction Energy in a Reaction Field ................ 74
The Self-Consistent Reaction Model .................... 75
The Configuration Interaction Calculation for Electronic Spectra ..... 78
Results . . . . 87

5 CONCLUSIONS ................................. 100

REFERENCES ...................................... 105

BIOGRAPHICAL SKETCH .............................. 111

Table 2-1:

Table 3-2:

Table 3-3:

Table 3-4:

Table 3-5:

Table 3-6:

Table 4-7:

Table 4-8:

Table 4-9:

Table 4-10:

Table 4-11:

Table 4-12:


Percentage of Cartesians needed to Compute . .

The Conversion Factors for 1 a.u. Multipole Moments .

Calculated Multipole Moments of H20 ............

Calculated Multipole Moments of Benzene . .

Calculated Multipole Moments of Hexafluorobenzene .

Calculated Multipole Moments for Pyrazine in Water Solution

H20 solvent cost energy convergence in H20 solution .

Calculated acetone n r* solvent shifts (cm-1): Theory A..

Calculated acetone n ir* solvent shifts (cm-1): Theory Al.

Calculated acetone n 7r* solvent shifts (cm-1): Theory B..

Calculated acetone n ir* solvent shifts (cm-1): Theory Bl.

Calculated acetone n 7r* solvent shifts (cm-1): Theory C..

Table 4-13: Calculated acetone n 7r* solvent shifts (cm-1): Theory Cl.

. 38

. 63

. 64

. .. 66

. .. 67

. .. 68

. .. 92

. .. 93

. 94

. .. 95

. .. 96

. .. 97

. .. 98

Table 4-14: Pyrazine n-7r* Solvent Shift (cm-1): From Acetonitrile to Diethyl Ether. 99


Figure 1: Rotation From Ellipsoidal Coordinate to Molecular Frame ........ 36

Figure 2: Flow chart of the HIMO computational program ............. 37

Figure 3: Input Molecular Orientations for Multipole Moment Calculations. ... 69

Figure 4: Dimer Structures of Hexafluorobenzene-Benzene and Benzene-Benzene
Based on Multipole Moment Interpretations . . 70

Figure 5: Spherical Cavity Model ........................... 90

Figure 6: Solvent Shift Dependence on Cavity Radii: Cl theory for acetone. 91

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




May, 1995

Chairman: Michael C. Zerner
Major Department: Chemistry

Molecular electric multiple moments are calculated to high orders using the semiem-

pirical Intermediate Differential Overlap (INDO) wavefunction at the Hartree-Fock level

of theory. A general approach is developed for multiple moment integrals arising in

computations using Slater-type orbitals (STO). A computational package implementing

this approach is produced in this study and in our experience it is fast and reliable.

Molecular multiple moments calculated with the existing INDO minium basis set agree

well with results in the available literature.

Calculated molecular multiple moments are included in the Self-Consistent Reaction

Field (SCRF) model through the first order multiple expansion of electrostatic potentials

between solute molecules and solvent dielectric media. Computations are implemented


in the Zemer's INDO (ZINDO) package to calculated solvent shifts of electronic spectra

with the parametrized Configuration Interaction Singles (CIS) method. Spectroscopic

shifts previously not accounted for by the Onsager term are produced in good agreement

with experiment by using the newly developed theories C and Cl. Multipole expansion

convergence, cavity radius dependence, and other aspects of the SCRF model are studied

and discussed.


Much chemistry is understood by simple interpretations based on dipole moments;

even more is understood by using the quadrupole moments. In fact, a study of multiple

moments is a popular, and arguably the most effective, language to address a variety

of phenomena at the molecular level. However, in practice this has been a language

with very limited vocabulary because a chemist can measure only dipole and quadrupole

moments of molecules with any accuracy. Understanding of more subtle molecular

interactions, structure, and physical and chemical behavior must resort to an examination

of higher order multiple moments that, regretfully, have not been obtainable from

experiment. In recent years, the need for high order molecular multiple moments

has become pressing due to the rapid development in molecular dynamics simulations

[1], computation of high order molecular polarizabilities [2-4], molecular optical and

magnetic properties [5-8], intermolecular force theory in terms of the London series

[4, 7, 9-13], and the quantum mechanical solvent effect in the formalism of the self-

consistent reaction field theory [14-16].

The experimentalist has not been able to meet this need, and it is reasonable to believe

this situation will remain unchanged for some time to come. To date quadrupole moments


are the highest-order multiple moments actually measured using a method that was

established thirty-five years ago [17]. These measurements still contain great experimental

errors, and in some cases even the sign is uncertain. Because multiple moments

traditionally bear named units, there has been a great incentive for the measurement

of high order moments! Yet even with this incentive there has been no breakthroughs.

On the theoretical "front" the challenge posed by calculating high-order multiple

moments is taken with enthusiasm and varied sophistication. After half a century of

development modem quantum mechanical calculation has finally come of age in attacking

a variety of problems that are unsolvable experimentally. With computational resources

readily available ranging from personal computer to massively parallel machines one can

nowadays use well established and tested procedures to solve the Schrodinger equation

for a given molecule or molecular system which gives an understanding of the underlying

electronic structure in terms of molecular orbitals, which in turn, can be used to calculate

desired electronic properties. By nature, molecular multiple moments are one-electron

properties of a molecule. In theory, they can be calculated to any order we want,

but current calculations have not exceeded the tenth order due to basis set reliability.

Modern electronic structure theory is formulated to solve the Schrodinger equation in

a basis set of finite dimension, which usually is expressed as linear combination of

atomic-like orbitals. There are two kinds of prevailing atomic orbitals in the current

stage of quantum mechanical calculations: Slater-type orbital (STO) with form e-ar

and Gaussian-type Orbital (GTO) with the form e -'r2. A molecular orbital description

based on an STO or GTO basis can often give molecular dipole moments with greater

precision than quadrupole moments because dipole moments are by and large caused by

electronic density near the molecular region. For moments higher than quadrupole, it

is known that basis set choice is very important [18], because high moments are very

sensitive to electronic distributions in the outer molecular region, regions which are not

important in the said variation of the energy calculated from the Schr6dinger equation.

Basis set effects are found to be far more important than correlation effect on these high

moments [19, 20]. However, most multipole-moment calculations are carried out using

Gaussian basis set; in fact we are not aware of any recent STO basis set calculations

on moments higher than quadrupole for polyatomic molecules. It is well understood

that the proper choice for molecular electric property calculations is STO, because STO

has the correct long range behavior. But the use of STO in solving the Schrodinger

equation is very difficult, and even today represents an unsolved problem for certain type

of integrals. In practice in approximate and semiempirical models not all STO integrals

are evaluated; rather some of them are parameterized by experimental data when they can

be tangibly related to observed physical properties. Today, however, all ab initio models

use GTOs in place of STOs to evaluate all the integrals necessary. The use of GTOs in

molecular orbital theory is a practical choice to facilitate integration, but not a desirable

one. Multipole moment integrals are conveniently worked out in GTOs by Matsuoka

[21], Augspurger and Dykstra [22], and Mikkelsen et al. [15]. A study of these moments

using numerical techniques can be found in the calculations of Pyykko [20]. High moment

integrals over STO are not well studied, although they would not appear to be particularly

intractable. Originally molecular orbital theory of polyatomic molecules was formulated

using STO basis sets and a series of now classical work on STO integrals was carried

out by Mulliken, Rieke, Orloff and Orloff [23], Roothaan [24, 25], Ruedenberg [26],

and Ruedenberg, Roothaan and Jaunzemis [27]. Local dipole moment integrals were

formulated by Hamilton using C functions [28]. Quadrupole moment integrals were

included in the study by Wahl, Cade, and Roothaan [29], and implemented by Bagus,

Liu, McLean and Yoshimine [30] for linear molecules. A specific and detailed treatment

of dipole moment integrals for polyatomic molecules was presented by Rein, Clarke and

Harris [31], which was extended to handle quadrupole moment integrals by Rabinowitz

and Rein [32], and octopole moment integrals by Swissler and Rein [33]. We have not

found studies of higher moment integrals over STO basis set for polyatomic molecules;

Rather, such studies are carried out with GTO basis set. The popular use of Gaussian

basis set here has an inherent drawback: multiple moment description of molecular

phenomena is a good approximation to the underlying electrostatic interaction only at

long range, but at long range Gaussian orbitals are known to fall too sharply. For this

reason there have been ab initio calculations using Gaussian basis set that are then fitted

to STOs to study long-range interactions [34].

An important interaction is the electrostatic interaction between a solute molecule and

its solvent. It is important as chemistry seldom takes place in a vacuum while quantum

chemistry usually considers only the vacuum environment. The theoretical treatment of

the thermodynamic and kinetic solvent effects can be undertaken by means of force-

field methods [35], wherein a classical Hamiltonian is used to treat both solute and

solvent molecules in the framework of molecular dynamics (MD) and Monte Carlo (MC)

simulations. However, these methods do not take account of the mutual solute-solvent

polarization, which is accounted for by quantum mechanical (QM) methods to allow

charge redistributions. There are three major approaches in the QM calculations of the

solvent effect: the supermolecular method [36, 37] the mixed QM- molecular mechanics

method (QM-MM) [38, 39], and the self-consistent reaction field (SCRF) method [40-47].

All three treat the solute molecule quantum mechanically, and the differences lie in the

description of the solvent. In the supermolecule model, the solvent molecules are also

described at the QM level, but only a few of them can be included explicitly. In the QM-

MM approach the solvent molecules are treated as classical entities with point dipoles or

partial atomic charges, and the dynamics of the solute-solvent system is accounted for by

MD or MC calculations. Finally, the SCRF model considers the solvent to be a continuum

dielectric medium, which reacts against the solute charge distribution generating a reaction

field, which is then included as a perturbation term in the vacuum Hamiltonian of a solute

molecule. The SCRF model, although it lacks microscopic representation of the solvent

molecules, is an inexpensive method to obtain solute molecular orbitals in solution, and

therefore is an excellent approach to examine spectroscopic shift and molecular property

changes of a chromophore in solution, and indeed, from the experience in this laboratory,

these examinations reveal a great deal about molecular electronic structure in solution.

However, the implementation of the SCRF model in our studies has not been complete

enough to study a broader range of molecules, particularly those that do not possess

permanent dipole and quadrupole moments. The impediment, of course, lies in the lack

of reliable computation of accurate high order multiple moments that describe the solute-

solvent interaction in the perturbation of the vacuum Hamiltonian and the conceptual

difficulty in the high order formalism.

Solvation has received a great deal of interest in the quantum chemistry community,

and some vigorous research to tackle the associated problems has been done in this group

over the years [48, 16, 49]. Hard pressed by the growing needs in the course of this

effort, the work presented in this thesis is carried out and it proves to be quite successful

in addressing the problem. Chapter 2 establishes the method of computing multiple

moment integrals necessary in this study [50]. Chapter 3 surveys the conventions and

presents calculational results of molecular multiple moments. Chapter 4 recaps the

essence of SCRF and electronic spectra theory, and gives the UV-Vis spectral shift in

solution as calculated by our method. And, finally, Chapter 5 gives the conclusion and

future outlook of this work.



As will be surveyed in Chapter 3, there exist definitions of electric multiple moment

in the literature. Nevertheless, and as confusing as it may seem, these different forms

are related to each other! Let us then choose one of them

MI,m = rt[ 4 1]/2Yt,m((0,q) (2.1)

where Yi,m(O, 0) is a spherical harmonics. Following Rabinowitz and Rein [32], it is

now convenient to define a primitive multiple moment operator as

Mi,j,k = xi jzk (2.2)

with its expectation values giving all the integrals necessary to compute the mutipole

moments defined in equation (2.1). The computation of the nuclear contribution to the

multiple moments is straightforward, simply replacing the electronic coordinates x, y,

z in Mi, j, k with the nuclear ones X, Y, Z. We therefore consider only the electronic

contribution in this chapter.

The integration of a primitive multiple moment operator of equation (2.2) over a

pair of chosen basis functions is in essence a one-electron, and at most, three-center


computation. A method of evaluating integrals of this nature is generally outlined by

Harris and Michels [51], Steinborn and Ruedenberg [52], and Guseinov and Sadichov

[53]. These traditional approaches all invoke translational, rotational and combinational

properties of spherical harmonics of various kinds that are nicely examined in the works

of Steinborn and Ruedenberg [54], and Steinborn [55]. This chapter presents a novel

method of evaluating multipole-moment integrals over STO basis functions. We will

show explicit formulas for up to hexadecapole moment in a basis set that includes up

to g orbitals. The method, although somewhat tedious, is completely generalizable and

results in very fast computer code.

Preliminary Considerations

Cartesians and Tesseral Harmonics

An STO is usually given by

(22n+1 1/2
Xn,t,m = [(2n)! rn-le-r7(9, 4) (2.3)

where Yi (0, 0) are the normalized real spherical harmonics [56], and the square root

factor is the normalization constant for the radial component. We regroup these terms as

( 2n+- 1/2
Xn,l,m = (2n)! rn-l-le-(r[r Y (0, )] (2.4)

which allows a convenient representation of the angular part in terms of Cartesian

components. The quantities in the second bracket are know as Tesseral harmonics [57],

which we note as

i [=r' (,)] (2.5)

where 1 has its usual meaning ( 1 = 0, 1, 2, 3 ... for s, p, d, f, ... ) and m enumerates

the real components obtained from adding Y/m -m.

We define a product of the powers of x, y and z to be a Cartesian,

(i, j,k) = xiy zk (2.6)

where i, j, and k are integers. In this notation the normalized Tesseral harmonics are

written below up to g type [58],

s harmonic:

s = v 0, ).

p harmonics:

P, = (1, 0, 0), py= (0,1, Pz = (0, 0,1).

d harmonics:

dz2 = J[2(0, 0, 2) (2,0,0) (0,2,0)],

d2-y2 = -[(2, 0,0) (0,2,0)],

d -= 4(1,1, 0), dxa = T (1,0,1), dyz = (0, 1,1).

f harmonics:


fz3 = [T ( 0,0,3) (2,0,1) (0,2,1)],

fxz2 = [4(1, 0, 2) (3,0,0) (1,2,0)],

fyz2 = [4(0, 1, 2) (2,1,0) (0,3,0)],

fz(2-y2) = -[(2, 0,1) (0,2,1)], fxyz = -(1, 1, 1),

f13-3xy3 = 2[(3,0,0) 3(1,2,0)],

f3x2y-Y3 = [3(2, 1, 0) (0, 3, 0)].

g harmonics:

gz = [(4, 0, 0) + 2(2, 2, 0) + (0, 4, 0) 8(2, 0,2) 8(0,2,2) + (0, 0, 4)],

gzXz = (1, 0, 3)- (3,0,1) (1,2,1)],

9y3 = (0, -1, 3) (2,1, 1) (0,3,1)],

gz2(SW-y2) = [-(4, 0, 0) + (0, 4, 0) + 6(2, 0, 2) 6(0, 2, 2)],

gxyZ = [6(1, 1,2) (3,1,0) (1,3, 0)],

z(X3-3y3) = 70[(3,0,1) -3(1,2,1)],

gz(32y-y,3) = j [3(2, 1, 1) (0,3, 1)],

2y2 = [(4, 0, 0) 6(2, 2,0) + (0,4,0)],

zy(X2-y2)= i/5[(3, 1,0) (1, 3,0)].

Coordinate Systems and The Rotational Matrix

The evaluation of STO integrals is usually first carried out locally regardless of the

orientation of molecular input that defines a molecular frame coordinate, and then through

rotational transform the integral value is obtained in the coordinate of the molecular frame.

The two atomic basis functions are usually defined in the individual atomic coordinate

systems. The relation between two individual atomic coordinate systems is established

through a third coordinate system centered in between the two atoms. The explicit

formulas to be established form the basis of the transformation from a local integral

value to that in the molecular frame.

The STO basis functions are generally associated with centers (usually on nuclei)

denoted as A, B, ... The center of a primitive moment operator of equation (2.6) is

denoted as 0. Lower-case subscripts a, b, and o specify the vector components measured

from the centers A, B and 0, respectively. Unlike Cartesian Gaussians, STO two center

integrals are generally evaluated over a local system (see Figure 1) with the axis between

the two centers defining the local Z, Za, and Zb axes, with the local X, Xa, Xb, and Y,

Ya, Yb axes parallel. Here the upper-case notations X, Y, Z represent axes in the local

system. The Za and Zb axes are either pointing at each other, or are directed in the same

direction. We adapt the latter, as shown in Figure 1.

With the origin of the Cartesian coordinate system [X, Y, Z] placed halfway between

the two basis function centers A and B, the local system described above is used to define

a prolated spheroidal, or, ellipsoidal coordinate system focused on A and B, respectively

[52]. Adapting the definition and notations of reference [52], we list some important

relations below for an ellipsoidal coordinate system,

ra = R +q), rb = R( v7),
R I t
2 2

Xa = b 2 1)(- 2)cos

Ya = Yb = 2 1)(1 -2)sinq, (6)

za -=-(l +-4), zb-- (1-- ),
2 2

dT = ( r2)d^dydo.

The integrals evaluated in a local ellipsoidal coordinate system then need to be back

transformed into a laboratory, or molecular based coordinate system, which we define

to be the molecular frame [x', y', z']. The rotation necessary for this transformation is

based on the directional cosine and sine, (see Figure 1)

cosO = --A
sine =
X/Ax'2 + Ay'2

and the rotation matrix R(0,0) can be obtained in two steps:

i) rotate about the y axis to align the z axis: R, (0),
X"1 (X cos6 0 sin0 X
Y"| =Ri(0)Y = 0 1 0 ) (8)
Z" Z \-sin0 0 cosO Z
ii) rotate about the z" axis to align the x" and y" axes: R2(0),

= R2() Y" = sing coso 00 (9)
z' ZI z" 0 0 z"
the final rotation matrix based on matrix equations ((8)) and ((9)) is therefore
(cosOcos4 -sine sinOcosq
R(0, ) -= R2())Ri(0) = cosOsin cosq sinOcosq (10)
k -sinO 0 cosO /
Translation of Centers

The most common way to proceed from here is, assuming the moment definition of

equation (2.1), to evaluate the integrals such as

I= dAYm, (Oa,a)Ym," (0o,4o)YI,' (Ob, b) (2.12)

and to translate the center 0 to the center A or B with the lowest I value, say, center B

with I"'. Making use of the translational property of spherical harmonics [44] and Gaunt

coefficients [45] this yields an integral with sigma ( m = 0 ), pi ( m = 1 ), delta ( m

= 2 ), phi ( m = 3 ), etc., components. The maximum m component required is then

the lesser of I" + 1'" or 1', i.e.,
min(l" +1'",l')
I= E Rm Im = RI + (R, + Rr,')I, + (R6 + R6)I6 +... (2.13)


where the Rm's are the factors due to rotating the integral back to the molecular frame,

and the Im's are the appropriate integrals evaluated in the local frame. Note that Rm =

Rm(A)m(B) each center must be rotated, and there are two pi components, two delta

components, etc.

In this study we will translate the angular components of both the operator and the

orbital centered on A to B. By doing so, center A will only have a sigma component and

equation (2.13) will have the simple appearance of

I = R,(B)1 (2.14)

We will show that the evaluation of I, and the rotation RU(A)Ra(B) = Ra(B) are very

simple. The use of equation (2.14) rather than equation (2.13), and its simplicity, are

at the heart of this work.

Considering a general laboratory coordinate system [x, y, z] and a particular one [xa,

ya, Za] centered at A(ax, ay, az), we have

Xa = x ax

Ya = y ay (12)

Za = Z -az

to translate the system [Xa, Ya, Za] to another one [Xb, Yb, Zb] centered at B(bx, by, bz)

we write

Xa = Xb + (bx ax) = Xb + bax

Xa = Yb + (by ay) = Yb + bay (2.16)

Ya = Zb + (bz az) = Zb + baz
An arbitrary Cartesian defined in the previous section is therefore translated from center

A to center B by way of

i j k
i = k E( D ) -rr j)I s (k) tt
r=0 s=0 t=0
SCr,s,t (bax, ay,baz) xyyz
The translations in terms of equation (2.17) enable a multiple moment integral, apart

from a normalization factor, to assume the following form,
na-la-1 ji J yk j .ra. I m n nb-lb-1 e f g -(brb
(MI,m,n) = (rXaxaya ze xO zO rb 6 b /

(i,j,k) (l,m,n)
= Cr,s,t(bax,bay,baz) Cu,v,w(box, boy,boz) (2.18)
(r,s,t)=(0,0,0) (u,v,w)=(0,0,0)

r (rna-I.-1 o -(ar |Ir. t-lb-1 .e+r+u /f+s+v e-(br,
"a b 1b b /
where center O(Ox, Oy, Oz) places the origin with which the multiple moment integral

is evaluated.

In this development it is clear that a moment integral is a combination of overlap

integrals, which, in terms of equation (2.18), are always between an s orbital and a higher

order orbital. The translation to achieve this is straightforward. However, the sixfold sum

that appears in equation (2.18), if directly implemented, is computationally inefficient,

as was first experienced by Taketa, Huzinaga and O-Ohata in their work on Gaussian

integrals [46]. For this reason we translate an STO centered at A(xa, ya, Za) and the

moment operator centered on O(xo, Yo, Zo) separately onto the center B(Xb, Yb, Zb) using

two sets of predetermined formulas. With the convention that all left-hand quantities in

a formula are centered on A, and those on the right-hand side centered on B we list

these formulas as follows:

Master formulas for STO translation:

Im = E (a,bay, baz)(i, j, k)B (2.19)

Px = [ba (O, 0, 0,) + (1,0,0)] (2.20)

py = [ay (0, 0, 0) + (0,1,0)] (2.21)

Pz = 4f[az(0, 0, 0) + (0, 0,1)] (2.22)
-2 2 -2
dz2 = [(2baz -b ba)(, 0, 0) 2ba(1, 0,0) 2ba(0,1,0)+
4ba(O, 0, 1) (2,0,0) (0,2,0) + 2(0,0,2)]

d2_y2 =Vi [(ba bay)(0, 0, 0) 2ay(0, 1, 0) + 2bax (1, 0, 0)+

(2,0,0) (0,2,0)]

d-y = [baxbay(0, 0, 0) + bay(1, 0,0) + bax(O, 1, 0) + (1,1, 0)]
V 47r
d= --[babaz(0, 0, 0)+baz(1, 0, 0) +bax(0,0,1) + (1,0,1)]

dyz = L [-ayba (0, 0, 0) + ba z(0, 1, 0) + bay(0, 0, 1) + (0, 1, 1)]
V 47r
f = [2baz 3ba baz 3bbaz) (0, 0,0) 6baxbaz(1, 0, 0)-

6baybaz(0, 1,0) + (6ba 3bTa 3ba) (0, 0,1) 3ba(2, 0, 0)-

6ba- (1, 0, 1) 3ba z(0, 2, 0) 6bay (0, 1, 1) + 6ba,(0, 0, 2)-

3(2,0,1)- 3(0,2,1) +2(0,0,3)]

fxz2 = -2 -[ (4baxba bax aTba y)(0, 0, 0)+

4ba 3ba2 bay) (1, 0, 0) 2babay(0, 1,0) + 8baxbaz(0, 0,1)-

3baz(2, 0,0) 2ba,(1, 1,0) bax(0, 2,0) + 8baz(1, 0, 1)+

4ax (0, 0, 2) (3,0,0) (1, 2,0) + 4(1,0,2)]
21 2 -2 F
fyz2= ~ [4bayba- ay ) (0, 0, 0) 2baxbay(1, 0, 0)+

-2 -2 -2
8baybaz(0, 0,1) + (4ba2 ba2 3ba2)(0, 1,0) bay (2, 0,0)-

2bax(1, 1,0) 3bay(0,2,0) + 8baz(0,1,1) + 4bay(0,0,2)-


(2,1,0) (0,3,0) + 4(0, 1,2)]






z(x2-2) = V[ba ba (0,0,0) + 2a, (1, 0, 0)-

247baba(0, 1, 0) + ba y bda (0, 0, 1) + b2(2, 0, 0)+

2bar (1, 0, 1) baz(0, 2, 0) 2bay(0, 1, 1) + (2, 0, 1) (0, 2, 1)]

f/yz = 4l [ya -ba a(0,0,0) + baybaz(1,0,0) + barbaz(0, 1,0)+
V 4r

baxbay(0, 0, 1) + baz(1, 1, 0) + bay(1, 0, 1) + ba (0, 1, 1) + (1, 1, 1)]

fx3_3-y= 3- [(ba3 b ab- )(0, 0, 0) + 3(ba2 a,) (1, 0, 0)-

6baxbay(0, 1, 0) + 3bax(2, 0, 0) 6bay(1, 1, 0)-

3ba-(0, 2, 0) + (3,0,0) 3(1, 2,0)]

35 2 3
f3x2y-y3 = 2[(3a ay a, )(0, 0, 0) + 6bax(1, 0, 0)+

3(bTa ay )(0,1,0)+3bay(2,0,0) + 6bax(1,1,0)-

3bay(0, 2, 0) + 3(2, 1, 0) (0, 3, 0)]





3 -4 -2-2 --4 -2-2 -4 -2-2
gz4 = 16 I {(3ba, + 6ba bay + 3bay 24baybaz + 8baz 24ba8baz)(0,0,0)

-3 -2 2 0,0) + 12(bay +a -- 2
+ (12ba, + 2bazba 48baxbaz)(1, 0, 0) + 12(baxbay + bay 4ba )0 1, 0)

+ 16 (2 3bba a (0,1)
16 (2ba 3ba baz 3ba baz (0, 0, 1)

(3b ba 4 (2,0,0)
+ 6 (3bax + bay 4bz (2, 0, 0)

+ 24baxbay(1, 1, 0) 96baxbaz(1, 0, 1) + 6 (ba + 3ba 4ba (0, 2, 0)

- 96bayba(0, 1, 1) 24( + ba 2a)(0, 0,2) + 12 [(3, 0, 0) + (1, 2, 0)

- 4(1, 0, 2)] + 12bay[(2, 1, 0) + (0, 3, 0) 4(0, 1, 2)] 16baz [3(2, 0, 1)

+ 3(0,2,1) 2(0,0,3)] + 3(4,0,0) + 6(2, 2,0) 24(2,0,2) + 3(0,4,0)

- 24(0, 2, 2) + 8(0,0,4)}

3 10 --3- -- 2 3 -2 2
xz3 = -{ (4baxbaz 3ba-baz 3baxbaaz) (0, 0, 0) + 4ba babaz ba

(1,0,0) 6baxbayba,(0, 1, 0) + 3 (4baxba ba -ba (0, 0,1)

-2 -2
+ 3(4ba, 3bax bay)(1, 0, 1) 3baxbaz[3(2, 0, 0) + (0, 2, 0) 4(0, 0, 2)]

6bayba,(1, 1,0) 6baxbay(0, 1, 1) + 3ba [4(1, 0, 2) (3,0,0) (1,2,0)]

+ bax[4(0, 0, 3) 3(0,2,1) 9(2,0,1)]- 6bay(1, 1, 1) 3(3,0,1)

-3(1,2,1) + 4(1,0,3)}

3 10,
9a = { LO

(baba baba babba (0,0,0) + ba baba 3ba2ba
(4ayba~z 3baybTaz 3ba x aybaz) (0, 0, 0) + (4ba z 9ba'ybaz 3Ta xTz)

(0,1,0) 6barbaybaz(1, 0,0) + 3 (4bayba ba ba (0,0,1)

-2 -2 2
+ 3(4baz 3bay ba)(0, 1, 1) 3baybaz[(2,0, 0) + 3(0,2,0) 4(0,0,2)]

- 6baxbaz(1, 1,0) 6baxbay(1, 0,1) + 3ba,[4(0, 1, 2) (0,3,0) (2, 1,0)]

+ bay[4(0, 0, 3) 3(2, 0, 1) 9(0, 2, 1)] 6bT (1, 1, 1) 3(0, 3, 1)

- 3(2,1, 1) + 4(0, 1,3)}

gz2(x2-y2) --

V + 6a 6-,a) (0, 0, 0) + 4 (3baz a)(1, 0,0)
(-- -37-2r2\2-

4(3baybaz ba (0, 1, 0) + 12 baba, babaz) (0, 0, 1) 6(ba2 baz)

* (2, 0, 0) + 24baxbaz(1, 0, 1) + 6(ba, ba)(0, 2, 0) 24baybaz(0, 1, 1)

+ 6(ba2 bay)(0, 0, 2) + 4bax[3(1, 0, 2) (3,0,0)] + 12baz[(2, 0, 1) (0,2,1)]

+ 4bay [(0, 3, 0) 3(0,1, 2)] (4, 0, 0) + 6(2, 0, 2) + (0, 4, 0) 6(0, 2, 2)}

3/5 --2 -3 -b-3
9gxyz2 = :- {(6baxbaybaz2 bax ba baxba)(0, 0, 0) +

yba 3ba2ba -
(Ty- 3baxbay Ty

*(1,0,0) + (6babaz ba3 3baxbaY (0,1,0) + 12barbayba,(0, 0,1)

+3(2ba a2 -ay)(1, 1, 0) 3baxbay[(2, 0, 0) + (0, 2, 0) 2(0, 0, 2)]

+ 12baybaz(1, 0, 1) + 12baxbaz(0, 1,1) bay[(3, 0,0) + 3(1, 2,0) 6(1,0, 2)]

- ax[3(2, 1,0) + (0, 3,0) 6(0,1,2)] + 12Fz(1, 1, 1) (3, 1,0) (1,3,0)


3 70 -3 -2- -2- -2
gz(x3-3xy2) = {(ba(baz 3baxbaybaz)(0,0,0) + 3 ba6baz ay az)(1,0,0)

6baxbaybaz(0, 1, 0) + (b 3baxbay)(0, 0,1) 6baybaz(1, 1, 0)

+ 3baxbaz [(2, 0, 0) (0, 2, 0)] + 3 (ba bay) (1, 0, 1) 6babay (0, 1, 1)

+ ba [(3, 0, 0) 3(1,2,0)] + 3ba [(2, 0, 1) (0, 2, 1)] 6ay(1, 1, 1)

+ (3,0,1) 3(1, 2, 1)}

gz(3x2y-y3) =

3 V0
8 7rr

{(3babaybaz ba3baz)(0, 0, 0) + 3 bab bayaz) (0, 1, 0)

+ 6baxbaybaz(1, 0, 0) + (3baxbay ba) (0,0, 1) + 6ba6baz(1, 1,0)

+ 3baybaz [(2, 0, 0) (0, 2, 0)] + 3 (ax ba) (0, 1, 1) 6baxbay(1, 0, 1)

+ baz [3(2, 1,0) -(0,3,0)] + 3bay [(2, 0, 1) -(0,2,1)] + 6-ax(1, 1,1)

- (0,3,1) + 3(2, 1, 1)}

3 35 -4 -2-2 --4 b 3)2
gx2 2 = -{(ba- 6bax ba+ bay)(0, 0, 0) + 4a 3baba)(1,0,0)

+ 4(ba 3aay(0, 1,0) + 6ba a[(2,0,0) (0,2,0)] 24baxbay(1, 1,0)

+ 4bax[(3, 0, 0) 3(1, 2, 0)] + 4ba [(0, 3, 0) 3(2, 1, 0)] + (4, 0, 0) 6(2, 2, 0)

+ (0,4,0)}
3 35 -3-- --2- -3 (1, 0O+
9zy(x -y=) {(baba {babba)(0, 0,0) + (3ba bay ba) (1,0,0)+

-a 3ba- b (0, 1,0) + 3baxbay[(2, 0, 0) (0, 2,0)] + 3 (a ay) (1,1, 0)

+ Tay[(3,0,0) 3(1, 2, 0)] + az[3(2, 1, 0) -

(0, 3,0)] + (3,1,0) (1,3,0)}

Master formulas for operator translation:

(n, l,m)0 = f(bo ,b,booz) (i, j, k)

(1,0,0) = o, (O0, 0,) + (1,0,0)

(0, 1, 0) = oy (0, 0, 0) + (0, 1, 0)

(0,0,1) = oz(O0,0,0) + (0,0,1)
(2,0,0) = o(0, 0,-2 ) + 2o(1,0,0) + (2,0,0)
(2, 0, 0) = box (O,0,0) + 2box (1,o0,o) (2, 0,0)

(1, 1, 0) = boXboy(0, 0, 0) + boy(1, 0, 0) + box(0, 1, 0) + (1, 1, 0)

(1,0,1) = boTbo,(0,0,0) + bo,(1, 0,0) + bo(0, 0,1) + (1,0,1)

(0,2,0) = bo,(0,0,0) + 2Yoy(0,1, 0) + (0, 2, 0)

(0,1,1) = boybo,(, 0,0) + YO(0,1,0) + boy(0,0,1) + (0,1,1)
(0,0,2) = o(O,(0,00) + 2Foz(0,0,1) + (0,0,2)
--3 --2
(3,0,0) = b(0,0,0) + 3bo(1,0,0) + 3bo(2,0,0) + (3,0,0)
(1,2,0) = b-obo (0, 0,0) + bo-(1, 0,0) + 2boboy(0, 1,0) + 2boy(1, 1,0)

+ bo(0,2,0)+ (1,2,0)
(1,0,2) =boboz(0, 0,0) + bo(1, 0,0) + 2bobo(0, 0,1) + 2boz(1, 0,1)

+ bo(0,0,2) + (1,0,2)
-2 2
(2,1,0) = boboy(0, 0,0) + 2boboy(1, 0,0) + bo(0, 1,0) + 2bo (1, 1,0)

+ boy(2,0,0) + (2,1,0)
(2,0,1) = bozboz(0, 0,0) + 2bomboz(1, 0,0) + bo2(0, 0, 1) + 2Lo,(1, 0,1)

+6oz(2,0,0) + (2,0,1)
(1, 1,1) = boboboy(0, 0,0) + boyboz(1, 0,0) + bo.boz(0, 1,0) + bo.boy(0, 0,1)

+ bo(1, 1,0) +boy(1,0,1) +bo(0, 1,1) + (1,1,1)
(0,3,0) = bo(0, 0,0) + 3o2(0,1,0) + 3oy(0, 2,0) + (0,3,0)
(0,1,2) = =oybobo(0, 0,0) + b-o(0, 1,0) + 2boyboz(0,0,1) + 2oz(0,1,1)

+6oy(0,0,2) + (0,1,2)
(0,2,1) = boyo(0, 0,0) + 2boyboz(0, 1,0) + b-o(0, 0, 1) + 2oy,(0,1,1)

+ oz(0,2,0) + (0,2,1)
--3 --2
(0,0,3) = bo3(0, 0,0) + 3bo-(0, 0, 1) + 3oz(0, 0, 2) + (0,0,3)

(4,0,0) = O (0,0,0) + 4o3(1,0,0) + 6o2(2,0,0) + 4o(3,0,0) + (4,0,0)

-3 -2- -3
(3,1,0) = bomboy(0,0,0) + 3boYboy(1, 0,0) + boz(0,1,0) + 3boboy(2,0,0)
+3Fbo(1,1,0) + oy,(3,0,0)+3boz(2,1,0)+ (3,1,0)
(3,0,1) = boYboz(0, 0,0) + 3boboz(1, 0,0) + bo(0, 0,1) + 3bo (1, 0,1)

+ 3bo.boz(2, 0,0) + (3, 0,0) + 3bo(2,0,1) + (3,0,1)
(2,2,0) = o 2 (0,0,0) + 2o., (1, 0,0) + 2bo2bo-(0,1,0) + -o2(2,0,0)
+4bobo(1,1,0)+bo(0, 2,0) + 2bo (1, 2,0) + 2boy (2, 1,0) + (2, 2,0)
-2 -2- -2
(2,1,1) = boxboyboz(0, 0,0) + 2bomboyboz(1, 0,0) + bobo(, 1,0) + boboy(0, 0,1)
+ boyboz(2,0,0) + 2boboz(1, 1,0) + 2boboy(1, 0,1) + bo(0,1,1)

+ bo(2,1,0) + boy(2,0, 1) + 2bo.(1,1, 1) + (2,1,1)
-2-2 -2 -2-
(2,0,2) = boboz(0, 0,0) + 2-bobo(1, 0,0) + 2bo2boz(0, 0,1) + 4boFbo,(1, 0,1)

+ boz(2, 0,0) + b-o(0, 0,2) + 2bo-,(1, 0,2)+ 2Fbo(2, 0,1) + (2,0,2)
(1,3,0) = bobo(0, 0,0) + 3bobo(0,1,0) + bo(1, 0,0) + 3boboy(0, 2,0)

+ 3o- (1,1,0)+ Fo.(0,3,0) + 3boy(1, 2,0) + (1,3,0)
(1,2,1) = =o2bo- bo(0, 0,0) + 2boboybo(0,1,0) bo 0) + bo) b (0,0,1)

+ bobo(0, 2,0) + 2boyboz(, 1,,) + 2boboy(0, 1,1) + Foy(1, 0,1)

+ boz(1, 2, 0) + bo(0, 2, 1) + 2boy(1, 1, 1) + (1,2, 1)
(,1 )-2 -- ,2 2- -2
(1,1,2) = boboboO, ) boyboboz(1, 0,0) + bo+bo(0, 1,0) + 2boxboyboz(0, 0,1)
+ bo,(1,1,0) + 2boyboz(1, 0,1) + 2bo.boz(0, 1,1) + boboy(0, 0, 2)

+ boy(1,0,2)+ 2bo(1,1,1) + bo,(0, 1,2)+ (1,1,2)
-3(1,0,3) = 0,0) (1, 0,0) 31) + 3 0,1)
(1,0,3) == boxbo(0, 0,0) bo(1,0,0) + 3bombo(0,0,1)+3bo(1,0,1)

+3bo(bo(0, 0,2) + 3bo(1,0,2)+ bo(0,0, 3) + (1,0,3)

(0, 4, 0) = o,(0, 0, 0) + 4bo(0, 1, 0) + 6bo(0, 2, 0) + 4b-oy (0, 3, 0) + (0, 4, 0)
(0,3,1) = oy (0, 0, 0) + 3boyoz (0, 1, 0) + boy (0, 0, 1) + 3boy(0,1,1)

+ 3boxboz(0, 2, 0) + boz,(0, 3, 0) + 3Foy(0, 2,1) + (0, 3,1)
(0,2,2) = boboz(O, 0, 0)1,02) + 2boybo(0, 0, 1) + 4boyboz(0, 1, 1)

+ bo(0, 2, 0) + bo(0, 0, 2) + 2boy(0, 1, 2) + 2Fo (0, 2, 1) + (0, 2, 2)
= 3 -3 -3 1, '-2-2
(0,1,3) = boybo(0, 0,0) + bo(0, ) + 3oboo(0, 0,1) + 3bo(0, 1,1)

+ 3boyboz(0, 0, 2) + 3boz(0, 1,2) + boy(0, 0, 3) + (0, 1, 3)
(0, 0, 4) = o (0, 0,0) + 4bo(1, 0,0) + 6oz (0, 0,2) + 4boz(0, 0,3) + (0,0,4)

With this accomplishment we can now write the sixfold sum of equation (2.18) in

terms of a twofold sum over the two Cartesians, both centered on B, in the following

ala-Va aaea--a XI oyo 1oh I b"ybfzbge I
(MI,m n) = ( zXIr"blb 04ybz e4rb

(rna-tl-1 ga(baxz,ay, baz)x '"tY'" t'" e-"r.

= ~g (bax, ba y, baz) g (box, boy, boz)

na-la-I (- a b Ab yb -b e
where the coefficients g~(bar, bay, baz)and fb (box, boy, boz) are taken from the master

formulas through two computational do loops.

Local Integral Evaluation

One-center Integrals

After the translations described by equation (2.44), all one-center integrals, apart

from a normalization factor, take the form

I n (rn ---.-b 2 ij yzfk e- (a+Cb)rb
c0 7

J fna.+-i.-al+i+j+k e-((a+(b)rdrb, Sini++l ObCosk bdOb
10 (2.45)

J Sini CosJdo
and we proceed to evaluate the non-vanishing angular part of Ii first.

Assuming m and n to be integers, we have, in general

7r,27r 7r,27r
Scosmsin d -m = m- f cosm-2 sijnnkd (2.46)
Sm+n I
0 0
applying equation (2.46) in a recursive manner, we get
7r,2r 7r,27r
Scosmsind = m- M mn -m--3 / cosm-4sinn"do
Ir m+n m+n-2
0 0
m-1 m-3 m-5 f
-- + 2 ------ 4 cost-6 sinn do = *
m+n m+n-2 m+n-4 J
(m- 1)!!n!! !f s.n d
(m + n)!! J ne


In equation (2.47) we need to examine two cases of the sine integral. First, if the

n is even,
7r,27r x,27 r,27r

in n n-2J
0 0 o (2.48)

(n 1)!!

Therefore for even m and even n, equations (2.47) and (2.48) give a closed formulae,

fcmOS sinndnO = (m 1)!!(n- 1)!!r
J (m + n)!!
c in (m 1)!!(n 1)!!
]cosmosinnd- = n (27) = F(m, n) (2.50)
J (m + n)!!
We then examine the sine integrals with odd n in equation (2.47),
i0 2 2-)- (n- I )2(nn-)!1 = -nn-l___n-l__
sin -d = -cos =- (2.51)
n!', n!
f sinnodb = 0 (2.52)
and finally equations (2.47) and (2.51) give, for even m and odd n, the closed formulae

c s (m 1)!!(n 1)!! 2"(T- )!(-1)!
cosmosmodo = (m 1)(n 2 2 (2.53)
J (m + n)!! n!2.53)
All the angular integrals necessary for Ii are given by equations (2.51) and (2.53).

The radial part of Ii is simply

rne-cr = (2.54)

Two-center Integrals

In an ellipsoidal coordinate system we use equation ((6)) noted previously to express

a two-center integral that evolves from the translation of equation (2.44),

12 = (rn-la-1le-(ar alr--l-1 xiy'zke brb

00 1 27r
= d f do ()(R)n+nb-1a-1b+i+J+k+1( +rl)na+1(_ q)lnb+l
1 -1 0

S[(2 )(1 )](i+)(1 ke- ossin (2.55)

00 1
= d dR( )n.+nb-la.-l +i+j+k+l ( + )a+1( )nb+l
1 -1

[(q2 1)(1 2]+)- )ke-P- F(i,j)
p = y(Ca + (b)
S(o (b)
( + C(b)
and F (i, j) is given by equation (2.50).

At this point we feel it is convenient if a new function is introduced as
00 1
Z,0,7,6 (p, ) = fd d7 (+ )a((- r)'[(2-1)(1 q2)(1 -r)6e-p-'P (2.57)
1 -1
which we shall call "Z functions". The recursion formulas of Z are simpler than those

derived by Ruedenberg, Roothaan and Jaunzemis [27]. In this work we adapt an algorithm

proposed by Stevens [59, 60] based on the binomial theorem,

Za,,,6(f,T)= -)
i,j,k,l,m (2.58)

Aa+3+27-i-j-2k+m (p) Bi+j+21+m (p, T)
where the A and B functions are as defined by Mulliken, Rieke, Orloff and Orloff [61],

and have been well studied by, among others, Miller, Gerhauser, Matsen, and Harris [62].

In this formalism, the two-center integral of equation (2.55) is evaluated as the

product of Z and F functions,
12 = (rna-le- ara Xr'-'- lxy ztebrb)
= Zn.-1-,nb-lb,(i+j)/2,k (P, T)F(i,j)
Rotational transform into Molecular Frame

As discussed in section 1.3, all integrals evaluated in section 111.2 need to be back

transformed into the molecular frame [x', y', z'] (see Figure 1) through a matrix rotation
(x) fn r12 7r13 \(Y'
R(90,) = yr21 r22 r23 = y (2.60)
z r3 0 r33 z
where the matrix elements are given by equation ((10)). We now let the rotation described

by equation (2.60) act on an arbitrary Cartesian in the ellipsoidal coordinate system,
Rxi jzk = (rll + r12Y + r13z) (r2x + r22y + r23z) (r31x + r33z)k

= Cr,s,t(rll, r12, 13)4,v,w(r21, r22, r23)Cp,q(r31, 33) (2.61)

*. r+u+P s+v t+w+q

t (r11, r12, r13) = r r8l2r 3

4,,,v'w( r21 72, r23) = !W! 21r'r23 (2.62)

ck, q(r31,r33)= r'3

Considering I2 of equation (2.55), we now have

RI2 = R(r'oe- ra Tn ixZ Yk e-rbCb
bb zb b eb)

= (rn^e-crajrn(Ray bz )e-r^" (.3

where we have applied equation (2.61) and used the rotational equivalence of coordinate

systems [x, y, z] and [xb, Yb, Zb] (see again, Figure 1).

In this approach it is notable that the rotation R(0,0) needs to act only on the right-

hand side of equation (2.55) to rotate the integral 12 into the molecular frame.

Computational Procedure

The method presented here proceeds as follows:

i) For two orbitals on the same center the operator is translated, the symmetry is

checked, and if non-zero, integral evaluation proceeds as outlined two subsections


ii) For orbitals on different centers:

(1) Choose the orbital with bigger I in equation (2.3) to be centered on B.

(2) Translation is made in two do loops of equation (2.44) by using the master


(3) Determine the order L = e + f + g + t,, + t, + t,, + t + t,, + t,v of

the Cartesian on the right hand side of equation (2.44).

(4) Sort out the angular symmetry "allowed" Cartesians of the order L determined

in step (3) by way of the F (i,j) functions shown two subsections ago.

(5) Compute the necessary Z functions [equation (2.58)], and then the local

integrals [equation (2.59)] required by step (4).

(6) Rotate the computed local integrals into the molecular frame by using equation


Based on this procedure, a program package HIMO is produced based on a flow

chart shown in Figure 1.

Discussion and Conclusions

A general method of evaluating electric moment integrals over Slater-type orbitals

is developed. A program that implements the procedure described in the last section has

been in use in this laboratory for some time. We have been using these moments for

studies on solvation using the self-consistent reaction field model, and for examining the

long-range interactions between molecules.

The method we outline is somewhat unusual in that the Cartesian components of one

orbital and the operator are both translated to the center of another orbital and thus only

a single sigma symmetry integration remains between an s orbital on one center and a

complex Cartesian orbital of the form xiyjzkrle-(r on the other. The local evaluation of

this integral is easy, as is the rotation back into the molecular frame. Although the method

is quite general, the price we pay for this computational simplicity is formal complexity.

Generating explicit formula for higher moments, or higher angular momentum atomic-

like orbitals, is tedious. We do not maintain the "symmetry" of the spherical harmonics

and we pay a penalty for this.

A definite advantage of the method described is that we need to consider z component

of one atomic orbital, and the Z function we need [equation (2.57)] is simpler than the

traditional L function [29],

00 1
L`` (p,T) = d d r( + r)( r)(1 + )'r1 )
1 -1 (2.64)

[(2 1 -2)])e-p(--

or the C function [27],

(1 ,a+/3+y+6+2e+1
C"^(P, ) = (bR L6', (p,T) (2.65)


Both L and C functions, if evaluated in the fashion of equation (2.58) require a six-fold

sum, while the evaluation of Z functions needs a five-fold sum, a saving of the innermost

do loop in computation. It would be worthwhile to further investigate the Z function

recursion relations.

Before concluding we further note that the "Cartesian STOs" are prescreened based

on angular symmetry before the integration takes place. Note that i and j in the Cartesian

(i, j, k) that occurs in equation (2.55) must both be even or the integral vanishes. For

example, a typical non-vanishing quadrupole integral between d orbitals contains on

average 18 different terms of 64 possible, but only 5 or 6 are non-zero. In Table 2-1 we

list the number of Cartesians for a given order, and the number of Cartesians that have

non-vanishing angular integral based on the prescreening of F (i,j). This table shows that

in the integration with Cartesian STO's, there is a saving of about 70%.


e Yb

Ya Yb

Za" Zb"

Xa" n"

Ya" Yb"

Za' Zb'

Xa' Xa'

Figure 1: Rotation From Ellipsoidal Coordinate to Molecular Frame






Symmetry zero


Figure 2: Flow chart of the HIMO computational program

Table 2-1: Percentage of Cartesians needed to Compute

Number of
Number of Cartesians with
order B/A
Cartesians (A) Nonzero Angular
Integral (B)

0 1 1 1.00

1 3 1 0.33
2 6 3 0.50
3 10 3 0.30
4 15 6 0.40
5 21 6 0.29
6 28 10 0.36
7 36 10 0.28
8 45 15 0.33
9 55 15 0.27
10 66 21 0.32
11 78 21 0.27
12 91 28 0.31


Definitions of Electric Multipole Moments

There exist many definitions of electric multiple moments in the literature. Very

often this situation causes a great deal of confusion in comparing results of different

workers, therefore caution must be taken here.

Physicists usually define multiple moments as complex quantities. The potential

outside the encompassing sphere of a charge density p (x) can be written as multiple

expansion of complex spherical harmonics:

E(X) = 47 qlm 0, (3.1)
l=0 m=-1

and on the other hand Coulomb's law has

(X) X) d3x' (3.2)

. Now applying the Von Neumann expansion to equation (3.2), we can rewrite equation

(3.1) as

(X) 4rZ 21r 1(0,'' x') '] 1 (3.3)

the coefficient of the multiple expansion of equation (3.1) is therefore defined as

multiple moments,

qi,m= fY*(0 ,' p (x' d 3x' (3.4)

with the first several multiple moments expressed as follows:

Monopole: 1 = 0

1 3X 1
o ,= r 3)r, = q (3.5)

where q is the charge of the density distribution; and

Dipole: 1 = 1

q9l, = --/ Jf ( iy' = (P- ipy)

q1,0 = Pz p( d = 3z' (3.6)

where pi is the i-th dipole component; and

Quadrupole: I = 2

2,2 = ( i x')d3x' = V (Q11 2iQ12 Q22)

q2,1 = 7/J- (z p TT (Q13 iQ23)

q2,0 = 3z'2 -_ r2 pX) d = 33 (3.7)

where Qij is the quadrupole moment tensor

Qij = (3x'.' r'26)p ')d3X' (3.8)

The above type of definition can be found in Jackson [63].

A much simpler definition is given by Bottcher [64], as


M x p (x)d3x' (3.9)


Q = xx p x dx (3.10)


U = X '" p(x d (3.11)

resulting in an expression of the potential that can be written as

= q -M.V +Q: VV- U:VVV1+... (3.12)
r r r r

Although this set of multiple moments may seem to be overly simple, it has

much mathematical convenience when organized in the formation of polytensor in

the work of Applequist [65, 66], and subsequently applied to molecular electrical

property calculations by Liu, et al. [67] and to the electrostatic interaction potential

formulation by Dykstra [68].

A third definition is put forth by Buckingham [69], Kielich [70], Stogryn and Stogryn

[71], and Krishnaji and Prakash [72]. This definition comes out very naturally from the

Taylor series expansion of electrostatic potential of a charge distribution, therefore it is

easy to use in applications and enjoys great popularity among chemists. For convenience

we choose Buckingham's definition of multiple moments, and we carefully establish

it as follows.

We consider a charge density p(r) at point (x, y,z) represented by the vector r from

a chosen origin 0. The electrostatic potential due to this charge density at an arbitrary

point P(X,Y,Z), denoted by a vector R from the origin 0, is written as

S(P) = pf d3r

= fp(r) d3r (3.13)
/(X x)2+ (Y y + (Z z)2

where r is the vector pointing from (x,y,z) to P(X,Y,Z). In the vicinity of 0, where R

> r, one can expand in equation (3.13) to achieve

a 1/a2/r 1 0(1/r') )
(P)6 rop(r)[o r+r + ( (r 99r1 o1 /r r +'" .d
+ 2 0 / 0

= Jp(r){+ (x + y + z)+

+ ) X + (2ZRI;X Y) z2+
-xy + Fxz + 7yz]+

9(5X2Y YR2) X2y + 9(5X2Z ZR2)X2z + 9(5XY2 XR2) xy2

9(5XZ2 XR2)xz2 + 9(5Y2Z ZR2)y2z + 9(5YZ2 YR2)yz2+


4[(105X4 90X2R2 + 9R4)X4 + (105Y4 90Y2R2 + 9R4)y4+

(105Z4 90Z2R2 + 9R4)z4 + 4(105X3y 45XYR12)x3y+

4(105X3Z 45XZR2)x3z + 4(105Xy3 45XYR2)xy 3

4(105XZ3 45XZR2)xz3 + 4(105Y3Z 45YZR2)y3z+

4(105YZ3 45YZR2)yz3 + 6(105X2Y2 15X2R2 15Y2R2 + 3R4)X2y2+

6(105X2Z2 15X2R2 15Z2R2 + 3R4) 2z2 +

6(105Y2Z2 15Y2R2 15Z2R2 + 3R4)y2z2 +

12(105X2YZ 15YZR2) x2yz+

12(105XY2Z 15XZR2)xy2z+

12(105XYZ2 15XYR2)xyz2] + -}d3r (3.14)

To extend equation (3.14) is obviously getting very cumbersome as we proceed to

higher orders. Therefore, to prevent the terms from getting intractable, it is convenient

to define a set of electric multiple moments in the following fashion,


q = p(r)d3r (3.15)


Ax = Jf xp(r)d3r

Py = f yp(r)d3r

z = zP(r)d3r (3.16)


Q = f f (22 y Z2)p()d3r

Qyy = f (2y2-2 2)p(r)d3r

Qzz = f (2z2 -2 y2)p(r)d3r

Qxy = f yp(r)d3r

Qxz = fxzp(r)d3r

Qyz = 3 yzp(r)d3r (3.17)


zxxx = f (2x 3-3xy2 -3 z2)d3r

yyy = f (2y3-3x2y- 3yz2)d3r

zzz = f (2y3 3x2y 3y2z)d3r

zXy f (4X2y y yz2)p(r)dr

Qz= f (4X2Z y2z z3)p(r)d3r

yy = J f (4xy2 X3 z2)p(r)d3r

xyz = fxyzp(r)d3r

zz = f (4xz2 3 xy2)p(r)d3r

yyz i f (4y2z x2z z3)p(r)d3r

=yzz = (4yz X2Y y3)p(r)d3r (3.18)
r x = f [x- 3(x2y2 x2z2) +2 (z4 + y4) + y22]p(r)d3r
rxy = f [ -x3y (3 + yz2)p(r)d3r

x = J [xz (y2z z3)]rdr
xxy= [ (9x2y2 x2z2 2z2) (x4 y4) + z4]p()d3p
1rxyz = [x2yz + y3z)]p(r)d3r
rzz = J [(9x2z2 xy2 yz2) z4+ ) + y]p(r)d

rxyy = j y (yz2 + x3y)]p(r)dZr
Fx!yz = f 1xy2z z + x3z)p(r)d r
S= J [yz2 -(3 3y)] p(r)d3r
rr z = f [x (xy2z + x)]p(r)d
ryyyy = f [y + (z4 + x4) 3(y2z2 x ) + 2]p(r)d3
ryyz = J [yz (yz3 + 2yz)] p(r)d3r
r.yzz = f[ f (9y2-2z 2- 2 2 I 4 ( y4) + lz4]p(r)d3r
rFzzz = f[~z3 15(y3z + 2yz)]p(r)d3

Fzzzz = z4 + 2 (y4 + X4) 3(x2z2 + 2z2) p(r)d3r (3.19)

so that the potential at any point P can now be rewritten as

,p q XI + Yty + Z1z
R R3

XXQXX + YYQyy + ZZQzz + 2(XYQxy + XZQZ + YZQyz) +

R-[XXXf2xx + YYYlyyy + ZZZQ.z + 3(XXYQzy+

XXZQXz, + XYYQxyy + XZZQzz, + YYZQ ,vz + YZZQwY)+


1[XXXxr"XZ + YYYYFyyyy + ZZZZPzzzz + 4(XXXYFTxy+


ZZPYzzz) + 6(XXYYFxxyy + XXZZTxxzz + YYZZvyyzz)+

12(XXYZ Py, + XYYZTFyy, + XYZZxy^zz)+


We shall not dwell on the significance of equation (3.20) until Chapter 4, suffice

it to say at this moment that it clearly reveals the electric multiple moments defined

in equations (3.15, 3.16,3.17, 3.18 and 3.19) to be the coefficients of a Taylor series

expansion. We will use this set of multiple moments throughout the remaining part

of this thesis.

Before closing this section, we note another related definition of electric moment that

surfaces in the literature from time to time, the "ordinal moments", defined as

Mi,j,k(O0) = J xiyJzkp(r)d3r (3.21)

commonly known as first moments when i+j+k = 1; second moments wheni+j+k = 2,

etc. The ordinal moments are the de facto primitive moments defined in Chapter 2.

Although they can clearly be assembled into multiple moments, they are not multiple

moments per se.

The Theory and Calculation of Molecular Orbitals

As a one-electron property, multiple moments are to be calculated for a given

molecular electronic state. However, there is no exact solution of the Schr6dinger

equation for a many-electron molecule. In this section we outline the basic theory and

technique that are used for calculating the electronic structure of a molecule. This begins

with the calculation of the approximate molecular orbitals that will be used to calculate

electronic properties like the molecular multiple moments under study.


1. The Hartree-Fock Theory

The electronic Hamiltonian of a molecule under Born-Oppenheimer approximation

can be written in atomic units (F = m = e = 1) as

nV2 1n n M Z
i i a
where the first term is the electronic kinetic energy operator, and the following sum

accounts for the Coulomb repulsion between pairs of electrons and the Coulomb attraction

between the electrons and nuclei. The Schrbdinger equation of the molecular electronic

system is therefore

HIb) = EIV) (3.23)

As mentioned, it is impossible to solve equation (3.23) for many-electron molecules

exactly. Yet approximations can be made to attain an approximate solution from which

a lot of the chemistry of a molecule can be uncovered. One such approximation is the

Hartree-Fock approximation where a Slater determinant made of independent molecular

spin orbitals {4i(i), i = 1, n} is used to approximate the molecular wavefunction

|I) = 101(1)02(2) ... (n)) (3.24)

For convenience, we restrict ourselves to a closed-shell formalism. There is no loss in

generality, but only, for the moment, in complexity. The energy of the molecular system

from equation (3.23) is
n -2 M Z.
E = (i(1)1 E I, i(1))+
i=1 a=1

1 E (0i(1) 0j(2)1 j[0i(1)0j(2) i(2)j(1)]) (3.25)

n n
= h + 2 E (JI Kii)
i=1 i,j=l
where ,i is a molecular orbital, hii is one-electron matrix element and Ji and Kij are

two-electron Coulomb and exchange integrals, respectively.

The optimal energy E of equation (3.25) can be attained through the variation of the

molecular orbitals while insuring their orthonormality. One minimizes the Lagrangian

L = E E ji[(i 1}j) Sij] (3.26)
where ij is a Lagrange multiplier,

6L = 6{E E Kji[(ilj) 6kil} = 0 (3.27)
to obtain the least E = (?^I//|)/( '|) > Eexact this leads to

n 1
E(j(2)|- 1(1 P12) j(2))}| i(1)) + Complex Conjugate (3.28)

n n
= (60i(1){ E Eji j(1)} + Complex Conjugate
i j

Since the variation is arbitrary, we must have
n n
{hi (1) + ij ki]} i4(1)) = EEji4j|(1))
j j

where Jj and Kj are called Coulomb and exchange operators, and f(1) is the Fock

operator. Since the Lagrange multiplier in this case is Hermitian, and the Fock operator

is invariant under unitary transform, we can choose a set of spin orbitals that enable a

diagonal form of Lagrangian multipliers, so that

f (1) i (1)) = Ei i (1)) (3.30)

Equation (3.30) is called the Hartee-Fock equation, from which molecular orbitals of the

given state [I) are obtained.

2. The LCAO-SCF Equations

For the Hartree-Fock equation (3.30) to have any practical use two problems remain

to be solved: (1) what are the explicit forms of the molecular orbitals? and (2) how is

the Fock operator constructed? The latter depends on the unknown molecular orbitals.

A successful solution to the first problem is the linear combination of atomic orbitals

(LCAO) method established by Roothaan [73], which can be written as

1i) = E ciilX/) (3.31)


where Ix,)'s are atomic orbitals usually centered at the constituent atoms of a molecular

system, and the coefficients of the combination is to be determined from the formalism

derived below.

The Hartree-Fock equation (3.30) in LCAO expansion takes the form

/ Cvlxv) = Ei ECviv) (3.32)
v=1 v=l
"bracketing" both sides of equation (3.32) by a bra (x, I (i.e. multiplying both sides by

X* and integrating), we get
SF,,Cvi = ei SICvi (3.33)
v=1 v=1
or in matrix form,

FC = SCe (3.34)

where F is the Fock matrix, S the overlap matrix, and C the solution of molecular

orbitals. Equation (3.34) is known as the Roothaan equation.

The overlap matrix S is Hermitian, and in the absence of basis set dependence can

be orthonormalized

XfSX = 1 (3.35)

through, for example, the Lowdin symmetric orthogonalization. Consider now a new

coefficient matrix C' related to the old one by

C' = X-1C, C = XC' (3.36)

% P

then Equation (3.34) can be written as

FXC' = SXC's (3.37)


XfFXC' = XtSXC'e (3.38)

or simply,

F'C' = C's (3.39)

This procedure transforms the Roothaan equation into a standard matrix eigenvalue

problem. Upon obtaining the solution of this matrix equation, electrons are assigned

to basis MO's defined by the linear combination coefficient contained in a column of the

C' matrix. Such columns define occupied molecular orbitals. is said to be occupying that

molecular orbital. In a calculation these occupations are of course assigned a priori. In

the case discussed above every has an occupation of two or zero electrons. The formalism

thus developed is called "restricted Hartree-Fock" (RHF) because every pair of a and 3

spin electrons are restricted to the same spatial molecular orbital.

The second problem posed in the beginning of this section leads to a strategy called

Self-Consistent Field (SCF). In an SCF method, a starting guess of the molecular orbitals

is given for computing the Fock matrix, then equation (3.39) is solved for a resulting set

of molecular orbitals, which is then used to estimate a new set of orbitals that are then


used to compute the new Fock matrix to solve for a new set of molecular orbitals, and

this iterative procedure continues cycle after cycle until self-consistency is achieved.

3. The INDO Technique

The Fock matrix elements of equation (3.39) for a restricted Hartree-Fock formalism

can be expressed in the atomic basis as follows,
fo = hJIV + ZECAaC,*[2(po7vA) (yoIAv)]
a aA

= h, + 2 CAaCa[({Lo7I vA) 2 (,iaAv)]
a O' (3.40)

= V + EP\-tK A) -{a\\v)]

= hyv + Gv
where the PA,'s are the first order density matrix elements, defined by

P = 2 E CAaCa (3.41)

and h,, is the one-electron Hamiltonian matrix element, and Gp,v is the two-electron

matrix element composed of Coulomb and exchange integrals. The Fock matrix may be

evaluated in several ways. Two major schools are: ab initio, in which all the integrals

in equation (3.40) are evaluated; and semiempirical, in which some of the integrals are

neglected and then the remaining ones are parametrized to reproduce experimental data.

The work in this thesis falls into the second school.

In our calculation we invoke a basic approximation, known as the Zero Differential

Overlap (ZDO) approximation, to neglect certain integrals,

X0(1)X,(1)dT(1) = 0, p # v (3.42)

whenever this differential is met in an integral over a symmetric operator. Thus the matrix

of atomic overlap integrals becomes the unit matrix. The effect of this approximation

on the one- and two-electron integrals, in terms of the number of such integrals and

the physics they contribute to the theory, is pronounced. In particular, all three- and

four-center Coulomb integrals are removed, and no exchange integrals remain, or

SPAVB = (3.43)

(IpAAIVBUBD) = 6pv~(ptAjI\A) (3.44)

Throughout the development of quantum chemistry the ZDO approximation has been

used to develop theory as well as to establish reliable models. Today there are numerous

ZDO methods in use. They are all, in one way or another, based on the following

approximations, the Complete Neglect of Differential Overlap (CNDO) approximation

[74], the Intermediate Neglect Differential Overlap (INDO) approximation [75], and

the Neglect of Diatomic Differential Overlap (NDDO) approximation [75]. The CNDO

method utilizes the ZDO approximation for every integral evaluation encountered, while

the INDO includes all one-center two-electron integrals and NDDO methods also include

all two-center hybrid integrals of the form (A VA 'AB B) where the subscript refers to

atomic center. All SCF calculations performed within this thesis were carried out at the

INDO level of approximation and so it is the INDO method we will outline below.

Within the INDO approximation the Roothaan equation, relation (3.39), is constructed

and iteratively solved as was previously discussed in the last section. The difference

between INDO and ab-initio Hartree-Fock lies in the type of integral included, and their

functional form, in the construction of the Fock matrix, with two caveats. The first being

that the INDO method described here explicitly uses only the valence atomic orbitals,

and the second being that in the LCAO-MO expansion only one basis function is used

for each atomic orbital. Thus we will be describing, and working with, a valence orbital

minimum basis set method. Below the INDO equations analogous to equation (3.40) are

presented, p, v, a and A again label the atomic basis functions, { I x)}, while A, B,C, ......

label the atomic center the AO basis functions are associated with. The imposition of

the INDO approximation on equation (3.40) leads to different forms for specific cases,

and these are listed below.

1. Case p = v; A = B:


2. Case I 0 v; A = B:



3. Case pI $ v; A # B:



Several new terms were used in expressing equations (3.45) through (3.47), these are

the terms used throughout the literature. We shall explain them below starting with the

one-center core integral UpA A of equation (3.45),

U ==AJ A A- 1 + VAIA})


the one-center electron-nuclear attraction,


the first order atomic density Pcc is defined as,

Pcc = E Ac Ac

the two-center two-electron 7yAc is defined as,



YAC = (IAAc A Ac)


and VAC the two-center electron-nuclear attraction as

VAC = (PAIVCIPA) (3.52)

The p13AVB term is explained as replacing

/hLAVB = (/AI V V V B) (3.53)

There are many established methods to approximate the various terms above. The

details of this subject can be found else where [76, 77]. It should be said the various

integrals are replaced by functions of experimentally obtained parameters such as ioniza-

tion potentials and electron affinities, with the exception of /3 in equation (3.53). The 0f

term is, in many implementations, a totally empirical parameter and is usually obtained

by varying it, while calculating a specific property or reproducing results of a higher

level of theory, until acceptable results are obtained. Such is the case in this work. We

shall obtain the SCF orbitals via two different INDO parametrization, one having been

parameterized to reproduce molecular structures and the second to reproduce the low

energy ultraviolet and visible absorption spectra of molecular systems.

Some Calculated Molecular Multipole Moments

As a one-electron property, a general multiple moment M can be calculated from

quantum mechanics as the expectation value for the state I'b),
M = (I\MIp)

= (01(1)02(2) .. 2N(2N)IM1|1(1)02(2) 02N(2N)) (3.54)


where the density matrix is obtained in this work from the SCF INDO calculation

described in the preceding section, and the moment integrals over the atomic orbital

basis are evaluated by using the method established in Chapter 2.

Multipole moments thus calculated are reported below for a variety of molecules. In

general, the calculated value of a multiple moment tensor depends on which center is

chosen to evaluate the moment integrals. To illustrate this, let us consider a system of

many particles with the i-th particle having mass mi and charge qj located along the z-axis

at a distance zi from the origin 0. The zeroth moments are mi = M and ej = q,

and are the total mass and total charge of the system. The first moments about the origin

0 are mizi = Cz and eizi = pz. If 0 is the center of gravity, Cz would vanish. To
i i
determine the effect on the moments of a change of origin, consider the dipole moment

I/ relative to a new origin 0' at the point -Z. Clearly

S= ei(zi + Z) = iz + qZ (3.55)
so that the dipole moment is independent of the origin only if q, the charge, is zero. The

choice of origin is ideally the center of charge. However, it is usually chosen to be the

center of mass (COM) defined as #c. The reason for this is that negative masses do not

exist, so that no real system can have M = 0 and a COM can always be found, while the

definition of center of charge varies. In molecular system, the center of mass is very close

to the center of charge, as charge is generally proportional to the mass of an atom. The

results reported below are all calculated from COM. Whenever possible, comparisons are

made with either experimental or existing theoretical values, or both. For convenience, a

conversion table is made as Table 3-2 for all the units involved in the following report.


Water is fundamental in chemistry and biology, and it has always received a great deal

of attention both theoretically and experimentally. Calculated electric multiple moments

in this study are listed in Table 3-3. The calculation uses experimental geometry (RO-H

= 0.958 A, LHOH = 104.450, and the oxygen is taken as the coordinate origin, and the

two hydrogens are placed in the XY plane, pointing towards the negative X direction,

where X is the C2 axis (Figure 3 (A)). We have tested the performance of two sets of

the two-electron integral INDO parametrization of 7, one which is fitted to produce good

UV-Vis spectroscopy, and the other which is ab initio calculated -7 used for obtaining

geometry. Within each set of 7, we further tested the effect of neglecting two-center

terms of the moment integrals. From this study we favor the inclusion of the two-center

and the spectroscopic 7. Our results compare well with those reported in the existing


literature. The quadrupole calculation favors the experimental value from the reference

c over that from the reference d. The lack of experimental higher moments makes other

conclusions difficult.


A favorite test to see if a theory can hold water (which our theory just did!) is to see if

it holds benzene. Besides intriguing chemists for many years, it helped establishing many

concepts and calculational models in molecular structure and spectra. Using HF-631G**

optimized geometry with GAMESS,(RC-C=l.386 A, RC-H=l.076 A), the orientation of

the input geometry is shown in Figure 3 (B). Benzene multiple moments are calculated

and tabulated in Table 3-4. Notice that the dipole and octopole moments for benzene both

vanish due to D6h symmetry. We realized that the basis set of Zerner and Gouterman

[78] used in ZINDO is not able to quantitatively reproduce the experimental moments.

Zerner basis set is designed to reproduce geometries, and a better basis set for molecular

moments is recommended by Rein, Nir and Swissler [79]. We use the Rein basis set for

carbon, and the results compares very well with the experimental ones (Table 3-4). The

most noticeable chemistry that can be explained from these results is the herring-shape

structure of the benzene dimer in gas phase. The dimer would be a perfect perpendicular

structure if only the quadrupoles are to be matched, but the hexadecapole interactions


exist as a perturbation to tilt the perpendicular dimer, giving the observed structure, as

shown in Figure 4 (A).


Calculated multiple moments of hexafluorobenzene, with its geometry optimized

using the same procedure as that for benzene (RC-c=1.379 A, Rc-F=1.314 A) are listed

in Table 3-5. The input orientation is shown in Figure 3 (C). The most salient electronic

feature of this molecule is the removal of electronic density along the ring. The calculated

moments reflect this feature, reversing the sign of the moments from that observed for

benzene. The Rein basis set does not have parameters for the fluorine atom. We therefore

extrapolate the fluorine exponents from the existing C, 0, and N exponents, and used

((2s)=2.7 and ((2p)=2.9, which, together with the Rein carbon exponents yields calculated

moments again in very good agreement with experiments (Table 3-5). With the same

argument for benzene dimers, we would predict the parallel docking structure for a

Hexafluorobenzene-Benzene dimer ( see Figure 4 (B)), as is observed in experiment.


Table 3-6 gives calculated quadrupole and hexadecapole moments as calculated for

pyrazine with the orientation shown in Figure 3 (D). The result is obtained for pyrazine

imbedded in water solution, by using Self-consistent Reaction Field Model to be described

in the next Chapter. The calculation uses the Zerner basis, default in ZINDO program.

It is marvelous that in solution the Zemer basis can adjust itself to give quadrupole

moments comparable to those obtained from the very expensive calculation by Zeng,

Woywod, Hush and Reimers [80]. This is very important because we will use default

basis in solution calculations in the next chapter.

Table 3-2: The Conversion Factors for 1 a.u.

Electric Charge (Zeroth Moment)

Dipole/First Moment

Quadrupole/Second Moment

Octopole/Third Moment

Hexadecapole/Fourth Moment

Multipole Moments

Conversion Factor
1.602189 x 10-19 C (Coulomb)
4.803242 x 10-10 esu
8.478418 x 10-30 m C
2.541765 x 10-18 cm esu
2.541765 D (Debyes)
4.486584 x 10-40 m2 C
1.345044 x 10-26 cm2 esu
1.345044 B (Buckinghams)
2.374197 x 10-50 m3 C
7.117664 x 10-35 cm3 esu
7.117664 U (unnamed)
1.256371 x 10-60 m4 C
3.766505 x 10-43 cm4 esu
3.766505 N (noname)

Table 3-3: Calculated Multipole Moments of H20

This work

Moments Spec. 7 Theor. 7 Rein exp

one-c two-c one-c two-c et al.

/p (D) 2.306 2.235 2.239 1.903 2.52a 1.85

Qxx -0.23 -0.01 -0.12 0.03 -0.32b -0.13c 0.4d
Qyy 1.35 1.57 1.24 1.43 1.00 2.63 1.6
Qzz -1.12 -1.56 -1.12 -1.46 -0.77 -2.50 -2.0

Qxxx 5.07 8.72 4.92 7.93 9.2e
Oxyy -8.46 -15.8 -7.95 -13.7 -16.8
2xzz 3.34 7.05 3.02 5.78 7.6

rxxxx -4.70 -3.42 -4.35 -8.93
Fxxyy 5.68 4.23 5.28 11.04
Fxxzz -0.98 -0.81 -0.93 -2.11
ryyyy -1.98 -1.64 -1.84 -4.14
ryyzz -3.70 -2.59 -3.44 -6.91
Pzzzz 4.68 3.41 4.37 9.02

Table 3-3 Continued:

a. Reference [79]
b. Reference [32]
c. From Magnetic Susceptibility, see reference [81]
d. Reference [82]
e. Reference [33]

Table 3-4: Calculated Multipole Moments of Benzene

Moments This Work Rabinowitz Exp.d
et al.c E

Zerner Basisa Rein Basisb

QXX 1.419 4.132 2.86 4.34
QYY 1.419 4.132 2.86 4.34
QZZ -2.838 -8.264 -8.69 -8.69

rXXXX 0.99 2.10
rXXYY 0.33 0.70
rXXZZ -1.31 -2.80
FYYYY 0.99 2.10
FYYZZ -1.31 -2.80
FZZZZ 2.63 5.60

a. Reference
b. Reference
c. Reference
d. Reference


Table 3-5: Calculated Multipole Moments of Hexafluorobenzene

Moments This Work Laidigc Expd

Zerner Basisa ERBb

QXX -8.766 -4.653 -4.718 -4.77
QYY -8.766 -4.653 -4.718 -4.77
QZZ 17.532 9.306 9.434 9.54

rXXXX -5.56 -1.98
FXXYY -1.85 -0.66
rXXZZ -7.41 -2.64
rYYYY -5.56 -1.98
rYYZZ 7.41 2.64
FZZZZ -14.82 -5.27

a. Reference [78]
b. Expanded Rein
c. Reference [84]
d. Reference [85]

Basis. See Reference [79] and text

Table 3-6: Calculated Multipole Moments for Pyrazine in Water Solution

Moments This Work Hush aet al.

Qxx 11.6 10.2

Qyy 1.14 2.2

Qzz -12.8 -12.4

rXXXX 1.53

FXXYY 2.94

rXXZZ 4.47

TYYYY 1.04

FYYZZ 1.90

rZZZZ -6.37


a. Reference [80]

(A) Water

(B) Benzene








(D) Pyrazine

Figure 3: Input Molecular Orientations for Multipole Moment Calculations.







Figure 4: Dimer Structures of Hexafluorobenzene-Benzene and Benzene-Benzene Based on
Multipole Moment Interpretations


As mentioned in Chapter 1, we are now going to examine the solvent effect on

molecular electronic structure. As we know, most chemistry reactions take place in

solution, and the solvation phenomenon is very challenging. A successful theoretical

method should be able to calculate two basic experimental measurements: the solvation

energy and the electronic spectroscopic shifts in solutions. We tackle the problem

by approximating the solvent media as a dielectric continuum surrounding a quantum

mechanical solute molecule, and by examining the electrostatic interactions between the

molecular multiple moments and the solvent dielectric continuum.

The Reaction Field

We consider a solute molecule in a spherical cavity of radius a inside a solvent media

with a uniform dielectric constant E, as shown in Figure 5. The electrostatic potential of

the solute molecule is expanded in multiple moments and allowed to interact with the

dielectric continuum. For now, let us examine the potential at any point P(r) outside the

spheric cavity caused by a unit charge s distance away from the center of the sphere 0.


The potential can then be expressed as
1 1
tout(r) = = 1
re Vs2 2srcosO + r2

1 s s2) 2
=-1 2T- +
r r r2

where we have

the CosO.

1 2) 1 3 4F 82)2
= [1 + -2rS-+ + 2 'T -2)T+r- + (4.1)

1 s 3T2 -_ 1 S2 57 3T s 3
= -[1+-+ + -..]
r r r \r+ 2 r
1 0 00 (A pr

n=O n=
used binomial expansion and Legendre polynomials, Pn(r), and r is

At a point inside the sphere equation (4.1) may not converge because s is not

guaranteed to be larger than r. However, the potential must satisfy the Laplace equation,

which, if the cavity is overall neutral, can be written as

V2in = 0


its general solution being


in = E rn + ) P(T)

We then require that the potential on the sphere be single valued,

in(r) = ou.t(r)jr=a


and due to linear independence of Legendre polynomials, this gives, for a particular term

n = 1 in equations (4.1) and (4.3),

a1+ = Bal + at+ (4.5)

and we further require continuity on the sphere,

in(r) = ,out(r)\r=a (4.6)

which gives, again for a particular term n = 1,

-e(l + 1)Ala-(+2) = lBal-1 (1 + 1)Dra-(l+2) (4.7)

Equations (4.5) and (4.7) are standard boundary conditions for cavity models, and they


(l+ 1)(1- E) 1
B + (l + 1) a2+1 D = fz(e)DI (4.8)

With BI determined by equation (4.8) we can write the potential inside the cavity as

bn+(r) = -+ f()rn DnPn(-) (4.9)

where the second term is caused by the solvent dielectric interaction with the solute

molecular multiple, and it is called the reaction field,

R = fn(E)DnraP(T) (4.10)

the factor fi(e) is called the reaction field factor, first derived by Kirkwood [86, 87].

The Interaction Energy in a Reaction Field

The interaction energy between the solute molecule with a quantum mechanically

described charge distribution p(r) and the dielectrically averaged solvent media is

Eint= p(r)Rd3 r

= fID fp(r)rl'P(r)d r (4.11)

1=0 m=-1
where Mtm's are multiple moments of the solute molecule in a general form. In many

applications the definition of the first moment differs from that of the second to gain

considerable computational advantages, a point that we shall address more later.

The first term of equation (4.11) is the well-known Born term

E = (1 )()q2 (4.12)

with q being the net charge of the solute molecule. The second term, first examined by

Onsager [88], is known as the Onsager term,

(2) 1 2(E 1) 1 9( (4.13)
Eint = 2(2E+1) a3 2 413)

with p being the dipole moment of the solute molecule. The higher terms are all derived

in this study, and they are

(3) 3(E -1) ( 2 =- 2 (4.14)
int 2(3E + 2) Qa=(


(4) 4 (E-1) (1\ ) 2 = -k(E)Q2 (4.15)

(5) 5 (E- ) (1 )2 -t(e)F2 (4.16)
t 2 (5 + 4) 9 0 -

In the development of this work conceptual difficulty lies in the uncertainty on how

to include higher terms into equation (4.11). Experience tells us these moments do not

need to be in any particular definition as long as they are consistent throughout the same

formulation, and we choose the form of equation (3.20) of Chapter 3 for two major

reasons. The first reason is that the multiple moments as defined by Buckingham are

already calculated and stored in our multiple moments analysis as described in Chapter

3. The second reason is that the coordinates appearing in equation (3.20) are primitive

multiple moment operators, with their integrals in STO basis evaluated as described

in Chapter 2 and stored as 'primitives'. It should also be added that equation (3.20)

has aesthetic appeal with the numerical coefficients equal to the number of the possible

permutations of the index, a feature we first mentioned in Chapter 3.

The Self-Consistent Reaction Model

When the solute molecule, the quantum mechanical motif that will be treated below,

is placed in the reaction field, the energy of the system becomes

Eu = Eo 1 M M
=Eo- ft(OJMrm))(OJMrJO)

where 1V> is the state of the molecule, and Eo the vacuum energy of the molecule. We

form the functional

L = E mfi (Mrhmb)(IMrk)
1,m (4.18)

W((1) 1)
where W is the Lagrangian multiplier. Variation of equation (4.18) proceeds as

6L = 6Eo W ( ) 1 6 f(IMrIm )(0 MIm)

= (6VIHo 4) Z fz(6 1Mm|0)(0 Mm1) W(6, I) + c.c. (4.19)


which further lead to the Schr6dinger equation of the molecular state

(Ho Ef (IMzm?)Mi)m I4) = WI\) (4.20)

It is clear now W is the quantum mechanical energy of the solute molecule, given by

W = (lHo |MI|Mmi
= Eo- E f I(Mrm |))(OIMrlm )

In a similar fashion as we did in Chapter 3, specifically, the equations (3.26) through

(3.30), we approximate I4> by a Slater determinant composed of molecular orbitals

{1i>} in the form of linear expansion of atomic basis functions. This procedure leads

to a Fock-like equation

(f ,ft(M |)Mm) 1i) = Ei\oi) (4.22)

where Ei is the solute molecular orbital energy, and fo is the vacuum Fock operator defined

in equation (3.29) of Chapter 3. In the usually way, we can define the new Fock operator

f = fo E f(|IM |k)Mm (4.23)
to give

f biM) = Ei) (4.24)

Equation (4.24) is solved iteratively because the construction of the new Fock operator

requires knowledge of the molecular orbitals both to form fo and to calculate ( IMjMI).

Upon an initial guess of the molecular orbitals, the iteration proceeds cycle after cycle

until a self-consistency in achieved. This treatment of the molecule in solution is therefore

known as the Self-Consistent Reaction Field (SCRF) model [40]. We have found later

that faster convergence results when the moments are updated each cycle, as opposed to

first converge the SCF and then the SCF and SCRF. Similarly we have found there is

no advantage to perform the SCRF at the dipole level first, and gradually the quadrupole

and higher moments after the initial calculation convergence.

Following the general approach of SCRF outlined above, there have been numerous

theories developed recently. One of them, called A theory", views that the W from


equation (4.21) is the energy of the quantum motif, i.e. that of vacuum solute molecule

and its interaction with the solvent. To counter the neglect of the solvent in the variation

process, the system energy of equation (4.17) should add to W an extra term from the

E, = W+ Ec,c

1 (4.25)
=W+M ( M }|)(|Mrlm)
which yields again the energy of equation (4.17) and the Ec,c herein is called the "solvent

cost", reflecting the fact that the solvent has lost energy to dissolve the solute. The second

theory, called "B theory", argues that the total energy of equation (4.17) is should be given

by the energy functional [89], and therefore the Fock operator equation (4.23) becomes,

f = fo Ef ( Mr ) M (4.26)

In B theory the solvent cost is included in equation (4.26) in the SCF energy variations.

The Configuration Interaction Calculation for Electronic Spectra

Solvents affect molecular electronic structure, causing sometimes sizable energy

shifts. The excited states of molecules are often calculated using the configuration

interaction (CI) method. Conceptually, CI is a straightforward method of linear variation

over the determinants of many-electron wave functions. Again we desire to solve the

Schr6dinger equation by assuming an approximate solution. In the CI approach the

approximate wave function I') is expanded as a linear combination of Slater determinants

I = C, 0 (4.27)

Where the Cs's, are the parameters which are varied to make E[T] stationary. N is the

number of expansion functions. We examine,

E[T] = (T 1Hj\I) / (F I') (4.28)

in which H in equation (4.28) above is the electronic Hamiltonian. This leads to the

general matrix eigenvalue equation

HC = SCe (4.29)


Hst = (OjHIlt) (4.30)


Sst = ('s|It) (4.31)

The N linearly independent eigenvectors c, the columns of C, can be chosen to be

orthonormal with respect to S,

c*S cq -= cStctq = 6pq (4.32)

leading to the simple eigenvalue equation

HcP = cE, (4.33)

We now need to find the basis over which CI is to be performed. A very natural

and intuitive choice of the determinantal wave functions {|J)} is constructed in the

following fashion,

occ vir occ vir occ vir
I)=Co10)+EEqI) + EECi + E E -bijklvi)jk +***
i a i (4.34)

where 10o), in the first term of equation (4.34), is the Hartree-Fock reference function

obtained in the SCF procedure previously described. The set of functions { ij ..')

are spin-adapted configurations generated by removing electrons from the occupied spin

molecular orbital {|ij), Ij|),...} and placing them in the virtual, or unoccupied, spin

orbital { 1a), |b), ...}. These are known collectively as the electronic excitations from

the molecular ground state described by the Hartree-Fock theory. If the summation in

equation (4.34) is carried to its limit, the calculation so performed is said to be a full

CI calculation. Computational resources usually limit the summation of equation (4.34)

to low "order", resulting in the CI Singles (CIS) procedure, CI Singles and Doubles

(CISD), and so on. Even with CI truncated in this fashion it is not unusual to reduce

it even further by limiting the size of the active space. That is, to limit the number of

active orbitals from which electrons are removed or into which they are placed, usually

centered around the highest energy occupied molecular orbital (HOMO) and the lowest

energy unoccupied molecular orbital (LUMO), or the "HOMO-LUMO gap". For the

work in this thesis the CI expansion, equation (4.34), is truncated at the second term,

and thus belongs to CIS model. The popular ZINDO code has been parametrized for

spectroscopy at this level of theory [90].

In these calculations the CI eigenvalues of equation (4.33)

El < E2 < E3 ..... < EN (4.35)

can be shown to be upper bounds to the corresponding true eigenvalues of the CI

Hamiltonian, and they approximate the energy of excited molecular states. Furthermore,

if additional configuration functions are added to the expansion, equation (4.27), each

eigenvalue EN+1 of the (N + 1) term expansion satisfies the bracketing theorem

E,, < EN+1) < E(N) (4.36)

and as {f N)}EN approaches completeness the {E"} must approach the exact eigen-

values of the CI Hamiltonian from above. It should be noted, however, that there is no

variational theorem for the transition energies between states.

We consider I0o), and one singly excited configuration function, 109), or

ICIS) = (C 0100) + CWa |i ))


and proceed to evaluate ('Pc1sIH cjs) for its lowest eigenvalue and associated eigen-

vector we obtain the matrix eigenvalue problem

((liH ) (lHlf) Co E, (C (4.38)
(V{W|H o) (OW\H| )c O \CC)

It is obvious that any mixing of the diagonal elements is through the off-diagonal element

()Of|,1j) = (ilha) + 1 [(irlar) (irl a)] (4.39)

or its adjoint. The ia-th element of the Fock operator, again in the basis of spin molecular

orbitals, is given by

(qif1j0a) = (ill>a) + Z[(irlar) (irjra>] (4.40)

Notice that the right-hand-side of both equation (4.39) and equation (4.40) are identical

and thus

(o01H\ 7) = (<|lfl|a) (4.41)

Now by definition solving the Hartree-Fock eigenvalue problem requires the off-diagonal

elements of the Fock matrix to vanish, or upon solution

(ill Wa) = 0 (4.42)

and then so must equation (4.39). The above exercise is an example of Brillouin's

Theorem [91]. Stated in words as, singly excited Slater determinates i\k) will not

directly interact with a closed-shell Hartree-Fock reference determinant I|'o) through

the electronic Hamiltonian, or

(b01HI07) = 0 (4.43)

{lf})} leaves the closed-shell Hartree-Fock reference unchanged. This is a desirable

feature for the INDO technique that has already included the experimental data in its

parametrization. This approach to the CI method has proven quite successful for the

calculation of the low energy ultraviolet, and visible, absorption spectra of molecular


Up to this stage the discussion of CI has been completely general in the sense that

it has not dealt with the specific forms of the Hamiltonian. When a molecule is placed

in the reaction field, the CI Hamiltonian matrix elements are

HIj = (IajHo h (0oMIM o)}MrI )
= (0jHo Ij) Efi(NoIMrtm o)(IM IM j)

where HO is the vacuum Hamiltonian, and in the formalism we follow the A theory

without losing generality. The rest of the CI implementation follows the discussion in

a straightforward fashion. There are, however, two important modifications to the CI

energy after the diagonalization.

The first deals with the absorption process in a dielectric media. According to

Classius-Mosotti, the bulk solvent dielectric constant c has the following relation

E-1 D' -1 n2 -1
I = +1 (4.45)
(l+1)e+l (l+ 1)D' +1 + (l+ 1)n2 +l

with D' being the contribution due to the slow motion of solvent reorientation, and n, the

solvent index of reflection, reflecting the electron polarization. In this spirit, we factor

equation (4.10) into two parts

f,(e) = ft (D') + fi(n2) (4.46)

We further assume that in an absorption process, generally of the order of 10-18 seconds,

there is no motion of the solvent molecules, but that the electronic polarization is

instantaneous. The energy transition is then modified and it becomes

AEI = fi(n)((?oiMtm'Io)('IlM I[i) (I(IMMi}\)2) (4.47)

where the first term removes the incorrect term included in using the SCRF orbitals

in forming the CI matrix, and the second adds the response of the electron polariza-

tion -fz(n)(li\Mimp\) to the excited state multipoles (OIIM Imi). Note that this

term does not allow for complete relaxation which would otherwise be described by

- fl(-)(I IM IM).

The second modification has to do with the solvent cost. Similar to equation (4.25),

we must add to the quantum mechanical result an extra term, that is,

WIA = ( H|I ) + ft(e) Mim o)(M I) (4.48)

to count for the solvent cost in the excited state.

With the above discussion, we are able to calculate the spectroscopic absorption

1o0) --+ ),

WIA WOA = [(',IHIi) (oIH[4o)]+

f (e) (<|01MMlo)]+
2 E fo] ( 1M 00 (,1(4.49)

1 fz(n)(VI(Mrj'iM[ l,m
where the first term is a direct result of the CI, the second term is due to the solvent

cost, and the third term is due to instantaneous electronic relaxation. And similarly,

discounting the solvent cost term, we can write for B theory

WfB Wo = [(H } (VoIHI o)]+

E 2 1 0 MI ) IM (4.50)
A second way to view the absorption process is to assume that the ground state and

excited states are intimately coupled to form a mean field, and the states feel the molecular

charge distribution through the instantaneous electron polarization of the solvent. In such

a case, the energy of a state \0k) is shifted by

f (n2) [I(oIM mlo) 2-
(OkkMrIOCk)((iolM71io) + (ziIMrlIi))}/2
and this leads to "Al theory",
W W Wo 1 = [ I|H ) ( oH 1o)]+

I f (e) f()o Mm 1 O)[(I MrI z) (4o MrmI'o)]- (4.52)

1 fl (n2)[(i oMrIoM (VIlMr bi)]2
and to B theory
W 1 W1 o= [(MilHIV) (0o8IHIo)]-

1 (4.53)
S:fi (n2)[(oiMrmo) (I|Mm )]2
Recently a new theory, called "C theory", is developed and tested in this study. This

theory, based on A theory, calculate the energy for the final state ]|J) as

Ef = Wf+ f(n2)(i IlMrIlf)[('olMImlo) (of jMjjIm )]
+ ft (n2)( IMm If)2/2 + ft (D')()oIMrm o)2/2
where the difference from A theory is that the final state assumes the ground state

geometrical contribution in the solvent cost energy. Its absorption is expressed as

Ef EgC = Wf Wo f(n 2)(( IMrIml) (_ oIjMYmWo))2 (4.55)

and in the mean field we have for "Cl theory",

E1 Ego = Wf Wo + f(n2((f) (oMI o))
1,m (4.56)

(3(oIMm 10o) (

In this section results of the performance of high moments in the SCRF-CI model

are reported. Our Hartree-Fock reference function is obtained by using the semiempirical

INDO/S technique. In the SCRF model, the convenient reference for multiple expansion

is at the center of mass (as discussed in the last chapter). The cavity redias is obtained

from the mass density, for which we have molecular cavity volume

4 3 M
V = -Tra =
3 dNA

a = 0.7346(M/d)I A

where M is the molecular weight in amu, d is the mass density of the solute (g/cm3)

and a, of course, the cavity radius. Unless specified, we assume the reference standard

in our calculations below.


Although H20 does not have interesting UV-Vis spectroscopic feature, it is a favorite

molecule to test the convergence of the SCRF model. Our result for H20 solvent cost

energy in H20 solution is shown in Table 4-7. The convergence this work achieved is

satisfactory from dipole to hexadecapole expansion, and it compares well with the results

from Multiconfigurational SCRF calculation by Mikkelsen, Argen, Jensen and Helgaker


We do not hesitate to point out here that there is clearly an error in reference [15],

where Table II and Table IV do not agree. And Judging from our values, their dipole

level calculation may be reported incorrectly.


The n 7r* UV-Vis spectroscopic solvent shifts of acetone is calculated with all

theories A, Al, B, Bl, C, Cl, at different cavity radii. Our starting point for the cavity

radii is the mass density radius, 2.84 A. The calculation is done first at the dipole level,

and then include higher moments up to hexdecapole. The results are listed in Tables 4-8

to 4-13 and the cavity radius dependence of calculated solvent shifts with Cl theory is

shown in Figure 6.

It is a general feature that solvent shift decreases as the cavity radius increases, as

shown in Figure 6 for the case of Cl theory. At the radius region examined, all previous

theories A, Al, B, and Bl underestimate the solvatochromore shift significantly for

acetone. Only theories C and Cl give results that compare well with experiment. Within

the same theory, dipole moment expansion usually underestimates the experimental shift,

and the inclusion of higher moments brings additional blue shift to yield better agreement

with experiment.


Solvent shifts of pyrazine are calculated at quadrupole and hexadecapole levels of

multiple expansion with all theories at cavity radius 3.61 A, a radius that is found in

our experience appropriate for diazene series of molecules. The results are reported in

Table 4-14.

Pyrazine is an interesting molecule, because it has not been able to be studied by

the SCRF methods that include only dipole moment expansions. Our conclusions made

earlier regarding the performance of multiple moments with difference theories for

acetone molecule is once again confirmed here by the results in Table 4-14.




Figure 5: Spherical Cavity Model

- -

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd