Overlap Integral – Obara-Saika Method –

Overlap integrals are one of the core parts of quantum computing chemistry such as Hartree Fock Method and DFT.

The general form of overlap integral is

∫χp(𝐫−𝐀)χq(𝐫−𝐁)d𝐫\int \chi_p(\mathbf {r -A} )\chi_q(\mathbf {r -B} ) d\mathbf r

where

A and B are the centers of atoms.

The orbital function is

χ(𝐫)=Nd(x−Ax)n(y−Ay)l(z−Az)mexp[−α|𝐫−𝐀|2]\chi(\mathbf r) = N\,d(x-A_x)^n(y-A_y)^l(z-A_z)^m\,exp[-\alpha\,|\mathbf {r-A}|^2]

So, overlap integral is described as

NpNqdpdq∫−∞∞(x−Ax)n(y−Ay)l(z−Az)m(x−Bx)n′(y−By)l′(z−Bz)m′exp[−(α|𝐫−𝐀|2+β|𝐫−𝐁|2)]dxdydz(1)\begin{aligned} N_p N_q\,d_p d_q\int_{-\infin}^{\infin} & (x-A_x)^n(y-A_y)^l(z-A_z)^m (x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’} \\ & exp[-(\alpha_\,|\mathbf {r-A}|^2 + \beta|\mathbf {r-B}|^2)]dxdydz \tag{1} \end{aligned}

Overlap Integral of s-orbitals

Overlap integral of s-orbitals can be very simple with (n, l, m) = (n‘, l‘, m‘) = (0, 0, 0)

(s|s)=NpNqdpdq∫−∞∞exp[−(α|𝐫−𝐀|2+β|𝐫−𝐁|2)]dxdydz(s|s)= N_p N_q\,d_p d_q\int_{-\infin}^{\infin} exp[-(\alpha_\,|\mathbf {r-A}|^2 + \beta|\mathbf {r-B}|^2)]dxdydz

This is easy to calculate the integral analyticaly

(s|s)=NpNqdpdq(πα+β)32exp[−αβα+β|𝐑AB|2](2)(s|s)=N_p N_q\,d_p d_q (\frac \pi {\alpha+\beta})^{\frac 32} exp[-\frac {\alpha \beta}{\alpha+\beta}|\mathbf R_{AB}|^2] \tag{2}

Where

𝐑𝐀𝐁=𝐀−𝐁\mathbf {R_{AB}} = \mathbf {A – B}

But, general cases like s-orbital and p-orbital, d-orbital and f-orbital etc. are not easy to calculate.

In addition, it is very bothering to code all combinations of overlap integrals. ( Neuclear attraction and electron repulsive integrals are much more boring!)

So many general solutions have been developed.

I selected obara-saika method among of great methods to code Hartree Fock method because it is clear to undarstand its algorism.

Obara-Saika Method for Overlap Integrals

There are two major Obara-Saika formulas.

Vartival Recurrence Relation (VRR)

(𝐚+𝟏i|𝐛)=(Pi−Ai)(𝐚|𝐛)+ai2p(𝐚−𝟏𝐢|𝐛)+bi2p(𝐚|𝐛−𝟏i)(3)(\mathbf {a+1}_i|\mathbf b) = ({P_i-A_i})(\mathbf {a}|\mathbf b)+\frac {a_i} {2p}(\mathbf {a-1_i}|\mathbf b)+\frac {b_i} {2p}(\mathbf {a}|\mathbf {b-1}_i) \tag {3}

or specifically

S(𝐚+𝟏i|𝐛)=(Pi−Ai) S(𝐚|𝐛)+ai2p S(𝐚−𝟏𝐢|𝐛)+bi2p S(𝐚|𝐛−𝟏i)(3’)S(\mathbf {a+1}_i|\mathbf b) = ({P_i-A_i}) \ S(\mathbf {a}|\mathbf b)+\frac {a_i} {2p} \ S(\mathbf {a-1_i}|\mathbf b)+\frac {b_i} {2p} \ S(\mathbf {a}|\mathbf {b-1}_i) \tag {3′}

