Semiempirical Methods¶

vibe-qc ships a self-contained semiempirical platform covering four method families, DFTB, GFN-xTB, NDDO, and INDO, for molecules and periodic solids. The platform is vibe-qc’s own implementation, not a wrapper around external programs (External QC codes (ORCA / Psi4 / others)). The MACE machine-learning interatomic potential is documented separately because it is an external pre-trained model, not a semiempirical electronic-structure method (Machine-learning interatomic potentials (MACE)).

Warning

Production readiness varies by method. DFTB is intended for screening/preoptimization; GFN2-xTB is gated experimental; PM6, OM2 and OM3 are production for molecular work (PES shape pinned by regression, published-benchmark validation for OM2, see the status table); OM1 is experimental until its analytic core-valence ECP lands; MSINDO is reference-parity validated for its current molecular scope. See the status table and the comparative production brief in Semiempirical and MACE method comparison.

Method families¶

Family	Methods	Best for	Cost
DFTB	DFTB0, SCC-DFTB, UDFTB0, USCC-DFTB	Screening, preoptimization	Fastest
GFN-xTB	GFN2-xTB	Organic / main-group	Fast
NDDO	PM6, OM1, OM2, OM3	Development benchmarking, pre-screening	Fast
INDO	MSINDO	Reference-parity molecular semiempirical runs inside the supported element/spin scope	Fast

Quick start¶

DFTB0, non-self-consistent tight-binding¶

from vibeqc.semiempirical import DFTB0Model
from vibeqc import Molecule, Atom

mol = Molecule([
    Atom(8, [0.00, 0.00, 0.00]),
    Atom(1, [1.55, 0.90, 0.00]),
    Atom(1, [-1.55, 0.90, 0.00]),
])

model = DFTB0Model(mol)
print(f"Energy: {model.energy():.6f} Ha")
print(f"Gradient shape: {model.gradient().shape}")  # (3, 3)

SCC-DFTB, self-consistent charges¶

from vibeqc.semiempirical import SCCDFTBModel

model = SCCDFTBModel(mol, charge_mixing=0.2)
print(f"Energy: {model.energy():.6f} Ha")

charge_mixing is the maximum fraction used by the molecular SCC charge update. At zero electronic temperature, bounded vector Aitken relaxation may reduce and subsequently recover that fraction when consecutive Mulliken residuals reveal charge sloshing; it never raises the fraction above the value requested by the caller. Molecular geometry optimization also seeds each new geometry from the preceding converged Mulliken charges, keeping the optimizer on a continuous SCC branch.

GFN2-xTB, published parameters (experimental)¶

from vibeqc.semiempirical import GFN2Model
from vibeqc.semiempirical.methods.gfn2_params import load_gfn2_params

params = load_gfn2_params()  # auto-fetches 86-element Grimme parameter set
model = GFN2Model(mol, params=params, warn=False)  # warn=False silences experimental gate
print(f"Energy: {model.energy():.6f} Ha")

Warning

GFN2-xTB is gated experimental. The H⁰/overlap deep-state bug was fixed 2026-06-01, and molecular AES (dipole + quadrupole), the GAM3 third-order term, and post-SCF native D4-style dispersion are now implemented. The remaining production gates are periodic AES image-cell multipole Ewald terms, stricter periodic molecular-limit parity, difficult periodic/polar fixtures, and a full external-parity matrix against xtb. See Method status. The bounded h-BN expansion sweep at lattice scales 1.02 through 1.06 converges on a smooth branch when given an explicit 1000-iteration total budget; its scale-1.06 point correctly fails under a 500-iteration cap. The seven-point bulk-Si sweep converges within 200 iterations. Periodic GFN2 honors max_iter as a hard total across its ordinary and automatic-stabilization attempts, and the public runner forwards both that budget and conv_tol to the native SCC loop. These bounded results do not relax the experimental GFN2 gate. Separately, the graphene periodic-PM6 sweep now has finite image exchange and converges. The one-carbon-chain periodic OMx sweep now cancels neutral image monopoles and crosses its 15-bohr image-shell boundary smoothly; diamond-Si periodic PM6 now uses one shared gamma kernel for electronic and core image monopoles and crosses its 15-bohr shell boundary smoothly. The h-BN PM6 and graphene OMx physics gaps remain open. The post-SCF native D4 term is included for H, He, B, C, N, O, F, and Ne; outside that reference-data set GFN2 returns zero D4 and emits GFN2D4UnsupportedWarning.

PM6, NDDO with published parameters¶

