QC

QC runs orbital-basis density-functional calculations through the native qc.rs engine. It accepts Atomli structures and ASE-provided Atoms objects.

Molecular calculations

With no settings, QC constructs a molecular calculation:

from atomli.build import molecule
from atomli.calculators.qc import QC
from atomli.optimize import BFGS

atoms = molecule("H2O")
atoms.calc = QC(xc="PBE", basis="def2-SVP")

energy = atoms.get_potential_energy()
forces = atoms.get_forces()
BFGS(atoms).run(fmax=1e-3)

Water relaxed with PBE/def2-QZVP: O-H 0.9685 A, H-O-H 104.18 degrees, against 0.9572 A and 104.52 degrees measured.

Basis-set comparison

The basis is the first decision and the one with the widest cost range, so here it is measured rather than described. Water relaxed with PBE at three basis sets, against the spectroscopic geometry (Benedict, Gailar and Plyler, J. Chem. Phys. 24, 1139 (1956)):

Basis	Energy	O-H	H-O-H	Single point
def2-SVP	-2075.48 eV	0.9746 Å (+0.0174)	102.05° (−2.47)	42 ms
def2-TZVPP	-2078.42 eV	0.9690 Å (+0.0118)	103.96° (−0.56)	61 ms
def2-QZVP	-2078.60 eV	0.9685 Å (+0.0113)	104.18° (−0.34)	144 ms

Deviations from experiment in brackets. The bond angle is the sensitive quantity: def2-SVP is 2.5° too narrow, which def2-TZVPP mostly repairs and def2-QZVP finishes, at 3× the cost of the smallest. The bond length converges much sooner and stays about 11 milliangstrom long, which is PBE itself rather than the basis: no further basis functions will remove it.

The total energies drop by 2.9 eV from def2-SVP to def2-TZVPP and only 0.18 eV for the next step. Absolute total energies from different bases are not comparable with each other; only differences computed in one basis are.

r2SCAN on the same molecule and basis returns -2079.58 eV in 67 ms, so the meta-GGA is not meaningfully more expensive per single point here.

Common molecular bases include def2-SVP, def2-TZVPP, and def2-QZVP. Molecular calculations use RI-J density fitting by default (density_fitting="default"). RI-J is the fast production path. Pass density_fitting="none" to select the conventional Coulomb build. Periodic mode requires "none" and uses it by default.

auxiliary_basis pins the fitting set that RI-J expands the density in, where the default lets the engine choose the set matching the functional:

atoms.calc = QC(xc="PBE", basis="def2-SVP",
                auxiliary_basis="def2-universal-jkfit")

It changes the Coulomb build, not just the record — the same functional and orbital basis fitted against a different auxiliary set is a different number. Because it is the fitting set, it is molecular-only, and passing it with density_fitting="none" or with periodic settings raises rather than being ignored: neither of those runs a fitted Coulomb build at all. Set or not, it reads back from calc.parameters["auxiliary_basis"], so a recorded run states that nothing was pinned instead of leaving the question open.

Set the total charge and number of unpaired electrons on Atoms, or pass charge and unpaired overrides to the calculator. qc.rs validates electron count and spin parity before SCF.

Exchange-correlation functionals

The xc argument selects the functional. Molecular and periodic mode accept PBE and r2SCAN. Periodic mode also accepts LDA/VWN with the matched GTH-PADE pseudopotential. SKALA 1.1 is molecular-only.

PBE

PBE is a generalized-gradient (GGA) functional and the QC default. Choose it for general molecular and periodic work when speed matters more than meta-GGA accuracy.

atoms.calc = QC(xc="PBE")

r2SCAN

r2SCAN is a meta-GGA functional. It adds the kinetic-energy density to the GGA form. Choose it when energetics matter more than the extra cost over PBE.

atoms.calc = QC(xc="r2SCAN")

Single-point energies are solid. Geometry optimization on this build is not: relaxing water at def2-TZVPP for 100 steps against an fmax of 0.001 eV/Å finishes at 0.64 eV/Å, three orders of magnitude above the requested tolerance, having wandered to an O-H length of 0.9500 Å. The same optimizer reaches 7.9e-4 eV/Å on PBE with the same basis. Check the residual force after any r2SCAN relaxation rather than trusting that the optimizer stopped because it converged:

