Evaluation of Gaussian Molecular Integrals

I. Overlap Integrals

Minhhuy Hô 
Julio Manuel Hernández-Pérez

This article discusses the evaluation of molecular overlap integrals for Gaussian-type functions with arbitrary angular dependence. As an example, we calculate the overlap matrix for the water molecule in the STO-3G basis set.

Introduction

Computational quantum chemistry makes extensive use of various integrals (and their derivatives) of the general form [, , ]

(∫_(-∞)^∞)(Subscript[χ, a])(r)O (Subscript[χ, b])(r)dr,

where (Subscript[χ, a])(r) is an unnormalized Cartesian Gaussian function centered at A={(Subscript[A, x]),(Subscript[A, y]),(Subscript[A, z])}:

({{(Subscript[χ, a])(r;α,A,a)=(((x-(Subscript[A, x])))^(Subscript[a, x]))(((y-(Subscript[A, y])))^(Subscript[a, y]))(((z-(Subscript[A, z])))^(Subscript[a, z]))(e^(-α|r-A((|)^2))),}})

where A is normally taken at the nucleus, α is the orbital exponent, and the polynomial represents the angular part, in that the sum of the Cartesian angular momenta (Subscript[a, x])+(Subscript[a, y])+(Subscript[a, z])=0,1,2,…  corresponds to functions of type s, p, d, f, …. When the operator O is 1, one simply has the overlap/density integral; otherwise it can be the energy operator for kinetic energy -(1/2)(∇^2), electron-nuclear attraction |r-R((|)^(-1)), or electron-electron repulsion |(Subscript[r, i])-(Subscript[r, j])((|)^(-1)) (which would involve double integrals). Other molecular properties involving external fields (response functions) or transition moments can also be computed from integrals of this form.

Gaussian-Type Functions

Gaussian-type functions are not the most natural choice for expanding the wavefunction. Slater-type functions, where the exponent is -α|r-A| instead, can describe atomic systems more realistically. However, complications quickly arise in molecular calculations, which has led to the use of Gaussian functions in the overwhelming majority of current computational programs. Gaussian functions possess several desirable computational properties [, , ]. (Much credit is due to S. F. Boys for the introduction of the Cartesian Gaussian function into computational chemistry and its early development in a series of 12 papers under the general title “Electronic Wavefunctions,” the first being [].) A given Slater function can be approximated by a linear combination of several Gaussians.

The first useful property is that the product of two Gaussian functions located at A and B is another Gaussian located at a point P somewhere between A and B. (The proof of this can be found in [].) The product of two Gaussians (Subscript[χ, a])(r;α,A,a) and (Subscript[χ, b])(r;β,B, b) is:

({{(Subscript[χ, a])(r;α,A,a)}})×(Subscript[χ, b])(r;β,B,b)=(((x-(Subscript[A, x])))^(Subscript[a, x]))(((y-(Subscript[A, y])))^(Subscript[a, y]))(((z-(Subscript[A, z])))^(Subscript[a, z]))(e^(-α|r-A((|)^2)))×(((x-(Subscript[B, x])))^(Subscript[b, x]))(((y-(Subscript[B, y])))^(Subscript[b, y]))(((z-(Subscript[B, z])))^(Subscript[b, z]))(e^(-β|r-B((|)^2)))=(Subscript[E, A,B])(e^(-(α+β)|r -P((|)^2))),

with (Subscript[E, A,B])=((e^(-((αβ)/(α+β))|A-B|))^2) and P=((α A+β B)/(α+β)).

A second desirable property is that a derivative of a Gaussian can be expressed as a sum of Gaussians of lower and higher Cartesian angular values.

Overlap Integrals

The simplest molecular integral is the overlap integral S=(∫_(-∞)^∞)({{(Subscript[χ, a])(r;α,A,a)}})×(Subscript[χ, b])(r;β,B, b)dr. We first separate the integral into its orthogonal components:

