In quantum computational chemistry, gaussian basis set is often used as wave functions of atomic orbitals.
Primitive Gaussian-Type Orbital (GTO)
One of the most basic atomic orbital is Slayter-type orbitals (STOs) like hydrogen like orbitals. However, analytical integrals of STOs are difficult. So Gaussian-type orbitals (GTOs) are often used instead of STOs. In addition, atomic orbitals are described as linear combination of primitive gaussian type orbital.
Primitive GTO is described as
where
- N is normalization constant.
- A is an atom center in cartesian coordinate.
- n, l, m are angular momentum powers (they are different from quantum number n, l, m ).
GTO is described using vector form. To clarify, GTO can be written as
Because
STO-nG
STO-nG is a basis set of immitating STO by “n primitive GTO”.
For example, in case of STO-3G which is well known basis set in quantum chemistry, atomic orbitals are described using 3 GTOs.
where
Note that normalization constants of gaussian basis sets have 2 factors.
- Primitive Normalization: The square of each GTO is needed to be 1.
Primitive normalization constant can be calculated analyticaly as
- Contracted Normalization: The square of total orbital φ(r) must be 1.
STO-3G Hydorogen 1s orbital
Let’s compare STO with STO-3G in case of Hydrogen 1s orbital. The figure is shown below.
Almost all orbital shape of STO-3G is close to STO. But change of STO-3G around zero is smaller than that of STO.
Code of Drawing STO and STO-3G
I show my code drawing orbitals.
"""
Created on Sat Jul 11 23:58:53 2026
@author: kaz2255pq
"""
import numpy as np
import scipy.special
import matplotlib.pyplot as plt
def basis_set(Z,qnl):
# STO-3G
n = qnl[0]
l = qnl[1]
m = qnl[2]
if Z == 1:
if n == 1:
a1 = 0.1543289673E+00
a2 = 0.5353281423E+00
a3 = 0.4446345422E+00
A1 = 0.3425250914E+01
A2 = 0.6239137298E+00
A3 = 0.1688554040E+00
if Z == 3:
if n == 1:
a1 = 0.1543289673E+00
a2 = 0.5353281423E+00
a3 = 0.4446345422E+00
A1 = 0.1611957475E+02
A2 = 0.2936200663E+01
A3 = 0.7946504870E+00
elif n == 2 and l == 0:
a1 = -0.9996722919E-01
a2 = 0.3995128261E+00
a3 = 0.7001154689E+00
A1 = 0.6362897469E+00
A2 = 0.1478600533E+00
A3 = 0.4808867840E-01
elif n == 2 and l == 1:
a1 = 0.1559162750E+00
a2 = 0.6076837186E+00
a3 = 0.3919573931E+00
A1 = 0.6362897469E+00
A2 = 0.1478600533E+00
A3 = 0.4808867840E-01
if Z == 6:
if n == 1:
a1 = 0.1543289673E+00
a2 = 0.5353281423E+00
a3 = 0.4446345422E+00
A1 = 0.7161683735E+02
A2 = 0.1304509632E+02
A3 = 0.3530512160E+01
elif n == 2 and l == 0:
a1 = -0.9996722919E-01
a2 = 0.3995128261E+00
a3 = 0.7001154689E+00
A1 = 0.2941249355E+01
A2 = 0.6834830964E+00
A3 = 0.2222899159E+00
elif n == 2 and l == 1:
a1 = 0.1559162750E+00
a2 = 0.6076837186E+00
a3 = 0.3919573931E+00
A1 = 0.2941249355E+01
A2 = 0.6834830964E+00
A3 = 0.2222899159E+00
if Z == 8:
if n == 1:
a1 = 0.1543289673E+00
a2 = 0.5353281423E+00
a3 = 0.4446345422E+00
A1 = 0.1307093214E+03
A2 = 0.2380886605E+02
A3 = 0.6443608313E+01
elif n == 2 and l == 0:
a1 = -0.9996722919E-01
a2 = 0.3995128261E+00
a3 = 0.7001154689E+00
A1 = 0.5033151319E+01
A2 = 0.1169596125E+01
A3 = 0.3803889600E+00
elif n == 2 and l == 1:
a1 = 0.1559162750E+00
a2 = 0.6076837186E+00
a3 = 0.3919573931E+00
A1 = 0.5033151319E+01
A2 = 0.1169596125E+01
A3 = 0.3803889600E+00
return np.array([a1,a2,a3,A1,A2,A3])
def angm_to_nlm(l,m):
if l == 0:
return (0,0,0)
if l == 1:
if m == -1:
return (1,0,0) #px
elif m == 1:
return (0,1,0) #py
elif m == 0:
return (0,0,1) #pz
def fact2(n):
if n <= 0:
return 1
return scipy.special.factorial2(n)
def primitive_gto(x,y,z,A,dA,aL):
l,m,n = aL
dx = x - A[0]
dy = y - A[1]
dz = z - A[2]
r2 = dx**2 + dy**2 + dz**2
poly = dx**l * dy**m * dz**n
val = poly * np.exp(-dA*r2)
#print(val.shape)
return val
def contracted_gto(x,y,z,A,dA,coef, aL):
val = 0
CN = contraction_norm(coef, dA, aL,tuple(A),tuple(A))
for beta,alpha in zip(dA,coef):
N = primitive_norm(beta,aL)
val += alpha *primitive_gto(x,y,z,A,beta,aL) * N
return val*CN
def primitive_norm(alpha,nL):
l = nL[0]
m = nL[1]
n = nL[2]
L = l+m+n
#num = (4*alpha)**L
num = (4*alpha)**L
den = fact2(2*l-1)*fact2(2*m-1)*fact2(2*n-1)
return (2*alpha/np.pi)**0.75 * np.sqrt(num/den)
def contraction_norm(coeff, d, aL,keyA,keyB):
S = 0.0
A = np.zeros(3)
for i in range(len(d)):
for j in range(len(d)):
Ni = primitive_norm(d[i],aL)
Nj = primitive_norm(d[j],aL)
Sij = overlap_os( aL,aL,A,A,d[i],d[j],keyA,keyA )
S += ( coeff[i] * coeff[j] * Ni * Nj * Sij )
return 1/np.sqrt(S)
def overlap_ss(aL,bL,A,B,da,db):
gamma = da + db
mu = da*db/gamma
Rab2 = np.dot(A-B,A-B)
return (np.pi/gamma)**1.5*np.exp(-mu*Rab2)
cache_ol = {}
def overlap_os(aL,bL,A,B,da,db,keyA,keyB):
key = (
aL,
bL,
keyA,
keyB,
da,
db
)
if key in cache_ol:
return cache_ol[key]
if aL==(0,0,0) and bL==(0,0,0):
val = overlap_ss(aL,bL,A,B,da,db)
cache_ol[key] = val
return val
gamma = da+db
P = (da*A + db*B)/gamma
for i in range(3):
if aL[i] > 0:
a_down = list(aL)
a_down[i] -= 1
a_down = tuple(a_down)
term1 = ( (P[i]-A[i]) * overlap_os(a_down,bL,A,B,da,db,keyA,keyB) )
term2 = 0.0
if a_down[i] >= 0:
a_down2 = list(a_down)
a_down2[i] -= 1
if a_down2[i] >= 0:
a_down2 = tuple(a_down2)
term2 = a_down[i]/(2*gamma) * overlap_os( a_down2,bL,A,B,da,db,keyA,keyB )
term3 = 0.0
if bL[i] > 0:
b_down = list(bL)
b_down[i] -= 1
b_down = tuple(b_down)
term3 = bL[i]/(2*gamma) * overlap_os( a_down,b_down,A,B,da,db,keyA,keyB )
val = (term1 + term2 + term3)
#if abs(val) < 1e-8:
# val = 0
cache_ol[key] = val
return val
for i in range(3):
if bL[i] > 0:
b_down = list(bL)
b_down[i] -= 1
b_down = tuple(b_down)
term1 = (P[i]-B[i]) * overlap_os(aL,b_down,A,B,da,db,keyA,keyB)
term2 = 0.0
if b_down[i] >= 0:
b_down2 = list(b_down)
b_down2[i] -= 1
if b_down2[i] >= 0:
b_down2 = tuple(b_down2)
term2 = b_down[i]/(2*gamma)* overlap_os(aL,b_down2,A,B,da,db,keyA,keyB)
term3 = 0.0
if aL[i] > 0:
a_down = list(aL)
a_down[i] -= 1
a_down = tuple(a_down)
term3 = aL[i]/(2*gamma)* overlap_os(a_down,b_down,A,B,da,db,keyA,keyB)
val = (term1 + term2 + term3)
#if abs(val) < 1e-8:
# val = 0
cache_ol[key] = val
return val
def sto(x,y,z,zeta,qnl):
#zeta = Z/n
n = qnl[0]
l = qnl[1]
m = qnl[2]
r = (x**2+y**2+z**2)**0.5
if n==1:
orb = 1/(np.pi)**0.5*zeta**1.5*np.exp(-zeta*r)
elif n==2 and l==0:
orb = 1/(3*np.pi)**0.5*zeta**2.5*r*np.exp(-zeta*r)
elif n==2 and l==1 and m==0:
orb = 1/(np.pi)**0.5*zeta**2.5*x*np.exp(-zeta*r)
elif n==2 and l==1 and m==-1:
orb = 1/(np.pi)**0.5*zeta**2.5*y*np.exp(-zeta*r)
elif n==2 and l==1 and m==1:
orb = 1/(np.pi)**0.5*zeta**2.5*z*np.exp(-zeta*r)
elif n==3 and l==0:
orb = 2/3/(10*np.pi)**0.5*zeta**2.7*r**2*np.exp(-zeta*r)
elif n==3 and l==1 and m==0:
orb = 1/(15*np.pi)**0.5*zeta**2.7*r*x*np.exp(-zeta*r)
elif n==3 and l==1 and m==-1:
orb = 1/(15*np.pi)**0.5*zeta**2.7*r*y*np.exp(-zeta*r)
elif n==3 and l==1 and m==1:
orb = 1/(15*np.pi)**0.5*zeta**2.7*r*z*np.exp(-zeta*r)
elif n==3 and l==2 and m==0:
orb = 1/(6*np.pi)**0.5*zeta**2.7*(3*z**2-r**2)/(3**0.5)*np.exp(-zeta*r)
elif n==3 and l==2 and m==-1:
orb = 1/(6*np.pi)**0.5*zeta**2.7*(2*x*z)*np.exp(-zeta*r)
elif n==3 and l==2 and m==1:
orb = 1/(6*np.pi)**0.5*zeta**2.7*(2*y*z)*np.exp(-zeta*r)
elif n==3 and l==2 and m==-2:
orb = 1/(6*np.pi)**0.5*zeta**2.7*(x**2-y**2)*np.exp(-zeta*r)
elif n==3 and l==2 and m==2:
orb = 1/(6*np.pi)**0.5*zeta**2.7*(2*x*y)*np.exp(-zeta*r)
return orb
############
### Main ###
qnl = [1,0,0]
basisA = basis_set(1,qnl)
aL = angm_to_nlm(0,0)
dA = basisA[3:]
coef = basisA[:3]
A = np.array([0,0,0])
x = np.linspace(0,5,100)
y = np.zeros(len(x))
z = np.zeros(len(x))
chi = contracted_gto(x,y,z,A,dA,coef,aL)
chi_sto = sto(x,y,z,1,qnl)
fig, ax = plt.subplots(1,1,figsize=(6,4))
plt.rcParams["font.family"] = "times new roman"
ax.set_title("Hydrogen")
ax.set_xlabel("r [atomic unit]")
ax.set_ylabel("χ(r)")
ax.plot(x,chi,label="STO-3G")
ax.plot(x,chi_sto, label ="STO")
plt.legend()
plt.show()
The basis set data included in this code were obtained from the
Basis Set Exchange (BSE).
Basis Set Exchange:
https://www.basissetexchange.org/
コメント