Semiempirical DFTB — fast preoptimization

vibe-qc ships a self-contained DFTB (density-functional tight-binding) implementation for rapid structure preoptimization and large-system screening. DFTB is 10²–10⁴× faster than full DFT while capturing qualitative bonding, making it ideal for:

  • geometry preoptimization before DFT refinement

  • screening hundreds of conformers

  • quick energy estimates for large systems

  • periodic Γ‑point relaxation

Methods available

Method

Description

Closed‑shell

Open‑shell

DFTB0

Non‑self‑consistent tight‑binding

DFTB0Model

UDFTB0Model

SCC‑DFTB

Self‑consistent‑charge DFTB

SCCDFTBModel

USCCDFTBModel

+D3(BJ)

Dispersion‑corrected

DispersionCorrectedModel

DispersionCorrectedModel

All methods have analytic gradients (DFTB0 and UDFTB0 are exact; SCC and USCC use an approximate fixed‑charge gradient with a CP‑SCC charge‑response correction).

Quick start — H₂O single‑point

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

# H₂O at experimental geometry (bohr)
import numpy as np
theta = np.deg2rad(104.5 / 2)
mol = Molecule(
    [
        Atom(8,  [0.0, 0.0, 0.0]),
        Atom(1,  [1.81 * np.sin(theta), 1.81 * np.cos(theta), 0.0]),
        Atom(1,  [-1.81 * np.sin(theta), 1.81 * np.cos(theta), 0.0]),
    ],
    charge=0,
    multiplicity=1,
)

# Non‑self‑consistent (fastest)
dftb0 = DFTB0Model(mol)
print(f"DFTB0:  {dftb0.energy():12.6f} Ha")
print(f"gradient:\n{dftb0.gradient()}")

# Self‑consistent‑charge (more accurate charges)
scc = SCCDFTBModel(mol)
print(f"SCC:    {scc.energy():12.6f} Ha")

Geometry optimization

Use the ASE calculator interface for geometry optimization:

from vibeqc import Atom, Molecule
from vibeqc.semiempirical import DFTB0Model, SemiempiricalParameters
from ase import Atoms, units
from ase.optimize import BFGSLineSearch
from ase.calculators.calculator import Calculator
import numpy as np

class DFTB0Calculator(Calculator):
    """ASE Calculator wrapping DFTB0."""
    implemented_properties = ["energy", "forces"]

    def __init__(self, params=None):
        super().__init__()
        self.params = params or SemiempiricalParameters.dftb0_default()

    def calculate(self, atoms, properties, system_changes):
        pos_bohr = atoms.positions / units.Bohr
        mol = Molecule(
            [Atom(int(z), list(xyz))
             for z, xyz in zip(atoms.numbers, pos_bohr)],
            charge=0, multiplicity=1,
        )
        model = DFTB0Model(mol, params=self.params)
        self.results["energy"] = model.energy() * units.Hartree
        self.results["forces"] = -model.gradient() * (
            units.Hartree / units.Bohr
        )

# Build water molecule
theta = np.deg2rad(104.5 / 2)
mol = Molecule(
    [Atom(8, [0, 0, 0]),
     Atom(1, [1.81 * np.sin(theta), 1.81 * np.cos(theta), 0]),
     Atom(1, [-1.81 * np.sin(theta), 1.81 * np.cos(theta), 0])],
    charge=0, multiplicity=1,
)

atoms = Atoms(
    numbers=[8, 1, 1],
    positions=np.array([a.xyz for a in mol.atoms]) * units.Bohr,
)
atoms.calc = DFTB0Calculator()
opt = BFGSLineSearch(atoms)
opt.run(fmax=0.05, steps=50)

final = atoms.positions / units.Bohr
for i in range(3):
    print(f"Atom {i}: {final[i]}")

Open‑shell — NO₂ radical

from vibeqc import Atom, Molecule
from vibeqc.semiempirical import UDFTB0Model, USCCDFTBModel

