MOPAC_RS

Modern High-Performance Data-Oriented Semi-Empirical Quantum Chemistry Engine in Rust

Executive Summary

mopac_rs is a ground-up architectural reimagining and rigorous modernization of the classic MOPAC (Molecular Orbital PACkage) quantum chemistry engine. Translated from legacy Fortran into idiomatic, zero-overhead Rust, it replaces decades of non-contiguous global arrays, static buffers, and triangular packing with a strictly Data-Oriented Programming (DOP) architecture.

Every single module, parameter table, and integral calculation is empirically verified against OpenMOPAC v23.2.5 references to rigorous double-precision tolerances.

Architectural Pillars

• Data-Oriented Memory Layout (SoA): Contiguous, 64-byte cache-line aligned Struct of Arrays (MolecularBatch) eliminating pointer-chasing and non-contiguous matrix indexing.
• Zero-Allocation Inner Loop Policy (0 malloc): Pre-allocated reusable workspaces (ScfWorkspace, UhfWorkspace, GradientWorkspace, CosmoState, IrcWorkspace, DrcWorkspace, MeciWorkspace, PbcWorkspace) ensure zero heap allocations during iterative Roothaan-Hall / Pople-Nesbet SCF cycles and dual Pulay DIIS extrapolations.
• Dual Compute Backend:
  • CPU SIMD: Vectorized AVX2 / FMA kernels with Rayon multi-threaded parallelism.
  • Universal GPU (Vulkan Compute): Cross-vendor hardware acceleration supporting NVIDIA RTX, AMD Radeon, Intel Arc, and Apple Silicon (via MoltenVK) with dedicated DMA host-to-device GDDR6 VRAM batch management.
• Axiomatic Verification: 89 automated scrutiny, unit, quantum theorem, and differential oracle tests validating physical invariance, rotation orthonormality, translational symmetry, and exact numerical parity against OpenMOPAC.

Scientific Capabilities

1. Semi-Empirical Hamiltonians

• MNDO (Modified Neglect of Diatomic Overlap; Dewar & Thiel 1977)
• AM1 (Austin Model 1; Dewar et al. 1985)
• PM3 (Parametric Method 3; Stewart 1989)
• RM1 (Recife Model 1; Rocha et al. 2006)
• PM6 (Parametric Method 6; Stewart 2007)
• PM7 (Parametric Method 7; Stewart 2013) with full diatomic parameters and d-orbital polarization for transition metals (e.g. Zn).
• NDDO 22-Multipole Integrals: Full diatomic charge separation multipoles ($dd, qq, am, ad, aq$) and 3D rotational coordinate transformations.
• Elemental Coverage: Complete authentic parameter sets for organic, heteroatomic, boron, organosilicon, and metal chemistry:
  • Hydrogen & Carbon-backbone: H (1), C (6)
  • Boron (Group 13): B (5) fully parameterized across MNDO, AM1, PM3, PM6, PM7 (including authentic diatomic pair parameters with H, C, F, Cl, Br, I) and isolated atom heats of formation ($135.700\text{ kcal/mol}$).
  • Pnictogens & Chalcogens: N (7), O (8), P (15), S (16)
  • Full Halogen Series: F (9), Cl (17), Br (35), I (53) across AM1, PM6, PM7, and RM1 with authentic diatomic pair parameters $(alpb, xfac)$.
  • Silicon (Group 14): Si (14) across MNDO, AM1, PM3, PM6, PM7 (including 11 authentic pair parameters with H, C, N, O, F, Al, Si, P, S, Cl, Br) exhibiting $< 2.4 \times 10^{-5}\text{ eV/\AA}$ analytical gradient parity.
  • Transition Metals: Zn (30) with $spd$ basis parameterization.

2. MOZYME $O(N)$ Linear Scaling Macromolecular Solver

• Lewis Chemical Topology Builder: Automatic bond order detection (single, double, triple, coordinate) and lone pair assignment.
• Directional Hybrid Atomic Orbitals (HAOs): $sp^3, sp^2, sp$ hybrid generation with symmetric Löwdin orthogonalization ($S^{-1/2}$).
• Localized Molecular Orbitals (LMOs): Orthogonal $\sigma, \pi$ bonding, $\sigma^, \pi^$ antibonding, and lone-pair initial basis.
• Pairwise 2x2 Jacobi Rotation Engine: Iterative localized orbital energy minimization with spatial distance cutoff screening ($R_{\text{cut}} = 8.5\text{ \AA}$).
• Linear-Scaling SCF Cycle: Smooth density damping ($\alpha = 0.65$) and exact density matrix reconstruction ($P = 2 \sum_i \phi_i \phi_i^T$) achieving $< 0.08%$ relative energy parity on macromolecules.