Horizontal Recurrence Relation (HRR)

(𝐚|𝐛+𝟏i)=(𝐚+𝟏i|𝐛)+(Ai−Bi)(𝐚|𝐛)(4)(\mathbf {a}|\mathbf {b+1}_i) = (\mathbf {a+1}_i|\mathbf b)+( {A_i-B_i})(\mathbf {a}|\mathbf b) \tag {4}
S(𝐚|𝐛+𝟏i)=S(𝐚+𝟏i|𝐛)+(Ai−Bi) S(𝐚|𝐛)(4’)S(\mathbf {a}|\mathbf {b+1}_i) = S(\mathbf {a+1}_i|\mathbf b)+( {A_i-B_i}) \ S(\mathbf {a}|\mathbf b) \tag {4′}

I selected VRR formula Eq.(3) for coding HF method because it is easier to code than HRR.

I explain the format of those fomulas.

a, b

They represent angular momentum number of each atom center like (n, l, m) and (n‘, l‘, m‘).

For example, Overlap integral is written as (without normalization coefficients)

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

This is the integral of mulitiplication of two orbitals (again without normalization coefficients)

χ1(𝐫)=(x−Ax)n(y−Ay)l(z−Az)mexp[−α|𝐫−𝐀|2]\chi_1(\mathbf r) = (x-A_x)^n(y-A_y)^l(z-A_z)^m\,exp[-\alpha\,|\mathbf {r-A}|^2]
χ2(𝐫)=(x−Bx)n′(y−By)l′(z−Bz)m′exp[−α|𝐫−𝐁|2]\chi_2(\mathbf r) = (x-B_x)^{n’}(y-B_y)^{l’}(z-B_z)^{m’}\,exp[-\alpha\,|\mathbf {r-B}|^2]

So, we reduce the description of the overlap integral as

( (n,l,m) | (n′,l′,m′) )(\space(n, l, m)\space|\space(n’, l’,m’)\space)

Simpler version is

( 𝐚 | 𝐛 )(\space\mathbf{a}\space|\space \mathbf{b}\space)

So what is subscription i ?

Subscription i means the index of x, y, z elements.

For example, x, y, z are usually assigned as

i=0:xi = 0:\\x
i=1:yi = 1:\\y
i=2:zi = 2:\\z

Let’s see i = 0 case of VRR.

((n+1,l,m)|(n′,l′,m′))=(\\(n+1,l,m)\\|\\(n’,l’,m’)\\) =
(Px−Ax)((n,l,m)|(n′,l′,m′))+n2p((n−1,l,m)|(n′,l′,m′))+n′2p((n,l,m)|(n′−1,l′,m′))({P_x-A_x})(\\(n,l,m)\\|\\(n’,l’,m’)\\)+\frac {n} {2p}(\\(n-1,l,m)\\|\\(n’,l’,m’)\\)+\frac {n’} {2p}(\\(n,l,m)\\|\\(n’-1,l’,m’)\\)

P

Here, new variant appeared. P is a single gaussian cemter brought by gaussian product theorem.

P=α𝐀+β𝐁α+βP=\frac{\alpha\mathbf A+\beta\mathbf B}{\alpha+\beta}

How to get VRR formula

I explain how to get VRR formula.

To simplify, I consider 1-D overlap integral. But the core procedure is not changed because x, y, z elements of exponential are independent.

Here, 1-D overlap integral

(n,n′)=∫−∞∞(x−Ax)n(x−Bx)n′exp[−(α(x−Ax)2+β(x−Bx)2)]dx(6)(n,n’)= \int_{-\infin}^{\infin} (x-A_x)^n(x-B_x)^{n’}exp[-(\alpha(x-A_x)^2 + \beta(x-B_x)^2)]dx \tag {6}

1. Transform Eq.(6) by gaussian product theorem

By gaussian product therem, we get

