4. Mcmurchie Davidson Scheme
4.1. Overlap
Let us frist consider the simple overlap integral \begin{equation} S_{ab} = \langle G_a|G_b \rangle \end{equation}
it can be factorized in the three Cartesian directons \begin{equation} S_{ab} = S_{ij}S_{kl}S_{mn} \end{equation}
The Obara-Saika recurrence relations for the Cartesian overlap integral over one direction is:
\begin{equation} \begin{aligned} S_{i+1,j} &= X_{PA}S_{ij} + \frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})\\ S_{i,j+1} &= X_{PB}S_{ij} + \frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})\\ \end{aligned} \end{equation}
With the boundary condtion
\begin{equation} S_{0,0} = \sqrt(\frac{\pi}{p}) exp(-\mu X_{AB}^2) \end{equation}
[1]:
class Overlap(object):
"""The Obara-Saika scheme for three-dimensional overlap integral over
primitive Gaussian orbitals.
Parameters
----------
a : float
Gaussian exponent facotr.
b : float
Gaussian exponent facotr.
i : int
Angular momentum quantum number.
j : int
Angular momentum quantum number.
A : float
Coordinate in on direction.
B : float
Coordinate in on direction.
Returns
-------
result : float
The non-normalizecd overlap interal.
"""
def __init__(self):
self.p = 0
self.mu = 0
self.P = ()
def __call__(self, pga, pgb):
"""Evaluates non-normalizecd overlap integral over two primitive gaussian orbitals.
Parameters
----------
pga: PrimitiveGaussian
The first primitive gaussian orbital.
pgb: PrimitiveGaussian
The second primitive gaussian orbital.
C: List[float,float,float]
Coordinate of nuclei.
Return
------
result : float
Integral value.
"""
result = 1
for r in range(3):
result *= self.S1d(r, pga, pgb)
return result
def S1d(self, r, pga, pgb):
a = pga.exponent
b = pgb.exponent
p = a + b
mu = (a*b)/(a+b)
A = np.array(pga.origin)
B = np.array(pgb.origin)
XAB = A-B
if pga.shell[r] > 0:
return self.recursive(r, *self.gaussian_factory(r, pga, pgb))
elif pgb.shell[r] > 0:
return self.recursive(r, *self.gaussian_factory(r, pgb, pga))
else:
# Starting from the spherical Gaussians.
S00 = np.power(np.pi/p,0.5)*np.exp(-mu*XAB[r]**2)
return S00
def recursive(self, r, pga, pgb, pga_1, pga_2, pgb_1):
term1 = term2 = term3 = 0
a = pga.exponent
b = pgb.exponent
p = a + b
mu = (a*b)/(a+b)
A = np.array(pga.origin)
B = np.array(pgb.origin)
P = (a*A+b*B)/p
XPA = P-A
if XPA[r] != 0:
term1 = XPA[r] * self.S1d(r, pga_1, pgb)
if pga_1.shell[r] >= 0:
term2 = pga_1.shell[r] * (1 / (2 * p)) * self.S1d(r, pga_2, pgb)
if pgb.shell[r] >= 0:
term3 = pgb.shell[r] * (1 / (2 * p)) * self.S1d(r, pga_1, pgb_1)
return term1 + term2 + term3
def gaussian_factory(self, r, pga, pgb):
ca = pga.coefficient
cb = pgb.coefficient
a = pga.exponent
b = pgb.exponent
A = pga.origin
B = pgb.origin
i,k,m = pga.shell
j,l,n = pgb.shell
if r == 0:
pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)
pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)
pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)
return pga, pgb, pga_i_1, pga_i_2, pgb_j_1
elif r == 1:
pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)
pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)
pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)
return pga, pgb, pga_k_1, pga_k_2, pgb_l_1
elif r == 2:
pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)
pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)
pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)
return pga, pgb, pga_m_1, pga_m_2, pgb_n_1
[4]:
# The overlap integral between two primitive Gaussian type orbital
S_integral = Overlap()
print(S_integral(pga,pgb))
-8.880194157791243e-05
4.2. Kinetic
Consider the integral of electrons kinetic energy
\begin{equation} \begin{aligned} T_{a,b} &= \langle G_a \vert -\sum^N_{i=1}\frac{\hbar^2}{2m_i}\boldsymbol{\nabla}_i^2 \vert G_b \rangle\\ T_{a,b} &= T_{ij}S_{kl}S_{mn} + S_{ij}T_{kl}S_{mn} + S_{ij}S_{kl}T_{mn} \end{aligned} \end{equation}
The Obara-Saika recurrence relations for the Cartesian kinetic integrals over one direction is:
\begin{equation} \begin{aligned} T_{i+1,j} &= X_{PA}T_{ij} + \frac{1}{2p}(iT_{i-1,j}+jT_{i,j-1}) + \frac{b}{p}(2aS_{i+1,j}-iS_{i-1,j})\\ T_{i,j+1} &= X_{PB}T_{ij} + \frac{1}{2p}(iT_{i-1,j}+jT_{i,j-1})+ \frac{a}{p}(2bS_{i,j+1}-jS_{i,j-1})\\ \end{aligned} \end{equation}
With the boundary condtion
\begin{equation} T_{0,0} = \left[a-2a^2(X_{PA}^2 + \frac{1}{2p})\right] S_{00} \end{equation}
[5]:
overlap = Overlap()
class Kinetic(object):
"""The Obara-Saika scheme for three-dimensional kinetic energy integral over
primitive Gaussian orbitals.
Parameters
----------
a : float
Gaussian exponent facotr.
b : float
Gaussian exponent facotr.
i : int
Angular momentum quantum number.
j : int
Angular momentum quantum number.
A : float
Coordinate in on direction.
B : float
Coordinate in on direction.
Returns
-------
result : float
The non-normalizecd kinetic interals in one dimension.
"""
def __init__(self):
self.p = 0
self.mu = 0
self.P = ()
def __call__(self, pga, pgb):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
pga: PrimitiveGaussian
The first primitive gaussian orbital.
pgb: PrimitiveGaussian
The second primitive gaussian orbital.
C: List[float,float,float]
Coordinate of nuclei.
Return
------
result : float
Integral value.
"""
Sij = overlap.S1d(0,pga,pgb)
Skl = overlap.S1d(1,pga,pgb)
Smn = overlap.S1d(2,pga,pgb)
Tij = self.T1d(0,pga,pgb)
Tkl = self.T1d(1,pga,pgb)
Tmn = self.T1d(2,pga,pgb)
Tab = Tij*Skl*Smn+Sij*Tkl*Smn+Sij*Skl*Tmn
return Tab
def T1d(self, r, pga, pgb):
a = pga.exponent
b = pgb.exponent
p = a + b
mu = (a*b)/(a+b)
A = np.array(pga.origin)
B = np.array(pgb.origin)
P = (a*A+b*B)/p
XAB = A-B
XPA = P-A
if pga.shell[r] > 0:
return self.recursive(r, *self.gaussian_factory(r, pga, pgb))
elif pgb.shell[r] > 0:
return self.recursive(r, *self.gaussian_factory(r, pgb, pga))
else:
# Starting from the spherical Gaussians.
S00 = np.power(np.pi/p,0.5)*np.exp(-mu*XAB[r]**2)
T00 = (a-2*a**2*(XPA[r]**2+1./(2*p)))*S00
return T00
def recursive(self, r, pga, pgb, pga_1, pga_2, pgb_1):
term1 = term2 = term3 = term4 = term5 = 0
a = pga.exponent
b = pgb.exponent
p = a + b
mu = (a*b)/(a+b)
A = np.array(pga.origin)
B = np.array(pgb.origin)
P = (a*A+b*B)/p
XPA = P-A
if XPA[r] != 0:
term1 = XPA[r] * self.T1d(r, pga_1, pgb)
if pga_1.shell[r] >= 0:
term2 = pga_1.shell[r] * (1 / (2 * p)) * self.T1d(r, pga_2, pgb)
if pgb.shell[r] >= 0:
term3 = pgb.shell[r] * (1 / (2 * p)) * self.T1d(r, pga_1, pgb_1)
term4 = (2*a*b) / p * overlap.S1d(r, pga, pgb)
if pga_1.shell[r] >= 0:
term5 = pgb.shell[r] * (b / p) * overlap.S1d(r, pga_2, pgb)
return term1 + term2 + term3 + term4 - term5
def gaussian_factory(self, r, pga, pgb):
ca = pga.coefficient
cb = pgb.coefficient
a = pga.exponent
b = pgb.exponent
A = pga.origin
B = pgb.origin
i,k,m = pga.shell
j,l,n = pgb.shell
if r == 0:
pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)
pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)
pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)
return pga, pgb, pga_i_1, pga_i_2, pgb_j_1
elif r == 1:
pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)
pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)
pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)
return pga, pgb, pga_k_1, pga_k_2, pgb_l_1
elif r == 2:
pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)
pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)
pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)
return pga, pgb, pga_m_1, pga_m_2, pgb_n_1
[6]:
# The kinetic integral between two primitive Gaussian type orbital
T_integral = Kinetic()
print(T_integral(pga,pgb))
0.0016734344532709933
4.3. Nuclear Attraction
\begin{equation} \begin{aligned} V_{a,b} &=\langle G_a \vert -\sum^N_{i=1}\sum^M_{\alpha=1} \frac{Z_\alpha e^2}{\textbf{r}_{i\alpha}} \vert G_b \rangle \\ V_{a,b} &= V_{ijklmn}^{000} = \Theta_{ijklmn}^{0} \end{aligned} \end{equation}
The Obara-Saika recurrence relations for the Cartesian nuclear attraction integrals:
\begin{equation} \begin{aligned} \Theta_{i+1,j,k,l,m,n}^{N} &= X_{PA}\Theta_{ijkln}^{N} +\frac{1}{2p}(i\Theta_{i-1,j,k,l,m,n}^{N} + j\Theta_{i,j-1,k,l,m,n}^{N})\\ &-X_{PC}\Theta_{ijklmn}^{N+1} -\frac{1}{2p}(i\Theta_{i-1,j,k,l,m,n}^{N+1} + j\Theta_{i,j-1,k,l,m,n}^{N+1})\\ \Theta_{i,j+1,k,l,m,n}^{N} &= X_{PB}\Theta_{ijkln}^{N} +\frac{1}{2p}(i\Theta_{i-1,j,k,l,m,n}^{N} + j\Theta_{i,j-1,k,l,m,n}^{N})\\ &-X_{PC}\Theta_{ijklmn}^{N+1} -\frac{1}{2p}(i\Theta_{i-1,j,k,l,m,n}^{N+1} + j\Theta_{i,j-1,k,l,m,n}^{N+1})\\ \end{aligned} \end{equation}
With the boundary condtion
\begin{equation} \Theta_{000000}^{N} = \frac{2\pi}{p}(-2p)^{-N}K_{ab}^{xyz}R_{000}^{N} = \frac{2\pi}{p} K_{ab}^{xyz}F_N(pR_{PC}^2) \end{equation}
[7]:
class NuclearAttraction(object):
"""The Obara-Saika scheme for three-dimensional nuclear attraction integral over
primitive Gaussian orbitals.
Attributes
----------
p : float
The total exponent.
mu : float
The reduced exponent.
P : List[float,float,float]
The centre of charge coordinate.
C : List[float,float,float]
The coordinate of given nuclei.
Kab : float
The pre-exponential factor.
Methods
-------
__init__(self)
Initialize the instance.
"""
def __init__(self):
"""Initialize the instance.
"""
self.p = 0
self.mu = 0
self.Kab = 0
self.P = []
self.C = []
self.boys_dict = {}
def __call__(self, pga, pgb, C):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
pga: PrimitiveGaussian
The first primitive gaussian orbital.
pgb: PrimitiveGaussian
The second primitive gaussian orbital.
C: List[float,float,float]
Coordinate of nuclei.
Return
------
result : float
Integral value.
"""
l_total = sum(pga.shell) + sum(pgb.shell)
a = pga.exponent
b = pgb.exponent
p = a + b
mu = (a*b)/(a+b)
A = np.array(pga.origin)
B = np.array(pgb.origin)
P = (a*A+b*B)/p
RAB = np.linalg.norm(A-B)
RPA = np.linalg.norm(P-A)
RPB = np.linalg.norm(P-B)
RPC = np.linalg.norm(P-C)
Kab = exp(-mu*RAB**2)
self.p = p
self.mu = mu
self.P = P
self.C = C
self.Kab = Kab
# Build boys function F_{N}(x)
N = l_total
x = p*RPC**2
boys_pre_factor = (2*np.pi)/p*Kab
boys_function = boys(l_total, x)
Theta_N_000000 = boys_pre_factor * boys_function
self.boys_dict = {l_total: Theta_N_000000}
while N >= 1:
boys_function = boys_recursion(N, x, boys_function)
N -= 1
Theta_N_000000 = boys_pre_factor * boys_function
self.boys_dict[N] = Theta_N_000000
result = self.V(0, pga, pgb)
return result
def V(self, N, pga, pgb):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
N : int
Order of the boys function F_{N}(x).
pga : PrimitiveGaussian
The first primitive gaussian orbital.
pgb : PrimitiveGaussian
The second primitive gaussian orbital.
Return
------
vlaue : float
Integral value.
"""
if pga.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pga, pgb))
elif pga.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pga, pgb))
elif pga.shell[2] > 0:
return self.recursive(2, N, *self.gaussian_factory(2, pga, pgb))
elif pgb.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pgb, pga))
elif pgb.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pgb, pga))
elif pgb.shell[2] > 0:
result = self.recursive(2, N, *self.gaussian_factory(2, pgb, pga))
return result
else:
return self.boys_dict[N]
def recursive(self, r, N, pga, pgb, pga_1, pga_2, pgb_1):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
r : int
Cartesian index 0, 1, 2.
N : int
Order of the boys function F_{N}(x).
pga_1 : PrimitiveGaussian
The primitive gaussian orbital.
pgb : PrimitiveGaussian
The primitive gaussian orbital.
pga_2 : PrimitiveGaussian
The primitive gaussian orbital.
pgb_1 : PrimitiveGaussian
The primitive gaussian orbital.
Return
------
result : float
Integral value.
"""
term1 = term2 = term3 = term4 = term5 = term6 = 0
a = pga.exponent
b = pgb.exponent
p = a+b
A = np.array(pga.origin)
B = np.array(pgb.origin)
P = (a*A+b*B)/p
C = self.C
XPA = np.array(P) - np.array(A)
XPC = np.array(P) - np.array(C)
if np.array_equal(P,A) is False:
term1 = XPA[r] * self.V(N, pga_1, pgb)
if pga_1.shell[r] > 0:
term2 = pga_1.shell[r] * (1 / (2 * p)) * self.V(N, pga_2, pgb)
if pgb.shell[r] > 0:
term3 = pgb.shell[r] * (1 / (2 * p)) * self.V(N, pga_1, pgb_1)
if np.array_equal(P,C) is False:
term4 = XPC[r] * self.V(N+1, pga_1, pgb)
if pga_1.shell[r] > 0:
term5 = pga_1.shell[r] * (1 / (2 * p)) * self.V(N+1, pga_2, pgb)
if pgb.shell[r] > 0:
term6 = pgb.shell[r] * (1 / (2 * p)) * self.V(N+1, pga_1, pgb_1)
result = term1+term2+term3-term4-term5-term6
return result
def gaussian_factory(self, r, pga, pgb):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
r : int
Cartesian index 0, 1, 2.
N : int
Order of the boys function F_{N}(x).
pga : PrimitiveGaussian
The primitive gaussian orbital.
pgb : PrimitiveGaussian
The primitive gaussian orbital.
Return
------
result : Tuple(pg,pg,pg,pg)
Tuple of 4 PrimitiveGaussian orbital instance.
"""
ca = pga.coefficient
cb = pgb.coefficient
a = pga.exponent
b = pgb.exponent
A = pga.origin
B = pgb.origin
i,k,m = pga.shell
j,l,n = pgb.shell
if r == 0:
pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)
pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)
pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)
return pga, pgb, pga_i_1, pga_i_2, pgb_j_1
elif r == 1:
pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)
pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)
pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)
return pga, pgb, pga_k_1, pga_k_2, pgb_l_1
elif r == 2:
pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)
pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)
pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)
return pga, pgb, pga_m_1, pga_m_2, pgb_n_1
[8]:
# The nuclear attraction integral between two primitive Gaussian type orbital
V_integral = NuclearAttraction()
print(V_integral(pga,pgb,FCenter[0]))
-8.543259485351415e-05
4.4. Electron Repulsion
\begin{equation} V_{abcd} = \langle G_a G_b\vert \sum^N_{i=1}\sum^N_{j>i} \frac{e^2}{\textbf{r}_{ij}}\vert G_c G_d \rangle \end{equation}
The source and target integrals: \begin{equation} \begin{aligned} \Theta_{0000;0000;0000}^{N} &= \frac{2\pi^{2.5}}{pq\sqrt{p+q}} K_{ab}^{xyz}K_{cd}^{xyz}F_N(\alpha R_{PQ}^2)\\ \Theta_{i_xj_xk_xl_x;i_yj_yk_yl_y;i_zj_zk_zl_z}^{0} &= g_{i_xj_xk_xl_x;i_yj_yk_yl_y;i_zj_zk_zl_z} \end{aligned} \end{equation}
The Obara-Saika two electron recurrence relation \begin{equation} \begin{aligned} \Theta_{i+1,j,k,l}^{N} &= X_{PA}\Theta_{ijkl}^{N} -\frac{\alpha}{p}X_{PQ}\Theta_{i,j,k,l}^{N+1} +\frac{i}{2p}\left(\Theta_{i-1,j,k,l}^{N}-\frac{\alpha}{p}\Theta_{i-1,j,k,l}^{N+1}\right)\\ &+\frac{j}{2p}\left(\Theta_{i,j-1,k,l}^{N}-\frac{\alpha}{p}\Theta_{i,j-1,k,l}^{N+1}\right) -\frac{k}{2(p+q)}\Theta_{i,j,k-1,l}^{N+1} -\frac{l}{2(p+q)}\Theta_{i,j,k,l-1}^{N+1} \end{aligned} \end{equation} Using the horizontal recurrence relation, a similar relation may be written down for increments in j, replaceing \(X_{PA}\) with \(X_{PB}\).
\begin{equation} \begin{aligned} \Theta_{i,j,k+1,l}^{N} &= X_{QC}\Theta_{ijkl}^{N} -\frac{\alpha}{q}X_{PQ}\Theta_{i,j,k,l}^{N+1} +\frac{k}{2q}\left(\Theta_{i,j,k-1,l}^{N}-\frac{\alpha}{q}\Theta_{i,j,k-1,l}^{N+1}\right)\\ &+\frac{l}{2q}\left(\Theta_{i,j,k,l-1}^{N}-\frac{\alpha}{q}\Theta_{i,j,k,l-1}^{N+1}\right) -\frac{i}{2(p+q)}\Theta_{i-1,j,k,l}^{N+1} -\frac{j}{2(p+q)}\Theta_{i,j-1,k,l}^{N+1} \end{aligned} \end{equation}
[9]:
class ElectronRepulsion:
"""The Obara-Saika scheme for three-dimensional nuclear attraction integral over
primitive Gaussian orbitals.
Attributes
----------
p : float
The total exponent.
mu : float
The reduced exponent.
P : List[float,float,float]
The centre of charge coordinate.
C : List[float,float,float]
The coordinate of given nuclei.
Kab : float
The pre-exponential factor.
Methods
-------
__init__(self)
Initialize the instance.
"""
def __init__(self):
"""Initialize the instance.
"""
self.alpha = 0
self.R = []
self.boys_dict = {}
def __call__(self, pga, pgb, pgc, pgd):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
pga: PrimitiveGaussian
The first primitive gaussian orbital.
pgb: PrimitiveGaussian
The second primitive gaussian orbital.
C: List[float,float,float]
Coordinate of nuclei.
Return
------
result : float
Integral value.
"""
l_total = sum(pga.shell) + sum(pgb.shell) + sum(pgc.shell) + sum(pgd.shell)
a = pga.exponent
b = pgb.exponent
c = pgc.exponent
d = pgd.exponent
p = a + b
q = c + d
mu = (a*b)/(a+b)
nu = (c*d)/(c+d)
alpha = (p * q) / (p + q)
self.alpha = alpha
A = np.array(pga.origin)
B = np.array(pgb.origin)
C = np.array(pgc.origin)
D = np.array(pgd.origin)
P = (a*A+b*B)/(a+b)
Q = (c*C+d*D)/(c+d)
R = (p*P+q*Q)/(p+q)
self.R = R
RAB = np.linalg.norm(A-B)
RCD = np.linalg.norm(C-D)
RPQ = np.linalg.norm(P-Q)
# Build boys function F_{N}(x)
Kab = np.exp(-mu*RAB**2)
Kcd = np.exp(-nu*RCD**2)
boys_pre_factor = (2*np.pi**(5/2))/(p*q*np.sqrt(p+q))*Kab*Kcd
N = l_total
x = alpha*RPQ**2
boys_function = boys(l_total, x)
Theta_N_0000_0000_0000 = boys_pre_factor * boys_function
self.boys_dict = {l_total: Theta_N_0000_0000_0000}
while N >= 1:
boys_function = boys_recursion(N, x, boys_function)
N -= 1
Theta_N_0000_0000_0000 = boys_pre_factor * boys_function
self.boys_dict[N] = Theta_N_0000_0000_0000
result = self.Eri(0, pga, pgb, pgc, pgd)
return result
def Eri(self, N, pga, pgb, pgc, pgd):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
N : int
Order of the boys function F_{N}(x).
pga : PrimitiveGaussian
The primitive gaussian orbital.
pgb : PrimitiveGaussian
The primitive gaussian orbital.
pgc : PrimitiveGaussian
The primitive gaussian orbital.
pgd : PrimitiveGaussian
The primitive gaussian orbital.
Return
------
vlaue : float
Integral value.
"""
if pga.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pga, pgb, pgc, pgd))
elif pga.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pga, pgb, pgc, pgd))
elif pga.shell[2] > 0:
return self.recursive(2, N, *self.gaussian_factory(2, pga, pgb, pgc, pgd))
elif pgb.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pgb, pga, pgd, pgc))
elif pgb.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pgb, pga, pgd, pgc))
elif pgb.shell[2] > 0:
return self.recursive(2, N, *self.gaussian_factory(2, pgb, pga, pgd, pgc))
elif pgc.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pgc, pgd, pga, pgb))
elif pgc.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pgc, pgd, pga, pgb))
elif pgc.shell[2] > 0:
return self.recursive(2, N, *self.gaussian_factory(2, pgc, pgd, pga, pgb))
elif pgd.shell[0] > 0:
return self.recursive(0, N, *self.gaussian_factory(0, pgd, pgc, pgb, pga))
elif pgd.shell[1] > 0:
return self.recursive(1, N, *self.gaussian_factory(1, pgd, pgc, pgb, pga))
elif pgd.shell[2] > 0:
return self.recursive(2, N, *self.gaussian_factory(2, pgd, pgc, pgb, pga))
else:
return self.boys_dict[N]
def recursive(self, r, N, pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
r : int
Cartesian index 0, 1, 2.
N : int
Order of the boys function F_{N}(x).
pga_1 : PrimitiveGaussian
The primitive gaussian orbital.
pgb : PrimitiveGaussian
The primitive gaussian orbital.
pga_2 : PrimitiveGaussian
The primitive gaussian orbital.
pgb_1 : PrimitiveGaussian
The primitive gaussian orbital.
Return
------
result : float
Integral value.
"""
term1 = term2 = term3 = term4 = term5 = term6 = term7 = term8 = 0
a = pga.exponent
b = pgb.exponent
c = pgc.exponent
d = pgd.exponent
p = a + b
q = c + d
alpha = (p*q)/(p+q)
#alpha = self.alpha
A = np.array(pga.origin)
B = np.array(pgb.origin)
C = np.array(pgc.origin)
D = np.array(pgd.origin)
P = (a*A+b*B)/(a+b)
Q = (c*C+d*D)/(c+d)
XPA = P - A
XPQ = P - Q
if XPA[r] != 0:
term1 = XPA[r] * self.Eri(N, pga_1, pgb, pgc, pgd)
if XPQ[r] != 0:
term2 = alpha/p*XPQ[r] * self.Eri(N+1, pga_1, pgb, pgc, pgd)
if pga_1.shell[r] > 0:
term3 = pga_1.shell[r] * (1 / (2 * p)) * self.Eri(N, pga_2, pgb, pgc, pgd)
term4 = pga_1.shell[r] * (alpha / (2 * p ** 2)) * self.Eri(N+1, pga_2, pgb, pgc, pgd)
if pgb.shell[r] > 0:
term5 = pgb.shell[r] * (1 / (2 * p)) * self.Eri(N, pga_1, pgb_1, pgc, pgd)
term6 = pgb.shell[r] * (alpha / (2 * p ** 2)) * self.Eri(N+1, pga_1, pgb_1, pgc, pgd)
if pgc.shell[r] > 0:
term7 = pgc.shell[r] * (1 / (2 * (p + q))) * self.Eri(N+1, pga_1, pgb, pgc_1, pgd)
if pgd.shell[r] > 0:
term8 = pgd.shell[r] * (1 / (2 * (p + q))) * self.Eri(N+1, pga_1, pgb, pgc, pgd_1)
return term1 - term2 + term3 - term4 + term5 - term6 + term7 + term8
def gaussian_factory(self, r, pga, pgb, pgc, pgd):
"""Evaluates nuclear attraction integral over two primitive gaussian orbitals.
Parameters
----------
r : int
Cartesian index 0, 1, 2.
N : int
Order of the boys function F_{N}(x).
pga : PrimitiveGaussian
The primitive gaussian orbital.
pgb : PrimitiveGaussian
The primitive gaussian orbital.
Return
------
result : Tuple(pg,pg,pg,pg)
Tuple of 4 PrimitiveGaussian orbital instance.
"""
ca = pga.coefficient
cb = pgb.coefficient
cc = pgc.coefficient
cd = pgd.coefficient
a = pga.exponent
b = pgb.exponent
c = pgc.exponent
d = pgd.exponent
A = pga.origin
B = pgb.origin
C = pgc.origin
D = pgd.origin
ix,iy,iz = pga.shell
jx,jy,jz = pgb.shell
kx,ky,kz = pgc.shell
lx,ly,lz = pgd.shell
if r == 0:
pga_1 = PrimitiveGaussian(ca, A, (ix - 1, iy, iz), a)
pga_2 = PrimitiveGaussian(ca, A, (ix - 2, iy, iz), a)
pgb_1 = PrimitiveGaussian(cb, B, (jx - 1, jy, jz), b)
pgc_1 = PrimitiveGaussian(cc, C, (kx - 1, ky, kz), c)
pgd_1 = PrimitiveGaussian(cd, D, (lx - 1, ly, lz), d)
return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1
elif r == 1:
pga_1 = PrimitiveGaussian(ca, A, (ix, iy-1, iz), a)
pga_2 = PrimitiveGaussian(ca, A, (ix, iy-2, iz), a)
pgb_1 = PrimitiveGaussian(cb, B, (jx, jy-1, jz), b)
pgc_1 = PrimitiveGaussian(cc, C, (kx, ky-1, kz), c)
pgd_1 = PrimitiveGaussian(cd, D, (lx, ly-1, lz), d)
return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1
elif r == 2:
pga_1 = PrimitiveGaussian(ca, A, (ix, iy, iz-1), a)
pga_2 = PrimitiveGaussian(ca, A, (ix, iy, iz-2), a)
pgb_1 = PrimitiveGaussian(cb, B, (jx, jy, jz-1), b)
pgc_1 = PrimitiveGaussian(cc, C, (kx, ky, kz-1), c)
pgd_1 = PrimitiveGaussian(cd, D, (lx, ly, lz-1), d)
return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1
[10]:
# The overlap integral between two primitive Gaussian type orbital
Eri_integral = ElectronRepulsion()
print(Eri_integral(pga,pgb,pga,pgb))
1.9060888184873294e-08