pydftb_torch

pydftb_torch is a scoped LibTorch-backed Python package for DFTB calculations. Slater-Koster parameter files are not bundled; pass a local SKF directory through skf_dir.

Supported calculation scope:

• DFTB1/2/3 molecular energy, forces, Hessian, and force constants.
• DFTB1/2/3 periodic energy, forces, Hessian, and stress.
• mDFTB2/3 molecular energy, forces, and Hessian through the molecular q/d/Q multipole correction.
• mDFTB2/3 periodic calculations are intentionally rejected in the current public API.

Current status

Version 0.1.0 is an initial development release focused on correctness, API stabilization, and validation coverage.

Performance optimization has not yet been a primary target. Some code paths are slower than expected. Future releases will focus on reducing memory use, improving batching, and optimizing critical C++ kernels.

This package should currently be treated as a prototype rather than a production-optimized DFTB engine.

Install

python -m pip install . --no-build-isolation --no-deps

For ASE examples, install the optional ASE dependency as well:

python -m pip install '.[ase]' --no-build-isolation

Minimal public API use

import numpy as np
from pydftb_torch import DftbCalculator

z = np.array([8, 1, 1], dtype=np.int64)
positions = np.array(
    [
        [0.000000, 0.000000, 0.000000],
        [0.957200, 0.000000, 0.000000],
        [-0.239987, 0.927297, 0.037500],
    ],
    dtype=np.float64,
)

calc = DftbCalculator(skf_dir="/path/to/skf", order=2)
energy_ev = calc.energy(z, positions)
forces_ev_angstrom = calc.forces(z, positions)
hessian_ev_angstrom2 = calc.hessian(z, positions)

Periodic public API use

import numpy as np
from pydftb_torch import DftbCalculator

z = np.array([8, 1, 1], dtype=np.int64)
positions = np.array(
    [
        [3.250000, 3.100000, 3.700000],
        [4.207200, 3.100000, 3.700000],
        [3.010013, 4.027297, 3.737500],
    ],
    dtype=np.float64,
)
cell = np.diag([8.0, 8.5, 9.0])

calc = DftbCalculator(skf_dir="/path/to/skf", order=3)
energy_ev = calc.energy(z, positions, cell_angstrom=cell)
forces_ev_angstrom = calc.forces(z, positions, cell_angstrom=cell)
hessian_ev_angstrom2 = calc.hessian(z, positions, cell_angstrom=cell)
stress_ev_angstrom3 = calc.stress(z, positions, cell_angstrom=cell)

mDFTB2/3 molecular use (Experimental)

The mDFTB implementation in this package is the molecular q/d/Q multipole correction on top of SCC-DFTB2 or DFTB3. Here q/d/Q means atomic charge, atomic dipole, and atomic quadrupole moments. The public helper mdftbn_options(model) prepares a DftbCalculator keyword dictionary with order=2 or order=3, multipole_model="mdftb2_moment_correction", multipole_max_order=2, traceless quadrupoles, and conservative SCC settings.

The standard mDFTB examples use self-consistent moment feedback: multipole_self_consistent_hamiltonian=True. In this mode the q/d/Q correction is inserted into the Hamiltonian during the SCC cycle, so the electronic state, energy, forces, and molecular Hessian are all evaluated on the same self-consistent mDFTB fixed point.

import numpy as np
from pydftb_torch import (
    DftbCalculator,
    build_multipole_integral_table_from_skf,
    mdftbn_options,
)

skf_dir = "/path/to/skf"
symbols = ("O", "H")
z = np.array([8, 1, 1], dtype=np.int64)
positions = np.array(
    [
        [0.000000, 0.000000, 0.000000],
        [0.957200, 0.000000, 0.000000],
        [-0.239987, 0.927297, 0.037500],
    ],
    dtype=np.float64,
)

# Build an auditable AO multipole-integral table from the homonuclear SKF files.
# The table supplies the on-site dipole and quadrupole factors required by the
# q/d/Q correction; no fitted table is bundled with the package.
table = build_multipole_integral_table_from_skf(skf_dir, symbols)

