Evaluation of Gaussian Molecular Integrals

II. Kinetic-Energy Integrals

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

This article carries out the evaluation of kinetic energy integrals using Gaussian-type functions with arbitrary Cartesian angular values. As an example, we calculate the kinetic matrix for the water molecule in the STO-3G basis set.

Introduction

In this article, the second of a series describing algorithms for evaluating molecular integrals, we detail the evaluation of the kinetic energy integrals. Detailed accounts of molecular integrals can be found in the references of []. The electronic kinetic energy in atomic units ((Subscript[m, e])=e=ℏ=1/4π (Subscript[ϵ, 0])=1) involves integrals of the type

K=(∫_(-∞)^∞)(Subscript[χ, a])(r;α,A,a)(-(1/2)(∇^2))(Subscript[χ, b])(r;β,B,b)dr,

in which (Subscript[χ, a])(r) is an unnormalized Cartesian Gaussian primitive centered at the nucleus 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)))=(((x-(Subscript[A, x])))^(Subscript[a, x]))(((e^(-α|x-(Subscript[A, x])((|)^2)))(y-(Subscript[A, y])))^(Subscript[a, y]))(((e^(-α|y-(Subscript[A, y])((|)^2)))(z-(Subscript[A, z])))^(Subscript[a, z]))(e^(-α|z-(Subscript[A, z])((|)^2)))= (Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))·(Subscript[χ, a])(y;α,(Subscript[A, y]),(Subscript[a, y]))·(Subscript[χ, a])(z;α,(Subscript[A, z]),(Subscript[a, z])).

Here α 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,3,… corresponds to functions of type s, p, d, f, …. One notable property of the Gaussian function, which will be used here, is that the derivative of a Gaussian function can be expressed as a sum containing Gaussians of lower and higher Cartesian angular values. In particular,

(∂/(∂x))(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))=(Subscript[a, x])(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])-1)-2α (Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])+1).

Similarly, (∂/(∂(Subscript[A, x])))(Subscript[χ, a])(x;α,(Subscript[A, x,])(Subscript[a, x]))=-(∂/(∂x))(Subscript[χ, a])(x;α,(Subscript[A, x,])(Subscript[a, x])). This property is useful in integrals involving differential operators, such as the kinetic energy; in calculations involving the gradient of the energy or Hamiltonian; and in deriving one of the most important algorithms of Gaussian function computation, the recurrence relation.

The kinetic operator, expressed in Cartesian coordinates as  -(1/2)(∇^2)=-(1/2)(((∂^2)/(∂(x^2)))+((∂^2)/(∂(y^2)))+((∂^2)/(∂(z^2)))), enables separation of the kinetic integrals into three components K=(Subscript[K, x])+(Subscript[K, y])+(Subscript[K, z]). For example,

(Subscript[K, x])=-(1/2)(∫_(-∞)^∞)(Subscript[χ, a])(r;α,A,a)((∂^2)/(∂(x^2)))(Subscript[χ, b])(r;β,B,b)dr.

Furthermore, if we also separate (Subscript[χ, a]) into its Cartesian components, (Subscript[K, x]) becomes a product of three integrals

(Subscript[K, x])=-(1/2)(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))((∂^2)/(∂(x^2)))(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x]))dx(∫_(-∞)^∞)(Subscript[χ, a])(y;α,(Subscript[A, y]),(Subscript[a, y]))·(Subscript[χ, b])(y;β,(Subscript[B, y]),(Subscript[b, y]))dy(∫_(-∞)^∞)(Subscript[χ, a])(z;α,(Subscript[A, z]),(Subscript[a, z]))·(Subscript[χ, b])(z; β,(Subscript[B, z]),(Subscript[b, z]))dz.

Derivative of Gaussian Functions

The last two integrals in (Subscript[K, x]) are simply overlap integrals and can be evaluated as outlined in []. Integrating by parts, the first integral is

-(1/2)[(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))(∂/(∂x))(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x]))(Subscript[|, x=-∞]^(x=∞))-(∫_(-∞)^∞)(∂/(∂x))(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))(∂/(∂x))(Subscript[χ, b])(x,β,(Subscript[B, x]),(Subscript[b, x]))dx]=(1/2)(∫_(-∞)^∞)(∂/(∂x))(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))(∂/(∂x))(Subscript[χ, b])(x,β,(Subscript[B, x]),(Subscript[b, x]))dx

where, according to the property mentioned above, the first term can be expanded into two Gaussian products, which vanish at the integral limits.

Next, we consider the first term of the integrand. Expressing equation () in terms of () and rearranging, we get