3. Open-Shell Unrestricted Hartree-Fock (UHF)

• Pople-Nesbet Spin Orbitals: Independent spin Fock operators ($F^\alpha, F^\beta$) and density matrices ($P^\alpha, P^\beta$).
• Dual DIIS Acceleration: Decoupled alpha and beta error subspace inversion buffers avoiding inter-spin oscillation damping traps.
• Spin Expectation $\langle S^2 \rangle$: Analytic computation of total spin angular momentum with spin contamination tracking ($\Delta \langle S^2 \rangle = 0.000178$ against OpenMOPAC on methyl radical).
• Closed-Shell Equivalence: Exact convergence to RHF energy for singlet states ($|\Delta E| < 10^{-10}\text{ eV}$).

4. Robust SCF Convergers

• Pulay DIIS Acceleration: Direct Inversion in the Iterative Subspace with B-matrix SVD stabilization and history pruning.
• Camp-King Unitary Interpolator: Monotonic electronic energy minimization for oscillating densities.
• Saunders-Hillier Virtual Orbital Level Shifting: Dynamic shift ($\sigma = 2.0\text{ -- }8.0\text{ eV}$) eliminating HOMO-LUMO degeneracy traps and limit-cycle density oscillations.
• Adaptive Multi-Tier Escalation: Automated converger pipeline achieving 100.0% convergence across all benchmark sets.

5. Non-Covalent Corrections

• Grimme D3-BJ Empirical Dispersion: Becke-Johnson rational damping ($s_6, s_8, a_1, a_2$) with exact analytical Cartesian gradients matching finite differences to $1.11 \times 10^{-11}\text{ kcal/(mol}\cdot\text{\AA)}$.
• PM6-DH+ and PM7 Dispersion: Empirically parameterized dispersion corrections with analytical derivatives ($\nabla E_{\text{disp}} < 1.88 \times 10^{-11}\text{ kcal}/(\text{mol}\cdot\text{\AA})$ error).
• H4 Hydrogen Bonding: Septic switching functions for covalent valence attenuation and 7th-order radial/angular polynomials ($D, A \in {N, O}$).
• Short-Range H-H Repulsion: Continuous piecewise potential with exact analytical derivatives.
• Composite PM6-D3H4 Method: Verbatim heat of formation parity on water dimer benchmark ($-71.99024\text{ kcal/mol}$).

6. COSMO Implicit Solvation & Analytical Nuclear Gradients

• Boundary Element Method (BEM): Regular icosahedral sphere tessellations ($N=12, 42, 1082$, dvfill).
• Solvent-Accessible Cavity: Klamt and Bondi van der Waals radii with analytical segment surface areas and volumes.
• Self-Consistent Reaction Field: In-place Cholesky decomposition of electrostatic boundary matrix $A$ and multipole coupling matrix $B$, modifying $H_{\text{core}}$ and Fock matrix $F$ self-consistently with 0 malloc per iteration.
• Analytical Nuclear Gradients ($\nabla E_{\text{diel}}$): Inter-segment screening forces and segment-solute net charge electrostatic derivatives matching OpenMOPAC diegrd with exact Newton's third law translational zero-sum invariance ($\sum_A \nabla_A E_{\text{diel}} < 10^{-14}\text{ eV/\AA}$).

7. Transition States & Reaction Paths

• Eigenvector Following (Baker P-RFO / TS): Search for first-order saddle points (transition states) with one negative Hessian eigenvalue.
• Two-Ended SADDLE Interpolation: Automated barrier search connecting Reactants and Products via geodesic coordinate relaxation.
• Intrinsic Reaction Coordinate (IRC / González-Schlegel): Exact mass-weighted steepest descent reaction path tracing forward and backward from transition state to minima.
• Dynamic Reaction Coordinate (DRC / Born-Oppenheimer MD): Microcanonical ($NVE$) and canonical ($NVT$) molecular dynamics trajectories with symplectic Velocity-Verlet integrator.

8. Photochemistry, MECI & UV-Vis Spectroscopy

• Multi-Electron Configuration Interaction (MECI): Full CI within active spaces ($N \le 10$ orbitals) computing multi-determinant ground and excited state roots ($S_0, S_1, T_1$, etc.).
• Analytical Excited-State Nuclear Gradients: State-specific Hellmann-Feynman and relaxed density matrix gradients for excited state geometry optimization.
• UV-Vis Spectral Simulation: Transition dipole moments, oscillator strengths ($f_{\text{osc}}$), and Lorentzian line broadening.

