import numpy as np
from ..auxiliary import gaussian_product_center, boys
[docs]def E(i,j,t,A,B,a,b):
''' Recursive definition of Hermite Gaussian coefficients.
Returns a float.
a: orbital exponent on Gaussian 'a' (e.g. alpha in the text)
b: orbital exponent on Gaussian 'b' (e.g. beta in the text)
i,j: orbital angular momentum number on Gaussian 'a' and 'b'
t: number nodes in Hermite (depends on type of integral,
e.g. always zero for overlap integrals)
Qx: distance between origins of Gaussian 'a' and 'b'
'''
p = a + b
q = a*b/p
AB = A-B
if (t < 0) or (t > (i + j)):
# out of bounds for t
return 0.0
elif i == j == t == 0:
# base case
return np.exp(-q*AB*AB) # K_AB
elif j == 0:
# decrement index i
return (1/(2*p))*E(i-1,j,t-1,A,B,a,b) - \
(q*AB/a)*E(i-1,j,t,A,B,a,b) + \
(t+1)*E(i-1,j,t+1,A,B,a,b)
else:
# decrement index j
return (1/(2*p))*E(i,j-1,t-1,A,B,a,b) + \
(q*AB/b)*E(i,j-1,t,A,B,a,b) + \
(t+1)*E(i,j-1,t+1,A,B,a,b)
[docs]def R(t,u,v,n,p,PCx,PCy,PCz,RPC):
''' Returns the Coulomb auxiliary Hermite integrals
Returns a float.
Arguments:
t,u,v: order of Coulomb Hermite derivative in x,y,z
(see defs in Helgaker and Taylor)
n: order of Boys function
PCx,y,z: Cartesian vector distance between Gaussian
composite center P and nuclear center C
RPC: Distance between P and C
'''
T = p*RPC*RPC
val = 0.0
if t == u == v == 0:
val += np.power(-2*p,n)*boys(n,T)
elif t == u == 0:
if v > 1:
val += (v-1)*R(t,u,v-2,n+1,p,PCx,PCy,PCz,RPC)
val += PCz*R(t,u,v-1,n+1,p,PCx,PCy,PCz,RPC)
elif t == 0:
if u > 1:
val += (u-1)*R(t,u-2,v,n+1,p,PCx,PCy,PCz,RPC)
val += PCy*R(t,u-1,v,n+1,p,PCx,PCy,PCz,RPC)
else:
if t > 1:
val += (t-1)*R(t-2,u,v,n+1,p,PCx,PCy,PCz,RPC)
val += PCx*R(t-1,u,v,n+1,p,PCx,PCy,PCz,RPC)
return val
[docs]def Vxyz(a,lmn1,A,b,lmn2,B,C):
''' Evaluates kinetic energy integral between two Gaussians
Returns a float.
a: orbital exponent on Gaussian 'a' (e.g. alpha in the text)
b: orbital exponent on Gaussian 'b' (e.g. beta in the text)
lmn1: int tuple containing orbital angular momentum (e.g. (1,0,0))
for Gaussian 'a'
lmn2: int tuple containing orbital angular momentum for Gaussian 'b'
A: list containing origin of Gaussian 'a', e.g. [1.0, 2.0, 0.0]
B: list containing origin of Gaussian 'b'
C: list containing origin of nuclear center 'C'
'''
i,k,m = lmn1
j,l,n = lmn2
p = a + b
P = gaussian_product_center(a,A,b,B) # Gaussian composite center
RPC = np.linalg.norm(P-C)
val = 0.0
for t in range(i+j+1):
for u in range(k+l+1):
for v in range(m+n+1):
val += E(i,j,t,A[0],B[0],a,b) * \
E(k,l,u,A[1],B[1],a,b) * \
E(m,n,v,A[2],B[2],a,b) * \
R(t,u,v,0,p,P[0]-C[0],P[1]-C[1],P[2]-C[2],RPC)
val *= 2*np.pi/p
return val
[docs]def nuclear_attraction(a,b,C):
'''Evaluates overlap between two contracted Gaussians
Returns float.
Arguments:
a: contracted Gaussian 'a', BasisFunction object
b: contracted Gaussian 'b', BasisFunction object
C: center of nucleus
'''
v = 0.0
for ia, ca in enumerate(a.coefs):
for ib, cb in enumerate(b.coefs):
v += a.norm[ia]*b.norm[ib]*ca*cb*\
Vxyz(a.exps[ia],a.shell,a.origin,
b.exps[ib],b.shell,b.origin,C)
return v