Fitting the parameters in STO-LG
In computational chemistry, the Slater-type orbital (STO) more accurately describes the qualitative features of the molecular orbitals than Gaussian functions (GF). However, calculating the two-electron integral using STO can be costly. On the other hand, integrating GFs is relatively cheap. One way to solve this problem is to use a linear combination of GFs to approximate a STO. Such linear combination of Gaussian functions is called contracted Gaussian functions (CGF).
$$\phi_\mu^{CGF}(\vec{r}-\vec{R}_A) = \sum_{p=1}^L d_{p\mu} \phi_p^{GF} (\alpha_{p\mu}, \vec{r} - \vec{R}_A)$$
where L is the length of the contraction, $d_{p\mu}$ and $\alpha_{p\mu}$ are contraction coefficients and contraction exponents, respectively. Hence, the so-called STO-LG strategy uses L Gaussian-type orbitals to approximate one STO function.
Approximating 1s Slater-type function using STO-LG
The expressions for 1s STO and GF are
\begin{align} \phi_{1s}^{STO} (\zeta, \vec{r}) = \left( \frac{\zeta^3}{\pi} \right)^{1/2} e^{-\zeta |\vec{r}-{\vec{R_A}}|} \notag \\ \phi_{1s}^{GF}(\alpha, \vec{r}) = \left(\frac{2\alpha}{\pi}\right)^{3/4} e^{-\alpha |\vec{r}-{\vec{R _A}}|^2} \notag \end{align}
where both orbitals have their corresponding parameters. The goal is to find the $d_{p}$ and $\alpha_{p}$ in the following equation
$$ \phi_{1s}^{STO} (\zeta, \vec{r}) = \sum_p^L d_{p}\phi_{1s}^{GF}(\alpha_p, \vec{r}) $$
STO-1G with $\zeta = 1.0$
In the first case, we will show the process of fitting the simplest function STO-1G by assuming $\zeta=1.0$ in the STO. Basically we will solve the following equation for $\alpha_{11}$
$$ \phi_{1s}^{STO} (\zeta=1.0, \vec{r}) = \phi_{1s}^{GF}(\alpha_{11}, \vec{r}) $$
from sympy import *
import numpy as np
zeta, r, alpha, a11 = symbols("zeta r alpha a11", positive=True) # define the variables
sto = (zeta **3 / pi) ** (1/2) * exp(-zeta * r) # general expression for one STO
gf = (2 * alpha /pi) ** (3/4) * exp(-alpha * r**2) # general expression for one GF
sto_1 = sto.subs(zeta, 1.0) # zeta = 1.0
gf_1 = gto.subs(alpha, a11) # alpha = a11
Instead of minimize the differences between the two functions, we will maximize the overlap between the GF and the STO following Szabo’s book.
S = integrate(sto_1 * gf_1 * r**2, (r, 0, oo)) * 4 * pi # the overlap between STO(1.0, r) and GF(alpha, r)
We will maximize this overlap in terms of $a_{11}$. Since we will use scipy, we turn the maximization problem into minimzation of the negative of the overlap S.
# function to minimize
def func(a):
res = S.subs(a11, a[0]).evalf()
return -res
from scipy.optimize import minimize
res = minimize(func, x0=[0.2], bounds=[(0.1, 1)])
res.x[0]
0.2709502078346618
We get the $\alpha_{11}$ value as 0.2709497296298158. This is almost identical to the result from Szabo’s book.
Show the STO-1G plot
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['font.size'] = 22
plt.rcParams['font.family'] = 'Arial'
gf_a1 = lambdify(r, gf_1.subs(a11, res.x[0]), "numpy")
sto_1_np = lambdify(r, sto_1, "numpy")
r_np = np.linspace(0, 10, 101)
plt.plot(r_np, gf_a1(r_np), label='STO-1G')
plt.plot(r_np, sto_1_np(r_np), label='STO')
plt.legend()
plt.xlabel("$r$")
<matplotlib.text.Text at 0x7fcc2089b898>
The results are reasonably good.
STO-LG
We will code the general procedure to calculate $L>1$ CGFs.
from IPython.display import display, Math
def get_gto(d, alpha):
return d * (2 * alpha /pi) ** (3/4) * exp(-alpha * r**2)
def get_symbols(L):
ds = symbols(f'd:{L}', positive=True)
alphas = symbols(f'alpha:{L}', positive=True)
return ds, alphas
class STOLG:
def __init__(self, L=3, zta=1.0):
self.L = L
self.ds, self.alphas = get_symbols(L)
self.GFs = [get_gto(d, a) for d, a in zip(self.ds, self.alphas)]
self.GF_sum = sum(self.GFs)
self.gg_int = integrate(self.GF_sum * self.GF_sum * r**2, (r, 0, oo)) * 4 * pi
self.gg_int = self.gg_int.evalf()
self.sto = sto.subs(zeta, zta)
self.S = integrate(self.GF_sum * self.sto * r**2, (r, 0, oo)) * 4 * pi
def fit(self):
def _func(x):
subs = {i: j for i, j in zip(self.ds[1:]+self.alphas, x)}
d0_val = solve(self.gg_int.subs(subs) - 1)
subs[self.ds[0]] = d0_val[0]
val = self.S.subs(subs).evalf()
# print(subs)
self.subs = subs
return -float(val)
# initial guesses
d_vals = np.linspace(0.1, 0.9, self.L-1).tolist()
alpha_vals = np.linspace(0.1, 0.9, self.L).tolist()
self.res = minimize(_func, x0=d_vals + alpha_vals, bounds=[(0, 10)] * (2 * self.L - 1))
return self.res
@property
def expression(self):
expr = []
for i, j in zip(self.ds, self.alphas):
expr.append(r"%.6f\phi^{GF}(%.6f)" % (self.subs[i], self.subs[j]))
return display(Math('+'.join(expr)))
@property
def func(self):
return lambdify(r, self.GF_sum.subs(self.subs), "numpy")
@property
def funcs(self):
return [lambdify(r, i.subs(self.subs), "numpy") for i in self.GFs]
def __call__(self, r):
return self.func(r)
def plot(self, r):
plt.plot(r, self(r), '-', label=f'STO-{self.L}G')
for i in range(self.L):
plt.plot(r, self.funcs[i](r), '--', label=f'GF-{i}')
plt.plot(r, sto_1_np(r), '-', label='STO')
plt.legend()
STO-2G
sto2g = STOLG(L=2)
sto2g.fit()
fun: -0.9984197028799346
hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
jac: array([ 0.00000000e+00, 2.40918396e-06, -4.10782519e-07])
message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
nfev: 96
nit: 19
status: 0
success: True
x: array([0.43013353, 0.15162213, 0.85180271])
The overlap has reached 0.998 with only two GFs.
sto2g.expression
$$0.678908\phi^{GF}(0.151622)+0.430134\phi^{GF}(0.851803)$$
The expression above matches with Equation (3.220) in Szabo’s book.
sto2g.plot(r_np)
STO-3G
sto3g = STOLG(3)
sto3g.fit()
fun: -0.9998347361981794
hess_inv: <5x5 LbfgsInvHessProduct with dtype=float64>
jac: array([5.09592368e-06, 1.04027897e-05, 1.17794663e-05, 5.55111512e-06,
2.22044605e-08])
message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
nfev: 276
nit: 39
status: 0
success: True
x: array([0.53532369, 0.15432918, 0.10982016, 0.40578573, 2.22784233])
sto3g.expression
$$0.444642\phi^{GF}(0.109820)+0.535324\phi^{GF}(0.405786)+0.154329\phi^{GF}(2.227842)$$
Again the STO-3G expression matches with Equation (3.221) in Szabo’s book.
sto3g.plot(r_np)
Summary
In this notebook, I show how we can fit the parameters in the contracted Gaussian functions. The results are relatively sensitive to the initial guesses given the optimizers. I believe it will be more so if L further increases.
I also see that with changing $\zeta$, the optimizer gives me results different from the scaling relationships. It will be interesting to further investigate the cause.