9. Molecular Properties & Population Analysis

• Electric Dipole Moments: Point-charge and intra-atomic $sp$ hybridization dipole moments in Debye.
• Polarizability & NLO Tensors (TD-CPHF / POLAR): Static and dynamic frequency-dependent polarizability tensor $\alpha(-\omega; \omega)$ and first hyperpolarizability $\beta$.
• Electrostatic Potential (ESP) Charges: Merz-Singh-Kollman grid fitting outside van der Waals envelope.
• Mayer Bond Orders & Valencies: Armstrong-Perkins-Stewart bond indices $B_{AB} = \sum_{\mu \in A, \nu \in B} (P S){\mu\nu} (P S){\nu\mu}$.
• Mulliken Population Analysis: Löwdin de-orthogonalization and gross atomic populations satisfying exact electron conservation ($\sum_A Pop_A \equiv N_{\text{elec}}$).

10. Geometry Optimization, Vibrations & Isotope Effects

• L-BFGS Optimizer: Quasi-Newton Cartesian minimization with two-loop history recursion and Armijo backtracking line search.
• Coordinate Pinning: Selective degree-of-freedom masking (frozen atoms/axes).
• Harmonic Vibrational Frequencies: Mass-weighted Cartesian Hessian with Eckart frame external projection (6 vanishing rotational/translational modes $< 10^{-5}\text{ cm}^{-1}$).
• Custom Isotopic Masses & Kinetic Isotope Effects (KIE): Direct mass substitution ($^2H, ^{13}C, ^{18}O$) for vibrational isotope shift analysis.
• Thermodynamic Properties: Zero-Point Vibrational Energy (ZPVE), thermal enthalpy ($H(T) - H(0)$), constant-pressure heat capacity ($C_p$), standard entropy ($S^\circ$), and Gibbs free energy correction ($G(T) - H(0)$).

11. Periodic Boundary Conditions (PBC)

• 1D, 2D, 3D Unit Cells: Lattice vector parameterization with Monkhorst-Pack reciprocal space k-point sampling and Bloch SCF band structures.

12. AM1-BCC Atomic Partial Charges

• AMBER / Open Force Field Parity: Clean-room implementation of the Jakalian et al. (2000, 2002) bond charge correction scheme directly atop semi-empirical AM1 Mulliken populations.
• Empirical BCC Parameterization: Includes standard bond charge corrections (BCC01–BCC28) covering aliphatic/aromatic C-H, carbonyl C=O, alcohols/ethers C-O, amines/imines C-N, and halogens C-X.
• Strict Charge Conservation: Enforces antisymmetric bond transfers $\delta_{ij} = -\delta_{ji}$, guaranteeing $\sum_i \Delta q_i \equiv 0$ strictly to machine precision ($10^{-16}$) and preserving total molecular charge.
• Direct CLI & Python Availability: Accessible via the --bcc / --am1-bcc CLI flag, input keyword AM1-BCC, and mopac_py.am1_bcc(atoms, coords).

13. QM/MM Electrostatic Embedding

• Classical Point Charge Coupling: Fully self-consistent electrostatic embedding of external point charges $Q_k$ at coordinates $\vec{R}k$ into the quantum core Hamiltonian: $$H{\mu\mu} \leftarrow H_{\mu\mu} - \sum_{k} \frac{Q_k}{\sqrt{R_{Ak}^2 + \rho_A^2}}$$
• Dewar-Klopman Damping: Regularized monopole interaction radius $\rho_A = \frac{\text{EV_ANGSTROM_FACTOR}}{2 g_{ss}^A}$ eliminating non-physical electron collapse / spillover catastrophes while converging exactly to Coulomb's law at long range ($R_{Ak} \gg \rho_A$).
• Nuclear Repulsion & Polarization: Includes external core-charge energy $E_{\text{core-ext}} = \sum_A \sum_k \frac{Z_A Q_k}{\sqrt{R_{Ak}^2 + \rho_A^2}}$ and captures explicit solvent polarization (e.g. water dipole moment increases from $1.854\text{ D}$ to $2.049\text{ D}$ in the presence of an external test charge).
• Exact Analytical Nuclear Gradients: Closed-form nuclear gradients $\nabla_A E_{\text{ext}}$ verified against central finite differences to within $< 10^{-5}\text{ eV/\AA}$ with zero translational drift.

14. Python Bindings & Ecosystem Bridge (mopac_py)

