STO-nG Fitting

Quantum Chemistry

STO-nG basis sets approximate Slater-type orbitals (STOs) as linear combinations of Gaussian-type orbitals (GTOs). This page demonstrates how to fit an STO using a Levenberg-Marquardt optimization algorithm.

Mathematical Formulation

A normalized 1s Slater orbital centered at the origin has the form:

We approximate this as a linear combination of Gaussian primitives:

where are contraction coefficients and are Gaussian exponents.

Optimization Problem

The fitting problem involves finding parameters that minimize the residual:

Note that we use in the parameterization to ensure positivity of the exponents.

Implementation

The following Python code implements STO-nG fitting using scipy's Levenberg-Marquardt algorithm:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.animation import FFMpegWriter
from scipy.optimize import least_squares

# ============================================================
# Orbitals
# ============================================================

def slater_1s(r, zeta=1.0):
    """Normalized Slater 1s orbital"""
    return (zeta**3 / np.pi)**0.5 * np.exp(-zeta * r)

def gaussian_1s(r, alpha):
    """Normalized Gaussian 1s primitive"""
    return (2 * alpha / np.pi)**0.75 * np.exp(-alpha * r**2)

def sto_ng_model(r, params):
    """
    params = [c1,...,cn, log(alpha1),...,log(alphan)]
    log(alpha) parametrization ensures positivity
    """
    n = len(params) // 2
    c = params[:n]
    alpha = np.exp(params[n:])

    psi = np.zeros_like(r)
    for i in range(n):
        psi += c[i] * gaussian_1s(r, alpha[i])
    return psi

def residuals(params, r, target):
    return sto_ng_model(r, params) - target

# ============================================================
# Grid
# ============================================================

r = np.linspace(0.0, 6.0, 800)
sto = slater_1s(r)

# ============================================================
# Plot styling (publication quality)
# ============================================================

def set_publication_style():
    """Set publication-quality matplotlib style."""
    mpl.rcParams.update({
        'font.family': 'serif',
        'font.size': 12,
        'axes.labelsize': 14,
        'axes.titlesize': 16,
        'axes.linewidth': 1.2,
        'axes.labelpad': 8,
        'axes.titlepad': 10,
        'xtick.labelsize': 12,
        'ytick.labelsize': 12,
        'xtick.direction': 'in',
        'ytick.direction': 'in',
        'xtick.top': True,
        'ytick.right': True,
        'xtick.major.size': 6,
        'ytick.major.size': 6,
        'xtick.major.width': 1.2,
        'ytick.major.width': 1.2,
        'legend.fontsize': 12,
        'legend.frameon': False,
        'lines.linewidth': 2,
        'lines.markersize': 6,
        'figure.dpi': 100,
        'savefig.dpi': 300,
        'savefig.bbox': 'tight'
    })

set_publication_style()

fig, ax = plt.subplots()
ax.set_xlim(0, 6)
ax.set_ylim(0, 1.1 * sto.max())
ax.set_xlabel(r"$r$ (a.u.)")
ax.set_ylabel(r"$\psi(r)$")

# Exact reference
ax.plot(r, sto, "k--", label="Slater 1s (exact)")

# Animated fit
fit_line, = ax.plot(r, np.zeros_like(r), label="STO-1G (LM fit)")

legend = ax.legend(loc="upper right", frameon=True, framealpha=1.0)
title = ax.set_title("STO-nG Levenberg–Marquardt Fit")

# ============================================================
# Animation (built-in LM)
# ============================================================

writer = FFMpegWriter(
    fps=1,
    metadata={"title": "STO-nG Built-in LM Fit"}
)

output_file = "sto_ng_lm.mp4"

with writer.saving(fig, output_file, dpi=400):

    for n in range(1, 7):

        # Initial guess
        c0 = np.ones(n) / n
        alpha0 = np.linspace(0.3, 2.5, n)
        p0 = np.concatenate([c0, np.log(alpha0)])

        # Built-in Levenberg–Marquardt
        res = least_squares(
            residuals,
            p0,
            args=(r, sto),
            method="lm",
            xtol=1e-12,
            ftol=1e-12,
            gtol=1e-12,
            max_nfev=2000
        )

        fit = sto_ng_model(r, res.x)

        fit_line.set_ydata(fit)
        fit_line.set_label(f"STO-{n}G (LM fit)")

        legend.remove()
        legend = ax.legend(loc="upper right", frameon=True, framealpha=1.0)

        title.set_text(f"STO-{n}G Levenberg–Marquardt Fit")

        writer.grab_frame()

plt.close(fig)
print(f"Saved animation to {output_file}")

Key Features