(∂/(∂x))(((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-((α(x-(Subscript[A, x])))^2)))=(((Subscript[a, x])(x-(Subscript[A, x])))^((Subscript[a, x])-1))(e^(-((α(x-(Subscript[A, x])))^2)))+[-2α(x-(Subscript[A, x]))](((x-(Subscript[A, x])))^(Subscript[a, x]))(e^(-((α(x-(Subscript[A, x])))^2)))=[(((Subscript[a, x])(x-(Subscript[A, x])))^((Subscript[a, x])-1))−2((α(x−A x))^((Subscript[a, x])+1))](e^(-((α(x-(Subscript[A, x])))^2))).

Substituting the result into equation (), we have

(1/2)(∫_(-∞)^∞)[(Subscript[a, x])(((x-(Subscript[A, x])))^((Subscript[a, x])-1))−2((α(x−A x))^((Subscript[a, x])+1))](e^(-((α(x-(Subscript[A, x])))^2)))[(((Subscript[b, x])(x-(Subscript[B, x])))^((Subscript[b, x])-1))−2((β(x−B x))^((Subscript[b, x])+1))](e^(-((β(x-(Subscript[B, x])))^2)))dx,

and expanding the polynomial part yields

(1/2)(∫_(-∞)^∞)[(Subscript[a, x])(((Subscript[b, x])(x-(Subscript[A, x])))^((Subscript[a, x])-1))(((x-(Subscript[B, x])))^((Subscript[b, x])-1))-2α (((Subscript[b, x])(x−A x))^((Subscript[a, x])+1))(((x-(Subscript[B, x])))^((Subscript[b, x])-1))-2(Subscript[a, x])((β(x-(Subscript[A, x])))^((Subscript[a, x])-1))(((x-(Subscript[B, x])))^((Subscript[b, x])+1))+4α ((β(x-(Subscript[A, x])))^((Subscript[a, x])+1))(((x-(Subscript[B, x])))^((Subscript[b, x])+1))](e^(-((α(x-(Subscript[A, x])))^2)))(e^(-((β(x-(Subscript[B, x])))^2)))dx.

We can now express the first integral of equation () in terms of one-dimensional Gaussian primitives:

(1/2)[(Subscript[a, x])(Subscript[b, x])(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])-1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])-1)dx-2α (Subscript[b, x])(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])+1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])-1)dx-2(Subscript[a, x])β(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])-1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])+1)dx+4α β(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])+1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])+1)dx].

Kinetic-Energy Integral

We have now shown that the kinetic-energy integral can be written as a product of overlap integrals. The x component, for example, is given by

(Subscript[K, x])=(1/2)[(Subscript[a, x])(Subscript[b, x])(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])-1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])-1)dx-2α (Subscript[b, x])(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])+1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])-1)dx-2(Subscript[a, x])β(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])-1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])+1)dx+4α β(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x])+1)·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x])+1)dx](∫_(-∞)^∞)(Subscript[χ, a])(y;α,(Subscript[A, y]),(Subscript[a, y]))·(Subscript[χ, b])(y;β,(Subscript[B, y]),(Subscript[b, y]))dy(∫_(-∞)^∞)(Subscript[χ, a])(z;α,(Subscript[A, z]),(Subscript[a, z]))·(Subscript[χ, b])(z; β,(Subscript[B, z]),(Subscript[b, z]))dz.

Recurrence Relations

Using the notation of [], where the overlap integral of two Gaussians is expressed in terms of its orthogonal components,

S=(∫_(-∞)^∞)(Subscript[χ, a])(x;α,(Subscript[A, x]),(Subscript[a, x]))·(Subscript[χ, b])(x;β,(Subscript[B, x]),(Subscript[b, x]))dx=(((Subscript[E, A,B])[(π/(α+β))])^(3/2))(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])),

the x component of the kinetic-energy integral, equation (), can be written as

