pm6-rs

Rust-native implementation of the PM6 semiempirical NDDO method (J. J. P. Stewart, J. Mol. Model. 13, 1173 (2007)). Heats of formation, fully analytic nuclear gradients and Hessians, L-BFGS geometry optimization, P-RFO transition-state search, and the PM6-D3 / PM6-D3H4 / PM6-D3H4X post-SCF corrections, with Rust, Python-native, ASE and command-line interfaces.

• Author: ss0832
• License: GPL-3.0-or-later
• Pure Rust linear algebra (faer + rayon); no LAPACK/BLAS.

Scope and validation status

Validated against the MOPAC v23.2.5 oracle (openmopac/mopac, Apache-2.0); parameters are extracted from the same release.

Capability	Status
SCF energy (RHF + UHF), heats of formation	✅ frozen molecular set; s/p ≤1e-6, hypervalent/TM ≤1e-4 kcal/mol vs MOPAC; see limitations below
d-orbital two-center integrals (S, P, Cl, transition metals, …)	✅ e.g. TiCl₄ ≤2e-5, HCl ≤1e-8 kcal/mol
Analytic gradient (dual-number, Hellmann–Feynman)	✅ s/p and d, matches MOPAC & FD
Analytic Hessian + frequencies	✅ closed-shell s/p/d CPHF; open-shell d uses FD-of-analytic-gradient
L-BFGS geometry optimization	✅ all elements, matches MOPAC minima
Mulliken charges	✅ (≤1e-4 e across s/p and d)
Dipole (full vector, s/p and d)	✅ vs MOPAC ≤3e-4 D; the p–d channel was missing before 0.2.1 and cost H₂S 1.39 D
External electric field: energy, analytic gradient, analytic Hessian	✅ RHF + UHF + DC, vs MOPAC FIELD= both signs; periodic refused
Dipole derivatives (atomic polar tensor), IR intensities	✅ three independent routes agree to 4e-8 e
Molecular orbitals, molden export, first-order orbital response	✅ round-tripped; the expensive products are lazy
SCF stability analysis (Seeger–Pople: singlet-real, triplet-real, UHF internal) + instability removal	✅ all three operators vs FD of the energy ≤1e-6; square H₄ saddle → minimum 0.58 eV lower; stretched H₂ singlet-stable / triplet-unstable, UHF switch 3.56 eV lower. Off by default; real→complex not implemented
⟨S²⟩ and spin contamination (UHF)	✅ CH₃· 0.7541 vs exact 0.75; broken-symmetry H₂ 0.9889 vs exact 0
Lanthanide sparkles (Ce…Yb)	⚠️ functional (0-orbital cores, oracle-derived ΔH_f); energy ~10 kcal/mol from MOPAC
D3 dispersion (PM6-D3)	✅ matches MOPAC ≤1e-5 kcal/mol (incl. 102 atoms)
H4 / X corrections (PM6-D3H4 / -D3H4X)	✅ matches MOPAC ≤1e-5 kcal/mol; see the H4 gate note below
PM6-D3/H4/X under PBC (1D/2D/3D lattice sums)	✅ supercell-exact; energy, gradient, stress, Hessian
Python native + ASE APIs (energy/grad/forces/Hessian/opt/freq)	✅
PBC — periodic SCF at Γ and k, RHF + UHF, 1D/2D/3D, s/p and d	✅ supercell Γ ⇔ k mesh ≤1e-8 eV/atom
PBC — analytic gradient and analytic stress (6 Voigt components)	✅ vs E(ε) finite difference ≤1e-6
PBC — analytic Hessian (periodic CPHF at Γ; complex CPHF at k, q = 0)	✅ vs FD of analytic forces ≤1e-6
PBC — phonons at any q of the k-mesh, band structure, ZPE, free energy	✅ supercell force constants; ASR imposed
PBC — perturbation theory (DFPT) at any q of the k-mesh, primitive cell	✅ vs 2× supercell spectrum 0.0018 cm⁻¹; LO–TO dipole term pinned
PBC — charged cells (background / Makov–Payne), 1D/2D/3D	✅ NaCl Madelung matches; E(L) follows q²α_M/2L
Divide-and-conquer (molecular, PBC Γ, PBC k): energy, gradient, stress	✅ exponent 1.13 vs dense 1.81; crossover ~98 atoms

