1s Gaussian Type Orbital Integrals

Quantum Chemistry

1s Gaussian-type orbitals are some of the simplest basis functions to work with from a computational perspective. This page shows you the key results that you need to get up and running with a 1s program. The derivation of these integrals can be at times tedious and challenging so they are omitted here, with the exception of the overlap integral. If you want to get computing quickly, you will have to just trust these integrals and when you would like to learn more you can always study the derivations later. Treat this like the quantum chemistry version of your integral table in your freshman calculus book.

Gaussian Product Theorem

This is the workhorse of these integrals

g_{1 s, α} (r - R_{A}) \times g_{1 s, β} (r - R_{B}) = K g_{1 s, γ} (r - R_{p})

K = exp [- \frac{α β}{α + β} ∣ R_{A} - R_{B} ∣^{2}]

γ = α + β

R_{p} = \frac{α R _{A} + β R _{B}}{α + β}

import numpy as np

def gaussian_product_theorem(alpha, RA, beta, RB):
    """
    Compute the Gaussian product theorem parameters.
    
    Parameters:
    -----------
    alpha : float
        Exponent of first Gaussian
    RA : array-like, shape (3,)
        Center of first Gaussian
    beta : float
        Exponent of second Gaussian
    RB : array-like, shape (3,)
        Center of second Gaussian
    
    Returns:
    --------
    K : float
        Pre-exponential factor
    gamma : float
        Combined exponent
    Rp : array, shape (3,)
        Product center
    """
    RA = np.array(RA)
    RB = np.array(RB)
    
    # Distance squared between centers
    Rab2 = np.dot(RA - RB, RA - RB)
    
    # Combined exponent
    gamma = alpha + beta
    
    # Pre-exponential factor
    K = np.exp(-alpha * beta / gamma * Rab2)
    
    # Product center
    Rp = (alpha * RA + beta * RB) / gamma
    
    return K, gamma, Rp

Overlap Integral

⟨ 1 s, α, A ∣1 s, β, B ⟩ = 4 π \int d r r^{2} K e^{γ r^{2}}

= (\frac{π}{α + β})^{3/2} exp [- \frac{α β}{α + β} ∣ R_{A} - R_{B} ∣^{2}]

import numpy as np

def overlap(alpha, RA, beta, RB):
    """
    Compute the overlap integral between two 1s Gaussian-type orbitals.
    
    Parameters:
    -----------
    alpha : float
        Exponent of first Gaussian
    RA : array-like, shape (3,)
        Center of first Gaussian
    beta : float
        Exponent of second Gaussian
    RB : array-like, shape (3,)
        Center of second Gaussian
    
    Returns:
    --------
    S : float
        Overlap integral value
    """
    RA = np.array(RA)
    RB = np.array(RB)
    
    # Distance squared between centers
    Rab2 = np.dot(RA - RB, RA - RB)
    
    # Combined exponent
    p = alpha + beta
    
    # Overlap integral
    S = (np.pi / p) ** (3.0 / 2.0) * np.exp(-alpha * beta / p * Rab2)
    
    return S

Kinetic Energy Integral

⟨ 1 s, α, A ∣ - \nabla^{2} ∣1 s, β, B ⟩

= \frac{α β}{α + β} [6 - \frac{4 α β}{α + β} ∣ R_{A} - R_{B} ∣^{2}] (\frac{π}{α + β})^{3/2} exp [- \frac{α β}{α + β} ∣ R_{A} - R_{B} ∣^{2}]

import numpy as np

def kinetic_energy(alpha, RA, beta, RB):
    """
    Compute the kinetic energy integral between two 1s Gaussian-type orbitals.
    
    Parameters:
    -----------
    alpha : float
        Exponent of first Gaussian
    RA : array-like, shape (3,)
        Center of first Gaussian
    beta : float
        Exponent of second Gaussian
    RB : array-like, shape (3,)
        Center of second Gaussian
    
    Returns:
    --------
    T : float
        Kinetic energy integral value
    """
    RA = np.array(RA)
    RB = np.array(RB)
    
    # Distance squared between centers
    Rab2 = np.dot(RA - RB, RA - RB)
    
    # Combined exponent
    p = alpha + beta
    
    # Reduced mass-like parameter
    mu = alpha * beta / p
    
    # Kinetic energy integral
    T = (
        mu
        * (6.0 - 4.0 * mu * Rab2)
        * (np.pi / p) ** (3.0 / 2.0)
        * np.exp(-mu * Rab2)
    )
    
    return T

Electron Nucleus Coulomb Integral

⟨ 1 s, α, A ∣ - \frac{Z}{r _{N}} ∣1 s, β, B ⟩