S=(Subscript[E, A,B])(Subscript[S, x])((Subscript[a, x]),(Subscript[b, x]))(Subscript[S, y])((Subscript[a, y]),(Subscript[b, y]))(Subscript[S, z])((Subscript[a, z]),(Subscript[b, z])),

where the notation (Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])) expresses its functional dependence on the Cartesian angular components. The x component, for instance, is:

(Subscript[S, x])((Subscript[a, x]),(Subscript[b, x]))=(∫_(-∞)^∞)(((x-(Subscript[A, x])))^(Subscript[a, x]))(((x-(Subscript[B, x])))^(Subscript[b, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx.

Using a binomial expansion in the polynomial part,

(((x-(Subscript[A, x])))^(Subscript[a, x]))(((x-(Subscript[B, x])))^(Subscript[b, x]))=(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))(((x-(Subscript[P, x])))^((Subscript[i, x])+(Subscript[j, x])))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))

yields:

(Subscript[S, x])((Subscript[a, x]),(Subscript[b, x]))=(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))(∫_(-∞)^∞)(((x-(Subscript[P, x])))^((Subscript[i, x])+(Subscript[j, x])))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx.

Odd values of (Subscript[i, x])+(Subscript[j, x]) result in odd functions whose integrals vanish. For even values of (Subscript[i, x])+(Subscript[j, x]), a solution [] for the integrals is given by ((((π/(α+β))))^(1/2))((((Subscript[i, x])+(Subscript[j, x])-1) !!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2))), and, in those cases:

(Subscript[S, x])((Subscript[a, x]),(Subscript[b, x]))=(Sqrt[(π/(α+β))])(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}})) ((((Subscript[i, x])+(Subscript[j, x])-1) !!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x]))).

We keep in mind that for the summations only terms of even values of (Subscript[i, x])+(Subscript[j, x]) survive. We have thus obtained an expression to evaluate the overlap of two Gaussians with arbitrary Cartesian angular factors. Equation () should be sufficient for programming or even for manual evaluation with small Cartesian angular values. One can, however, further reduce the number of operations by using a recurrence relation, one of the most useful techniques of computational chemistry. Recurrence relations let us efficiently calculate molecular integrations of higher angular values using previously obtained results with lower angular values. Recurrence relations are used in most computational chemistry programs and their application to other molecular integrations will be shown here.

Recurrence Relations

The derivation of the following relations involves straightforward algebraic manipulations, but is rather lengthy. Its omission does not impede our understanding of the recurrence relation, but helps maintain a reasonable continuity in our discussion. Interested readers can find a detailed derivation in the Appendix.

We start by defining the function

(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=(Sqrt[((α+β)/π)])(Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])),

where we have removed the factor (Sqrt[π/(α+β)]) from (Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])) in (). If we take the derivative of (Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) with respect to the nuclear coordinate (Subscript[A, x]) using the definition of (Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])) in equation (),

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[B, x])))^(bx))(∂/(∂(Subscript[A, x])))(((x-(Subscript[A, x])))^(ax))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx=-(Subscript[a, x])(Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x]))+2α (Subscript[s, x])((Subscript[a, x])+1,(Subscript[b, x]))+2α((Subscript[A, x])-(Subscript[P, x]))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])),

and combining the result with the same derivative of (Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])), but using the definition of (Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])) in equation () instead,

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}})) ((((Subscript[i, x])+(Subscript[j, x])-1) !!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))(∂/(∂(Subscript[A, x])))[((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))]=(Subscript[a, x])((α/(α+β))+1) (Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x]))+(Subscript[b, x])((α/(α+β)))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])-1),

we obtain

(Subscript[s, x])((Subscript[a, x])+1,(Subscript[b, x]))=-((Subscript[A, x])-(Subscript[P, x]))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) +((Subscript[a, x])/(2(α+β)))(Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x])) +((Subscript[b, x])/(2(α+β)))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])-1).

The same approach in which we combine the derivatives of () and () with respect to the nuclear coordinate (Subscript[B, x]) gives