from vibeqc.semiempirical import PM6Model

model = PM6Model(mol)
print(f"Energy: {model.energy():.6f} Ha")

The high-level closed-shell run_job(method="pm6") route first attempts the ordinary zero-temperature SCF. If hard occupations cycle across near-degenerate fragment orbitals, it uses a bounded finite-temperature occupation homotopy and then cools back to an idempotent zero-temperature density. Only the cooled PM6 energy is accepted; the finite-temperature intermediate is never reported as a successful result.

Gamma-periodic PM6 instead uses Pulay extrapolation of the physical commutator residual. Convergence requires both the density change and [F, D] residual to meet conv_tol; this prevents a repeated extrapolated Fock matrix from being accepted when its density is not stationary under the physical Fock operator. It uses the shared guarded Pulay history, which shortens a linearly dependent history before extrapolating and falls back to the current physical Fock matrix when no stable history remains. This avoids the cancellation-dominated extrapolation that previously made the archived Pa-3 dry-ice scale-0.96 point fail after 1200 iterations and made its scale-0.98 neighbor unusually slow. The archived proton-ordered Ice-Ih and bounded dry-ice expansion slices now converge on smooth solver branches.

Periodic PM6 remains experimental until its energies receive independent external-reference validation. In particular, dry-ice energies remain sensitive to direct-space image-shell cutoff; solver convergence does not establish a physical molecular-crystal equation of state. The periodic diatomic core branch also cancels the neutral-pair monopole retained by the electronic decomposition. A closed-shell fcc-Ar regression pins that cancellation and the repulsive compressed branch. Bare PM6 has no validated rare-gas dispersion model here, however, so its dissociative Ar curve is not a physical rare-gas equation of state.

For elements such as silicon, whose electronic block uses the multi-term NDDO gamma expansion, the neutral core baseline uses that identical kernel rather than the separate PM6 short-range core kernel. Damping and Gaussian pair terms remain corrections on top. This keeps a neutral image shell neutral before truncation: the diamond-Si 15-bohr sphere changes from 55 to 43 cells between lattice scales 1.02 and 1.04, but the total-energy step is about 0.00195 Ha rather than the former 1.01917 Ha. This is a bounded implementation regression, not independent MOPAC validation.

OMx, orthogonalization-corrected NDDO¶

from vibeqc.semiempirical import OMxModel

model = OMxModel(mol, variant="om2")  # "om1", "om2", or "om3"
print(f"Energy: {model.energy():.6f} Ha")

Through `run_job`¶

All seven methods are available via vibeqc.runner.run_job():

from vibeqc import run_job

run_job(mol, method="dftb0", optimize=True, output="h2o_dftb0")
run_job(mol, method="pm6", output="h2o_pm6")
run_job(mol, method="gfn2_xtb", output="h2o_gfn2")   # emits experimental warning

run_job treats Molden and population sidecars independently. DFTB0, SCC-DFTB, and GFN2-xTB emit .population.txt and .population.json by default from the engine’s native net atomic Mulliken charges. DFTB0 forms these charges from its final one-shot density and overlap; UDFTB0 uses the sum of its alpha and beta densities. Analyses that require the Gaussian-AO property stack are present as explicit unsupported: entries. Molden remains unavailable for every semiempirical route because the current writer serializes Gaussian GTO shells, not the methods’ minimal Slater-orbital bases. PM6/OMx and MSINDO population sidecars remain unavailable until their native result adapters expose a validated atomic-population contract. An explicit unsupported write_molden_file=True or write_population_file=True fails before calculation; None (the default) selects only the sidecars the chosen route can produce truthfully.

DFTB0-SECCM (experimental)¶

The first non-MSINDO SECCM adapter is available through the direct run_dftb0_seccm API. It consumes a separately built and finite-group-bound SECCM topology; it does not reinterpret the ordinary Gamma-periodic DFTB0 cutoff domain. The current T3a gate accepts neutral, closed-shell, insulating one-dimensional H/C cyclic clusters with the built-in in-house screening parameter set:

import numpy as np

from vibeqc import Atom, Molecule
from vibeqc.semiempirical import run_dftb0_seccm
from vibeqc.semiempirical.seccm import (
    bind_finite_group,
    build_seccm_topology,
)

