UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
MOLECULAR MULTIPLE MOMENTS AND THEIR PERFORMANCE IN THE SELFCONSISTENT REACTION FIELD MODEL By XUEHE ZHENG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1995 To the Clbhs ACKNOWLEDGMENTS 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 indepth 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! iii I certainly thank Dr. Marshall G. Cory (QTP, FL) for his help in constructing this thesis. 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. TABLE OF CONTENTS ACKNOWLEDGMENTS ...... ... ... ....... ............ iii LIST OF TABLES ......... ... .................... vii LIST OF FIGURES ....... .. ............. ........... viii ABSTRACT ........... ... ............ .... .......... ix CHAPTERS 1 INTRODUCTION .................................. 1 2 THE EVALUATION OF ELECTRIC MULTIPLE MOMENT INTEGRALS OVER SLATERTYPE ORBITALS ........................ 8 Introduction ........... ... ........ .. .. ........ 8 Preliminary Considerations.............................. 9 M ethod .. .. . .. .. .. 14 Discussion and Conclusions .......................... 33 3 THE CALCULATION OF MOLECULAR ELECTRIC MULTIPLE M OM ENTS ..................................... 39 Definitions of Electric Multipole Moments ................ 39 The Theory and Calculation of Molecular Orbitals ............. 47 Some Calculated Molecular Multipole Moments ............. 57 4 A SELFCONSISTENT REACTION FIELD MODEL INCLUDING HIGH MOLECULAR ELECTRIC MULTIPLE MOMENTS APPLIED TO THE CALCULATION OF MOLECULAR ELECTRONIC SPECTRA ...... 71 The Reaction Field ............................. 71 The Interaction Energy in a Reaction Field ................ 74 The SelfConsistent Reaction Model .................... 75 The Configuration Interaction Calculation for Electronic Spectra ..... 78 Results . . . . 87 5 CONCLUSIONS ................................. 100 REFERENCES ...................................... 105 BIOGRAPHICAL SKETCH .............................. 111 Table 21: Table 32: Table 33: Table 34: Table 35: Table 36: Table 47: Table 48: Table 49: Table 410: Table 411: Table 412: LIST OF TABLES 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 (cm1): Theory A.. Calculated acetone n ir* solvent shifts (cm1): Theory Al. Calculated acetone n 7r* solvent shifts (cm1): Theory B.. Calculated acetone n ir* solvent shifts (cm1): Theory Bl. Calculated acetone n 7r* solvent shifts (cm1): Theory C.. Table 413: Calculated acetone n 7r* solvent shifts (cm1): Theory Cl. . 38 . 63 . 64 . .. 66 . .. 67 . .. 68 . .. 92 . .. 93 . 94 . .. 95 . .. 96 . .. 97 . .. 98 Table 414: Pyrazine n7r* Solvent Shift (cm1): From Acetonitrile to Diethyl Ether. 99 LIST OF FIGURES 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 HexafluorobenzeneBenzene and BenzeneBenzene 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 MOLECULAR MULTIPLE MOMENTS AND THEIR PERFORMANCE IN THE SELFCONSISTENT REACTION FIELD MODEL By XUEHE ZHENG 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 HartreeFock level of theory. A general approach is developed for multiple moment integrals arising in computations using Slatertype 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 SelfConsistent Reaction Field (SCRF) model through the first order multiple expansion of electrostatic potentials between solute molecules and solvent dielectric media. Computations are implemented ix 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. CHAPTER 1 INTRODUCTION 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 [24], molecular optical and magnetic properties [58], intermolecular force theory in terms of the London series [4, 7, 913], and the quantum mechanical solvent effect in the formalism of the self consistent reaction field theory [1416]. 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 1 2 are the highestorder multiple moments actually measured using a method that was established thirtyfive 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 highorder 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 oneelectron 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 atomiclike orbitals. There are two kinds of prevailing atomic orbitals in the current 3 stage of quantum mechanical calculations: Slatertype orbital (STO) with form ear and Gaussiantype 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 multipolemoment 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 4 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 5 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 longrange 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 solutesolvent 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 (QMMM) [38, 39], and the selfconsistent reaction field (SCRF) method [4047]. 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 solutesolvent system is accounted for by 6 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 7 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 UVVis spectral shift in solution as calculated by our method. And, finally, Chapter 5 gives the conclusion and future outlook of this work. CHAPTER 2 THE EVALUATION OF ELECTRIC MULTIPLE MOMENT INTEGRALS OVER SLATERTYPE ORBITALS Introduction 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) 2(1+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 oneelectron, and at most, threecenter 8 9 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 multipolemoment 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)! rnler7(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)! rnlle(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)], d2y2 = [(2, 0,0) (0,2,0)], d = 4(1,1, 0), dxa = T (1,0,1), dyz = (0, 1,1). f harmonics: 11 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(2y2) = [(2, 0,1) (0,2,1)], fxyz = (1, 1, 1), f133xy3 = 2[(3,0,0) 3(1,2,0)], f3x2yY3 = [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(SWy2) = [(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(X33y3) = 70[(3,0,1) 3(1,2,1)], gz(32yy,3) = j [3(2, 1, 1) (0,3, 1)], 2y2 = [(4, 0, 0) 6(2, 2,0) + (0,4,0)], zy(X2y2)= i/5[(3, 1,0) (1, 3,0)]. 12 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. Lowercase 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 uppercase 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 13 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) R R 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 R' (7) sine = X/Ax'2 + Ay'2 and the rotation matrix R(0,0) can be obtained in two steps: 14 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 / Method 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) m=O 15 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 (2.17) (i,j,k) SCr,s,t (bax, ay,baz) xyyz (r,s,t)=(0,0,0) The translations in terms of equation (2.17) enable a multiple moment integral, apart from a normalization factor, to assume the following form, nala1 ji J yk j .ra. I m n nblb1 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 (rnaI.1 o (ar Ir. tlb1 .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 17 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 OOhata 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 lefthand quantities in a formula are centered on A, and those on the righthand 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) 4w 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)+ (2.23) 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.24) (2,0,0) (0,2,0)] dy = [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.30) (2,1,0) (0,3,0) + 4(0, 1,2)] (2.25) (2.26) (2.27) (2.28) (2.29) z(x22) = 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_3y= 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 f3x2yy3 = 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)] (2.31) (2.32) (2.33) (2.34) 3 4 22 4 22 4 22 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)} (2.35) 3 10 3  2 3 2 2 xz3 = { (4baxbaz 3babaz 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)} (2.36) 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)} (2.37) gz2(x2y2)  V + 6a 6,a) (0, 0, 0) + 4 (3baz a)(1, 0,0) ( 372r2\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)} (2.38) 3/5 2 3 b3 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) +6(1,1,2)} (2.39) 3 70 3 2 2 2 gz(x33xy2) = {(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)} (2.40) gz(3x2yy3) = 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)] + 6ax(1, 1,1)  (0,3,1) + 3(2, 1, 1)} (2.41) 3 35 4 22 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)} (2.42) 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)} (2.43) 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) 2 (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) = bobo (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) + bo(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) + bo(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) 2 +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) 2 +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) 2 + 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) 22 2 2 (2,0,2) = boboz(0, 0,0) + 2bobo(1, 0,0) + 2bo2boz(0, 0,1) + 4boFbo,(1, 0,1) + boz(2, 0,0) + bo(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) 2 + 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) 27 (0, 4, 0) = o,(0, 0, 0) + 4bo(0, 1, 0) + 6bo(0, 2, 0) + 4boy (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, '22 (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 fashion, alaVa aaeaa XI oyo 1oh I b"ybfzbge I (MI,m n) = (r.yt zXIr"blb 04ybz e4rb (rnatl1 ga(baxz,ay, baz)x '"tY'" t'" e"r. = ~g (bax, ba y, baz) g (box, boy, boz) nalaI ( 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 Onecenter Integrals After the translations described by equation (2.44), all onecenter 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) 2r J Sini CosJdo 0 and we proceed to evaluate the nonvanishing 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 cosm2 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 m3 / cosm4sinn"do Ir m+n m+n2 0 0 r,27 m1 m3 m5 f  + 2  4 cost6 sinn do = * m+n m+n2 m+n4 J 0 r,27r (m 1)!!n!! !f s.n d (m + n)!! J ne 0 (2.47) 29 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 n2J 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)!! 0 27r c in (m 1)!!(n 1)!! ]cosmosinnd = n (27) = F(m, n) (2.50) J (m + n)!! 0 We then examine the sine integrals with odd n in equation (2.47), i0 2 2) (n I )2(nn)!1 = nnl___nl__ sin d = cos = (2.51) n!', n! 0 2r f sinnodb = 0 (2.52) 0 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) 0 All the angular integrals necessary for Ii are given by equations (2.51) and (2.53). The radial part of Ii is simply rnecr = (2.54) 0 Twocenter Integrals In an ellipsoidal coordinate system we use equation ((6)) noted previously to express a twocenter integral that evolves from the translation of equation (2.44), 12 = (rnla1le(ar alrl1 xiy'zke brb 00 1 27r = d f do ()(R)n+nb1a1b+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.+nbla.l +i+j+k+l ( + )a+1( )nb+l 1 1 [(q2 1)(1 2]+) )keP F(i,j) where p = y(Ca + (b) (2.56) 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)'[(21)(1 q2)(1 r)6ep'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 31 proposed by Stevens [59, 60] based on the binomial theorem, Za,,,6(f,T)= ) i,j,k,l,m (2.58) Aa+3+27ij2k+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 twocenter integral of equation (2.55) is evaluated as the product of Z and F functions, 12 = (rnale ara Xr'' lxy ztebrb) (2.59) = Zn.1,nblb,(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,s,t)(u,v,w)(p,q) *. r+u+P s+v t+w+q where t (r11, r12, r13) = r r8l2r 3 4,,,v'w( r21 72, r23) = !W! 21r'r23 (2.62) k! ck, q(r31,r33)= r'3 Considering I2 of equation (2.55), we now have RI2 = R(r'oe ra Tn ixZ Yk erbCb bb zb b eb) = (rn^ecrajrn(Ray bz )er^" (.3 (r,s,t)(u,v,w)(p,q) 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 nonzero, integral evaluation proceeds as outlined two subsections ago. 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 formulas. (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 (2.63). 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 Slatertype 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 34 studies on solvation using the selfconsistent reaction field model, and for examining the longrange 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)])ep( or the C function [27], (1 ,a+/3+y+6+2e+1 C"^(P, ) = (bR L6', (p,T) (2.65) 35 Both L and C functions, if evaluated in the fashion of equation (2.58) require a sixfold sum, while the evaluation of Z functions needs a fivefold 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 nonvanishing quadrupole integral between d orbitals contains on average 18 different terms of 64 possible, but only 5 or 6 are nonzero. In Table 21 we list the number of Cartesians for a given order, and the number of Cartesians that have nonvanishing 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%. Zb e Yb Theta Ya Yb Za" Zb" Xa" n" Ya" Yb" Za' Zb' Xa' Xa' Figure 1: Rotation From Ellipsoidal Coordinate to Molecular Frame Za Ya X'i X' A=B Symmetry zero Value Figure 2: Flow chart of the HIMO computational program Table 21: 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 CHAPTER 3 THE CALCULATION OF MOLECULAR ELECTRIC MULTIPLE MOMENTS 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) l,m 40 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 ith 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) 41 The above type of definition can be found in Jackson [63]. A much simpler definition is given by Bottcher [64], as Dipole: M x p (x)d3x' (3.9) Quadrupole: Q = xx p x dx (3.10) Octopole: 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 42 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+ 6XY 6XZ 6YZ 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+ 6(XYZ)xyz]+ 43 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, Monopole: q = p(r)d3r (3.15) Dipole: Ax = Jf xp(r)d3r 44 Py = f yp(r)d3r z = zP(r)d3r (3.16) Quadrupole: Q = f f (22 y Z2)p()d3r Qyy = f (2y22 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) Octopole: zxxx = f (2x 33xy2 3 z2)d3r yyy = f (2y33x2y 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 45 yyz i f (4y2z x2z z3)p(r)d3r =yzz = (4yz X2Y y3)p(r)d3r (3.18) Hexadecapole: 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 (9y22z 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) 46 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) + R5 R[XXXf2xx + YYYlyyy + ZZZQ.z + 3(XXYQzy+ XXZQXz, + XYYQxyy + XZZQzz, + YYZQ ,vz + YZZQwY)+ 6XYZQYz]+ 1[XXXxr"XZ + YYYYFyyyy + ZZZZPzzzz + 4(XXXYFTxy+ XXXZFTx, + XYYYr,Yy +XZZZXzzz +YYYZVyyyz+ ZZPYzzz) + 6(XXYYFxxyy + XXZZTxxzz + YYZZvyyzz)+ 12(XXYZ Py, + XYYZTFyy, + XYZZxy^zz)+ (3.20) 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 47 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 oneelectron 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 manyelectron 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. 48 1. The HartreeFock Theory The electronic Hamiltonian of a molecule under BornOppenheimer approximation can be written in atomic units (F = m = e = 1) as nV2 1n n M Z (3.22) 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 manyelectron 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 HartreeFock 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 closedshell 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) i,j=l n n = h + 2 E (JI Kii) i=1 i,j=l where ,i is a molecular orbital, hii is oneelectron matrix element and Ji and Kij are twoelectron 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) i,j where ij is a Lagrange multiplier, 6L = 6{E E Kji[(ilj) 6kil} = 0 (3.27) ij to obtain the least E = (?^I//)/( ') > Eexact this leads to n (6i(1){hii(1)+ 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 50 Since the variation is arbitrary, we must have n n {hi (1) + ij ki]} i4(1)) = EEji4j(1)) j j (3.29) n 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 HarteeFock equation, from which molecular orbitals of the given state [I) are obtained. 2. The LCAOSCF Equations For the HartreeFock 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 N 1i) = E ciilX/) (3.31) uL=l 51 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 HartreeFock equation (3.30) in LCAO expansion takes the form N N / 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 N N 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' = X1C, C = XC' (3.36) % P then Equation (3.34) can be written as FXC' = SXC's (3.37) and 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 HartreeFock" (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 SelfConsistent 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 53 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 selfconsistency is achieved. 3. The INDO Technique The Fock matrix elements of equation (3.39) for a restricted HartreeFock formalism can be expressed in the atomic basis as follows, N/2 fo = hJIV + ZECAaC,*[2(po7vA) (yoIAv)] a aA N/2 = 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 N/2 P = 2 E CAaCa (3.41) a and h,, is the oneelectron Hamiltonian matrix element, and Gp,v is the twoelectron 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. 54 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 twoelectron integrals, in terms of the number of such integrals and the physics they contribute to the theory, is pronounced. In particular, all three and fourcenter 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 onecenter twoelectron integrals and NDDO methods also include 55 all twocenter 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 abinitio HartreeFock 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 LCAOMO 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: FIAIA = UIAIA + E [PCC 7AC VAC] C#A (3.45) +" Pt'AAA [(AAAAIIAO`A) I(IAAAAI0'AILA)] AA UA 2. Case I 0 v; A = B: (3.46) FIPAA P= Z P A [A A A IVAO'A) A (AAAIAA)] AA 'A 3. Case pI $ v; A # B: F AVB AVB 2 PAVB TAB (3.47) 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 onecenter core integral UpA A of equation (3.45), U ==AJ A A 1 + VAIA}) (3.48) the onecenter electronnuclear attraction, VA = ZA TiA the first order atomic density Pcc is defined as, Pcc = E Ac Ac Ac the twocenter twoelectron 7yAc is defined as, (3.49) (3.50) YAC = (IAAc A Ac) (3.51) 57 and VAC the twocenter electronnuclear 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 oneelectron property, a general multiple moment M can be calculated from 58 quantum mechanics as the expectation value for the state I'b), M = (I\MIp) = (01(1)02(2) .. 2N(2N)IM11(1)02(2) 02N(2N)) (3.54) N 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 ith particle having mass mi and charge qj located along the zaxis 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) i 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 59 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 32 for all the units involved in the following report. H20 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 33. The calculation uses experimental geometry (ROH = 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 twoelectron integral INDO parametrization of 7, one which is fitted to produce good UVVis 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 twocenter terms of the moment integrals. From this study we favor the inclusion of the twocenter and the spectroscopic 7. Our results compare well with those reported in the existing 60 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. Benzene 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 HF631G** optimized geometry with GAMESS,(RCC=l.386 A, RCH=l.076 A), the orientation of the input geometry is shown in Figure 3 (B). Benzene multiple moments are calculated and tabulated in Table 34. 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 34). The most noticeable chemistry that can be explained from these results is the herringshape 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 61 exist as a perturbation to tilt the perpendicular dimer, giving the observed structure, as shown in Figure 4 (A). Hexafluorobenzene Calculated multiple moments of hexafluorobenzene, with its geometry optimized using the same procedure as that for benzene (RCc=1.379 A, RcF=1.314 A) are listed in Table 35. 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 35). With the same argument for benzene dimers, we would predict the parallel docking structure for a HexafluorobenzeneBenzene dimer ( see Figure 4 (B)), as is observed in experiment. Pyrazine Table 36 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 Selfconsistent Reaction Field Model to be described in the next Chapter. The calculation uses the Zerner basis, default in ZINDO program. 62 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 32: The Conversion Factors for 1 a.u. Property Electric Charge (Zeroth Moment) Dipole/First Moment Quadrupole/Second Moment Octopole/Third Moment Hexadecapole/Fourth Moment 63 Multipole Moments Conversion Factor 1.602189 x 1019 C (Coulomb) 4.803242 x 1010 esu 8.478418 x 1030 m C 2.541765 x 1018 cm esu 2.541765 D (Debyes) 4.486584 x 1040 m2 C 1.345044 x 1026 cm2 esu 1.345044 B (Buckinghams) 2.374197 x 1050 m3 C 7.117664 x 1035 cm3 esu 7.117664 U (unnamed) 1.256371 x 1060 m4 C 3.766505 x 1043 cm4 esu 3.766505 N (noname) Table 33: Calculated Multipole Moments of H20 This work Moments Spec. 7 Theor. 7 Rein exp onec twoc onec twoc 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 (B) 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 (U) 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 (N) 65 Table 33 Continued: a. Reference [79] b. Reference [32] c. From Magnetic Susceptibility, see reference [81] d. Reference [82] e. Reference [33] Table 34: 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 (B) 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 (NX100) a. Reference b. Reference c. Reference d. Reference [78] [79] [32] [83] Table 35: 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 (B) 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 (NX100) a. Reference [78] b. Expanded Rein c. Reference [84] d. Reference [85] Basis. See Reference [79] and text 68 Table 36: 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 (B) rXXXX 1.53 FXXYY 2.94 rXXZZ 4.47 TYYYY 1.04 FYYZZ 1.90 rZZZZ 6.37 (NX100) a. Reference [80] (A) Water (B) Benzene J ^'N F ,,/F lurobenzene K ) (D) Pyrazine Figure 3: Input Molecular Orientations for Multipole Moment Calculations. 9.55 8.26 8.26 \4.13 (B) (A) Figure 4: Dimer Structures of HexafluorobenzeneBenzene and BenzeneBenzene Based on Multipole Moment Interpretations CHAPTER 4 A SELFCONSISTENT REACTION FIELD MODEL INCLUDING HIGH MOLECULAR ELECTRIC MULTIPLE MOMENTS APPLIED TO THE CALCULATION OF MOLECULAR ELECTRONIC SPECTRA 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. 71 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 (4.2) its general solution being (4.3) in = E rn + ) P(T) n=0 We then require that the potential on the sphere be single valued, in(r) = ou.t(r)jr=a (4.4) 73 and due to linear independence of Legendre polynomials, this gives, for a particular term n = 1 in equations (4.1) and (4.3), At DI 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) = lBal1 (1 + 1)Dra(l+2) (4.7) Equations (4.5) and (4.7) are standard boundary conditions for cavity models, and they give (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 00 bn+(r) = + f()rn DnPn() (4.9) n=o where the second term is caused by the solvent dielectric interaction with the solute molecular multiple, and it is called the reaction field, 00 R = fn(E)DnraP(T) (4.10) n=O the factor fi(e) is called the reaction field factor, first derived by Kirkwood [86, 87]. 74 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 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 wellknown 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=( 75 (4) 4 (E1) (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 SelfConsistent 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 l,m (4.17) =Eo ft(OJMrm))(OJMrJO) l,m 76 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) l,m = (6VIHo 4) Z fz(6 1Mm0)(0 Mm1) W(6, I) + c.c. (4.19) l,rm =0 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 MIMmi 1.m (4.21) = Eo E f I(Mrm ))(OIMrlm ) l.m 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 Focklike 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) l,m 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 selfconsistency in achieved. This treatment of the molecule in solution is therefore known as the SelfConsistent 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 78 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 solvent, E, = W+ Ec,c 1 (4.25) =W+M ( M })(Mrlm) 1,m 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) 1,m 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 manyelectron wave functions. Again we desire to solve the Schr6dinger equation by assuming an approximate solution. In the CI approach the 79 approximate wave function I') is expanded as a linear combination of Slater determinants I') N 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) where Hst = (OjHIlt) (4.30) and Sst = ('sIt) (4.31) The N linearly independent eigenvectors c, the columns of C, can be chosen to be orthonormal with respect to S, N c*S cq = cStctq = 6pq (4.32) s,t 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 where 10o), in the first term of equation (4.34), is the HartreeFock reference function obtained in the SCF procedure previously described. The set of functions { ij ..') are spinadapted 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 HartreeFock 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 81 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 "HOMOLUMO 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 )) (4.37) 82 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{WH o) (OW\H )c O \CC) It is obvious that any mixing of the diagonal elements is through the offdiagonal element ()Of,1j) = (ilha) + 1 [(irlar) (irl a)] (4.39) r or its adjoint. The iath element of the Fock operator, again in the basis of spin molecular orbitals, is given by (qif1j0a) = (ill>a) + Z[(irlar) (irjra>] (4.40) r Notice that the righthandside of both equation (4.39) and equation (4.40) are identical and thus (o01H\ 7) = (<lfla) (4.41) Now by definition solving the HartreeFock eigenvalue problem requires the offdiagonal 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 83 directly interact with a closedshell HartreeFock reference determinant I'o) through the electronic Hamiltonian, or (b01HI07) = 0 (4.43) {lf})} leaves the closedshell HartreeFock 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 systems. 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 ) l,m (4.44) = (0jHo Ij) Efi(NoIMrtm o)(IM IM j) l,m 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. 84 The first deals with the absorption process in a dielectric media. According to ClassiusMosotti, the bulk solvent dielectric constant c has the following relation E1 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 1018 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) l,m 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). 85 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 = ( HI ) + ft(e) Mim o)(M I) (4.48) l,m 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) 2 E fo] ( 1M 00 (,1(4.49) l,m 1 fz(n)(VI(Mrj'iM[ 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) 1,m 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 86 a case, the energy of a state \0k) is shifted by f (n2) [I(oIM mlo) 2 (4.51) (OkkMrIOCk)((iolM71io) + (ziIMrlIi))}/2 and this leads to "Al theory", W W Wo 1 = [ IH ) ( oH 1o)]+ I f (e) f()o Mm 1 O)[(I MrI z) (4o MrmI'o)] (4.52) l,m 1 fl (n2)[(i oMrIoM (VIlMr bi)]2 l,m and to B theory W 1 W1 o= [(MilHIV) (0o8IHIo)] 1 (4.53) S:fi (n2)[(oiMrmo) (IMm )]2 l,m 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 )] 1,m (4.54) + 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) l,m 87 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) ( Results In this section results of the performance of high moments in the SCRFCI model are reported. Our HartreeFock 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 (4.57) 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. H20 Although H20 does not have interesting UVVis 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 47. The convergence this work achieved is 88 satisfactory from dipole to hexadecapole expansion, and it compares well with the results from Multiconfigurational SCRF calculation by Mikkelsen, Argen, Jensen and Helgaker [15]. 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. Acetone The n 7r* UVVis 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 48 to 413 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, 89 and the inclusion of higher moments brings additional blue shift to yield better agreement with experiment. Pyrazine 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 414. 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 414. / E ( Figure 5: Spherical Cavity Model   
Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID E9XZTFKWN_IZOQN3 INGEST_TIME 20120220T21:33:12Z PACKAGE AA00009043_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 