Kinetic Energy Integral Obara-Saika Method

Kinetic Energy Integral

In Hartree Fock equation and Kohn-Sham equation, kinetic energy term is written as

Kpq=−12 ∫−∞∞χp(𝐫−𝐀)∇2χq(𝐫−𝐁)d𝐫(1)\begin{aligned} K_{pq}=- \frac 12 \ \int_{-\infin}^{\infin} \chi_p(\mathbf {r-A}) \nabla^2 \chi_q(\mathbf {r-B})d\mathbf r \tag{1} \end{aligned}

In this page, I explain how to compute kinetic energy integral analytically using Obara-Saika Method.

Kinetic Energy Integral with gaussian basis

We use gaussian basis sets as atomic orbitals.

Then, kinetic energy integral is descraibed as (without normalization coefficients)

− 12∫−∞∞(x−Ax)n(y−Ay)l(z−Az)mexp[−α|𝐫−𝐀|2)]\begin{aligned} – \ \frac 12 \int_{-\infin}^{\infin} & (x-A_x)^n(y-A_y)^l(z-A_z)^m exp[-\alpha_\,|\mathbf {r-A}|^2 )] \end{aligned}
×∇2(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]dxdydz(2)\times \nabla^2(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]dxdydz \tag{2}

You can calculate Eq.(2) analitically. But it’s quite hard to compute all combinations of orbitals including higer angular momontum orbitals like d and f orbitals.

So you can use the great method called Obara-Saika method.

Obara-Saika method of Kinetic Energy Integral

General formula is shown below

T(𝐚,𝐛)=β∑i=02(2bi+1) S(𝐚,𝐛)−12∑i=02bi(bi−1)S(𝐚,𝐛i−2)−2β2∑i=02S(𝐚,𝐛i+2)(3)T(\mathbf a,\mathbf b) = \beta\sum_{i=0} ^{2} (2b_i+1) \ S(\mathbf a,\mathbf b) – \frac 12 \sum_{i=0} ^{2} b_i(b_i-1)S(\mathbf a,\mathbf b_i-2)- 2\beta^2 \sum_{i=0} ^{2} S(\mathbf a,\mathbf b_i+2) \tag {3}

Where

𝐚=(n,l,m)\mathbf a =(n, l, m)
𝐛=(n′,l′,m′)\mathbf b =(n’, l’, m’)

Kinetic Energy Integral is described using overlap integral S(a,b).

So it’s easy to code kinetic energy.

But my curiosity is “how can I get equation (3) ?”.

Let’s dive into obara-saika formula.

How to get Obara-Saika fomula of Kinetic Energy

First, we look at Eq.(2) again.

− 12∫−∞∞(x−Ax)n(y−Ay)l(z−Az)mexp[−α|𝐫−𝐀|2)]\begin{aligned} – \ \frac 12 \int_{-\infin}^{\infin} & (x-A_x)^n(y-A_y)^l(z-A_z)^m exp[-\alpha_\,|\mathbf {r-A}|^2 )] \end{aligned}
×∇2(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]dxdydz(2)\times \nabla^2(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]dxdydz \tag {2}

Derivation is taken to second row part of Eq.(2).

So we calulate second row part

∇2(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)](3)\nabla^2(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)] \tag {3}

To clarify, we rewrite nabula to derivation forms.

(3)=(d2dx2+d2dy2+d2dz2)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)](4)(3) = (\frac {d^2}{dx^2}+\frac {d^2}{dy^2}+\frac {d^2}{dz^2})(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)] \tag {4}

We need to calculate second order derivatives. So let’s dive into them.

1. Fisrt Order Derivative

To simplify, we calculate derivative with 1D equation.

Take derivative the atomic orbital on B center with respect to x, we get

ddx[ (x−Bx)n′exp{−β(x−Bx)2} ]\frac{d} {dx}[ \ (x-B_x)^{n’} exp\{- \beta(x-B_x)^2\} \ ]
=n′(x−Bx)n′−1exp{−β(x−Bx)2}=n'(x-B_x)^{n’-1} exp\{- \beta(x-B_x)^2\}
−2β(x−Bx)n′+1exp{−β(x−Bx)2} (5)-2\beta(x-B_x)^{n’+1} exp\{- \beta(x-B_x)^2\} \ \tag {5}

2. Second Order Derivative

Take derivative Eq.(4) with respect to x, we get