High-performance PyO3 bridge exposing the quantum chemical engine directly to Python for RDKit, ASE (Atomic Simulation Environment), and PyTorch / TorchMD:

• Zero-copy analytical gradients directly returned as NumPy arrays or torch tensors.
• Vibrational Frequencies & Thermochemistry: mopac_py.frequencies(atoms, coords, method="AM1") computing full normal modes, ZPVE, $H(T)$, $G(T)$, $C_p$, and $S^\circ$.
• AM1-BCC Charge Calculator: mopac_py.am1_bcc(atoms, coords) returning raw Mulliken and refined BCC partial charges.
• MOZYME Linear Scaling: mopac_py.mozyme(atoms, coords, method="PM6") for macromolecules.
• Direct RDKit Interoperability: mopac_py.from_rdkit(mol) extracts atomic numbers and 3D conformer coordinates directly.
• Native ASE Calculator: mopac_py.MopacASECalculator plugs directly into ASE dynamics (BFGS, VelocityVerlet, etc.).
• Differentiable PyTorch Potential (TorchMD / $\Delta$-ML): mopac_py.MopacPotential as a native torch.nn.Module with custom autograd function (mopac_py.MopacEnergyFunction) for hybrid QM/MM and machine-learned force field training.
• Full parameterization control (method="PM6", dispersion="D3-BJ", cosmo_eps=78.4, use_nddo=True).

import mopac_py

# Water single-point calculation (PM6 + COSMO solvation + Grimme D3-BJ)
atoms = [8, 1, 1]
coords = [[0.0, 0.0, 0.0655], [0.0, 0.7571, -0.5205], [0.0, -0.7571, -0.5205]]
res = mopac_py.calculate(atoms, coords, method="PM6", dispersion="D3-BJ", cosmo_eps=78.4)

print(f"Total Energy: {res.total_energy_ev:.6f} eV")
print(f"Heat of Formation: {res.heat_of_formation_kcal:.3f} kcal/mol")
print(f"Dipole: {res.dipole_debye[3]:.3f} Debye")
print(f"Mulliken Charges: {res.mulliken_charges}")
print(f"Gradients (eV/A): {res.gradients_ev_angstrom}")

# AM1-BCC atomic partial charges
bcc = mopac_py.am1_bcc(atoms, coords)
print(f"AM1-BCC Charges: {bcc.bcc_charges}")

# Vibrational frequencies and thermochemistry (298.15 K, 1 atm)
vib = mopac_py.frequencies(atoms, coords, method="AM1")
print(f"Vibrational Frequencies (cm^-1): {vib.vibrational_frequencies_cm1}")
print(f"ZPVE: {vib.zpve_kcal_mol:.3f} kcal/mol")
print(f"Standard Entropy S^o: {vib.thermo.entropy_total_cal_k_mol:.3f} cal/(mol*K)")

# Geometry optimization (L-BFGS)
opt = mopac_py.optimize(atoms, coords, method="PM6", max_cycles=50)
print(f"Optimized Energy: {opt.final_energy_ev:.6f} eV (Converged: {opt.converged})")

15. UCSF ChimeraX Extension (integrations/chimerax)

Interactive visual molecular modeling plugin for UCSF ChimeraX 1.8+:

• Seamless interactive semi-empirical calculations inside the 3D viewport.
• Surface electrostatic coloring from AM1-BCC / Mulliken charges.
• One-click L-BFGS geometry optimization and MOZYME macromolecular calculations via mopac calculate, mopac bcc, mopac optimize, and mopac mozyme.

Hardware Execution & GPU Dispatch Policy

A key architectural insight in quantum chemistry is the trade-off between SIMD latency and discrete GPU dispatch overhead:

• Small-to-Medium Molecules ($N < 150$ atoms):
  • Evaluated on the CPU SIMD (AVX2 / FMA + Rayon) backend.
  • With zero PCIe/DMA overhead and sub-millisecond execution times ($< 1\text{ -- }10\text{ ms}$ per SCF cycle), CPU execution strictly outperforms GPU dispatch for single molecules.
• Macromolecules & Batched Screening ($N > 200$ atoms or $M \ge 300$ conformations):
  • Evaluated on the Universal Vulkan GPU backend.
  • Pipelined host-to-device asynchronous GDDR6 memory transfers saturate thousands of FP32/FP64 GPU shader compute cores, amortizing Vulkan command buffer submission latency and delivering order-of-magnitude throughput scaling for high-throughput virtual screening (HTVS) and macromolecular assemblies.

CLI Usage