options = mdftbn_options(
    2,  # use 3 for mDFTB3
    dtype="float64",
    multipole_self_consistent_hamiltonian=True,
    multipole_auto_integral_table=False,
    **table.as_calculator_kwargs(),
)
calc = DftbCalculator(skf_dir=skf_dir, **options)
energy_ev = calc.energy(z, positions)
forces_ev_angstrom = calc.forces(z, positions)
hessian_ev_angstrom2 = calc.hessian(z, positions)

Important options for mDFTB are:

• model in mdftbn_options(model): 2 gives mDFTB2, i.e. the q/d/Q correction on an SCC-DFTB2 reference. 3 gives mDFTB3, i.e. the same multipole correction on a DFTB3 reference.
• multipole_max_order: 0 keeps only q-like terms, 1 includes q+d, and 2 includes q+d+Q. The default from mdftbn_options is 2.
• multipole_self_consistent_hamiltonian=True: feed the q/d/Q correction back into the Hamiltonian during SCC. The mDFTB examples use this self-consistent setting by default. Because this is a harder nonlinear SCC problem than ordinary DFTB2/3, conservative mixer settings from mdftbn_options() are kept.
• multipole_auto_integral_table: when True, the calculator tries to infer the required table from available homonuclear SKF files. Passing an explicit MultipoleIntegralTable with multipole_auto_integral_table=False is more reproducible and is preferred for examples, tests, and publications.

mDFTB2/3 support is molecular only in the current public API. Passing cell_angstrom to energy, forces, hessian, or the ASE calculator raises NotImplementedError. Molecular mDFTB Hessians are assembled as H = -dF/dR through PyTorch autograd over the molecular force route, not through the closed-form DFTB1/2/3 Hessian kernel. Periodic mDFTB stress and periodic mDFTB Hessians are therefore out of scope at present.

For ASE optimization, pass the same option dictionary to ase_calculator:

from ase import Atoms
from ase.optimize import BFGS
from pydftb_torch import ase_calculator

atoms = Atoms(
    "OH2",
    positions=[
        [0.000000, 0.000000, 0.000000],
        [0.957200, 0.000000, 0.000000],
        [-0.239987, 0.927297, 0.037500],
    ],
)
atoms.calc = ase_calculator(skf_dir=skf_dir, **options)
BFGS(atoms).run(fmax=0.05, steps=20)

ASE structure optimization

from ase import Atoms
from ase.optimize import BFGS
from pydftb_torch import ase_calculator

atoms = Atoms(
    "OH2",
    positions=[
        [0.000000, 0.000000, 0.000000],
        [0.957200, 0.000000, 0.000000],
        [-0.239987, 0.927297, 0.037500],
    ],
)
atoms.calc = ase_calculator(skf_dir="/path/to/skf", order=2)
BFGS(atoms).run(fmax=0.05, steps=20)

More complete examples are in examples/:

export PYDFTB_TORCH_SKF_DIR=/path/to/3ob-3-1
python examples/api_usage_matrix.py
python examples/ase_structure_optimization_matrix.py

References

• D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, “Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon,” Phys. Rev. B 51, 12947 (1995). DOI: 10.1103/PhysRevB.51.12947.
• M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, “Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties,” Phys. Rev. B 58, 7260 (1998). DOI: 10.1103/PhysRevB.58.7260.
• M. Gaus, Q. Cui, and M. Elstner, “DFTB3: Extension of the self-consistent-charge density-functional tight-binding method (SCC-DFTB),” J. Chem. Theory Comput. 7, 931-948 (2011). DOI: 10.1021/ct100684s.
• V.-Q. Vuong, T. W. Williams, J. Li, Q. Cui, and I. S. Ufimtsev, “Multipole expansion of atomic electron density fluctuation interactions in the density-functional tight-binding method,” J. Chem. Theory Comput. 19, 8152-8171 (2023). DOI: 10.1021/acs.jctc.3c00778.
• A. H. Larsen et al., “The atomic simulation environment — a Python library for working with atoms,” J. Phys.: Condens. Matter 29, 273002 (2017). DOI: 10.1088/1361-648X/aa680e.