d2dx2[ (x−Bx)n′exp{−β(x−Bx)2} ]\frac{d^2} {dx^2}[ \ (x-B_x)^{n’} exp\{- \beta(x-B_x)^2\} \ ]
=n′(n′−1)(x−Bx)n′−2exp{−β(x−Bx)2}=n'(n’-1)(x-B_x)^{n’-2} exp\{- \beta(x-B_x)^2\}
−2βn′(x−Bx)n′exp{−β(x−Bx)2} -2\beta n'(x-B_x)^{n’} exp\{- \beta(x-B_x)^2\} \
−2β(n′+1)(x−Bx)n′exp{−β(x−Bx)2} -2\beta(n’+1) (x-B_x)^{n’} exp\{- \beta(x-B_x)^2\} \
+4β2(x−Bx)n′+2exp{−β(x−Bx)2} +4\beta^2(x-B_x)^{n’+2} exp\{- \beta(x-B_x)^2\} \

Further more

=n′(n′−1)(x−Bx)n′−2exp{−β(x−Bx)2}=n'(n’-1)(x-B_x)^{n’-2} exp\{- \beta(x-B_x)^2\}
−2β(2n′+1)(x−Bx)n′exp{−β(x−Bx)2} -2\beta (2n’+1)(x-B_x)^{n’} exp\{- \beta(x-B_x)^2\} \
+4β2(x−Bx)n′+2exp{−β(x−Bx)2} (6)+4\beta^2(x-B_x)^{n’+2} exp\{- \beta(x-B_x)^2\} \ \tag {6}

We get back to 3D equation by adding y and z components

d2dx2[(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]]\frac{d^2} {dx^2}[(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]]
=n′(n′−1)(x−Bx)n′−2(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]=n'(n’-1)(x-B_x)^{n’-2} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
−2β(2n′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2n’+1)(x-B_x)^{n’} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+4β2(x−Bx)n′+2(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)](7)+4\beta^2(x-B_x)^{n’+2} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)] \tag {7}

In the same way of second order derivative with respect to x, we can get y and z version

d2dy2[(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]]\frac{d^2} {dy^2}[(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]]
=l′(l′−1)(x−Bx)n′(y−By)l′−2(z−Bz)m′exp[−β|𝐫−𝐁|2)]=l'(l’-1)(x-B_x)^{n’}(y-B_y)^{l’-2}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
−2β(2l′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2l’+1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+4β2(x−Bx)n′(y−By)l′+2(z−Bz)m′exp[−β|𝐫−𝐁|2)](7)+4\beta^2(x-B_x)^{n’}(y-B_y)^{l’+2} (z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)] \tag {7}

d2dz2[(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]]\frac{d^2} {dz^2}[(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]]
=m′(m′−1)(x−Bx)n′(y−By)l′(z−Bz)m′−2exp[−β|𝐫−𝐁|2)]=m'(m’-1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’-2} exp[- \beta|\mathbf {r-B}|^2)]
−2β(2m′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2m’+1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+4β2(x−Bx)n′(y−By)l′(z−Bz)m′+2exp[−β|𝐫−𝐁|2)](8)+4\beta^2(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’+2} exp[- \beta|\mathbf {r-B}|^2)] \tag {8}

3. Insert Eq.(6), (7), (8) into Eq.(4)

Before Inserting, we consider sum up of Eq.(6) ~ (8) one term by one term.

(d2dx2+d2dy2+d2dz2)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)] (\frac {d^2}{dx^2}+\frac {d^2}{dy^2}+\frac {d^2}{dz^2})(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]

First, look at second terms of right hand in Eq.(6) ~ (8).

We notice that aungular momentum exponents and exponential terms are the same.

So second terms of right hand can be described as (using ket form)

−2β(2n′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2n’+1)(x-B_x)^{n’} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
−2β(2l′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2l’+1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
−2β(2m′+1)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]-2\beta (2m’+1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
=−2β∑i=02(2bi+1) |𝐛)(9)=-2\beta\sum_{i=0} ^{2} (2b_i +1) \ |\mathbf b) \tag {9}

Next, take a look at first term, x, y, z component of angular momentum exponents are different.