import numpy as np
optimizer = BFGS(atoms)
optimizer.run(fmax=1e-3)
assert np.abs(atoms.get_forces()).max() < 1e-3

SKALA 1.1

The SKALA 1.1 neural exchange-correlation functional runs through the same calculator. Select it, and the pinned checkpoint (2.3 MiB) downloads automatically on first use. The checkpoint is then cached locally:

atoms.calc = QC(xc="SKALA-1.1", basis="def2-SVP")

The pinned SKALA 1.1 rev1 checkpoint is downloaded from the mlip-models GitHub Release. Atomli verifies the file against its pinned SHA-256 before use. Pass skala_checkpoint="skala-1.1.fun" to use a local copy. Pass download=False to require a cache hit.

Dispersion correction

Combine QC with DFTD3 through ASE’s SumCalculator to add a dispersion correction.

GPU calculations

device="wgpu" evaluates the exchange-correlation potential, and the RI-J Coulomb build, on the GPU through one portable path: Metal on Apple silicon, Vulkan or DX12 elsewhere. When device is omitted, QC runs a lightweight readiness probe and selects GPU only when the packaged runtime opens on real hardware; otherwise it falls back to CPU. Use device="cpu" to pin host execution or for host-only capabilities.

atoms.calc = QC(xc="PBE", basis="def2-SVP", device="wgpu")

Supported on PBE, r2SCAN, and SKALA-1.1. The energy is the same calculation, not an approximation of it.

Measured on the machine that built this page, water at def2-SVP returns -2075.4815 eV on both paths, agreeing to 2.9e-10 eV, which is floating-point ordering rather than method. That single point took 155 ms on the GPU against 42 ms on the host: a three-atom molecule is far too small to pay back the fixed device setup cost.

Three limits are worth knowing before you reach for it.

Closed-shell only. The semilocal device backend implements RKS. An open-shell system on device="wgpu" raises rather than silently falling back to the host:

>>> QC(xc="PBE", device="wgpu", unpaired=2)
ValueError: ... the existing semilocal device backend implements RKS only

Run those on device="cpu", or use SKALA-1.1, which has its own device route.

Nuclear gradients stay on the host. get_forces() works on device="wgpu" and returns host-computed derivatives, so a forces calculation is only partly accelerated. Single-point energies see the whole benefit.

Some capabilities stay host-only. Periodic SKALA, periodic DFT+U, and analytical molecular Hessians require device="cpu". Periodic PBE/r2SCAN have their native GPU route; unsupported combinations raise explicitly rather than silently falling back.

A GPU build is required. The published wheels are built with the GPU path included, so device="wgpu" is available in a normal pip install atomli.

Whether the GPU is faster depends on the system. It has a fixed setup cost per calculation, so small molecules can be slower on the GPU than on eight CPU threads — the benchmarks show both cases measured.

Periodic calculations and stress

Periodic mode is explicit and requires a nondegenerate cell with PBC on all three axes:

from ase.build import bulk
from atomli.calculators.qc import QC

atoms = bulk("Si", "diamond", a=5.43)
atoms.calc = QC(settings={"periodic": True})

energy = atoms.get_potential_energy()
forces = atoms.get_forces()
stress = atoms.get_stress()
stress_tensor = atoms.calc.get_stress(atoms, voigt=False)

A primitive silicon cell at gth-szv-molopt-sr with 2×2×2 k-points returns -210.30 eV in 0.5 s on this machine, which is the cheapest honest periodic calculation available here and a useful floor to scale from. Single-zeta is a starting basis, not a production one, so treat that number as a cost measurement rather than an accuracy claim.

Stress is returned in eV/ų. ASE’s default Voigt order is [xx, yy, zz, yz, xz, xy]; voigt=False returns a symmetric 3 × 3 tensor. Molecular QC does not advertise stress.

Periodic settings