(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])+1)=-((Subscript[B, x])-(Subscript[P, x]))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) +((Subscript[b, x])/(2(α+β)))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])-1) +((Subscript[a, x])/(2(α+β)))(Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x])).

For (Subscript[b, x])=0, we have the index recursion relation

(Subscript[s, x])((Subscript[a, x])+1,0)=-((Subscript[A, x])-(Subscript[P, x]))(Subscript[s, x])((Subscript[a, x]),0)+((Subscript[a, x])/(2(α+β)))(Subscript[s, x])((Subscript[a, x])-1,0),

and combining equations () and () yields the transfer equation

(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])+1)=(Subscript[s, x])((Subscript[a, x])+1,(Subscript[b, x])) +((Subscript[A, x])-(Subscript[B, x]))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) .

Starting with the initial values (Subscript[s, x])(0,0)=1 and (Subscript[s, x])(1,0)=-((Subscript[A, x])-(Subscript[P, x])), the recurrence relation and the transfer equation, we can build up the overlap of functions of higher Cartesian angular values from lower ones. This is particularly useful for contracted Gaussian basis primitives with different functions sharing the same exponent.

Implementation

The function Ov evaluates the overlap integral of two Gaussian functions; here alpha, beta, RA, RB, LA, and LB are α, β, A, B, a, and b as defined earlier.

Ov[alpha_,beta_,RA_,RB_,LA_,LB_]:=Module[{EAB,Overlap},
Do[
(* Initial Conditions *)

s[i_,0,0]:=1;
s[i_,1,0]:=-(RA[[i]]-((alpha*RA[[i]]+beta*RB[[i]])/(alpha+beta)));

(* Recurrence Index *)

s[i_,a_,0]:=-(RA[[i]]- (alpha*RA[[i]]+beta*RB[[i]])/(alpha+beta))*s[i,a-1,0]+((a-1)/(2*(alpha+beta)))*s[i,a-2,0];

(* Transfer Equation *)

s[i_,a_,b_]:= s[i,a+1,b-1]+(RA[[i]]-RB[[i]])*s[i,a,b-1]

,{i,1,3}];

EAB=Exp[-(alpha*beta/(alpha+beta))*(RA-RB).(RA-RB)];

Overlap=EAB*(Pi/(alpha+beta))^(3/2)*s[1,LA[[1]],LB[[1]]]*
s[2,LA[[2]],LB[[2]]]*
s[3,LA[[3]],LB[[3]]]
];

Examples

The published contracted Gaussian basis sets (see, for example, []) are usually not normalized; in our first example, we will calculate the normalization factor of the s and p Cartesian Gaussian functions that we will need later on. The normalization factor is simply the inverse square root of the overlap integral. Here we calculate the overlap between two s functions with the Cartesian angular vector a={0,0,0},

(Subscript[S, s])=(((Subscript[E, A,A])((π/(α+α))))^(3/2))(Subscript[s, x])(0,0)(Subscript[s, y])(0,0)(Subscript[s, z])(0,0)=((((π/(2 α))))^(3/2)),

and observe that (Subscript[E, AA])= (Subscript[s, x])(0,0)=(Subscript[s, y])(0,0)=(Subscript[s, z])(0,0)=1. The normalization factor for 1s functions is 1/(Sqrt[(Subscript[S, s])])=(((2α/π))^(3/4)). Similarly, the overlap between two (Subscript[p, x]) functions (a={1,0,0}) is

(Subscript[S, (Subscript[p, x])])=(((Subscript[E, A,A])((π/(α+α))))^(3/2))(Subscript[s, x])(1,1)(Subscript[s, y])(0,0)(Subscript[s, z])(0,0)=((((π/(α+α))))^(3/2))(Subscript[s, x])(2,0)=((((π/(α+α))))^(3/2))(1/(4α)),