n′(n′−1)(x−Bx)n′−2(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]n'(n’-1)(x-B_x)^{n’-2} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+l′(l′−1)(x−Bx)n′(y−By)l′−2(z−Bz)m′exp[−β|𝐫−𝐁|2)]+l'(l’-1)(x-B_x)^{n’}(y-B_y)^{l’-2}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+m′(m′−1)(x−Bx)n′(y−By)l′(z−Bz)m′−2exp[−β|𝐫−𝐁|2)]+m'(m’-1)(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’-2} exp[- \beta|\mathbf {r-B}|^2)]
=∑i=02bi(bi−1) |𝐛−𝟐𝐢)(10)=\sum_{i=0}^2 b_i(b_i-1) \ |\mathbf {b-2_i}) \tag {10}

Then, third terms of right hand can be described as

4β2(x−Bx)n′+2(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)]4\beta^2(x-B_x)^{n’+2} (y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
+4β2(x−Bx)n′(y−By)l′+2(z−Bz)m′exp[−β|𝐫−𝐁|2)]+4\beta^2(x-B_x)^{n’}(y-B_y)^{l’+2} (z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
4β2(x−Bx)n′(y−By)l′(z−Bz)m′+2exp[−β|𝐫−𝐁|2)]4\beta^2(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’+2} exp[- \beta|\mathbf {r-B}|^2)]
=4β2∑i=02 |𝐛+𝟐𝐢)(11)=4\beta^2\sum_{i=0}^2 \ |\mathbf {b+2_i}) \tag {11}

Them, we get

(d2dx2+d2dy2+d2dz2)(x−Bx)n′(y−By)l′(z−Bz)m′exp[−β|𝐫−𝐁|2)] (\frac {d^2}{dx^2}+\frac {d^2}{dy^2}+\frac {d^2}{dz^2})(x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} exp[- \beta|\mathbf {r-B}|^2)]
=−2β∑i=02(2bi+1) |𝐛)+∑i=02bi(bi−1) |𝐛−𝟐𝐢)+4β2∑i=02 |𝐛+𝟐𝐢)(12)=-2\beta\sum_{i=0} ^{2} (2b_i +1) \ |\mathbf b) + \sum_{i=0}^2 b_i(b_i-1) \ |\mathbf {b-2_i}) +4\beta^2\sum_{i=0}^2 \ |\mathbf {b+2_i}) \tag {12}

4. Insert Eq.(12) into Eq.(2)

Now we are ready to get Obara-Saika formula of kinetic energy integral

Orbital function on center A with angular momentum a is multiplied to left side.

Using bra form, we get

−12[−2β∑i=02(2bi+1) |𝐛)+∑i=02bi(bi−1) (𝐚|𝐛−𝟐𝐢)+4β2∑i=02 (𝐚|𝐛+𝟐𝐢)]-\frac 12[-2\beta\sum_{i=0} ^{2} (2b_i +1) \ |\mathbf b) + \sum_{i=0}^2 b_i(b_i-1) \ (\mathbf a|\mathbf {b-2_i}) + 4\beta^2\sum_{i=0}^2 \ (\mathbf a|\mathbf {b+2_i})]
=β∑i=02(2bi+1)  (𝐚|𝐛)−12∑i=02bi(bi−1) (𝐚|𝐛−𝟐𝐢)−2β2∑i=02 (𝐚|𝐛+𝟐𝐢)=\beta\sum_{i=0} ^{2} (2b_i +1) \ \ (\mathbf a|\mathbf b) -\frac 12 \sum_{i=0}^2 b_i(b_i-1) \ (\mathbf a|\mathbf {b-2_i}) – 2\beta^2\sum_{i=0}^2 \ (\mathbf a|\mathbf {b+2_i})

Recall the relatinship

(𝐚|𝐛)=S(𝐚|𝐛)(\mathbf a|\mathbf {b})=S(\mathbf a|\mathbf {b})

This is exactly overlap integral.

So using overlap ingegral form, we finally get

T(𝐚,𝐛)=β∑i=02(2bi+1)  S(𝐚|𝐛)−12∑i=02bi(bi−1) S(𝐚|𝐛−𝟐𝐢)−2β2∑i=02 S(𝐚|𝐛+𝟐𝐢)]T(\mathbf a,\mathbf b)=\beta\sum_{i=0} ^{2} (2b_i +1) \ \ S (\mathbf a|\mathbf b) -\frac 12 \sum_{i=0}^2 b_i(b_i-1) \ S(\mathbf a|\mathbf {b-2_i}) – 2\beta^2\sum_{i=0}^2 \ S(\mathbf a|\mathbf {b+2_i})]

コメント

タイトルとURLをコピーしました