a = np.array([8.0, 0.0, 0.0])
coords = np.array([[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]])
topology = build_seccm_topology(
    coords, [a], length_unit="bohr", geometry_quantum=1.0e-10
)
topology = bind_finite_group(
    topology,
    primitive_vectors=[a],
    replicas=(1, 1, 1),
    geometry_tolerance=1.0e-9,
    length_unit="bohr",
)
mol = Molecule([Atom(1, xyz.tolist()) for xyz in coords])
result = run_dftb0_seccm(mol, topology, compute_gradient=True)

The reported energy is per primitive cell. Its complete T3a closure is

energy = electronic_energy + repulsive_energy
       + long_range_energy + dispersion_energy + specific_energy

where the last three terms are explicitly zero. The adapter evaluates the full finite cyclic cluster, so a valid group-only defect topology is allowed; translation-orbit reduction is not used. With compute_gradient=True, result.gradient is the analytic derivative of the same per-cell electronic and repulsive closure over the frozen record set. It does not differentiate a topology switch or the cyclic translations.

Periodic reaction paths can select this adapter explicitly with run_neb(..., method="dftb0", seccm_topology=topology). The endpoint cell must match the topology translations. A topology containing exact ownership ties also requires a positive seccm_max_tie_score_excursion trust bound. SCC charges, open shells, charged cells, stress, smearing, broader parameter sets, and two- or three-dimensional public calls fail closed.

Periodic systems¶

Periodic support is route-specific. DFTB0, SCC-DFTB, GFN2-xTB, PM6, and OMx have Gamma-point periodic energy routes; their periodic gradients and stress remain experimental and route-labeled as analytic or finite-difference stopgaps in the status table below. Full k-point public semiempirical routes are currently limited to closed-shell DFTB0 and SCC-DFTB energy/band jobs plus the experimental lower-level DFTB derivative and NEB facades. MSINDO periodic work uses the SECCM cyclic-cluster route, not the Gamma/k-point runner.

DFTB periodic¶

from vibeqc._vibeqc_core import PeriodicSystem, Atom
from vibeqc._vibeqc_core import semiempirical as _se
import numpy as np

# 1D carbon chain
atoms = [Atom(6, [0.0, 0.0, 0.0])]
cell = np.diag([2.5, 30.0, 30.0])
system = PeriodicSystem(1, cell, atoms)

params = _se.SemiempiricalParameters.dftb0_default()
result = _se.run_dftb0_gamma(system, params)
print(f"Energy: {result.energy:.6f} Ha")

GFN2 periodic¶

from vibeqc._vibeqc_core.semiempirical import xtb as _xtb
from vibeqc.semiempirical.methods.gfn2_params import load_gfn2_params

params = load_gfn2_params()
result = _xtb.run_gfn2_xtb_gamma(system, params)

The native Gamma driver also accepts explicit finite-temperature occupations for metallic or near-degenerate cells:

opts = _xtb.XTBSccOptions()
opts.electronic_temperature = 0.001  # electronic k_B T in Hartree
result = _xtb.run_gfn2_xtb_gamma(system, params, opts)

This applies particle-number-conserving Fermi-Dirac occupations. The default remains exact zero-temperature Aufbau occupation; vibe-qc does not silently select a smearing temperature. The reported energy is the GFN2 internal energy at the fractional-occupation density, not a Helmholtz free energy.

PM6 / OMx periodic¶

from vibeqc.semiempirical import PeriodicPM6Model, PeriodicOMxModel

pm6 = PeriodicPM6Model(system)
print(f"PM6: {pm6.energy():.6f} Ha")

omx = PeriodicOMxModel(system, variant="om2")
print(f"OM2: {omx.energy():.6f} Ha")

Periodic OMx constructs charge-conserving nonzero-image potentials from Mulliken populations, using Tr(P S) rather than Tr(P) as the electron count in its non-orthogonal basis. Directed image overlap, resonance, and exchange blocks are counted once, while image core pairs carry their required one-half energy weight. The one-carbon-chain cutoff-shell regression is stable within 0.7e-6 Ha when its direct-space cutoff doubles from 15 to 30 bohr. This closes the neutral-monopole discontinuity, not the method’s production gate: independent MOPAC parity and broader polar/two-dimensional validation remain open.

Preoptimization workflows¶

Use semiempirical methods for fast structure preoptimization before an expensive DFT calculation:

from vibeqc.semiempirical import preoptimize_periodic

# Preoptimize a periodic system with DFTB0, then run DFT
preoptimize_periodic(
    system,
    method="dftb0",
    fmax=0.01,
)

For molecular systems, use optimize=True with run_job:

from vibeqc import run_job

# Preoptimize with DFTB0, then refine with DFT
run_job(mol, method="dftb0", optimize=True)
run_job(mol, method="rks", functional="PBE", basis="def2-svp", optimize=True)