where we have analogously used the recurrence relations to obtain (Subscript[s, x])(2,0) and afterward, (Subscript[s, x])(1,1). Similar results are obtained for (Subscript[S, (Subscript[p, y])]) and (Subscript[S, (Subscript[p, z])]). The normalization factor for the p function is then (((128(α^5)/(π^3)))^(1/4)). These results are special cases of the more general formula of the normalization factor

N(α,a)=(((((2α)/π)))^(3/4))((((4α))^(((Subscript[a, x])+(Subscript[a, y])+(Subscript[a, z]))/2))/((((2(Subscript[a, x])-1)!!(2(Subscript[a, y])-1)!!(2(Subscript[a, x])-1)!!))^(1/2))),

and we note that this value depends only on the orbital exponent and the Cartesian angular values.

NormCoeff[α_,{ax_,ay_,az_}] :=(((((2 α)/π)))^(3/4))((((4 α))^((ax+ay+az)/2))/((((2 ax-1)!!(2 ay-1)!!(2 az-1)!!))^(1/2))) ;

In the second example, we will calculate the overlap matrix of the water molecule ((Subscript[r, OH])=1.86942bohr, (Subscript[∡, HOH])=100.0269°, geometry optimized at the HF/STO-3G level). The molecule lies in the y-z plane with Cartesian coordinates in atomic units.

R={{0.,1.43233673,-0.96104039},
{0.,-1.43233673,-0.96104039},
{0.,0.,0.24026010}};

In the STO-3G basis set, each atomic orbital is described by a sum of three Gaussians; here are their primitive contraction coefficients and orbital exponents (taken from []).

PrimCoeff={{0.1543289673,0.5353281423,0.4446345422},
{0.1543289673,0.5353281423,0.4446345422},
{0.1543289673,0.5353281423,0.4446345422},
{-0.09996722919,0.3995128261,0.7001154689},
{0.155916275,0.6076837186,0.3919573931},
{0.155916275,0.6076837186,0.3919573931},
{0.155916275,0.6076837186,0.3919573931}};

OrbCoeff ={{3.425250914,0.6239137298,0.168855404},
{3.425250914,0.6239137298,0.168855404} ,
{130.7093214,23.80886605,6.443608313},
{5.033151319,1.169596125,0.38038896} ,
{5.033151319,1.169596125,0.38038896},
{5.033151319,1.169596125,0.38038896},
{5.033151319,1.169596125,0.38038896}};

Here are the centers and Cartesian angular values of the orbitals, in the following order: (Subscript[H, 1s]), (Subscript[H, 2s]), (Subscript[O, 1s]), (Subscript[O, 2s]), (Subscript[O, 2(Subscript[p, x])]), (Subscript[O, 2(Subscript[p, y])]), and (Subscript[O, 2(Subscript[p, z])]).

FCenter ={R[[1]],R[[2]],R[[3]],R[[3]],R[[3]],R[[3]],R[[3]]};

CartAng={{0,0,0},{0,0,0},{0,0,0},{0,0,0},
{1,0,0},{0,1,0},{0,0,1}};

For instance, the s atomic orbital of the hydrogen atom 1 is described as

(Subscript[ϕ, 1])(r)=(∑_(i=1)^3)N((Subscript[α, 1i]),a)×(Subscript[c, 1i])×(Subscript[χ, 1i])(r;(Subscript[α, 1i]),(Subscript[R, H]),a)=1.7944418×0.1543289673×(e^(-3.425250914 |r -(Subscript[R, H])((|)^2)))+0.50032649221111580306597943490487523377`8.× 0.5353281423×(e^(-0.6239137298 |r -(Subscript[R, H])((|)^2)))+0.18773546×0.4446345422×(e^(-0.168855404 |r -(Subscript[R, H])((|)^2))).

Similarly, the (Subscript[p, z]) orbital of the oxygen atom is

(Subscript[ϕ, 7])(r)=10.745833×0.155916275× (z-0.24026010)×(e^(-5.033151319 |r -(Subscript[R, O])((|)^2)))+1.7337441×0.6076837186×(z-0.24026010)×(e^(-1.169596125 |r -(Subscript[R, O])((|)^2)))+0.4258189334467701314679288771003484726`8.×0.391957393×(z-0.24026010)×(e^(-0.38038896 |r -(Subscript[R, O])((|)^2))).