✅ = validated to the stated tolerance; ⚠️ = implemented and functional, not yet bit-exact. Because MOPAC has no PM6 periodic path, the periodic code is validated against internal identities — supercell Γ ⇔ k mesh, analytic ⇔ finite difference, D(q) ⇔ supercell Γ — each a statement that two independently-coded routes must agree. See docs/scope.md for the full table of capabilities and limitations, and THIRD_PARTY_NOTICES.md.

The H4 covalent gate is discontinuous

PM6-D3H4 selects a donor–H–acceptor triple only while the donor–H bond is covalent (rDH < 1.15 Å), and drops the whole triple the instant that cutoff is crossed. That is a genuine jump discontinuity in the potential energy surface — 2.3–3.7 kcal/mol on an ammonium/water pair, a water dimer and a formic acid dimer — and MOPAC behaves identically, so it belongs to the H4 model rather than to this implementation. Where the energy jumps the gradient contains a delta function and the Hessian does not exist, but automatic differentiation cannot see a branch it did not take and returns a plausible number regardless.

• Default (H4Smoothing::Mopac): bit-exact MOPAC energies, discontinuous surface.
• H4Smoothing::C2 / CLI --h4-smoothing c2: tapers the gate over the last 0.05 Å with the C² switch the model already uses elsewhere, giving a differentiable surface for optimization, MD, Hessians and phonons. Outside that window the energy is bit-identical to MOPAC.

pm6_rs::h4_discontinuity_warning reports proximity to the gate; the CLI warns automatically for gradient, optimize, saddle, frequencies and hessian.

Units

Internal computation uses eV and Bohr with MOPAC's 2018-CODATA model constants (src/constants.rs). At the boundaries:

• Rust and Python-native APIs return atomic units (Hartree, Bohr), plus convenience eV fields and ΔH_f in kcal/mol.
• The ASE calculator uses eV / Å (PM6 is natively eV, so no round-trip).

Build

cargo build --release        # library + pm6_rs_cli
cargo test                   # unit + MOPAC-oracle regression tests

CLI

pm6_rs_cli energy      water.xyz
pm6_rs_cli gradient    water.xyz
pm6_rs_cli optimize    water.xyz            # writes water.pm6opt.xyz
pm6_rs_cli frequencies water.pm6opt.xyz
pm6_rs_cli charges     water.xyz --charge 0 --multiplicity 1
pm6_rs_cli energy      dimer.xyz --method PM6-D3H4X   # PM6 | PM6-D3 | PM6-D3H4 | PM6-D3H4X

pm6_rs_cli orbitals    water.xyz --coefficients --window 3
pm6_rs_cli molden      h2s.xyz                       # [GTO]; d elements included
pm6_rs_cli molden      water.xyz --molden-basis sto  # the exact Slater primitives

pip install pm6-rs-python installs the same interface as pm6-rs. Periodic runs read the cell from an extended-XYZ Lattice="..." header, or from --cell:

pm6-rs energy   examples/polyethylene.xyz --kpts 1,1,8
pm6-rs stress   examples/polyethylene.xyz --kpts 1,1,8
pm6-rs bands    examples/polyethylene.xyz --kpts 1,1,8 --path "0,0,0;0,0,0.5"
pm6-rs phonons  examples/polyethylene.xyz --kpts 1,1,8 --supercell 1,1,7
pm6-rs phonons-dfpt examples/polyethylene.xyz --kpts 1,1,8 --qpoints "0,0,0"
pm6-rs response examples/nacl.xyz --properties born,dielectric
pm6-rs energy   big.xyz --dc --dc-buffer 12.0        # divide-and-conquer