Method status¶

Method	Status	Energy accuracy	Gradient	Periodic	Open-shell	Elements
DFTB0 / UDFTB0	Screening/preopt	In-house parameters, not DFTB+ parity	Analytic	Gamma; full-k DFTB0 closed-shell	molecular/Gamma yes; full-k closed-shell only	91 in-house
DFTB0-SECCM	Gated experimental	In-house H/C screening parameters; finite-cluster T3a closure	Analytic, fixed topology	1-D direct finite torus	no	H, C
SCC-DFTB / USCC	Screening/preopt	In-house parameters, not DFTB+ parity	Analytic	Gamma; full-k SCC-DFTB closed-shell	molecular/Gamma yes; full-k closed-shell only	91 in-house
GFN2-xTB	Gated experimental	External xTB parity matrix still open	Analytic (FD-consistent to <1e-5 Ha/bohr on H2O/CH4/NH3: full H0 shape + shell-ES + 3rd-order + AES derivatives; periodic Gamma is fixed-charge)	Gamma experimental; full-k gated	no	86 fetched (LGPL)
PM6 / UPM6	Production (molecular)	Physical PES pinned by regression; spherical Klopman-Ohno fallback where MOPAC diatomic data is absent (wells ~0.1-0.3 bohr long)	FD	Gamma experimental, closed-shell	molecular yes; periodic no	82 bundled MOPAC
OM2 / OM3	Production (molecular)	Published OMx Hamiltonian (Dral 2016); relative energetics match published benchmarks (H3- bend 0.1 kcal/mol, ethane barrier ±0.6); bond minima ~0.1-0.2 Å long (documented integral stand-ins)	FD	Gamma, legacy-model parity only, closed-shell	molecular yes; periodic no	5 published (H,C,N,O,F)
OM1	Experimental (warns)	Analytic core-valence ECP (Kolb-Thiel 1993) not implemented; X-H bonds ~0.3 A short, close contacts can collapse	FD	Gamma, legacy-model parity only, closed-shell	molecular yes; periodic no	5 published (H,C,N,O,F)
MSINDO	Production within scope	Reference MSINDO parity <=1 uHa (INDO + NDDO)	FD and analytic molecular/CCM paths	SECCM 1-D/2-D/3-D + Ewald	molecular UHF s/p/d within validated fixtures; SECCM closed-shell	H-Xe (Z=1-54); NDDO H,Li-F,Na-Cl

The same implementation labels are available from Python through vibeqc.semiempirical.semiempirical_route_status(route). Routes such as msindo-cis, msindo-cis-gradient, msindo-ovgf, msindo-md, and msindo-metadynamics are explicitly marked python-reference until their hot loops move to native kernels or are declared intentionally orchestration-only. periodic-pm6 and periodic-omx are marked mixed-native: Gamma energy uses native NDDO kernels, while the current gradient, stress, and cell-optimization helpers still rely on finite-difference Python orchestration. The lookup also labels molecular pm6-gradient-fd and omx-gradient-fd as native-fd: their displacement loops are C++-backed, but they remain finite-difference stopgaps until analytic molecular NDDO gradients land. The lookup accepts public method aliases such as dftb0, scc-dftb, gfn2xtb, om2, om2-gradient-fd, gfn2, msindo, and ccm.

Direct run_semiempirical(...) results compute gradients lazily when result.gradient() is called. DFTB0/SCC-DFTB, GFN2-xTB, PM6/UPM6, and OMx expose their existing gradient surfaces through that adapter, while closed-shell MSINDO INDO uses the native analytic-gradient route there. MSINDO NDDO and open-shell MSINDO keep gradient() unset at the unified runner layer until those gradient surfaces are promoted. Geometry-optimizer SemiempiricalProvider calls the same unified runner, preserving the documented MSINDO finite-difference force fallback when a MSINDO result is energy-only.

Element coverage¶

DFTB, 91 elements (H-U except Po/Z=84), including 3d/4d/5d transition metals, lanthanides (La-Lu), and early actinides (Ac-U). All parameters are in-house estimates; DFT-fitted production repulsive potentials are deferred.

GFN2-xTB, 86 elements from the published Grimme-group parameter set. Parameters are fetched on demand at first use (LGPL-3.0 licensed, not bundled, see ADR-002).

PM6, 82 elements from the bundled MOPAC PM6 parameter cache (Apache-2.0 provenance in the TOML header), with the Stewart 2007 H/C/N/O/F subset still available. The public wrapper auto-selects the MOPAC-derived cache for elements outside H/C/N/O/F.