The overlap integral between the two orbitals expands to nine integrals involving the primitives. Using the function Ov, for example, the overlap integral between the first two primitives of equations () and () is

(∫_(-∞)^∞)(Subscript[χ, 11])(r)×(Subscript[χ, 71])(r)dr=(∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(e^(-3.425250914 |r -(Subscript[R, H])((|)^2)))×(z-0.24026010)×(e^(-5.033151319 |r -(Subscript[R, O])((|)^2)))dx dy dz.

Ov[OrbCoeff[[1,1]],OrbCoeff[[7,1]],
FCenter[[1]],FCenter[[7]],CartAng[[1]],CartAng[[7]]]

And the overlap between (Subscript[ϕ, 1])(r) and (Subscript[ϕ, 7])(r) is

(Subscript[S, 17])=(∫_(-∞)^∞)(Subscript[ϕ, 1])(r)×(Subscript[ϕ, 7])(r)dr.

Sum[
NormCoeff[OrbCoeff[[1,i]],CartAng[[1]]]*
NormCoeff[OrbCoeff[[7,j]],CartAng[[7]]]*
PrimCoeff[[1,i]]*
PrimCoeff[[7,j]]*
Ov[
OrbCoeff[[1,i]],OrbCoeff[[7,j]],
FCenter[[1]],FCenter[[7]],
CartAng[[1]],CartAng[[7]]
]
,{i,3},{j,3}]

The overlap matrix S for the entire water molecule in the STO-3G basis set can be calculated in a similar manner.

Table[Sum[
NormCoeff[OrbCoeff[[p,i]],CartAng[[p]]]*
NormCoeff[OrbCoeff[[q,j]],CartAng[[q]]]*
PrimCoeff[[p,i]]*
PrimCoeff[[q,j]]*
 Ov[
OrbCoeff[[p,i]],OrbCoeff[[q,j]],
FCenter[[p]],FCenter[[q]],
CartAng[[p]],CartAng[[q]]
]
,{i,3},{j,3}],{p,7},{q,7}]//MatrixForm//Chop

Since the overlap matrix is symmetrical, we need to calculate only the elements above the main diagonal. The basis functions are normalized, as indicated by the unit diagonal elements. We note that (Subscript[S, 15]) equals zero. This is the overlap between the 1s orbital of the first hydrogen atom and the 2(Subscript[p, x]) of the oxygen atom. The molecule lies in the y-z plane, so this overlap vanishes due to symmetry. Similar cases occur in the overlaps between the s and p orbitals of O ((Subscript[S, 35]),(Subscript[S, 36]),(Subscript[S, 37]),…) and among the 2p orbitals ((Subscript[S, 56]),(Subscript[S, 57]),(Subscript[S, 67])), etc. The strongest overlaps are those between the 1s hydrogen orbitals and the 2s oxygen orbital.

We next plot the atomic-orbital overlap between 1s of (Subscript[H, 1]) and 2(Subscript[p, z]) of O,

(∑_(m=1)^3)(∑_(n=1)^3)(Subscript[c, 1m])(Subscript[c, 7n])(Subscript[χ, m])(r;(Subscript[α, m]),(Subscript[R, H]), (Subscript[a, m]))×(Subscript[χ, n])(r;(Subscript[β, n]),(Subscript[R, O]), (Subscript[b, n])),

in the y-z plane, superimposing the molecule structure.

spz=Sum[PrimCoeff[[1,i]]PrimCoeff[[3,j]]
Exp[-OrbCoeff[[1,i]]
({y,z}-{R[[1,2]],R[[1,3]]}).
({y,z}-{R[[1,2]],R[[1,3]]})]
(z-R[[3,3]])Exp[-OrbCoeff[[1,j]]
({y,z}-{R[[3,2]],R[[3,3]]}).
({y,z}-{R[[3,2]],R[[3,3]]})]
,{i,3},{j,3}];

spzFig=ContourPlot[spz,{y,-3,3},{z,-3,3},Contours->25,PlotPoints->100,PlotRange->All];

Molecule=Graphics[{Thickness[0.02],Opacity[0.3],Pink,Line[{{R[[1,2]],R[[1,3]]},{R[[3,2]],R[[3,3]]}}],Line[{{R[[2,2]],R[[2,3]]},{R[[3,2]],R[[3,3]]}}]}];

Show[spzFig,Molecule]

We observe a strong distortion of the positive (lower) lobe of the 2(Subscript[p, z]) function toward the hydrogen atom. The negative lobe shows less deformation, whereas the node remains precisely at the atomic position. Note that this is simply the orbital overlap between noninteracting atoms, such as in the case of the “promolecule.”

Conclusion

We have provided an introduction to the evaluation of molecular integrals involving Gaussian-type basis functions both analytically and by use of recurrence relations. The results are general and relatively straightforward; the simple algorithm makes it suitable for implementation in a number of programming languages. Together with the kinetic, nuclear-electron attraction, and electron-electron repulsion energies, this is the first step toward the calculation of molecular energies and electronic properties.

Appendix

Here we provide the derivation of equations () and (). Differentiate (Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) using equation () with respect to (Subscript[A, x]):

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[B, x])))^(Subscript[b, x]))(∂/(∂(Subscript[A, x])))(((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx.

Consider the derivative term

(∂/(∂(Subscript[A, x])))(((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))=(e^(-(α+β)(((x-(Subscript[P, x])))^2)))(∂/(∂(Subscript[A, x])))(((x-(Subscript[A, x])))^(Subscript[a, x]))+(((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))(∂/(∂(Subscript[A, x])))[-(α+β)(((x-(Subscript[P, x])))^2)];

the first term is simply

-((Subscript[a, x]))(((x-(Subscript[A, x])))^((Subscript[a, x])-1))(e^(-(α+β)(((x-(Subscript[P, x])))^2))),

and, using the chain rule, the derivative of the second term is

-(α+β)(∂/(∂(Subscript[A, x])))[(((x-(Subscript[P, x])))^2)]=(α+β)2(x-(Subscript[P, x]))(∂/(∂(Subscript[A, x])))[((α (Subscript[A, x])+β (Subscript[B, x]))/(α+β))]=(α+β)2(x-(Subscript[P, x]))(α/(α+β))=2α(x-(Subscript[P, x])),

recalling that (Subscript[P, x])=(α (Subscript[A, x])+β (Subscript[B, x]))/(α+β). Substitute the results into the expression for ∂(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))/∂(Subscript[A, x]),

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[B, x])))^(Subscript[b, x]))[-((Subscript[a, x]))(((x-(Subscript[A, x])))^((Subscript[a, x])-1))(e^(-α+β)(((x-(Subscript[P, x])))^2)))+2α(x-(Subscript[P, x]))(((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))]dx.