(n,n′)=exp[αβα+β(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′exp[−(α+β)(x−αAx+βBxα+β)2]dx(7)(n,n’)= exp[\frac {\alpha \beta}{\alpha +\beta}(A_x-B_x)^2]\int_{-\infin}^{\infin} (x-A_x)^n(x-B_x)^{n’}exp[-(\alpha+\beta)(x-\frac {\alpha A_x + \beta B_x}{\alpha +\beta})^2]dx \tag {7}
Let μ=αβα+β ,    α+β=p,    Px=αAx+βBxα+β\begin{aligned} Let &\space\mu = \frac {\alpha \beta}{\alpha +\beta} \ ,\ \ \ \ \alpha+\beta=p ,\ \ \ \ P_x= \frac {\alpha A_x + \beta B_x}{\alpha +\beta} \end{aligned}

we get a simpler form

(n,n′)=exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′exp[−p(x−Px)2]dx(8)\begin{equation} (n,n’)= exp[ \ \mu(A_x-B_x)^2]\int_{-\infin}^{\infin} (x-A_x)^n(x-B_x)^{n’}exp[-p(x-P_x)^2]dx \tag{8} \end{equation}

2. Transform Eq.(9) by integral by parts

First, we consider the case of (n+1,n’)

(n+1,n′)=exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n+1(x−Bx)n′exp[−p(x−Px)2)]dx(14)(n+1,n’) = exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n+1}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx \tag {14}

Then, change the form of (x-Ax) term into

(n+1,n′)=exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)(x−Ax)n(x−Bx)n′exp[−p(x−Px)2)]dx(15)(n+1,n’) = exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)(x-A_x)^{n}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx \tag {15}

Here, let

(x−Ax)=(x−Px)+(Px−Ax)(x-A_x) = (x-P_x)+(P_x-A_x)

we get

(n+1,n′)=exp[ μ(Ax−Bx)2]∫−∞∞{(x−Px)+(Px−Ax)}(x−Ax)n(x−Bx)n′exp[−p(x−Px)2)]dx(n+1,n’) =exp[ \ \mu(A_x-B_x)^2] \int_{-\infin} ^{\infin} \{(x-P_x)+(P_x-A_x)\}(x-A_x)^{n}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx

Then

(n+1,n′)=exp[ μ(Ax−Bx)2]∫−∞∞(x−Px)(x−Ax)n(x−Bx)n′exp[−p(x−Px)2)]dx(n+1,n’) = exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-P_x)(x-A_x)^{n}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx
+(P−Ax)exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′exp[−p(x−Px)2)]dx(16)+ (P-A_x)exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx \tag {16}

Take a look at the frist term of right hand inEq.(16), we notice that

(x−Px)exp[−p(x−Px)2)]=− 12p ddx  exp{−p(x−Px)2)} (17)\begin{aligned} (x-P_x)exp[-p(x-P_x)^2)] = – \ \frac {1}{2p} \ \frac {d}{dx} \ \ exp\{-p(x-P_x)^2)\} \ \end{aligned} \tag {17}

So we can rewrite Eq.(16) using Eq.(17) as

(n+1,n′)=− 12p exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′ddxexp[−p(x−Px)2)]dx(n+1,n’) = – \ \frac {1}{2p} \ exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} \frac {d}{dx} exp[-p(x-P_x)^2)] dx
+(P−Ax)exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′  exp{−p(x−Px)2)}dx(18)+ (P-A_x)exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} \ \ exp \{-p(x-P_x)^2) \} dx \tag {18}

Now, let’s recall integration by parts formula

ddx[f(x)g(x)]=ddxf(x)×g(x)+f(x)×ddxg(x)\frac {d}{dx} [ f(x)g(x) ] = \frac {d}{dx} f(x) \cross g(x) + f(x) \cross \frac {d}{dx} g(x)
f(x)×ddxg(x)=ddx[f(x)g(x)]−ddxf(x)×g(x) f(x) \cross \frac {d}{dx} g(x)=\frac {d}{dx} [ f(x)g(x) ] -\frac {d}{dx} f(x) \cross g(x)

then, we get the fomula of integration by parts

