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.

✅ = 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

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   chain.xyz --pbc FFT --kpts 1,1,8
pm6-rs stress   chain.xyz --pbc FFT --kpts 1,1,8
pm6-rs bands    chain.xyz --pbc FFT --kpts 1,1,8 --path "0,0,0;0,0,0.5"
pm6-rs phonons  chain.xyz --pbc FFT --kpts 1,1,8 --supercell 1,1,4
pm6-rs phonons-dfpt chain.xyz --pbc FFT --kpts 1,1,8 --qpoints "0,0,0;0.5,0,0"
pm6-rs energy   big.xyz   --dc --dc-buffer 12.0     # divide-and-conquer

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_ev_per_angstrom3"])   # Voigt: xx yy zz yz xz xy

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

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.

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. J. P. Moussa, J. J. P. Stewart, J. Open Source Softw. 11(119), 8025 (2026) — MOPAC.