Writing x-(Subscript[P, x])=(x-(Subscript[A, x]))+((Subscript[A, x])-(Subscript[P, x])) in the second term inside the bracket, after expanding we have

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))=-(Subscript[a, x])(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[A, x])))^((Subscript[a, x])-1))(((x-(Subscript[B, x])))^(Subscript[b, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx+2α((Subscript[A, x])-(Subscript[P, x]))(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[A, x])))^(Subscript[a, x]))(((x-(Subscript[B, x])))^(Subscript[b, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx+2α(Sqrt[((α+β)/π)])(∫_(-∞)^∞)(((x-(Subscript[A, x])))^((Subscript[a, x])+1))(((x-(Subscript[B, x])))^(Subscript[b, x]))(e^(-(α+β)(((x-(Subscript[P, x])))^2)))dx.

Comparing the three integrals with the definitions of (Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) and (Subscript[S, x])((Subscript[a, x]),(Subscript[b, x])) we have the desired equation ().

To derive equation (), we differentiate (Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])) with respect to (Subscript[A, x]) using equation () instead,

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))
=(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))(∂/(∂(Subscript[A, x])))[((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))]
=(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))[((Subscript[a, x])-(Subscript[i, x]))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])-1))((α/(α+β))-1)((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))+((Subscript[b, x])-(Subscript[j, x]))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])-1))((α/(α+β)))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))]
=((α/(α+β))-1)(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))[((Subscript[a, x])-(Subscript[i, x]))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])-1))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))]+((α/(α+β)))(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))[((Subscript[b, x])-(Subscript[j, x]))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])-1))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))].