∫f(x)×ddxg(x)dx=|f(x)g(x)|−∫ddxf(x)×g(x)dx\int f(x) \cross \frac {d}{dx} g(x) dx= | f(x)g(x) | – \int \frac {d}{dx} f(x) \cross g(x) dx

Applying integration by parts to the second term of right hand in Eq.(18), we get

(n+1,n′)= 12p exp[ μ(Ax−Bx)2] |(x−Ax)n(x−Bx)n′]exp[−p(x−Px)2]|−∞∞(n+1,n’) = \ \frac {1}{2p} \ exp[ \ \mu(A_x-B_x)^2] \ | (x-A_x)^{n}(x-B_x)^{n’}] exp[-p(x-P_x)^2]|_{-\infin} ^{\infin}
+ 12p exp[ μ(Ax−Bx)2]∫−∞∞ddx[(x−Ax)n(x−Bx)n′]exp[−p(x−Px)2)]dx+ \ \frac {1}{2p} \ exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} \frac {d}{dx} [ (x-A_x)^{n}(x-B_x)^{n’}] exp[-p(x-P_x)^2)] dx
+(Px−Ax)exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′  exp{−p(x−Px)2)}dx+ (P_x-A_x)exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} \ \ exp \{-p(x-P_x)^2) \} dx

Because

limx→±∞⁡ exp[−p(x−Px)2]=0\lim_{x \to \pm \infin} \ exp[-p(x-P_x)^2] =0

the first term of the right hand become zero.

So, we get

(n+1,n′)= 12p exp[ μ(Ax−Bx)2]∫−∞∞ddx[(x−Ax)n(x−Bx)n′]exp[−p(x−Px)2)]dx(n+1,n’) = \ \frac {1}{2p} \ exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} \frac {d}{dx} [ (x-A_x)^{n}(x-B_x)^{n’}] exp[-p(x-P_x)^2)] dx
+(Px−Ax)exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′  exp{−p(x−Px)2)}dx+ (P_x-A_x)exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} \ \ exp \{-p(x-P_x)^2) \} dx

Then, we continue calculation

(n+1,n′)= n2pexp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n−1(x−Bx)n′  exp{−p(x−Px)2)}dx(n+1,n’) = \ \frac {n}{2p}exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n-1}(x-B_x)^{n’} \ \ exp \{-p(x-P_x)^2) \} dx
+ n′2pexp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′−1  exp{−p(x−Px)2)}dx + \ \frac {n’}{2p}exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’-1} \ \ exp \{-p(x-P_x)^2) \} dx
+(Px−Ax)exp[ μ(Ax−Bx)2]∫−∞∞(x−Ax)n(x−Bx)n′exp[−p(x−Px)2)]dx+(P_x-A_x)exp[ \ \mu(A_x-B_x)^2]\int_{-\infin} ^{\infin} (x-A_x)^{n}(x-B_x)^{n’} exp[-p(x-P_x)^2)] dx

Finally, we get

(n+1,n′)=(Px−Ax)(n,n′)+n2p(n−1,n′)+ n′2p(n,n′−1)(19)(n+1,n’) = (P_x-A_x)(n,n’) + \frac {n}{2p} (n-1,n’)+ \ \frac {n’}{2p} (n,n’-1) \tag {19}

y and z components of VRR

If you start with y or z components, you get with the same procedure as that of x components

(l+1,l′)=(Py−Ay)(l,l′)+ l2p(l−1,l′)+ l′2p(l,l′−1)(20)(l+1,l’) =(P_y-A_y)(l,l’)+ \ \frac {l}{2p} (l-1,l’)+ \ \frac {l’}{2p} (l,l’-1) \tag {20}
(m+1,m′)=(Pz−Az)(m,m′)+ m2p(m−1,m′)+ m′2p(m,m′−1)(21)(m+1,m’) =(P_z-A_z)(m,m’)+ \ \frac {m}{2p} (m-1,m’)+ \ \frac {m’}{2p} (m,m’-1) \tag {21}

Reference

コメント

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