Relaxing a cell (new in 0.2.6) — ions, lattice, or both:

pm6-rs optimize examples/polyethylene.xyz --kpts 1,1,8 --relax ions
pm6-rs optimize examples/polyethylene.xyz --kpts 1,1,8 \
       --relax all --cell-mask 0,0,1,0,0,0     # the c axis only

Which axes are periodic can be said in whichever way reads best, and a single --kpts N means N points per periodic axis:

pm6-rs energy chain.xyz --cell "20 20 2.55" --pbc z       --kpts 8
pm6-rs energy chain.xyz --cell "20 20 2.55" --pbc FFT     --kpts 1,1,8
pm6-rs energy slab.xyz  --cell "4 4 30"     --pbc-z false --kpts 6

Every SCF control the library has is now reachable too — --max-scf, --e-tol, --p-tol, --no-diis, --level-shift, --damping, --h4-smoothing, --projected, --strict, --sanity.

Rust API

use pm6_rs::{Molecule, Pm6Parameters, Pm6Options, run_pm6};

let mol = Molecule::from_xyz_file("water.xyz", 0.0)?;
let params = Pm6Parameters::standard()?;
let r = run_pm6(&mol, &params, &Pm6Options::default())?;
println!("ΔHf = {} kcal/mol", r.heat_of_formation_kcal);

Python

pip install maturin
maturin develop --release --features python   # or: pip install pm6-rs-python

import pm6_rs
import numpy as np

numbers = [8, 1, 1]
positions = np.array([[0.0, 0.0, 0.0], [0.9584, 0.0, 0.0], [-0.24, 0.9278, 0.0]])

# charge, multiplicity and reference ("auto"/"rhf"/"uhf") are per-call arguments.
out = pm6_rs.single_point(numbers, positions, charge=0.0, multiplicity=1, reference="auto")
print(out["heat_of_formation_kcal"], out["charges"])

g = pm6_rs.gradient(numbers, positions)          # dE/dR (Hartree/Bohr and eV/Å)
f = pm6_rs.forces(numbers, positions)            # −dE/dR
h = pm6_rs.hessian(numbers, positions)           # Hartree/Bohr², 3N×3N
w = pm6_rs.frequencies(numbers, positions)       # cm⁻¹ (evaluate at a minimum)

Periodic

cell = [[20.0, 0.0, 0.0], [0.0, 20.0, 0.0], [0.0, 0.0, 2.55]]   # Å
pbc = (False, False, True)                                       # a 1D chain

r = pm6_rs.calculate(numbers, positions, cell=cell, pbc=pbc, kpts=(1, 1, 8))
print(r["energy_ev"], r["stress_voigt_ev_per_angstrom3"])   # Voigt: xx yy zz yz xz xy

ph = pm6_rs.phonons(numbers, positions, cell, pbc=pbc,
                    supercell=(1, 1, 7), kpts=(1, 1, 8))
print(ph["frequencies_cm"][0], ph["decay_ratio"])

# Relax the ions, or the lattice, or both (new in 0.2.6).
out = pm6_rs.relax_cell(numbers, positions, cell, pbc=pbc, kpts=(1, 1, 8),
                        relax="all", cell_mask=[0, 0, 1, 0, 0, 0])
print(out["cell_angstrom"][2][2], out["max_stress_ev_per_angstrom3"])
print(out["cutoff_shell_crossed"])   # did it walk over the lattice-parameter step?

An explicitly-given kpts is now checked the way the Γ path has always been: a mesh whose Born–von Kármán supercell does not clear the integral cutoff is refused rather than answered. That is not conservatism — methane in a 6 Å cell at kpts=(2,2,2) returns −205.86 eV against a converged −177.17. allow_undersampled_cell=True overrides it, which is the right choice when the mesh is deliberately coarse (comparing an analytic derivative against a finite difference of its own energy, say).

charge is the net charge per cell, and heat_of_formation_kcal is None — ΔH_f is not defined per cell. kpts=None picks a mesh from the cell rather than defaulting to Γ, which is wrong (not merely coarse) on a small cell.

