H2+ Molecular Orbital (C++)

Quantum Chemistry

This page presents a C++ implementation for computing the H2+ molecular orbital (gerade/ungerade) and related matrix elements using numerical integration on a 3D grid.

H2+ Molecular Orbital

For H2+, the bonding (gerade) orbital is constructed as a linear combination of atomic orbitals:

ψ_{g} (r) = N_{g} [ϕ_{1} (r - R_{a}) + ϕ_{1} (r - R_{b})]

where $ϕ_{1}$ is a 1s hydrogen orbital and the normalization constant is:

N_{g} = \frac{1}{2 ( 1 + S )}

with overlap integral $S = ⟨ ϕ_{1} (r - R_{a}) ∣ ϕ_{1} (r - R_{b})⟩$ .

Matrix Elements

The Hamiltonian matrix elements are computed numerically on a 3D grid:

Kinetic Energy: $T = - \frac{1}{2} \int ψ \nabla^{2} ψ d V$
Potential Energy: $V = \int ψ (- \frac{1}{r _{a}} - \frac{1}{r _{b}} + \frac{1}{R}) ψ d V$
Total Energy: $E = T + V$

C++ Implementation

The following code implements the H2+ orbital and matrix element calculations with all helper functions included:

#include <vector>
#include <functional>
#include <cmath>
#include <iostream>
#include <numbers>
#include <algorithm>

// Helper function: linspace
std::vector<double> linspace(double start, double end, size_t num) {
    std::vector<double> result(num);
    double dx = (end - start) / (num - 1);
    for (size_t i = 0; i < num; ++i) {
        result[i] = start + i * dx;
    }
    return result;
}

// Euclidean distance
double distance(const std::vector<double>& x, const std::vector<double>& xp) {
    double result = 0.0;
    for (size_t i = 0; i < x.size(); i++) {
        result += std::pow(x[i] - xp[i], 2);
    }
    return std::sqrt(result);
}

// Kahan summation for numerical stability
enum class SumType {
    Naive,
    Kahan,
    Neumaier
};

class Summation {
public:
    explicit Summation(SumType t) : type(t), sum(0.0), c(0.0) {}

    void add(double x) {
        if (type == SumType::Kahan) {
            double y = x - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        } else {
            sum += x;
        }
    }

    double value() const {
        return sum + c;
    }

private:
    SumType type;
    double sum;
    double c;
};

// Hydrogen 1s orbital
double hydrogenOrbital(double r) {
    double a0 = 1.0;  // Bohr radius in atomic units
    return std::exp(-r/a0);
}

// H2+ gerade orbital class
struct h2orbital {
    const std::vector<double>& a;  // Position of nucleus A
    const std::vector<double>& b;  // Position of nucleus B

    h2orbital(const std::vector<double>& a_,
              const std::vector<double>& b_)
        : a(a_), b(b_) {}

    double operator()(double x, double y, double z) {
        std::vector<double> r = {x, y, z};
        double pi = std::numbers::pi_v<double>; 
        double R = distance(a, b);
        
        // Overlap integral (analytical form for 1s orbitals)
        double S = (1.0 + R + (1.0/3.0)*R*R) * std::exp(-R);
        
        // Normalization constant for gerade orbital
        double N_g = 1.0 / std::sqrt(2.0 * (1.0 + S));
        
        // 1s orbital normalization
        double oneSNorm = 1.0 / std::sqrt(pi);
        
        // Construct gerade orbital
        double result = oneSNorm * N_g * 
                       (hydrogenOrbital(distance(r, a)) + 
                        hydrogenOrbital(distance(r, b)));
        return result;
    }
};

// 3D Laplacian using finite differences
double templaplace3dcart(double x, double y, double z, 
                         std::function<double(double, double, double)> f) {
    double h = 0.001;
    double result = f(x + h, y, z) + f(x - h, y, z) 
                  + f(x, y + h, z) + f(x, y - h, z) 
                  + f(x, y, z + h) + f(x, y, z - h) 
                  - 6.0 * f(x, y, z);
    result = result / (h * h);
    return result;
}

// H2+ Hamiltonian calculation
void computeH2Hamiltonian(double R, double box, int points) {
    std::vector<double> R_a = {0, 0, 0.5*R};
    std::vector<double> R_b = {0, 0, -0.5*R};
    
    auto xs = linspace(-box, box, points);
    auto ys = linspace(-box, box, points);
    auto zs = linspace(-box, box, points);
    
    double dx = std::abs(xs[1] - xs[0]); 
    double dy = std::abs(ys[1] - ys[0]); 
    double dz = std::abs(zs[1] - zs[0]); 
    
    h2orbital orbital(R_a, R_b);
    std::vector<double> r = {0, 0, 0};
    
    Summation kahan_k(SumType::Kahan);
    Summation kahan_v(SumType::Kahan);
    
    for (size_t i = 0; i < xs.size(); i++) {
        for (size_t j = 0; j < ys.size(); j++) {
            for (size_t k = 0; k < zs.size(); k++) {
                r[0] = xs[i];
                r[1] = ys[j];
                r[2] = zs[k];
                
                double ra = distance(r, R_a);
                double rb = distance(r, R_b);
                
                double orb = orbital(xs[i], ys[j], zs[k]);
                double laplace = templaplace3dcart(xs[i], ys[j], zs[k], 
                                                   [&](double x, double y, double z) {
                                                       return orbital(x, y, z);
                                                   });
                
                // Kinetic energy density
                double T_density = (-0.5 * laplace) * orb * dx * dy * dz;
                kahan_k.add(T_density);
                
                // Potential energy density
                double V_density = (-1.0/std::max(rb, 1e-3) + 
                                   -1.0/std::max(ra, 1e-3) + 
                                   1.0/R) * orb * orb * dx * dy * dz;
                kahan_v.add(V_density);
            }
        }
    }
    
    double kinetic = kahan_k.value();
    double potential = kahan_v.value();
    double total_energy = kinetic + potential;
    
    std::cout << "Kinetic energy: " << kinetic << std::endl;
    std::cout << "Potential energy: " << potential << std::endl;
    std::cout << "Total energy: " << total_energy << std::endl;
}

// Example usage
int main() {
    double R = 2.0;  // Bond length in atomic units
    double box = 4.0;
    int points = 50;  // Reduced for faster execution
    
    computeH2Hamiltonian(R, box, points);
    
    return 0;
}

Key Features

Functor-based orbital class for easy evaluation
3D finite difference Laplacian
Kahan summation for stable integration
Handles singularities at nuclear positions

Numerical Considerations

Grid spacing must be small enough to resolve orbital features
Box size must be large enough to capture orbital decay
Singularities at $r_{a} = 0$ and $r_{b} = 0$ require careful handling
Kahan summation essential for accurate 3D integration