Baer M., Billing G.D. (eds.) - The role of degenerate states in chemistry (Adv.Chem.Phys. special issue, Wiley, 2002)

514	k. k. liang et al.

Although all of the nuclear coordinates participate in this kinetic energy operator, and in our previous discussions, all of the nuclear coordinates are expanded, with respect to an equivalent position, in power series of the parameter k, here in the specific case of a diatomic molecule, we found that only the R coordinate seems to have an equilibrium position in the molecular fixed coordinates. This means that actually we only have to, or we can only, expand the R coordinate, but not the other coordinates, in the way that

R ¼ R0 þ kR

ð53Þ

By replacing q=qR with q=kqR we have

2 h2 m q

2 q

4 h2

2

m

2

T^N ¼ k

R

k

nr0

þ

rXo

ð54Þ

2m R2 qR

qR

2m

R2

In other words, for calculating the second-order energy (the vibrational energy), we only have to keep the term to do with the interatomic distance. The other terms, then, will enter the total Schro¨dinger equation in higher orders.

Now let us use a more specific coordinate system shown in Figure 2 for H2 molecule. The z axis is taken to be along the internuclei vector R. Adapting this coordinate-system, the consequence is that only the z coordinate of the nucleus are to be expanded in powers of k. That is, we define

Z ¼ R0 þ kR

ð55Þ

Figure 2. Molecular-fixed coordinates.

the crude born–oppenheimer adiabatic approximation

515

Due to our choice of the coordinates, Z is actually identical with R. Therefore, we rewrite the following definitions:

2

Þ

h2m 1

q

2

q

^ ð

¼

R

H1

2m

R2

qR

qR

^ ð1Þ

qV

H0

¼ qR R¼R0 R

2

Þ

1 q2V

2

^ ð

¼

R¼R0 R

H0

2

qR2

1

q2V

Unð2Þ

¼ (

hfnð0Þj

R¼R0

2

qR2

				jhfnð0Þj	qV	R¼R0 jfmð0Þij2
	0		X		qR				2
	fð	Þ						R
j	n		i þ m n	Unð0Þ			Umð0Þ	)
			¼6

ð56Þ

ð57Þ

ð58Þ

ð59Þ

Thus, the equation we are going to solve for the zeroth-order nuclear motion is

m n

h2m 1 q 2 q 1 0

q2V

jhfnð0Þj

qV

0jfmð0Þij2

0

X

qZ

2 0