OMx, 5 elements (H, C, N, O, F) from Dral 2016 Tables 1-3.

# Check element coverage
from vibeqc.semiempirical import SemiempiricalParameters

params = SemiempiricalParameters.dftb0_default()
elements = [Z for Z in range(1, 93) if params.has_element(Z)]
print(f"DFTB covers {len(elements)} elements")

Parameter customisation¶

DFTB custom parameters¶

from vibeqc.semiempirical import SemiempiricalParameters

custom = SemiempiricalParameters()
custom.add_element(
    Z=1, on_site=[-0.21], zeta=[1.24],
    hubbard_u=0.42, valence_electrons=1,
)
custom.add_element(
    Z=8, on_site=[-0.89, -0.33], zeta=[2.25, 2.25],
    hubbard_u=0.45, valence_electrons=6,
)
# Set repulsive pair (R⁻¹² form)
custom.set_repulsive_pair_analytic(1, 1, A=5.0)
custom.set_repulsive_pair_analytic(1, 8, A=15.0)
custom.set_repulsive_pair_analytic(8, 8, A=40.0)

model = DFTB0Model(mol, params=custom)

GFN2 parameters¶

GFN2-xTB parameters are fetched automatically from the Grimme group’s GitHub repository. To force a refresh:

from vibeqc.semiempirical.methods.gfn2_params import load_gfn2_params

params = load_gfn2_params(force_refetch=True)

For offline batch planning, probe the local cache without opening the network:

from vibeqc.semiempirical import semiempirical_route_runtime_available

if not semiempirical_route_runtime_available("gfn2_xtb"):
    # Mark GFN2-xTB rows unavailable before submitting the batch.
    ...

PM6 parameters¶

from vibeqc.semiempirical.methods.pm6_params import load_pm6_params

params = load_pm6_params()
model = PM6Model(mol, params=params)

Comparing against external programs¶

Reference energies from external programs can be obtained via out-of-process subprocess runners (External QC codes (ORCA / Psi4 / others)):

from examples.regression.core.runner_xtb import energy as xtb_energy
from examples.regression.core.runner_mopac import energy as mopac_energy
from examples.regression.core.runner_dftbp import energy as dftbp_energy

print(f"xTB GFN2 H2O:  {xtb_energy('H2O'):.6f} Eh")
print(f"MOPAC PM6 H2O: {mopac_energy('H2O'):.6f} Ha")
print(f"DFTB+ H2O:     {dftbp_energy('H2O'):.6f} Ha")

These runners require the external program to be installed on $PATH (see each runner’s docstring for install instructions).

Performance tips¶

DFTB0 is 3-5× faster than SCC‑DFTB (no SCF loop). Use it for preoptimization where charge self-consistency is less important.
DFTB gradients (DFTB0 and SCC-DFTB) are analytic and match finite differences tightly; the SCC energy is variational in the density, so its fixed-charge gradient is exact at SCC convergence.
Periodic systems support Gamma-point energy routes for DFTB0/SCC-DFTB, GFN2-xTB, PM6, and OMx. Public full-k semiempirical support is DFTB0 and SCC-DFTB only, closed-shell only. Increase the lattice cutoff (cutoff_bohr) for tight cells.
Memory is negligible, the basis is minimal (one function per valence shell).

Known limitations¶

DFTB repulsive potentials are in-house R⁻¹² estimates; DFT-fitted production repulsives are deferred.
GFN2-xTB still lacks periodic AES image-cell multipole Ewald terms and a closed external xtb parity matrix (Method status).
PM6 reports a PM6-like total, not a MOPAC heat of formation; use MOPAC out-of-process when exact MOPAC convention parity is required.
Periodic PM6 cancels neutral-pair lattice monopoles, but bare PM6 has no validated rare-gas dispersion model; do not use its dissociative fcc-Ar curve as a physical equation of state.
OM2/OM3 are production molecular paths within their documented H/C/N/O/F scope; OM1 remains experimental until the analytic core-valence ECP lands.
Periodic GFN2/NDDO gradients are finite-difference only; analytic periodic NDDO gradients are deferred.
MSINDO molecular closed-shell analytic gradients and closed-shell CCM analytic / finite-difference gradients are available through the native route inside their documented scopes. NDDO, odd-electron analytic gradients, and excited-state/root-tracking gradients remain on their documented fallback or reference paths. See MSINDO (semiempirical INDO).