Key	Default	Meaning
periodic	required	True selects periodic QC
basis	"gth-dzvp-molopt-sr"	periodic orbital basis
pseudopotential	"gth-pbe"	use "gth-pade" for LDA/VWN
energy_cutoff	100.0	FFT-grid cutoff in Hartree
mesh	automatic	positive [nx, ny, nz] override
kpoint_density	15.0	minimum kpt × length product in Å
kpoints	automatic	positive [kx, ky, kz], or a full mesh identity
precision	1e-8	integral and grid target
derivative_mode	"automatic"	"analytical" or "finite_difference" override
displacement	1e-3	centered force-difference step in Bohr
strain	1e-3	dimensionless centered stress-difference step
smearing	absent	occupation broadening, {"sigma": …} in Hartree
numerical_integrator	"knumint"	"multigrid" selects PySCF’s MultiGrid path
dft_u	absent	Hubbard +U as {element: {shell: U_eV}}

Automatic k-point counts use max(1, ceil(kpoint_density / lattice_length_angstrom)) independently for each reciprocal dimension. mesh overrides the cutoff-derived grid; kpoints overrides the density-derived sampling. Invalid or unknown settings raise errors instead of selecting a fallback.

k-point sampling

kpts= is the primary spelling, the same keyword ASE calculators such as GPAW, VASP, and Espresso use. A mesh tuple selects gamma-centered sampling and implies periodic=True:

atoms.calc = QC(xc="PBE", kpts=(2, 2, 2), settings={
    "periodic": True,
    "basis": "gth-szv-molopt-sr",
    "pseudopotential": "gth-pbe",
})

kpts={"size": (2, 2, 2), "gamma": False} is ASE’s mesh dict, where gamma=False asks for the shifted Monkhorst-Pack coordinates instead. kpts={"density": d} targets d k-points per inverse angstrom, the ase.calculators.calculator.kpts2sizeandoffsets density convention (the engine stores the equivalent length target 2π·d as kpoint_density). The remaining ASE forms — even, explicit k-point lists — are rejected with an error naming what is supported.

Mesh centering and symmetry

Gamma-centered and shifted Monkhorst-Pack coordinates are different point sets at any even dimension, so the dimensions alone do not say what was sampled. KMesh carries the whole identity, and is accepted by kpts= and settings["kpoints"] alike:

from atomli.calculators.qc import KMesh

atoms.calc = QC(xc="PBE", kpts=KMesh((4, 4, 4), centering="monkhorst_pack"),
                settings={"periodic": True, "basis": "gth-szv-molopt-sr",
                          "pseudopotential": "gth-pbe"})

Its fields are size, centering ("gamma" by default, matching PySCF’s make_kpts), symmetry ("none", "time_reversal", "space_group", "space_group_and_time_reversal"), wrap_around, and ao_symmetry; KMesh.fields() prints them with their documentation. The equivalent dict ({"size": [4, 4, 4], "centering": "monkhorst_pack"}) works wherever the class does.

The settings-dict keys above are the override spelling of the same thing: exactly one of kpts= and settings["kpoints"]/settings["kpoint_density"] may declare the sampling, and declaring both raises an error naming the two spellings. The resolved values read back from calc.parameters["settings"], where kpoints is always the full mesh identity — a plain [2, 2, 2] echoes back with the centering it was resolved to, so a recorded run never has to infer which Brillouin-zone points it sampled.

Metals and magnetic structure

A metal has no gap for integer occupations to fill against, so its SCF oscillates between degenerate states and never settles. settings["smearing"] broadens the occupations, which is what makes such a cell solvable. It is absent by default, because integer occupations are the right description of an insulator.

sigma is the width in Hartree and is required; method is "fermi" (the default) or "gaussian". Two further fields have no ASE spelling and live only here: chemical_potential fixes the Hartree chemical potential instead of optimizing it against the electron count each cycle, and fixed_spin=True holds the alpha and beta counts apart so the requested moment survives rather than being relaxed away by one common chemical potential. With smearing on, the reported energy is the Mermin free energy E − σS, which is the quantity the returned forces and stress differentiate.

