Metadata-Version: 2.4
Name: pm6-rs-python
Version: 0.2.7
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Rust
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: Operating System :: OS Independent
Classifier: Typing :: Typed
Requires-Dist: numpy>=1.21
Requires-Dist: ase>=3.22 ; extra == 'ase'
Provides-Extra: ase
License-File: LICENSE
License-File: THIRD_PARTY_NOTICES.md
License-File: third_party/README.md
License-File: third_party/mopac/LICENSE
License-File: third_party/pyseqm/LICENSE
License-File: third_party/dftd3/NOTICE
License-File: third_party/h_bonds4/NOTICE
License-File: third_party/mopac/NOTICE
License-File: third_party/pyseqm/NOTICE
License-File: third_party/rust-crates/NOTICE
License-File: third_party/rust-crates/LICENSES.txt
Summary: PM6 semiempirical method (pm6-rs) with a Python-native API and an ASE calculator: periodic boundary conditions, analytic gradients, stress and Hessians, divide-and-conquer, and the PM6-D3/D3H4/D3H4X corrections
Keywords: chemistry,computational-chemistry,semiempirical,pm6,nddo,quantum-chemistry,ase,periodic-boundary-conditions,phonons
Home-Page: https://github.com/ss0832/pm6-rs
Author: ss0832
License-Expression: GPL-3.0-or-later
Requires-Python: >=3.9
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Changelog, https://github.com/ss0832/pm6-rs/blob/main/CHANGELOG.md
Project-URL: Homepage, https://github.com/ss0832/pm6-rs
Project-URL: Issues, https://github.com/ss0832/pm6-rs/issues
Project-URL: Source, https://github.com/ss0832/pm6-rs

# 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

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

## CLI

```sh
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
```

```sh
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`:

```sh
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:

```sh
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:

```sh
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

```rust
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

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

```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

```python
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

```python
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`:

```python
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](https://github.com/openmopac/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`](third_party/mopac/NOTICE) |
| [PySEQM](https://github.com/lanl/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`](third_party/pyseqm/NOTICE) |
| `dftd3` (through MOPAC) | see notice | the D3 model and its `C6` / `r0ab` / `r2r4` reference tables | [`third_party/dftd3/NOTICE`](third_party/dftd3/NOTICE) |
| `h_bonds4` (through MOPAC) | see notice | the H4 hydrogen-bond and X halogen-bond corrections | [`third_party/h_bonds4/NOTICE`](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`](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`](third_party/README.md) for the index and
[`THIRD_PARTY_NOTICES.md`](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.