= \frac{2 π Z}{α + β} exp [\frac{- α β}{α + β} ∣ R_{A} - R_{B} ∣^{2}] \times F_{0} [(α + β) ∣ R_{P} - R_{N} ∣^{2}]

F_{m} (t) = \int_{0}^{1} exp (- t s^{2}) s^{2 m} d s

F_{0} (t) = \int_{0}^{1} exp (- t s^{2}) d s

erf (x) = \frac{2}{π} \int_{0}^{x} d s exp (- s^{2})

One can perform a change of variables and show that

F_{0} (x) = \frac{π}{4 x} erf (x)

import numpy as np
from scipy.special import erf

def F0(t):
    """
    Compute the F0 function used in electron-nucleus Coulomb integrals.
    
    Parameters:
    -----------
    t : float
        Argument of F0 function
    
    Returns:
    --------
    F0_value : float
        Value of F0 function
    """
    if t < 1e-10:  # Avoid division by zero
        return 1.0
    
    sqrt_t = np.sqrt(t)
    F0_value = np.sqrt(np.pi / (4 * t)) * erf(sqrt_t)
    
    return F0_value

def electron_nucleus_coulomb(alpha, RA, beta, RB, Z, RN):
    """
    Compute the electron-nucleus Coulomb integral.
    
    Parameters:
    -----------
    alpha : float
        Exponent of first Gaussian
    RA : array-like, shape (3,)
        Center of first Gaussian
    beta : float
        Exponent of second Gaussian
    RB : array-like, shape (3,)
        Center of second Gaussian
    Z : float
        Nuclear charge
    RN : array-like, shape (3,)
        Nuclear position
    
    Returns:
    --------
    Vne : float
        Electron-nucleus Coulomb integral value
    """
    RA = np.array(RA)
    RB = np.array(RB)
    RN = np.array(RN)
    
    # Distance squared between Gaussian centers
    Rab2 = np.dot(RA - RB, RA - RB)
    
    # Combined exponent
    p = alpha + beta
    
    # Product center
    Rp = (alpha * RA + beta * RB) / p
    
    # Distance squared from product center to nucleus
    Rpn2 = np.dot(Rp - RN, Rp - RN)
    
    # Argument for F0 function
    t = p * Rpn2
    
    # Electron-nucleus Coulomb integral
    Vne = (
        2 * np.pi * Z / p
        * np.exp(-alpha * beta / p * Rab2)
        * F0(t)
    )
    
    return Vne

Two Electron Coulomb Integral

⟨ 1 s, α, A; 1 s, β, B ∣ \frac{1}{r _{12}} ∣1 s, γ, C; 1 s, δ, D ⟩

= \frac{2 π ^{(5/2)}}{( α + γ ) ( β + δ ) ( α + β + γ + δ ) ^{1/2}}

\times exp [- \frac{α γ}{α + γ} ∣ R_{A} - R_{C} ∣^{2} - \frac{β δ}{β + δ} ∣ R_{B} - R_{D} ∣^{2}] \times F_{0} [\frac{( α + γ ) ( β + δ )}{( α + β + γ + δ )} ∣ R_{P} - R_{Q} ∣^{2}]

R_{p} = \frac{α R _{A} + β R _{B}}{α + β}

R_{q} = \frac{γ R _{C} + δ R _{D}}{γ + δ}

import numpy as np
from scipy.special import erf

def F0(t):
    """
    Compute the F0 function used in two-electron Coulomb integrals.
    
    Parameters:
    -----------
    t : float
        Argument of F0 function
    
    Returns:
    --------
    F0_value : float
        Value of F0 function
    """
    if t < 1e-10:  # Avoid division by zero
        return 1.0
    
    sqrt_t = np.sqrt(t)
    F0_value = np.sqrt(np.pi / (4 * t)) * erf(sqrt_t)
    
    return F0_value