Per-site initial moments follow the ASE convention, and the calculator reads them off the structure:

atoms = bulk("Fe", "bcc", a=2.87)
atoms.set_initial_magnetic_moments([2.0])

atoms.calc = QC(xc="PBE", settings={
    "periodic": True,
    "basis": "gth-szv-molopt-sr",
    "pseudopotential": "gth-pbe",
    "kpoints": [1, 1, 1],
    "mesh": [49, 49, 49],
    "smearing": {"sigma": 0.02, "method": "fermi"},
    "dft_u": {"Fe": {"3d": 2.0}},
})

Broadening is not always the whole answer. That cell needs the explicit mesh too: on the cutoff-derived grid the same SCF oscillates past 200 cycles, and it settles in 13 at [49, 49, 49]. An under-resolved real-space grid and a metallic spectrum fail the same way from the outside, so try both.

The declared moments seed the unrestricted SCF with the matching per-site alpha/beta split, which is the only way to reach a ferrimagnetic or antiferromagnetic solution; any nonzero declaration selects spin-polarized treatment even at zero net moment. Given alone, their sum sets the net spin and so has to be an integer — the physical 2.2 μB of bcc iron is a converged result, not something to seed with. Given alongside unpaired, the two must agree or the calculation raises naming both numbers.

SCF settings

scf= carries the numerical controls the self-consistent field runs under, on the molecular and the periodic path alike. It takes an Scf instance, an equivalent dict, or a preset name:

from atomli.calculators.qc import Scf

QC(xc="PBE", scf=Scf(max_cycles=200, level_shift=0.3))
QC(xc="PBE", scf={"max_cycles": 200, "level_shift": 0.3})   # the same block
QC(xc="PBE", scf="robust")                                  # a preset

"default" is the engine’s own block — 50 cycles, 1e-9 Hartree, no damping or level shift — and it is what an unspecified scf= resolves to. "robust" is the one to reach for when a cell will not settle: 256 cycles, Fock damping at 0.2, a 0.25 Hartree level shift, DIIS held back to cycle 3, and a core Hamiltonian initial guess. Scf.default() and Scf.robust() are the same two blocks as objects, so a preset can be taken as a starting point and adjusted.

An unknown key raises naming the fields that exist rather than being dropped, and every field is documented on the class itself:

Scf.fields()          # [(name, documentation), ...] for all 16 fields
Scf(level_shift=0.25).to_dict()

The sixteen are max_cycles, energy_tolerance, density_tolerance, fock_damping, level_shift, diis_start, diis_size, diis_damp, initial_guess, grid_schedule, coulomb_workers, density_fitting_metric_regularization, pyscf_nelec_compat, pyscf_cycle_compat, extrapolation, and routing_policy. The resolved block — every field, not only the ones named — reads back from calc.parameters["scf"].

ASE parameter aliases

kpts= is one of four keywords spelled the way ASE and GPAW spell them. Each is sugar for a native key, and each converts units on the way in, because the ASE vocabulary is in eV and the engine’s is in Hartree:

Keyword	Native key	Unit conversion
kpts=(2, 2, 2)	settings["kpoints"]	none; {"density": d} stores 2π·d as kpoint_density
smearing={"name": …, "width": w}	settings["smearing"]	w eV → sigma Hartree
convergence={"energy": e}	scf["energy_tolerance"]	e eV → Hartree
maxiter=n	scf["max_cycles"]	none

Exactly one spelling may declare each of these. Writing both raises a ValueError naming the two, rather than picking a winner: the ASE keyword is sugar for the native key, never an override of it. The resolved value reads back from calc.parameters under the native key either way, so a recorded run states what ran and not which vocabulary asked for it.

atoms.calc = QC(xc="PBE", kpts=(4, 4, 4), maxiter=200,
                convergence={"energy": 1e-6},              # eV
                smearing={"name": "fermi-dirac", "width": 0.1},  # eV
                settings={"periodic": True,
                          "basis": "gth-szv-molopt-sr",
                          "pseudopotential": "gth-pbe"})