ASE

from ase.build import molecule
from pm6_rs.ase import PM6

atoms = molecule("H2O")
atoms.calc = PM6(charge=0, multiplicity=1, reference="auto")
print(atoms.get_potential_energy())   # eV
print(atoms.get_forces())             # eV/Å
print(atoms.calc.get_gradient())      # eV/Å  (= −forces)
print(atoms.calc.get_hessian())       # eV/Ų, 3N×3N

The calculator reads atoms.cell and atoms.pbc, so a periodic system needs no extra arguments, and get_stress() returns ASE's six-component Voigt vector in eV/ų — enough for StrainFilter, FrechetCellFilter and ase.md.npt.NPT:

from ase import Atoms
from ase.filters import FrechetCellFilter
from ase.optimize import BFGS

# Polyethylene: one C2H4 unit per cell, periodic along z only.
atoms = Atoms("C2H4",
              positions=[[0.000, 0.635, 0.000], [0.000, -0.635, 1.275],
                         [0.876, 1.235, 0.000], [-0.876, 1.235, 0.000],
                         [0.876, -1.235, 1.275], [-0.876, -1.235, 1.275]],
              cell=[20.0, 20.0, 2.55], pbc=(False, False, True))
atoms.calc = PM6(kpts=(1, 1, 8))

# mask relaxes the c axis only, which is the one that is periodic.
BFGS(FrechetCellFilter(atoms, mask=[0, 0, 1, 0, 0, 0])).run(fmax=0.05)
print(atoms.cell[2, 2])          # 2.529 A  (experiment: 2.53 A)
print(atoms.get_stress()[2])     # -2.6e-4 eV/A^3, i.e. relaxed

Ten BFGS steps, 7 s. The relaxed repeat lands on the experimental polyethylene c-axis, which nothing in the calculation was fitted to.

Attribution

pm6-rs is an independent Rust implementation, but the science in it is not original to this project. Four upstreams are involved, each with its license text and a notice naming the exact files here that came from it:

upstream	license	what came from it	notice
MOPAC v23.2.5	Apache-2.0	the PM6 parameter tables in src/data/, the reference algorithms reimplemented across src/, and the oracle every numeric claim on this page is measured against	third_party/mopac/NOTICE
PySEQM	BSD-3-Clause	the closed forms and frame convention for the s/p two-centre integrals in src/integrals.rs	third_party/pyseqm/NOTICE
dftd3 (through MOPAC)	see notice	the D3 model and its C6 / r0ab / r2r4 reference tables	third_party/dftd3/NOTICE
h_bonds4 (through MOPAC)	see notice	the H4 hydrogen-bond and X halogen-bond corrections	third_party/h_bonds4/NOTICE

MOPAC is not vendored here and is not required to run pm6-rs; it is required to re-run the oracle in tools/oracle/. Its copyright notice — Molecular Orbital PACkage (MOPAC), Copyright 2021 Virginia Polytechnic Institute and State University — is retained in full at third_party/mopac/NOTICE, and every licence text and notice under third_party/ ships inside the wheel as well as the source distribution. See third_party/README.md for the index and THIRD_PARTY_NOTICES.md for the same material organised by what it is rather than where it came from.

References

• J. J. P. Stewart, J. Mol. Model. 13, 1173 (2007) — PM6.
• W. Thiel, A. A. Voityuk, J. Phys. Chem. 100, 616 (1996) — MNDO-d.
• S. Grimme et al., J. Chem. Phys. 132, 154104 (2010) — D3.
• J. Řezáč, P. Hobza, J. Chem. Theory Comput. 8, 141 (2012) — H4.
• J. Řezáč, P. Hobza, Chem. Phys. Lett. 506, 286 (2011) — X (halogen).
• G. Zhou et al., J. Chem. Theory Comput. 16, 4951 (2020) — PySEQM.
• J. E. Moussa, J. J. P. Stewart, J. Open Source Softw. 11(119), 8025 (2026) — MOPAC.