no2 = Molecule(
    [
        Atom(7, [0.0, 0.0, 0.0]),
        Atom(8, [2.2, 0.0, 0.0]),
        Atom(8, [-1.1, 1.9, 0.0]),
    ],
    charge=0,
    multiplicity=2,  # doublet
)

udftb0 = UDFTB0Model(no2)
print(f"UDFTB0: {udftb0.energy():.6f} Ha")

uscc = USCCDFTBModel(no2)
print(f"USCC:   {uscc.energy():.6f} Ha")

Periodic systems

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

# 1D carbon chain at 3.0 bohr spacing
atoms = [Atom(6, [0.0, 0.0, 0.0])]
cell = np.eye(3) * 3.0
system = PeriodicSystem(3, cell, atoms, charge=0, multiplicity=1)

params = _se.SemiempiricalParameters.dftb0_default()

# Gamma‑point DFTB0
result = _se.run_dftb0_gamma(system, params)
print(f"Periodic DFTB0: {result.energy:.6f} Ha/cell")

# SCC variant
opts = _se.PeriodicSCCOptions()
opts.max_iter = 50
scc_result = _se.run_scc_dftb_gamma(system, params, opts)
print(f"Periodic SCC:   {scc_result.energy:.6f} Ha/cell")
print(f"  converged: {scc_result.converged}")

Dispersion corrections

from vibeqc.semiempirical import (
    SCCDFTBModel, DispersionCorrectedModel,
    DFTB_D3_DEFAULTS, d3bj_energy,
)

# Build a dispersion‑corrected model
scc = SCCDFTBModel(mol)
disp_model = DispersionCorrectedModel(scc, DFTB_D3_DEFAULTS)
print(f"SCC + D3(BJ): {disp_model.energy():.6f} Ha")

Preoptimization → DFT handoff

The primary use case: preoptimize with DFTB, then refine with DFT.

from vibeqc.semiempirical.preoptimize import preoptimize_molecule
from vibeqc import Molecule, Atom

# Starting geometry (rough guess)
mol = Molecule(
    [Atom(8, [0, 0, 0]),
     Atom(1, [1.5, 0, 0]),
     Atom(1, [-1.0, 1.2, 0])],
    charge=0, multiplicity=1,
)

# Preoptimize with DFTB0
opt_mol = preoptimize_molecule(mol, method="dftb0")
print("Preoptimized geometry:")
for a in opt_mol.atoms:
    print(f"  Z={a.Z}: {a.xyz}")

# Now feed into run_job for DFT refinement:
# from vibeqc import run_job
# run_job(opt_mol, basis="6-31G*", method="rks", functional="pbe")

Parameter customisation

from vibeqc.semiempirical import SemiempiricalParameters

# Built‑in parameter sets
default = SemiempiricalParameters.dftb0_default()
production = SemiempiricalParameters.dftb0_production()

# Query element data
print(f"O on‑site (s): {default.on_site_energy(8, 0):.4f} Ha")
print(f"O Hubbard U:   {default.hubbard_u(8):.4f} Ha")
print(f"kappa:          {default.kappa}")

# Build custom parameters
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)

Element coverage

The built‑in parameter sets cover 85 elements (H–U), including the 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.

Check coverage:

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

Performance tips

  • DFTB0 is 3–5× faster than SCC‑DFTB (no SCF loop). Use it for preoptimization where the SCC charge correction is less important.

  • Gradients are analytic for all methods. DFTB0 and UDFTB0 gradients match finite‑differences to machine precision.

  • Periodic systems use a Γ‑point lattice sum with a 15 bohr cutoff. For very tight cells, increase cutoff_bohr.

  • Memory is negligible — the basis is minimal (one function per valence shell).

Known limitations

  • SCC and USCC gradients use an approximate fixed‑charge formula (see roadmap). For quantitative geometry optimisation, prefer DFTB0.

  • Only Γ‑point periodic calculations are supported; general k‑points are planned.

  • Repulsive pair potentials for elements beyond the core set (H, C, N, O, F, P, S, Cl) use a default R⁻¹² estimator.