Build and Install

# Clone the repository
git clone https://github.com/zephsystems/mopac_rs.git
cd mopac_rs

# Build optimized release binary
cargo build --release --bin mopac

# Run entire test suite
cargo test --workspace

Command-Line Arguments

Usage: mopac [OPTIONS] <INPUT>

Arguments:
  <INPUT>  Input file path (.mop)

Options:
  -m, --mode <MODE>      Force calculation mode (1SCF or OPT)
      --method <METHOD>  Semi-empirical method override (PM6, PM7, PM3, RM1, AM1, MNDO)
      --nddo             Enable full NDDO diatomic 22-multipole integrals & 3D rotation
      --opt              Enable geometry optimization (L-BFGS)
      --ts               Enable transition state optimization via Eigenvector Following (P-RFO Baker)
      --irc              Enable Intrinsic Reaction Coordinate path tracing (González-Schlegel)
      --drc              Enable Dynamic Reaction Coordinate molecular dynamics (Velocity-Verlet)
      --force            Enable Cartesian Hessian & vibrational frequency analysis
      --bonds            Enable Mayer bond orders and atomic valencies calculation
      --mullik           Enable Mulliken population analysis
      --bcc, --am1-bcc   Compute AM1-BCC atomic partial charges
      --ci <N>           Enable Multi-Electron Configuration Interaction (MECI) with active space size N
      --uv-vis           Simulate UV-Vis electronic absorption spectrum
      --static, --polar  Calculate finite-field polarizability and NLO hyperpolarizability tensors
      --mozyme           Enable MOZYME localized molecular orbital linear-scaling SCF (O(N))
      --pbc              Enable Periodic Boundary Conditions (PBC) Bloch SCF & band structure
      --gpu              Enable Vulkan GPU compute acceleration
      --fp32             Use FP32 single-precision GPU pipeline
      --eps <EPS>        Solvent dielectric constant for COSMO implicit solvation (e.g. 78.4)
      --disp <DISP>      Empirical dispersion model (pm6-dh+, pm7)
      --d3h4             Enable D3H4 composite correction (dispersion, H4, H-H repulsion)
      --1scf             Force single-point calculation (1SCF)
      --threads <N>      Number of Rayon worker threads
  -o, --output <FILE>    Custom output report file path (.out)
      --arc <FILE>       Custom archive file path (.arc)
  -h, --help             Print help
  -V, --version          Print version

Sample Input File (water.mop)

PM6 EPS=78.4 BONDS MULLIK 1SCF
Water in aqueous solution COSMO implicit solvation test

O   0.000000 0   0.000000 0   0.000000 0
H   0.757000 0   0.586000 0   0.000000 0
H  -0.757000 0   0.586000 0   0.000000 0

Execute:

./target/release/mopac water.mop

Outputs generated:

• water.out: Comprehensive human-readable report with energetic breakdown, dipole contributions, charges, and bond orders.
• water.arc: Standard archive containing converged geometry and final heat of formation.

Physical Invariants & Quantum Theorems

Empirical verification of fundamental quantum mechanical and systems invariants in crates/mopac_core/tests/quantum_theorems.rs:

Canonical Oracle Differential Parity (vs OpenMOPAC v23.2.5)

Direct, double-precision differential validation executed by running the compiled native OpenMOPAC v23.2.5 binary side-by-side against mopac_rs (crates/mopac_core/tests/golden_parity.rs):

Criterion Microbenchmark Suite

Empirically measured single-thread performance benchmarks (crates/mopac_core/benches/quantum_benchmarks.rs):

Technical Documentation

Detailed mathematical derivations, Fortran audits, and GPU specifications:

• Translation Manifesto
• Fortran Codebase Audit & DOP Feasibility
• Universal GPU Acceleration (Vulkan Compute / WGPU)
• Strict Testing & Verification Policy
• Translation Devlog & Physical Constants Audit
• Vulkan GPU Acceleration & Level Shifting

Citation

If you use mopac_rs in your academic research, benchmarks, or software, please cite it as:

@software{mopac_rs_2026,
  author = {zeph.sys},
  title = {{MOPAC\_RS: Modern High-Performance Semi-Empirical Quantum Chemistry Engine in Rust}},
  year = {2026},
  publisher = {Zenodo},
  version = {0.1.1},
  doi = {10.5281/zenodo.22731253},
  url = {https://github.com/zephsystems/mopac_rs}
}

Direct metadata is also available in CITATION.cff.

License

Distributed under the Apache License Version 2.0 (Apache-2.0). See LICENSE for complete details.