gfn2-rs

A Rust implementation of GFN2-xTB-style energies, analytical nuclear gradients, stress, Hessians, and periodic (Gamma / Monkhorst-Pack k-point) electrostatics, with an optional Python/ASE interface. Dense linear algebra uses faer; no BLAS/LAPACK required. Unofficial — use with caution.

Names: Cargo package gfn2-rs; CLI + Python package gfn2_rs; Rust library target gfn2_rs_core (alias to gfn2_rs downstream if preferred).

Scope

Area	Included
SCC energy	Gamma (real path) and explicit Monkhorst-Pack k-points (Bloch path).
Gradients / stress	Analytical, Gamma and Bloch; stress for periodic cells (opt-in).
Hessian	Fully analytic CPHF response (default since 0.5.0), exact to the finite-difference truncation floor on molecular, Gamma-cell and k-point references (5.9e-8 / 7.6e-9 Eh/bohr² vs. the FD of the analytic gradient, all terms enabled). --hessian-mode analytic|semi-numerical|fd selects the strategy. Since 0.5.0 open-shell and spin-polarized molecular Gamma references (integer channel occupations) also take the analytic route, via a coupled two-channel spin CPHF; LC, EEQ-BC D4, and open-shell references carrying a lattice / k-point mesh / fractional occupations fall back to fd automatically, which HessianResult::mode_used reports.
Periodic AES/Ewald	3D + 2D-slab multipolar electrostatics (shell charge, dipole, quadrupole). Charged cells need --charged-pbc background.
D4 dispersion	Pair + ATM from in-tree DFT-D4 reference data.
TD-GFN2 / spin	spGFN2 spin projection, TD-GFN2/TDA excitations + analytical excited-state gradients (closed-shell + spin-conserving open-shell, Gamma + Bloch).
Python / ASE	PyO3 bindings + one ASE calculator GFN2RSCalculator (energy/forces/stress/charges + TD/Hessian methods).

Documentation

Start here for this release: docs/whats-new-0.6.1.md — what changed, what it costs, and an explicit account of what is still open. Read it first if you have periodic results from an earlier version.

CHANGELOG.md is the full release record, with the measurement behind every entry. The rest of docs/:

• docs/rust-api.md — Rust library API (run_electronic_core, native_analytic_gradient, compute_hessian, CPHF/TD-GFN2/spin entry points, divide-and-conquer SCC, analytic third derivatives, NN-xTB, parameter refitting).
• docs/python-api.md — Python (Gfn2NativeCalculator) and ASE (GFN2RSCalculator) API, including k-points, stress, TD-GFN2, Hessian, and spin controls.
• docs/measurements.md — measured scaling of the divide-and-conquer SCC on a water-cluster size series: fitted log-log exponents for the exact and DC paths, which stage dominates at which size, the DC-vs-exact accuracy, and a plain answer to "is this linear scaling?" (short version: DC removes the cubic eigensolve, but the calculation as a whole is not linear).
• CHANGELOG.md — the complete inventory, including the internal fixes and the per-test measured residuals.

Build

cargo build --release --bins          # CLI (target/release/gfn2_rs, gradient_check, stress_check, hessian_check, cubic_check)

Python extension (Python 3.9–3.14):

python -m pip install maturin "ase>=3.22" "numpy>=1.22"
maturin develop --release --features python

For a Python newer than PyO3's support window, set PYO3_USE_ABI3_FORWARD_COMPATIBILITY=1 before cargo/maturin.

Parameters

The GFN2-xTB parameter file is bundled: data/param_gfn2-xtb.txt is copied verbatim from the xtb repository (tag bleed, commit 2b5cd4829290775e575807daee21560f851ff7e1; LGPL-3.0-or-later, dual GPL-3.0-or-later — see THIRD_PARTY_NOTICES.md) and embedded into the binary with include_str!. Nothing has to be downloaded or configured: it is used by default. To override it with a different parameter file, pass --param PATH (which wins) or set GFN2_XTB_PARAM:

export GFN2_XTB_PARAM=/path/to/param_gfn2-xtb.txt   # (PowerShell: $env:GFN2_XTB_PARAM = "...")

Spin-polarization constants

Open-shell collinear spin-polarization (--spinpol) needs per-element spin constants W_ll' (Hartree). Since v0.5.0 these are bundled: the spGFN2 table for Z = 1..86 (all six shell pairs ss/sp/pp/sd/pd/dd per element) ships as data/param_spin_gfn2.txt and is embedded with include_str!, so --spinpol works with nothing else supplied. The values are transcribed verbatim (no unit conversion — they are already Hartree) from the tblite library (src/tblite/data/spin.f90, commit 6617517, LGPL-3.0-or-later); the file's own header and THIRD_PARTY_NOTICES.md section 3 carry the full provenance, and reference 4 below is the method citation. These atomic spin constants are not fitted to a specific GFNn-xTB method and are usable with GFN1/GFN2-xTB.

To override them, pass --spin-param FILE (which wins) or set GFN2_SPIN_PARAM — the same explicit > env > bundled order the GFN2-xTB parameter file uses. Format (one element per line, comments with #/!, - for shell pairs absent from the AO basis); see examples/spin_params_template.txt:

• H ss <W_ss> — a single named shell-pair value, or
• <El> <ss> <sp> <pp> <sd> <pd> <dd> — packed symmetric s/p/d matrix, or
• <El> <ss> <sp> <sd> <ps> <pp> <pd> <ds> <dp> <dd> — full row-major 3x3 matrix.

Use it from the CLI or either API:

# bundled constants - no extra files needed
gfn2_rs --xyz radical.xyz --electronic --charge 0 -u 1 --spinpol
# ...or override them (add --kgrid a b c [--kshift] for k-points)
gfn2_rs --xyz radical.xyz --electronic --charge 0 -u 1 \
  --spinpol --spin-param spin_params.txt

// Rust: explicit > GFN2_SPIN_PARAM > bundled, in one call
let params = gfn2_rs::SpinPolarizationParameters::load(None)?;   // bundled
// let params = gfn2_rs::SpinPolarizationParameters::load(Some("spin_params.txt"))?;
let spin = gfn2_rs::SpinPolarizationOptions { enabled: true, parameters: params };
// electronic_options.spin_polarization = spin;

# Python / ASE (low-level Gfn2NativeCalculator takes the same two kwargs)
calc = GFN2RSCalculator(charge=0.0, uhf=1, spin_polarized=True)
# spin_param_path="spin_params.txt" overrides the bundled constants

Lanthanides (Ln-xTB overlay)

--ln-xtb-param FILE (Python/ASE ln_xtb=True with ln_xtb_param=, Rust Gfn2Parameters::apply_lnxtb_overlay()) swaps the La-Lu (Z=57-71) element blocks for the re-optimized Ln-xTB blocks of reference 5 below (SI Part H).

The Ln-xTB parameters are not distributed with this project — see THIRD_PARTY_NOTICES.md section 2 for the reason.

python scripts/extract_lnxtb_params.py --si jcc70321-sup-0002-datas2.pdf -o param_lnxtb.txt

Then pass --ln-xtb-param param_lnxtb.txt, or point GFN2_LNXTB_PARAM at it and use plain --ln-xtb. Without one of those, --ln-xtb fails with a message saying so; it never falls back to plain GFN2 silently, because that would return a converged result that looks like Ln-xTB and is not. If you use the overlay, cite reference 5.

The 4f electrons are f-in-core. The SI reference protocol sets UHF to the Ln(III) 4f unpaired-electron count, parity-floored for even-electron systems — La/Ce/Yb/Lu 0, Pr/Nd/Er/Tm 2, Pm/Sm/Dy/Ho 4, Eu/Gd/Tb 6:

gfn2_rs --xyz NdF3.xyz --electronic --ln-xtb-param param_lnxtb.txt --charge 0 -u 2

SCC convergence

Since v0.5.0 the self-consistency loop is driven by a Fock-space ADIIS -> CDIIS accelerator (src/scc_accel.rs) by default. CDIIS is Pulay commutator DIIS on the orthonormal-basis error X^T(FPS-SPF)X; ADIIS is the Hu & Yang energy-model DIIS (J. Chem. Phys. 132, 054109 (2010)) minimized over the coefficient simplex. The driver runs ADIIS until the commutator infinity-norm falls below 1e-2, then switches to CDIIS.

gfn2_rs --xyz mol.xyz --electronic --scf-accel adiis-cdiis   # default
gfn2_rs --xyz mol.xyz --electronic --scf-accel charge        # legacy v0.4.x path

Python/ASE: scf_accelerator="adiis-cdiis" (default) or "charge"; Rust: ElectronicOptions::scf_accelerator.

The fixed points are identical to the legacy charge-space path (energies agree to ~1e-15..1e-13 Eh) and the iteration counts drop sharply — water dimer 29 -> 12, ferrocene 32 -> 15, and PrF3 under the Ln-xTB overlay, which previously stalled at 300 iterations, converges tightly in 64. Wall time of the full xtb-oracle suite is roughly halved.

A divergence guard inside the SCF loop watches the residual from iteration 8 onwards; if it exceeds 100x the best value seen so far, the run is handed back to the legacy charge-space DIIS/Broyden/mixing automatically, so pathological cases behave no worse than in v0.4.5. In v0.5.0 the accelerator drives the Gamma real path and the Bloch k-point path (per-k realified blocks with shared extrapolation coefficients; a 2x2x2 periodic water dimer drops from 23 to 12 iterations at fixed-point agreement 5e-15 Eh) — the spin-polarized and LC paths ignore the setting and use the legacy charge-space accelerators.

One caveat: on surfaces with multiple SCC solutions the energy-seeking ADIIS step can converge to a lower fixed point than damped iteration reaches. That is arguably the correct variational behaviour, but it can differ from loosely converged reference implementations (see the Pr/Sm note above).

Since v0.6.1: reading the iteration, not just its exit code

scf_converged true says the residual reached its tolerance. It does not say whether the iteration got there in three steps or nearly failed to. Every run now also prints scc_verdict — converging, sloshing (over-correcting) or stiff (converging at a rate that will not arrive) — and scc_class (insulator or metallic, from the frontier gap and the occupations). --scc-trace FILE writes the whole residual sequence with its ratios.

The two failing verdicts want opposite remedies, which is the reason for telling them apart: damping a stiff iteration makes it creep more slowly. On a stiff verdict the CLI converges at a raised electronic temperature and anneals back down to the requested one before it tries the convergence ladder. --scc-preset auto|insulator|metal selects which regime's settings to start from; auto decides by measuring the first solve.

docs/measurements.md is the census behind all of this, including the one thing it does not fix: a metallic cell at zero electronic temperature on a k-mesh. That failure is diagnosed there (it is the self-consistent multipole response, not charge sloshing) and its remedy is the finite electronic temperature the physics wants anyway — the shipped 300 K default is enough. --kerker-beta implements the charge-space analogue of plane-wave Kerker mixing; it is off by default because it was measured not to help.

CLI

Coordinates are XYZ/extXYZ in Angstrom (use --bohr for Bohr). See --help for the full option set; common workflows:

# single-point energy
gfn2_rs --xyz examples/water_dimer_distorted.xyz --electronic --charge 0 --uhf 0

# energy + analytical gradient (add --stress for periodic stress)
gfn2_rs --xyz examples/water_dimer_pbc_16a.xyz --electronic --gradients --stress

# L-BFGS optimization (analytical gradient; add --kgrid a b c for Bloch k-points)
gfn2_rs --xyz examples/water_dimer_distorted.xyz --electronic --optimize \
  --opt-output opt.xyz --opt-trace opt.tsv --opt-grad-tol 5e-4

# k-point mesh folded onto its time-reversal-irreducible half (opt-in)
gfn2_rs --xyz slab.xyz --electronic --kgrid 4 4 1 --kpoint-trs

# export the converged wavefunction for Multiwfn/Avogadro/Molden (Gamma point)
gfn2_rs --xyz examples/water_dimer_distorted.xyz --electronic --molden water.molden

# analytic Hessian / harmonic frequencies (analytic is the default since 0.5.0)
gfn2_rs --xyz examples/water_dimer_distorted.xyz --electronic --freq --hessian-mode analytic

# orbital energies, occupations and the frontier (0.6.1). --orbital-coefficients
# adds the MO matrix. On a k-point run the gap is the fundamental one, across
# the whole Brillouin zone, and is flagged when it is indirect.
gfn2_rs --xyz examples/water_dimer_distorted.xyz --orbitals

# relax the cell as well as the atoms (0.6.1): full | volume | fixed.
# Only strains between periodic axes move, so a slab relaxes in its plane and
# never along the vacuum.
gfn2_rs --xyz examples/h2_pbc_12a.xyz --opt-cell full --opt-stress-tol 1e-6

# is this SCC converging, oscillating or stalled? (0.6.1)
gfn2_rs --xyz slab.xyz --electronic --scc-trace scc.tsv   # prints scc_verdict, scc_class

# EXPERIMENTAL divide-and-conquer SCC: only the eigen-solve/density build is
# partitioned; H, S, SCC, AES, D4 and the Ewald electrostatics stay exact.
# Closed-shell only (molecular, Gamma-periodic or --kgrid meshes via a per-k DC
# step with one global Fermi level), needs --etemp > 0 (300 K default).
# Gradients and stress are fixed-partition analytic.
gfn2_rs --xyz cluster.xyz --dc --dc-core-atoms 8 --dc-buffer-radius 10 --gradients

--dc-fragments "1-3;4-6" supplies explicit 1-based cores instead of spatial boxes; Python/ASE take the same four knobs (dc=True, dc_core_atoms=, dc_buffer_radius=, dc_fragments=). See CHANGELOG.md for the accuracy, the quantified truncated-buffer residual and the scope errors, and docs/measurements.md for the measured scaling — DC turns the N^3 eigensolve into an N^1.2-1.4 partitioned solve, but the calculation as a whole stays ~N^2 because H/S, AES, D4 and the Ewald electrostatics are still exact and full-system every iteration.

--kpoint-trs (Python/ASE kpoint_trs=True, Rust ElectronicOptions::kpoint_time_reversal) merges k with -k and accumulates the weights. It is exact for the real SCC ground state, validated against the full mesh (energy to 1e-9 Eh, charges to 1e-8), and roughly halves the diagonalization work on non-Gamma meshes. It is off by default in v0.5.0. All k-point paths (SCF, gradients, band reconstruction, response) share one k-list builder, so they can never disagree about the mesh.

TD-GFN2/spin (--td-gfn2, --td-spin multiplicity:N, --sp-gfn2) and parameter finite-difference derivatives (param_deriv --active-targets ...) are also available; see --help.

Console commands in Python environments

Installing the package (maturin develop --release --features python, or pip install) puts the CLI drivers on PATH as console scripts: gfn2-rs (= gfn2_rs), gfn2-nn-train (= nn_train) and gfn2-nn-eval (= nn_eval). They run the same Rust code as the cargo-built binaries, with identical flags, output and exit codes:

gfn2-rs --xyz examples/h2.xyz --electronic
gfn2-nn-eval --weights model.gfn2nn --xyz examples/h2.xyz

Python / ASE

from ase import Atoms
from ase.optimize import LBFGS
from gfn2_rs.ase import GFN2RSCalculator      # reads GFN2_XTB_PARAM, or pass param_path=...

atoms = Atoms("H2", positions=[[0, 0, 0], [0.74, 0, 0]])
atoms.calc = GFN2RSCalculator(charge=0.0, uhf=0)
print(atoms.get_potential_energy())            # eV
print(atoms.get_forces())                      # eV/Angstrom

# periodic + k-points
atoms.set_cell([16, 16, 16]); atoms.set_pbc(True)
atoms.calc = GFN2RSCalculator(kgrid=[2, 1, 1], kshift=True)
print(atoms.get_potential_energy(), atoms.get_stress())

# explicit GFN2-RS workflows
tda = atoms.calc.get_tda(n_states=3, spin_adaptation="multiplicity:1")
hess = atoms.calc.get_hessian(relaxed=True)    # eV/Angstrom^2

ASE units throughout: Angstrom, eV, eV/A, eV/A^3. The low-level gfn2_rs.Gfn2NativeCalculator exposes the same controls without ASE.

Validation

cargo test --quiet --offline
GFN2_XTB_PARAM=/path/to/param_gfn2-xtb.txt maturin develop --release && python -m pytest -q

Finite-difference checkers (gradient_check, stress_check, hessian_check, cubic_check) compare analytical vs FD derivatives; GFN2_PROFILE=1 enables timing scopes. Example:

gfn2_rs/target/release/hessian_check --xyz examples/water_dimer_distorted.xyz --d4 \
  --step 1e-4 --out hessian_check.tsv

cubic_check is the experimental third-derivative (E^xyz) counterpart and front end. It prints max |analytic E^xyz - FD(analytic Hessian)| per term class (repulsion, D4 pair, D4 pair+ATM, electronic) together with the translation sum rule and permutation symmetry of each block, and --form assembles the tensor in one of the three public output forms. Molecular only by design — periodic (lattice-image / Ewald) third derivatives are out of scope, and a lattice is rejected with an explicit error.

gfn2_rs/target/release/cubic_check --xyz examples/water_dimer_distorted.xyz --step 1e-4
# optional: --charges FILE supplies the frozen D4 charges (default: all zero)

# --form dense: the whole tensor (symmetric-packed; capped at 3N <= 150)
# --form block: symmetric 3N x 3N slices M^(c)_ab = E_abc at --block-coords,
#               i.e. dH_ab/dR_c -- no cap, never materializes 3N^3
# --form vec:   symmetric 3N x 3N contractions sum_c E_abc v_c from --vectors FILE
gfn2_rs/target/release/cubic_check --xyz water.xyz --compare-forms
gfn2_rs/target/release/cubic_check --xyz big.xyz --no-fd --form block \
  --block-coords 0,1,2 --out fc3.tsv

# NN-xTB surface (v0.6.1, EXPERIMENTAL). Needs --no-fd: this binary's FD
# oracles displace a PLAIN-GFN2 closure, so comparing them against an NN
# tensor would measure the gap between two surfaces, not a derivative error.
gfn2_rs/target/release/cubic_check --xyz water.xyz --no-fd --form block \
  --block-coords 0 --nn-weights models/nnxtb-v0.6.1-provisional.gfn2nn

--nn-weights takes a plain path, resolved against the working directory, so the models/... spelling above assumes a checkout or an unpacked release archive. From a pip install the same model ships inside the package instead — python -c "import gfn2_rs; print(gfn2_rs.bundled_model_path())" prints its absolute path. See docs/pip_reachability.md.

On the NN path this is the frozen-parameter tensor: the exact E^xyz of the GFN2 surface pinned at the parameters the network produced for that geometry. Because those parameters move with the geometry, that is not the whole third derivative — it captures ~91% of it on water with the shipped model. The full chain-rule tensor is refused by name rather than approximated; see docs/nn-xtb.md.

The same three forms are on the Rust API (compute_third_derivatives with ThirdDerivativeRequest::{dense, block, vec}) and on the native Python calculator (calc.third_derivatives(..., form="dense"|"block"|"vec")); see docs/rust-api.md and docs/python-api.md. The assembly is complete for molecular systems — nothing is staged.

xtb oracle suite

tests/oracle_gfn2.rs (reference data in tests/oracle/xtb_data.rs, notes in tests/ORACLE_NOTES.md) checks this implementation against 18 geometries / 170 reference values transcribed from the xtb unit tests (tag bleed, commit 2b5cd48...). After fixing the radii-table conversion constant — xtb converts its atomicRad / covalentRadD3 tables with the legacy Bohr radius autoaa = 0.52917726, not the CODATA-2018 value — the suite agrees to:

quantity	agreement
water energy components	1e-11..1e-13 Eh
caffeine repulsion	3e-13 Eh
periodic GFN coordination numbers	<= 1.9e-9
analytic gradient vs finite difference	2.9e-10 (FD-exact)

Cross-implementation tolerances here are floored by the references' own generation settings (acc = 1.0: scfconv 1e-6 Eh, qconv 1e-4 e), not by this implementation; ORACLE_NOTES.md has the per-quantity details.

References

1. V. Alizadeh, J. Pototschnig, S. Ehlert, T. D. Kuehne, "Periodic GFN2-xTB in CP2K," ChemRxiv (2026). DOI: 10.26434/chemrxiv.15003939/v1.
2. C. Bannwarth, S. Ehlert, S. Grimme, "GFN2-xTB," J. Chem. Theory Comput. 15, 1652 (2019). DOI: 10.1021/acs.jctc.8b01176.
3. A. Hjorth Larsen et al., "The Atomic Simulation Environment," J. Phys.: Condens. Matter 29, 273002 (2017). DOI: 10.1088/1361-648X/aa680e.
4. H. Neugebauer, B. Bädorf, S. Ehlert, A. Hansen, S. Grimme, "High-throughput screening of spin states for transition metal complexes with spin-polarized extended tight-binding methods" (spGFNn-xTB), J. Comput. Chem. 44(27), 2120 (2023). DOI: 10.1002/jcc.27185. The spin constants W_ll' are distributed with the tblite library (https://github.com/tblite/tblite); the bundled data/param_spin_gfn2.txt is transcribed from tblite's src/tblite/data/spin.f90 at commit 66175172f1e80b899fd67f79128ee27b2be8dafa, see THIRD_PARTY_NOTICES.md.
5. Y.-Y. Zhang, "New Optimized Parameters of Extended Tight-binding Method Ln-xTB (Ln = La-Lu) to Explore Lanthanide Molecular Chemistry," J. Comput. Chem. 47(4), e70321 (2026). DOI: 10.1002/jcc.70321. Nothing from this paper is distributed with the project. Earlier releases bundled the parameter overlay, a reference bond-length table and 50 SI geometries; all were removed, because the publisher's terms of use forbid incorporating downloaded material into a repository and this project publishes to three of them. scripts/extract_lnxtb_params.py regenerates the parameter file from your own copy of the Supporting Information (Table S62 already prints it in the required syntax). The article is subscription-access and must be obtained from the publisher. See THIRD_PARTY_NOTICES.md section 2 and docs/lnxtb_licensing.md.
6. J. Xia et al., "Neural-network-parametrized extended tight-binding," Nat. Commun. 17, 7302 (2026). DOI: 10.1038/s41467-026-73184-z. The experimental NN-xTB module (src/nn/) follows this architecture; the implementation-level caveats and deliberate divergences are documented in docs/nn-xtb.md (section "Divergence from the paper and the official capsule"), and the training-data provenance in THIRD_PARTY_NOTICES.md. The paper's reference implementation (CodeOcean capsule 1479726, MIT, (c) 2025 Yufan Xia) was consulted for two published training settings; no capsule code or data is bundled here — see THIRD_PARTY_NOTICES.md § 4.
7. I. Batatia, D. P. Kovacs, G. N. C. Simm, C. Ortner, G. Csanyi, "MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields," NeurIPS 35 (2022). arXiv: 2206.07697. The NN-xTB descriptor encoder is a from-scratch MACE-type implementation (no e3nn parity claimed).
8. M. Nishida, K. Fujiwara, T. Taketsugu, M. Kobayashi, "Divide-And-Conquer Extended Tight-Binding Molecular Dynamics: A General-Purpose, Very Large-Scale Quantum Molecular Dynamics Method," J. Comput. Chem. 46, e70255 (2025). DOI: 10.1002/jcc.70255. The experimental divide-and-conquer SCC (--dc, src/dc.rs) is an independent implementation of the DC-xTB idea; see docs/measurements.md for what it does and does not deliver.

License

GPL-3.0-or-later.