smearing takes "fermi-dirac" or "gaussian" and nothing else; another name raises instead of being approximated with the nearer of the two. It needs a periodic mode to broaden — from settings["periodic"] or from kpts=, which implies it — and a molecular call raises. It reaches only sigma and method; the native block’s other two fields have no ASE spelling and stay in settings["smearing"].

convergence honors energy; any other key raises naming it, because a quietly dropped threshold would read back as a tighter calculation than the one that ran. maxiter is at least 1 — scf={"max_cycles": 0} is the spelling for evaluating the initial density without an SCF.

Both of those resolve into the scf block, and an scf= dict conflicts only on the key it names, so disjoint keys compose:

QC(xc="PBE", maxiter=99, scf={"level_shift": 0.25})   # both take effect
QC(xc="PBE", maxiter=99, scf={"max_cycles": 10})      # raises

An Scf instance or a preset name specifies the whole block, so either one conflicts with maxiter and convergence outright.

DFT+U

settings["dft_u"] names the Hubbard correction the way the literature quotes it — element, shell, U in eV — not as the AO indices the engine works in:

atoms.calc = QC(xc="PBE", kpts=(2, 2, 2), settings={
    "periodic": True,
    "basis": "gth-dzvp-molopt-sr",
    "pseudopotential": "gth-pbe",
    "dft_u": {"Fe": {"3d": 4.0}},
})

Every iron site in the cell takes the same 4 eV correction on its 3d shell. Several elements and several shells per element are allowed ({"Ni": {"3d": 6.0}, "O": {"2p": 1.0}}), and the mapping echoes back from calc.parameters["settings"]["dft_u"] exactly as written. DftU is the typed spelling of the same mapping, validated on construction and accepted wherever the dict is:

from atomli.calculators.qc import DftU

correction = DftU({"Fe": {"3d": 4.0}})
atoms.calc = QC(xc="PBE", kpts=(2, 2, 2), settings={
    "periodic": True,
    "basis": "gth-dzvp-molopt-sr",
    "pseudopotential": "gth-pbe",
    "dft_u": correction,
})

Each entry becomes a PySCF AO-label pattern ("Fe 3d") matched against the MINAO reference basis, which is what the projectors are built from — not the periodic basis the SCF runs in. So the shell has to exist in MINAO for that element. A pattern that selects nothing there is refused when the energy is requested, naming what it tried:

DFT+U selection "Fe 3q" matched no AOs in the "MINAO" reference

That check needs the cell’s own elements, so it cannot run at construction; a misspelt shell is an error rather than a run that quietly drops the +U term and reports plain DFT under a +U label. dft_u cannot be combined with numerical_integrator="multigrid", which the constructor refuses: PySCF’s MultiGrid path is wired for KRKS/KUKS, not KRKSpU.

SCF convergence troubleshooting

An unconverged SCF raises by default, on both the molecular and the periodic path. The error names the cycles spent, the threshold missed, and the settings to change — an unconverged energy is a point on an oscillating trajectory, not an observable.

QC(..., require_convergence=False) is the deliberate opt-out: it emits a UserWarning, returns the last iterate, and sets results["converged"] to False. It covers the energy only; forces and stress are derivatives of the converged state itself, so they are refused either way.

calc.converged and calc.scf_iterations report the last SCF programmatically (None before anything has run). Note that converged is necessary but not sufficient: an under-resolved grid can converge both engines to an agreeing, physically wrong state, so validate a physical observable when a fixture or workflow depends on the result.

Second derivatives

Molecular QC other than SKALA carries an analytical Hessian, shaped (n_atoms, 3, n_atoms, 3) in eV/Ų:

hessian = atoms.calc.get_hessian(atoms)

atomli.vibrations.Vibrations is the route to frequencies and normal modes. It takes that Hessian in one shot where it exists and displaces atoms where it does not, so SKALA and periodic cells work through the same call.

State reuse

The backend and available results remain attached across repeated property calls, ASE optimizers, and molecular dynamics. Atomic numbers, positions, cell, PBC, charge, and unpaired electrons participate in exact state comparison. Changing any of them invalidates the cached result; reset() also discards the retained backend.