(Subscript[K, x])=(1/2)(((Subscript[E, AB])[(π/(α+β))])^(3/2)){(Subscript[a, x])(Subscript[b, x])(Subscript[s, x])((Subscript[a, x])-1,(Subscript[b, x])-1)-2α (Subscript[b, x])(Subscript[s, x])((Subscript[a, x])+1,(Subscript[b, x])-1)-2(Subscript[a, x])(Subscript[βs, x])((Subscript[a, x])-1,(Subscript[b, x])+1)+4α β (Subscript[s, x])((Subscript[a, x])+1,(Subscript[b, x])+1)}(Subscript[s, y])((Subscript[a, y]),(Subscript[b, y]))(Subscript[s, z])((Subscript[a, z]),(Subscript[b, z]))=(((Subscript[E, AB])[(π/(α+β))])^(3/2))(Subscript[k, 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 factor (1/2) has been absorbed into the definition of the kinetic integral (Subscript[k, x])((Subscript[a, x]),(Subscript[b, x])). We will derive the recurrence relation for this term (Subscript[k, x])((Subscript[a, x]),(Subscript[b, x])). For (Subscript[a, x])=(Subscript[b, x])=0, the first three terms of equation () are zero, and we are left with (Subscript[k, x])(0,0)=2α β (Subscript[s, x])(1,1). Analogously, for (Subscript[a, x])=0 or (Subscript[b, x])=0, we have

(Subscript[k, x])((Subscript[a, x]),0)=-(Subscript[a, x]) β (Subscript[s, x])((Subscript[a, x])-1,1)+2α β (Subscript[s, x])((Subscript[a, x])+1,1),
(Subscript[k, x])(0,(Subscript[b, x]))=-α (Subscript[b, x])(Subscript[s, x])(1,(Subscript[b, x])-1)+2α β (Subscript[s, x])(1,(Subscript[b, x])+1),

and the recurrence relation for the kinetic integral function with any two Cartesian angular momenta generalizes to

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

Implementation

The function Kin evaluates the kinetic integral of two Gaussian primitives; here alpha, beta, RA, RB, LA, and LB are α, β, A, B, a, and b as defined earlier. The first step is the evaluation of the overlap integral as described in [].

Kin[alpha_,beta_,RA_,RB_,LA_,LB_]:=Module[{EAB,Kinetic},
Do[
(* Initial Conditions for overlap recurrence relation *)

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];

(* Initial Conditions for kinetic integral function *)

k[i,0,0]=  2*alpha*beta*s[i,1,1];
k[i_,a_,0]:=-a*beta*s[i,a-1,1]+2*alpha*beta*s[i,a+1,1];
k[i_,0,b_]:=-alpha*b*s[i,1,b-1]+2*alpha*beta*s[i,1,b+1];

(* Kinetic Energy Integral *)

k[i_,a_,b_]:=(a*b*s[i,a-1,b-1]
             -2*a*beta*s[i,a-1,b+1]
             -2*alpha*b*s[i,a+1,b-1]
             +4*alpha*beta*s[i,a+1,b+1])/2;

,{i,1,3}]; (* Loop over three Cartesian components *)

(* Exponential Overlap Factor *)

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

(* Kinetic Energy Integral *)

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

We describe in detail the evaluation of the kinetic-energy matrix for the water molecule ((Subscript[r, OH])=1.86942bohr, (Subscript[∡, HOH])=100.0269°, the 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 approximated by a sum of three Gaussians; their unnormalized primitive contraction coefficients and orbital exponents (taken from []) are as follows.

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 basis function origins and Cartesian angular values of the orbitals, listed in the order (Subscript[H, 1])(1s), (Subscript[H, 2])(1s), (Subscript[O, 3])(1s), (Subscript[O, 3])(2s), (Subscript[O, 3])(2(Subscript[p, x])), (Subscript[O, 3])(2(Subscript[p, y])), and (Subscript[O, 3])(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}};

Specifically, for the kinetic-energy integral of the first primitive of the 1s orbital of hydrogen atom 1, (Subscript[χ, 11])(r), and the first primitive of the 2(Subscript[p, z]) orbital of the oxygen atom, (Subscript[χ, 71])(r),

(∫_(-∞)^∞)(Subscript[χ, 11])(r)(-(1/2)(∇^2))(Subscript[χ, 71])(r)dr=(∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞)(e^(-3.425250914 |r -(Subscript[R, H1])((|)^2)))(-(1/2)(∇^2))(z-0.24026010)×(e^(-5.033151319 |r -(Subscript[R, O3])((|)^2)))dx dy dz,

where the indices (Subscript[χ, i,j]) indicate primitive j of basis function i. The kinetic integral in terms of the kinetic and overlap function is then

(((Subscript[E, AB])[(π/(α+β))])^(3/2)){(Subscript[k, x])(0,0)(Subscript[s, y])(0,0)(Subscript[s, z])(0,1)+(Subscript[s, x])(0,0)(Subscript[k, y])(0,0)(Subscript[s, z])(0,1)+(Subscript[s, x])(0,0)(Subscript[s, y])(0,0)(Subscript[k, z])(0,1)}.

We start with a={0,0,0} and b={0,0,1}, and from the first four equations of the module, we use the recurrence scheme to build up the overlap values needed for (Subscript[k, z])(0,1). We will need values up to (Subscript[s, z])(1,2),

(Subscript[s, z])(0,0)=1,
(Subscript[s, z])(1,0)=-((Subscript[H1, z])-((α (Subscript[H1, z])+β (Subscript[O, z]))/(α+β)))=0.7148308840373454`,
(Subscript[s, z])(2,0)=-((Subscript[H1, z])-((α (Subscript[H1, z])+β (Subscript[O, z]))/(α+β)))(Subscript[s, z])(1,0)+((2-1)/(2(α+β)))(Subscript[s, z])(0,0)=0.5700960946402821`,
(Subscript[s, z])(3,0)=-((Subscript[H1, z])-((α (Subscript[H1, z])+β (Subscript[O, z]))/(α+β)))(Subscript[s, z])(2,0)+((3-1)/(2(α+β)))(Subscript[s, z])(1,0)=0.4920335786647066`,
(Subscript[s, z])(1,1)=(Subscript[s, z])(2,0)+((Subscript[H1, z])-(Subscript[O, z]))(Subscript[s, z])(1,0)=-0.2886296292666999`,
(Subscript[s, z])(2,1)=(Subscript[s, z])(3,0)+((Subscript[H1, z])-(Subscript[O, z]))(Subscript[s, z])(2,0)=-0.19282260414703994`,
(Subscript[s, z])(1,2)=(Subscript[s, z])(2,1)+((Subscript[H1, z])-(Subscript[O, z]))(Subscript[s, z])(1,1)=0.15390836042870001`,

leading to (Subscript[k, z])(0,1)=-α·1·(Subscript[s, z])(1,0)+2 α·β·(Subscript[s, z])(1,2)=2.85821.

Similarly, for (Subscript[k, x])(0,0), we will need (Subscript[s, x])(1,1) which, in turn, comes from (Subscript[s, x])(1,0) and (Subscript[s, x])(2,0). We have

(Subscript[s, x])(1,0)=0,

(Subscript[s, x])(2,0)=0.05911281897298251`,

(Subscript[s, x])(1,1)=0.05911281897298251`,

leading to (Subscript[k, x])(0,0)=2α β (Subscript[s, x])(1,1)=2.038187080823005`.

For the y component,

(Subscript[s, y])(1,0)=-0.8523084269657298`,

(Subscript[s, y])(2,0)=0.7855424736497794`,

(Subscript[s, y])(1,1)=-0.43525019158175793`,

leading to (Subscript[k, y])(0,0)=15.0073.

With (Subscript[E, AB])=0.0008064444210737256`, the kinetic integral of equation () is K=0.0016734344532709955`, which we can also obtain from the module Kin.

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

Three Gaussian primitives (Subscript[χ, a])(r) for each atomic orbital (Subscript[ϕ, i])(r) result in nine integrals of the type that we have just evaluated. For example, the (Subscript[K, 17]) element of the kinetic-energy matrix,

(Subscript[K, 17])=(∑_(i=1)^3)(∑_(j=1)^3)(Subscript[N, 1i])(Subscript[N, 7j])(Subscript[d, 1i])(Subscript[d, 7j])(∫_(-∞)^∞)(Subscript[χ, 1i])(r)(-(1/2)(∇^2))(Subscript[χ, 7j])(r)dr,

is (∫_(-∞)^∞)(Subscript[ϕ, 1])(r)(-(1/2)(∇^2))(Subscript[ϕ, 7])(r)dr, derived via a contraction scheme that requires the following normalization factor.

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

This calculates (Subscript[K, 17]).

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

For larger basis sets, one needs only to replace the summation upper limit 3 in equation () with the appropriate number of primitives belonging to a particular basis function and an additional summation for the basis functions.

Finally, here is the resulting kinetic energy matrix.

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

Since the kinetic-energy matrix is symmetrical, we need only to calculate the upper elements. The value of the elements (Subscript[K, 11])=(Subscript[K, 22])=0.760032 hartree is the electronic kinetic energy of the hydrogen atom described by the STO-3G basis set (compared to the exact value of 0.5 hartree). Analogously, the element (Subscript[K, 33]) is the corresponding kinetic energy of an electron in the 1s orbital of the oxygen atom, and so forth.

Conclusion

We have provided an introduction to the evaluation of kinetic-energy integrals involving Gaussian-type basis functions both analytically and by use of recurrence relations. The results are sufficiently general so that no modification of the algorithm is needed when larger basis sets with more Gaussian primitives or primitives with larger angular momenta are employed.

References

M. Hô and J. M. Hernández-Pérez, “Evaluation of Gaussian Molecular Integrals I,” The Mathematica Journal, 2012. doi:10.3888/tmj.14-3.

“Basis Set Exchange.” (April 12, 2012) bse.pnl.gov/bse/portal.