( 2m R2 qR R qR þ "2 hfnð Þj qZ2

Unð0Þ

Umð0Þ

#R )jwnð Þi

0jfnð Þi þ

¼6

¼ Enð2Þjwnð0Þi

ð60Þ

Equation (60) is a standard equation of motion of a harmonic oscillator. It can be easily solved as long as the term corresponding to the force constant can be evaluated. To do that with Eq. (66), we need to know how to calculate the

matrix elements

hfn

j q0

=q

Z2

0jfn

i

and

hfn

jðq

V

=q

Z

Þ0jfm

i

for the given

ð0Þ

2V

ð0Þ

ð0Þ

ð0Þ

molecular basis set

fj

n

ig

.

fð Þ

In modern quantum chemistry packages, one can obtain molecular basis set at the optimized geometry, in which the wave functions of the molecular basis are expanded in terms of a set of orthogonal Gaussian basis set. Therefore, we need to derive efficient formulas for calculating the above-mentioned matrix elements, between Gaussian functions of the first and second derivatives of the Coulomb potential terms, especially the second derivative term that is not available in quantum chemistry packages. Section IV is devoted to the evaluation of these matrix elements.

In the work of King, Dupuis, and Rys [15,16], the matrix elements of the Coulomb interaction term in Gaussian basis set were evaluated by solving the differential equations satisfied by these matrix elements. Thus, the Coulomb matrix elements are expressed in the form of the Rys’ polynomials. The potential problem of this method is that to obtain the matrix elements of the higher derivatives of Coulomb interactions, we need to solve more complicated differential equations numerically. Great effort has to be taken to ensure that the differential equation solver can solve such differential equations stably, and to

516	k. k. liang et al.

make sure that the result is accurate. In this work, we have carried out the integrals explicitly up to the second-order derivatives of the Coulomb interaction. After the lengthy derivations, it should be clear that the results we obtain, in the form of a finite series (with number of terms <12 for the extreme case in which h orbitals are taken into the basis set) of error functions and exponential functions, are simple to calculate, and the numerical accuracy is high, thanks to the existing numerical libraries providing accurate computation of error function.

The starting point of our technique is to make use of the equality

1	1		1	2	2	ð61Þ
r	¼ pp		ð 1 e t		r dt

By using the geometry defined in Figure 2, the coordinates of the electrons with respect to the nucleus will be written as

1

1

1

2

2

2

2

¼ ðxi; yi; zi þ R=2Þ ¼

pp

ð 1 e t

½xi

þyi

þðzi

þR=2Þ

&dt

rAi

1

1

1

t2

xi2

yi2

zi

R=2 2

rBi

¼ ðxi; yi; zi R=2Þ ¼ pp ð 1 e

½

þ

þð

Þ

&dt

Thus, it can be seen that the first-derivatives w. r. t. R are

q

1

1

1

2

2

2

2

¼

pp

ð 1 t2zAie t

½xi þyi

þðziþR=2Þ

&dt

qR

rAi

q

1

1

1

2

t2

xi2

yi2

zi

R=2

2

qR rBi

¼ pp ð 1 t

zBie

½

þ

þð

Þ

&dt

and the second derivatives are

q2 1

1

1

t2

2

e

t2

x2

y2

z R 2

2

¼

pp

ð 1

þ t4zAi

½ i

þ

i

þð iþ

=

Þ

&dt

qR2 rAi

2

q2

1

1

1

t2

2

e

t2

x2

y2

z

R 2

2

¼

pp

ð 1

þ t4zBi

½ i

þ

i

þð i

=

Þ

&dt

qR2 rBi

2

ð62Þ

ð63Þ

ð64Þ

ð65Þ

ð66Þ

ð67Þ

The matrix elements of these derivatives are to be evaluated with R equal to its equilibrium value R0. However, to keep the notation simple, we shall still write R in place of R0 in later text unless ambiguity may occur.

With the following discussions, it can also be seen that higher order derivatives can be evaluated with similar technique.

the crude born–oppenheimer adiabatic approximation			517
IV. MATRIX ELEMENTS OF ANGULAR-MOMENTUM-ADOPTED
GAUSSIAN FUNCTIONS
In this work, we shall use the notation
jgaðiÞi ¼	2a	3=4e aðri RaÞ2	ð68Þ
	p

to represent a normalized Gaussian wave function. Here, a is the label to specify the nuclei, which will be either A or B in our cases, and i is the label to specify the electron 1 or 2. In the context in which the label of electron is immaterial, we shall drop the i index.

To incorporate the angular dependence of a basis function into Gaussian orbitals, either spherical harmonics or integer powers of the Cartesian coordinates have to be included. We shall discuss the latter case, in which a primitive basis function takes the form

jZi ¼ Nn;l;m;a xnylzme ar2

ð69Þ

This type of basis functions is frequently used in popular quantum chemistry packages. We shall discuss the way to evaluate different kinds of matrix elements in this basis set that are often used in quantum chemistry calculation.

A.Normalization Factor

hZjZi ¼ Nn2;l;m;a ð dt x2ny2lz2me 2ar2
1	2	1	2	1	2
¼ Nn2;l;m;a ð 1 dx x2ne 2ax		ð 1 dy y2le 2ay		ð 1 dz z2me 2az
¼ Nn2;l;m;aIxIyIz					ð70Þ

By using the integral formula in Appendix A, we can expressions for Ix, Iy, and Iz, so that

N2

ð2nÞ!ð2lÞ!ð2mÞ!

p

3=2

hZjZi ¼

2a

n;l;m;a n!l!m!ð8aÞðnþlþmÞ

Therefore, the normalization factor is

"

l m 8

ðnþlþmÞ 2

3=2

#

1=2

Nn;l;m;a ¼

n! !

!ð aÞ

a

2n

! 2l ! 2m

!

p

ð

Þ ð

Þ ð Þ

obtain the proper

ð71Þ

ð72Þ

518	k. k. liang et al.

B.The Overlap Integrals

Next, we consider the simple overlap integral of two such basis functions with different powers of Cartesian coordinates and different Gaussian width, centered at different points. Let nuclei 1 locate at the origin, and let nuclei 2 locate atR, then

Z1

i ¼

Nn1;l1;m1;a1 xn1 yl1 zm1 e a1ðx2þy2þz2

Þ

ð

73

Þ

j

n2

l2

m2 e a2½ðxþXÞ2

þðyþYÞ2þðzþZÞ2

&

Z2

i ¼

Nn2;l2;m2;a2

ð

x

þ

X

Þ

ð

y

þ

Y

Þ

ð

z

þ

Z

Þ

ð

74

Þ

j

It can be observed that we can separate the overlap integral into the product of three independent spatial integrals

hZ1jZ2i ¼ Nn1;l1;m1;a1 Nn2;l2;m2;a2 IxIyIz					ð75Þ
1	2			2
Ix ¼ ð 1 dx xn1 ðx þ XÞn2 e a1x e a2			ðxþXÞ		ð76Þ
1	2		2
Iy ¼ ð 1 dy yl1 ðy þ YÞl2 e a1y e a2		ðyþYÞ			ð77Þ
1	2		2
Iz ¼ ð 1 dz zm1 ðz þ ZÞm2 e a1z e a2			ðzþZÞ		ð78Þ

Obviously, we do not need to calculate all three integrals. We shall calculate Ix and apply the formula to the other integrals.

1

2

2

Ix ¼ ð 1 dx xn1 ðx þ XÞn2 e ða1þa2Þx

2a2Xx a2X

2

1

2

2

2

¼ e a2X

ð 1 dx xn1 ðx þ XÞn2 e ða1þa2Þ½xþa2X=ða1þa2Þ& þa2X

=ða1þa2

Þ

2

1

2

¼ e a1a2X

=ða1

þa2Þ

ð 1 dx xn1 ðx þ XÞn2 e ða1þa2Þx0

ð79Þ

where x0 ¼ x þ a2X=ða1 þ a2Þ. Therefore x ¼ x0

a2X=ða1 þ a2Þ, and x þ X ¼

x0 þ a1X=ða1 þ a2Þ. Thus, by changing the dummy variables, we have

2

1

2X

n1

1X

n2

2

Ix ¼ e a1a2X

=ða1

þa2

Þ ð 1 dx x

a

x þ

a

e ða1þa2Þx

a1 þ a2

a1 þ a2

2

X X

n !

n !

n1 n

n2 m

n1

n2

a2

a1

¼ e a1a2X

=ða1

þa2

Þ n

1

2

0 m

0

n1

n !n!

n2

m

!m!

a1 a2

a1 a2

¼

¼

ð

Þ

ð

Þ

þ

þ

1

2

Xn1þn2 n m ð 1 dx xnþme ða1þa2Þx

ð80Þ

the crude born–oppenheimer adiabatic approximation

519

According to the results shown in Appendix A, the integral in Eq. (80) is nonvanishing only when ðn þ mÞ is even. Consequently, Ix is a polynomial of order n1 þ n2 in X. If n1 þ n2 is even, Ix is an even polynomial; If n1 þ n2 is odd, Ix is an odd polynomial. Therefore, we can write the result as

rþ

p

e

a1a2

2

I

x ¼ a1

a1þa2

X

Xn1þn2

a2

X X

n1!n2!ðn þ mÞ!

1 þ ð 1Þnþm

n1

n2

n 0 m

¼

0 ðn1 nÞ!n!ðn2 mÞ!m!fn þ m=2g!

2

¼

½4ða1 þ a2Þ& ðnþmÞ=2X n m

ð81Þ

Iy and Iz are obtained in exactly the same way.

C.Interaction Terms with the Nuclei

Next, we shall consider four kinds of integrals. The first is the expectation value of the Coulomb potential by one nucleus for one of the primitive basis function centered at that nucleus. The second is the expectation value of the Coulomb potential by one nucleus for one of the primitive basis function centered at a different point (usually another nucleus). Then, we will consider the matrix element of a Coulomb term between two primitive basis functions at different centers. The third case is when one basis function is centered at the nucleus considered. The fourth case is when both basis functions are not centered at that nucleus. By that we mean, for two Gaussian basis functions defined in Eqs. (73) and (74), we are calculating

1

N2

1

2

2

2

hZ1j

jZ1i ¼ p1p

ð 1 dt

ð dt x2n1 y2l1 z2m1 e 2a1r e t

r

r

1

N2

1

2n

2l

2m

2

2

2

hZ2j

jZ2i ¼ p2p

ð 1 dt

ð dt ðx þ XÞ

1 ðy þ YÞ

1 ðz þ ZÞ

1 e 2a2jrþRj

e t

r

r

1

N1N2

1

n1

n2 l1

l2 m1

m2

a2

r

R

t2r2

þ YÞ z ðz þ ZÞ

a1r2

2

hZ1jrjZ2i ¼ pp

ð 1 dt

ð dt x

ðx þ XÞ y ðy

e

e

j þ

j e

1

N1N2

1

hZ1j

jZ2i ¼

pp

ð 1 dt ð dt xn1 ðx þ XÞn2 yl1 ðy þ YÞl2 zm1 ðz þ ZÞm2

r

þ

R0

j

j

a1r2

a2

r

R

2

t2 r

R0

2

e

e

j þ

j e

j þ

j

520

k. k. liang et al.

N1 and N2 are the normalization constants. Let us consider

1

2

2

2

I1 ¼ ð1 dt ð dt x2ny2lz2me 2ar e 2t

r

1

1

2

2

1

2

2

1

2

2

¼ ð1 dt ð1 dx x2ne ð2aþt

Þx

ð1 dy y2le ð2aþt

Þy

ð1 dz z2me ð2aþt

Þz

¼ p

ðn!4n

l!4l

m!4m

ð 1

a þ

3=2

2nÞ!

ð2lÞ!

ð2mÞ!

1

dt 2

t2

ðnþlþmþ3=2Þ

ðn!4n

ð aÞ

ðnþlþmþ3=2Þ

¼ p

l!4l

m!4m

ð 1

þ 2a

2n

!

2l

!

2m

!

1

t2

3=2

Þ

ð Þ

ð

Þ

2 ðnþlþmþ3=2Þ

dt 1

ðnþlþmþ3=2Þ

ðn!4n

l!4l

m!4m

þ 2a

¼ p

ð aÞ

ð 1

p2a

2n

!

2l

!

2m

!

1

t

t2

3=2

Þ

ð Þ

ð

Þ

2 ðnþlþmþ1Þ

1

d

1

2n !

3=2

3=2

2l ! 2m !

n

l m

1

2

n

l

m

¼ p

ð

ð aÞ

ð 1

þ

t

ð þ þ þ Þ

ð Þ

n!4n

l!4l

m!4m

Þ

ð Þ

ð

Þ

2

ð þ þ þ Þ

dt 1

82

By making use of Eq. (A.12), we find

1 ¼ p

ðn!4n l!4l

m!4m

ð aÞ

ð2n þ 2l þ 2m þ 1Þ!

2n

!

2l

!

2m

!

2

n

l

m

!

2

4nþlþm

I

3=2

Þ

ð

Þ

ð

Þ

2 ðnþlþmþ1Þ

½ð

þ

þ Þ

&

p

2

þ mÞ!&24nþlþm

p2

ð2nÞ!ð2lÞ!ð2mÞ!

3=2

2½ðn

þ l

¼

a n!l!m!ð8aÞnþlþm

2a

ð2n þ 2l þ 2m þ 1Þ!

ð

Þ

¼

n;l;m;a

ð2n þ 2l þ 2m þ 1Þ!

p2aN 2

2½ðn

þ l

þ mÞ!&

4nþlþm

83

Therefore

ð

þ

þ

þ

Þ

r

1

2½ðn1

þ l1 þ m1

Þ!&24n1

þl1

þm1

2a1

84

hZ1jrjZ1i ¼

1 !

p

ð

Þ

2n1

2l1

2m1

There are two kinds of integrals involving two different centers. The first case is

1

N2

1

2

2

2

hZ2j

jZ2i ¼ p2p

ð 1 dt

ð dtðx þ XÞ2n2 ðy þ YÞ2l2 ðz þ ZÞ2m2 e 2a2jrþRj e t

r

r

N2

1

1

2n2

2a2

x X 2 t2x2

2

ð 1 dt

ð 1 dxðx þ XÞ e ð þ Þ e

¼ pp

1

2l2

2a2

y

Y 2

t2y2

1

2m2

2a2

z Z

2

t2z2

ð 1dyðy þ YÞ

e

ð

þ Þ e

ð 1dzðz þ ZÞ

e

ð þ

Þ e

ð85Þ

the crude born–oppenheimer adiabatic approximation

521

Consider the integral

1

1

2

2

2

2

I2 ¼ ð 1 dt

ð 1 dx x2ne 2ax

e t

ðx

2XxþX

Þ

1

2

2

2

2

1

2

2

2

2

ð 1 dy y2le 2ay

e t

ðy

2YyþY

Þ

ð 1 dz z2me 2az

e t

ðz

2ZzþZ

Þ

1

2

2

1

2

2

2

¼ ð 1 dt e t

R

ð 1 dx x2ne ð2aþt

Þx

þ2t

Xx

1

2

2

2

1

2

2

2

ð 1 dy y2le ð2aþt

Þy

þ2t

Yy ð 1 dz z2me ð2aþt

Þz

þ2t

Zz

1

2

2

1

2

2

2

2

4

2

2

¼ ð 1 dt e t

R

ð 1 dx x2ne ð2aþt

Þ½x t

X=ð2aþt

Þ&

et

X

=ð2aþt

Þ

1

2

2

2

2

4

2

2

ð 1 dy y2le ð2aþt

Þ½y t

Y

=ð2aþt

Þ&

et

Y

=ð2aþt

Þ

1

2

2

2

2

4

2

2

ð 1 dz z2me ð2aþt

Þ½z t

Z=ð2aþt

Þ&

et

Z

=ð2aþt

Þ

1

1

t2X

2n

2

2

2

2

2

¼ ð 1 dt e 2at

=ð2aþt

ÞR

ð 1 dx x þ

e ð2aþt

Þx

2a þ t2

1

t2Y

2l

1

t2Z

2m

2

2

2

2

ð 1 dy y þ

e ð2aþt

Þy

ð 1 dz

z þ

e ð2aþt

Þz

2a þ t2

2a þ t2

1

2

2

2

"

2n

2n !

t2X

2n n 1

2

2

¼

dt e 2at

R

=ð2aþt

Þ

2nð

Þ

xne 2aþt

x dx

n !n! 2a

t2

ð 1

n¼

ð

Þ

þ

ð 1

"

2l

X

2l l 1

2l !

t2Y

2 2

ð

Þ

yle

2aþt

Þ

y dy

2l

l !l! 2a

t2

l¼0

ð

Þ

þ

ð 1

"

X

2m m 1

2m

2m !

t2Z

2

2

ð

Þ

zme

ð

2aþt

Þ

z

dz

2m

m !m! 2a

t2

m¼

ð

Þ

þ

ð 1

1

X

n

2n 2n

1

2

2

2

2n

t2X

2

2

¼

dt e 2at

R

=ð2aþt

Þ

"

ð

Þ!

x2ne ð2aþt

Þx dx#

n

0

2n

2n ! 2n !

2a t2

ð 1

¼

ð

Þ ð

Þ

þ

ð 1

l

X

t2Y

2l 2l

1

2l

2

2

"

ð Þ!

y2le ð2aþt

Þy

dy#

l¼0

ð

Þ ð Þ

þ

ð 1

X

2m 2m

m

m

2

Z

1

2

2

"

ð2 Þ!

t

z2me ð2aþt

Þz

dz#

m¼

0

ð

2m

2m ! 2m !

2a

t2

Þ ð Þ

þ

ð 1

ð1

X

¼dt e 2at2R2=ð2aþt2Þ

1

522

k. k. liang et al.

r

n

n

2n

t2X

2n 2n

2

"

X

ð

Þ!

ð nÞ!

p

#

0

2n

2n ! 2n !

2a

þ

t2

4nn!

2a

t2 n

2a

t2

"

¼

ð

Þ ð Þ

2l 2l

ð

þ

Þ

þ

#

l¼0

ð

Þ ð

Þ

þ

rþ

l

2l !

t2Y

2l !

p

X

2l

ð Þ

4ll!

ð Þ

t2

l

2l

! 2l !

2a t2

2a

2a t2

2m 2m

ð

þ

Þ

r

m

m

2m !

t2Z

2m !

"

X

2m

ð

Þ

4mm!

ð Þ

t2 m

#

0

¼

ð

Þ ð Þ

þ

ð þ Þ

þ

X X X

n

l

m

ð2nÞ!ð2lÞ!ð2mÞ!p3=2X2n 2nY2l 2lZ2m 2m

¼

2

2n

2n !n! 2l

2l !l! 2m

2m !m!4nþlþm

n¼0 l¼0 m¼0 ð

Þ ð

Þ ð

Þ

1

2

2

2

1

t2

2ðnþlþm n l mÞ

ð0

dt e 2at

R

=ð2aþt

Þ

ð2a þ t2Þnþlþmþ3=2

2a þ t2

Now, let

x p

t

2

Since

2a þ t

dx

¼

2a

dt

ð2a þ t2Þ3=2

we find

dt

¼

ð2a þ t2Þ3=2

d

x

2a

Also, noticing that when t ¼ 0, x ¼ 0, and x ¼ 1 when t ¼ 1, and

2a

¼ 1 x

2

we have

2a þ t2

1

X X X

2n ! 2l ! 2m !

3=2X2n 2nY2l 2lZ2m 2m

n

l

m

I

ð

Þ ð

Þ ð

Þ p

2

¼ a n

¼

0

l

¼

0

m

¼

0

ð2n 2nÞ!n!ð2l 2lÞ!l!ð2m 2mÞ!m!4nþlþm

2a

nþlþm

ð

0

ð

Þ

1 d

e 2aR2x2 1 x2

2½2ðnþlþm n l mÞ&

x

nþlþm

x

1 n

l

m

ð

2n !

2l !

2m !

p

3=2X2n 2nY2l 2lZ2m 2m

¼

X X X

Þ ð

Þ ð

Þ

a

2n

2n

!n!

2l

2l

!l! 2m

2m

!m! 8a nþlþm

ð

Þ

ð

Þ

ð

Þ

ð Þ

n¼0 l¼0 m¼0

X

ðn

þ l þ mÞ!

1

1 d

2½2ðnþlþm n l mÞþk&e 2aR2x2

nþlþm

k

k¼0

ðn þ l þ m kÞ!k!

ð Þ

ð0

x x

ð86Þ

ð87Þ

ð88Þ

ð89Þ

ð90Þ

ð91Þ

the crude born–oppenheimer adiabatic approximation

According to Eq. (A.16),

r

ð

0

ð

a

Þ s 1

1

2se 2aR2x2

ð2sÞ!

p

p2

d

x x

¼ s! 8 R2 s

1

erf

a

R2

2 2aR2

X

j!

j#

e 2aR2

j

0

ð

2j

þ

1

Þ

!

8aR2

¼

where s ¼ 2ðn þ l þ m n l mÞ þ k. Also,

1 n l

m

2n ! 2l ! 2m !

p

3=2X2n 2nY2l 2lZ2m 2m

I

X X X

ð

Þ ð

Þ ð

Þ

2

¼ a

2n

!l! 2m

2m

!m! 8a

nþlþm

ð

2n

!n! 2l

2l

Þ

Þ

ð

Þ

ð

Þ

ð

n¼0 l¼0 m¼0

X

ðn

þ l mÞ!

1

1 d

2½2ðnþlþm n l mÞþk&e 2aR2x2

nþlþm

k

k¼0

þ

ð

Þ

ð0

x x

ðn þ l þ m kÞ!k!

N 2

X X X

nþlþm n l m

2n

2l

2m 2m

n

l

m

2n

2l

¼ ð

2

aÞ

3=2

n;l;m;a

n!l!m!ð8aÞ

X

Y

Z

a

n 0 l 0 m 0

ð

2n

2

nÞ

n ð

2

lÞ

l ð

2m

2

mÞ

m

!

¼

¼

¼

!

!

2l

!

!

!

523

ð92Þ

nþXlþm

k¼0

Therefore

1

hZ2jrjZ2i ¼

ðn l mÞ!

1

k

1

2½2ðnþlþm n l mÞþk&e 2aR2x2

d

þ þ

ð

Þ

ð0

x x

ðn þ l þ m kÞ!k!

r

n

0 l 0

m 0

8a2

X X X

n2

l2

m2

p

¼

¼

¼

2m nþlþm

n !l !m ! 8

n2þl2þm2 n l mX2n2 2nY2l2 2lZ2m2

2 2

2 ð a2Þ

X

ð2n2 2nÞ!n!ð2l2 2lÞ!l!ð2m2 2mÞ!m!

k

¼

0

ðn þ l þ mÞ!

k

1

2½2ðn2þl2þm2 n l mÞþk&e 2a2

1

d

ð

Þ

ð0

x x

ðn þ l þ m kÞ!k!

ð93Þ

R2x2

ð94Þ

Next, we consider

1 hZ1jrjZ2i

¼	N1N2	1	2	2 2	r	2
	pp	ð1 dt	ð dt xn1 ðx þ XÞn2 yl1 ðy þ YÞl2 zm1 ðz þ ZÞm2 e a1r	a2jrþRj t