def two_electron_coulomb(alpha, RA, beta, RB, gamma, RC, delta, RD):
    """
    Compute the two-electron Coulomb integral (electron-electron repulsion).
    
    Parameters:
    -----------
    alpha : float
        Exponent of first Gaussian
    RA : array-like, shape (3,)
        Center of first Gaussian
    beta : float
        Exponent of second Gaussian
    RB : array-like, shape (3,)
        Center of second Gaussian
    gamma : float
        Exponent of third Gaussian
    RC : array-like, shape (3,)
        Center of third Gaussian
    delta : float
        Exponent of fourth Gaussian
    RD : array-like, shape (3,)
        Center of fourth Gaussian
    
    Returns:
    --------
    Vee : float
        Two-electron Coulomb integral value
    """
    RA = np.array(RA)
    RB = np.array(RB)
    RC = np.array(RC)
    RD = np.array(RD)
    
    # Combined exponents
    p = alpha + beta
    q = gamma + delta
    pq = p + q
    
    # Product centers
    Rp = (alpha * RA + beta * RB) / p
    Rq = (gamma * RC + delta * RD) / q
    
    # Distance squared between product centers
    Rpq2 = np.dot(Rp - Rq, Rp - Rq)
    
    # Distance squared terms for exponential
    Rac2 = np.dot(RA - RC, RA - RC)
    Rbd2 = np.dot(RB - RD, RB - RD)
    
    # Arguments for exponential and F0
    exp_arg1 = -alpha * gamma / p * Rac2
    exp_arg2 = -beta * delta / q * Rbd2
    t = p * q / pq * Rpq2
    
    # Two-electron Coulomb integral
    Vee = (
        2 * np.pi ** (5.0 / 2.0)
        / (p * q * np.sqrt(pq))
        * np.exp(exp_arg1 + exp_arg2)
        * F0(t)
    )
    
    return Vee

Complete Example

Here's a complete example showing how to use all the integral functions together:

import numpy as np
from scipy.special import erf

def F0(t):
    """F0 function for Coulomb integrals."""
    if t < 1e-10:
        return 1.0
    sqrt_t = np.sqrt(t)
    return np.sqrt(np.pi / (4 * t)) * erf(sqrt_t)

def overlap(alpha, RA, beta, RB):
    """Overlap integral."""
    RA = np.array(RA)
    RB = np.array(RB)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    return (np.pi / p) ** (3.0 / 2.0) * np.exp(-alpha * beta / p * Rab2)

def kinetic_energy(alpha, RA, beta, RB):
    """Kinetic energy integral."""
    RA = np.array(RA)
    RB = np.array(RB)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    mu = alpha * beta / p
    return mu * (6.0 - 4.0 * mu * Rab2) * (np.pi / p) ** (3.0 / 2.0) * np.exp(-mu * Rab2)

def electron_nucleus_coulomb(alpha, RA, beta, RB, Z, RN):
    """Electron-nucleus Coulomb integral."""
    RA = np.array(RA)
    RB = np.array(RB)
    RN = np.array(RN)
    Rab2 = np.dot(RA - RB, RA - RB)
    p = alpha + beta
    Rp = (alpha * RA + beta * RB) / p
    Rpn2 = np.dot(Rp - RN, Rp - RN)
    t = p * Rpn2
    return 2 * np.pi * Z / p * np.exp(-alpha * beta / p * Rab2) * F0(t)

def two_electron_coulomb(alpha, RA, beta, RB, gamma, RC, delta, RD):
    """Two-electron Coulomb integral."""
    RA = np.array(RA)
    RB = np.array(RB)
    RC = np.array(RC)
    RD = np.array(RD)
    p = alpha + beta
    q = gamma + delta
    pq = p + q
    Rp = (alpha * RA + beta * RB) / p
    Rq = (gamma * RC + delta * RD) / q
    Rpq2 = np.dot(Rp - Rq, Rp - Rq)
    Rac2 = np.dot(RA - RC, RA - RC)
    Rbd2 = np.dot(RB - RD, RB - RD)
    exp_arg1 = -alpha * gamma / p * Rac2
    exp_arg2 = -beta * delta / q * Rbd2
    t = p * q / pq * Rpq2
    return 2 * np.pi ** (5.0 / 2.0) / (p * q * np.sqrt(pq)) * np.exp(exp_arg1 + exp_arg2) * F0(t)

# Example: Compute integrals for H2+ molecule
# Two 1s Gaussians centered at different positions
alpha = 0.5
beta = 0.5
RA = np.array([0.0, 0.0, 0.0])
RB = np.array([1.4, 0.0, 0.0])  # Bond length ~1.4 Bohr

# Overlap integral
S = overlap(alpha, RA, beta, RB)
print(f"Overlap integral: {S:.6f}")

# Kinetic energy integral
T = kinetic_energy(alpha, RA, beta, RB)
print(f"Kinetic energy integral: {T:.6f}")

# Electron-nucleus Coulomb (nucleus at origin)
Z = 1.0
RN = np.array([0.0, 0.0, 0.0])
Vne = electron_nucleus_coulomb(alpha, RA, beta, RB, Z, RN)
print(f"Electron-nucleus Coulomb: {Vne:.6f}")

# Two-electron Coulomb (for H2, not H2+)
# Using same Gaussians for both electron pairs
Vee = two_electron_coulomb(alpha, RA, beta, RB, alpha, RA, beta, RB)
print(f"Two-electron Coulomb: {Vee:.6f}")