We make use of the relationship

(({{(Subscript[a, x])}, {(Subscript[i, x])}}))((Subscript[a, x])-(Subscript[i, x]))=((((Subscript[a, x])-1)!(Subscript[a, x])((Subscript[a, x])-(Subscript[i, x])))/((Subscript[i, x])!1×2×3×…×((Subscript[a, x])-(Subscript[i, x])-1)((Subscript[a, x])-(Subscript[i, x]))))=(Subscript[a, x])(({{(Subscript[a, x])-1}, {(Subscript[i, x])}}))

to reduce the derivative above to

(∂/(∂(Subscript[A, x])))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x]))
=(Subscript[a, x])((α/(α+β))-1)(∑_((Subscript[i, x])=0)^((Subscript[a, x])-1))(∑_((Subscript[j, x])=0)^(Subscript[b, x]))(({{(Subscript[a, x])-1}, {(Subscript[i, x])}}))(({{(Subscript[b, x])}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))[((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])-1))((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])))]+(Subscript[b, x])((α/(α+β)))(∑_((Subscript[i, x])=0)^(Subscript[a, x]))(∑_((Subscript[j, x])=0)^((Subscript[b, x])-1))(({{(Subscript[a, x])}, {(Subscript[i, x])}}))(({{(Subscript[b, x])-1}, {(Subscript[j, x])}}))((((Subscript[i, x])+(Subscript[j, x])-1)!!)/(([2(α+β)])^(((Subscript[i, x])+(Subscript[j, x]))/2)))[((((Subscript[P, x])-(Subscript[B, x])))^((Subscript[b, x])-(Subscript[j, x])-1))((((Subscript[P, x])-(Subscript[A, x])))^((Subscript[a, x])-(Subscript[i, x])))]
=(Subscript[a, x])((α/(α+β))-1)(Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x]))+(Subscript[b, x])((α/(α+β)))(Subscript[s, x])((Subscript[a, x]),(Subscript[b, x])-1).

This is equation (), which is what was needed to prove.

References

A. Szabo and N. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, New York: Dover, 1996.

T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory, New York: Wiley, 2000.

F. Jensen, Introduction to Computational Chemistry, New York: Wiley, 1999.

I. Shavitt, “The Gaussian Function in Calculations of Quantum Mechanics,” in Methods in Computational Physics, Vol. 2, B. Alder, S. Fernbach, and M. Rotenberg, eds., New York: Academic Press, 1963 pp. 1–45.

S. F. Boys, “A General Method of Calculation for Stationary States of Any Molecular System,” Proceedings of the Royal Society A, 200, 1950 pp. 542–554. rspa.royalsocietypublishing.org/content/200/1063/542.abstract.

H. Taketa, S. Huzinaga, and K. O-ohata, “Gaussian-Expansion Methods for Molecular Integrals,” Journal of the Physical Society of Japan, 21, 1966 pp. 2313–2324.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Dover, 1965, Eq. (7.4.4) p. 302.

“Basis Set Exchange.” (Sep 13, 2011) bse.pnl.gov/bse/portal.