Constrained DFT

Charge- and spin-constrained DFT, and the electron-transfer couplings that follow from it.

Run it

Python and the Rust library; there is no CLI section.

  • run_cdft(mol, basis_set, constraints, functional=None, ...) returns a CdftResult: a constrained UHF solve, or UKS when functional names a libxc functional other than "HF" (None and "HF", any case, give UHF).
  • CdftConstraint(atoms, target, kind="charge") defines one fragment constraint. atoms are 0-based atom indices. target is the electron population on the fragment, the Becke-weighted trace \( \mathrm{Tr}[W D] \), not a net charge: \( N_\alpha + N_\beta \) for kind="charge" and \( N_\alpha - N_\beta \) for kind="spin". A neutral He atom has a charge population of 2.0; He⁺ has 1.0.
  • cdft_coupling(state_a, state_b) returns a CdftCouplingResult with the Wu–Van Voorhis coupling h_ab, the determinant overlap s_ab and the two diabat energies e_a, e_b. The raw element uses each state's free energy F = E + λN, which makes the coupling independent of a constant shift of the constraint operator.

In Rust the entry points are ferric_scf::cdft_driver::solve_cdft_uhf and ferric_scf::cdft_coupling::coupling_hab; see the Rust API.

This example reproduces the HeNe⁺ constrained solution that crates/ferric-scf/tests/cdft_outer_loop.rs pins (E = −130.4021906 Ha, λ = −2.754 Ha per electron), with the He fragment held at 2.0 electrons:

import ferric

# HeNe+ doublet at 2.0 Å; hold the He atom (index 0) at 2.0 electrons.
mol = ferric.Molecule.from_xyz_string(
    "2\nHeNe+\nHe 0.0 0.0 0.0\nNe 0.0 0.0 2.0\n", 1, 2
)
r = ferric.run_cdft(
    mol,
    ferric.BasisSet.bundled("def2-svp"),
    [ferric.CdftConstraint([0], 2.0, kind="charge")],
    guess="hcore",
    level_shift=0.5,
    scf_max_iter=400,
    lambda_tol=1e-5,
    max_outer=40,
    stability_descent=False,
    grid_radial=99,
    grid_angular=302,
)
print(r.converged)
print(f"{r.energy:.6f}")
print(f"{r.populations[0]:.5f}")

True
-130.402190
2.00000

r.energy is the ordinary UHF/UKS energy at the constrained density, without the constraint term. r.lambdas holds the multipliers (Ha per electron) and r.weight_matrix(i) the AO weight operator of constraint i.

What to know before using it:

  • A returned result has a converged λ loop. When the outer loop exceeds max_outer, run_cdft raises RuntimeError. The inner SCF at the final λ can still be unconverged, so check r.converged: it is true only when the inner SCF converged and every constraint is met to lambda_tol.
  • A constraint is applied at least once before it counts as met. The loop starts at λ = 0, where the state is the unconstrained one. If its population is within lambda_tol of the target but more than 1e-3 electrons away, the loop still takes a Newton step rather than returning the unconstrained state as if it were constrained. Within 1e-3 the unconstrained state is returned with λ = 0: that is the measured spread of one state's population across grids, level shifts and guesses (up to 8.9e-5), with a margin. A loose tolerance can therefore return a state that is only partly constrained, so check r.max_constraint_error (or the per-constraint r.residuals) to see how far from the target it is. The stability-descent restarts follow the same rule.
  • stability_descent defaults to True here (in run_uhf it defaults to False). After the λ loop converges, a constrained saddle point is followed downhill to a lower state that still meets the constraint. The example turns it off to reproduce the pinned solution; with it on, this HeNe⁺ case ends about 0.0245 Ha lower. The descent is skipped for a KS reference.
  • The weight grid defaults to 99 × 302. It must resolve populations below lambda_tol; 75 × 110 resolves them only to about 1e-4. Setting grid_radial or grid_angular uses that grid for the XC quadrature too.
  • The UKS path is the one compared against another code. Constrained UKS/PBE matches NWChem 7.2.2 (see Accuracy). The UHF path (functional=None or "HF") has no comparison for energies along a target scan or for couplings: NWChem's standalone SCF (Hartree–Fock) module has no cDFT, because the Becke weight operator is built on the XC grid of its DFT module. Constrained UHF state selection on HeNe⁺ and LiH⁺ is compared through that DFT module, with Hartree–Fock exchange and a vanishing Slater term extrapolated to zero (ferric-scf/tests/validation_cdft_state.rs).
  • One constraint is the tested case. For one constraint the outer loop is a Newton step kept inside a sign-change bracket. It backs off from a λ whose inner SCF does not converge, and it discards a finite-difference derivative whose two inner solves landed in different SCF states. With several constraints the outer loop is a plain k × k Newton step without these safeguards.
  • cdft_coupling has strict preconditions and raises ValueError when one fails. Each state carries exactly one kind="charge" constraint, both are converged, and both come from the same molecule, geometry, basis, charge and multiplicity, and the same Hamiltonian: functional, df_j_aux/df_k_aux, k_builder, XC grid, point charges and external field. Two states that are the same determinant (\( |S_{ab}| \to 1 \)) also raise. The sign of h_ab is a determinant-phase convention; compare \( |H_{ab}| \).

Accuracy

Constrained UKS/PBE energies and multipliers are compared against NWChem 7.2.2 cdft ... pop becke on LiH, HF and H2O⁺ (charge and spin constraints) at 6-31G and def2-SVP (ferric-scf/tests/validation_cdft.rs): E(N) − E_unconstrained to 2.0e-7 Ha and λ to 1.1e-6 against NWChem's grid limit, and dE/dN = −λ to 2.3e-12 Ha (this identity is checked at def2-SVP, on the first target of each constraint kind); numbers are on Capabilities and validation. No external reference value is stated for UHF-cDFT energies apart from the HeNe⁺ and LiH⁺ states in validation_cdft_state.rs. On HeNe⁺/def2-SVP at R = 2.0 Å the hole-on-Ne state (N_He = 2) is NWChem's lower state once NWChem's Becke partition is given ferric's He/Ne radius ratio (0.667; NWChem's own is 0.700). Against NWChem's native partition it sits 0.563 eV higher. The offset comes from the partition: with NWChem's radii patched into ferric's radius table (a one-off measurement, not a test), ferric reaches −130.447404 Ha against NWChem's native −130.447401 Ha. The He₂⁺ coupling ingredients (determinant overlap, one- and two-electron transition elements, |V|) match NWChem's et module to ≤ 6e-11 Ha on NWChem's own determinants, and the KS diabats and couplings match end to end when started from NWChem's λ, except def2-SVP at 3.50 Å, where inner SCF solves within 1e-9 in λ of the root do not converge in 100 iterations (see Capabilities and validation). These diabats are validated only from NWChem's λ, set through the Rust RhfConfig::cdft_lambda_init; run_cdft starts from λ = 0, where the symmetric pair is delocalized and the inner SCF is bistable. The tests also check the coupling kernel on synthetic matrices and He₂⁺ identities (ferric-scf/tests/cdft_coupling.rs). On HeNe⁺, an asymmetric pair, the same file checks that both diabats carry the requested charge, that the hole-on-He state is reached through either fragment (3.9e-6 Ha apart), that the hole prefers Ne, the lower-IP atom (gap 1.7646 eV with the default stability descent; 1.0975 eV is a saddle point), and that H_ab is invariant under W → W − cS. These run at R = 2.0 Å only: from the hcore start, at 2.5, 3.5 and 4.0 Å the outer loop does not converge or returns the unconstrained state, and from MINAO at 2.5 Å the two routes to the hole-on-Ne state disagree by 6e-4 Ha. |H_ab| is not asserted on HeNe⁺; cdft_coupling_hene.rs probes it over a distance series. The tests also check exact identities on LiH/def2-SVP (cdft_uhf.rs): the constraint is satisfied, λ = 0 reproduces plain UHF, and the constraint composes with an external point charge. The Python tests (crates/ferric-python/tests/test_cdft.py) rerun those configurations through the bindings and check cdft_coupling against an independent transition-density construction. Apart from the HeNe⁺ and LiH⁺ states above, treat UHF-cDFT energies as unvalidated against other codes; see Capabilities and validation.

The response connection

A cDFT constraint couples a Lagrange multiplier \( \lambda \) to a fragment-weighted density operator. The derivative

\[ \frac{\partial N}{\partial \lambda} \]

is the charge moved per unit constraint potential, a susceptibility. cDFT therefore probes the same object as RPA and GW and attenuated MP2, through a different coupling.

Implementation

  • Fragment charge and spin constraints via a grid-Becke weight operator
  • A nested Lagrange-multiplier solve (Wu–Van Voorhis): an inner SCF at fixed \( \lambda \), an outer Newton iteration on \( \lambda \) itself

Each outer step is a full converged inner solve, so cDFT costs several plain SCFs.

Electron-transfer coupling

The coupling \( H_{ab} \) between two charge-localized diabatic states follows from a non-orthogonal determinant overlap, computed via Löwdin biorthogonalization. It is the matrix element governing electron-transfer rates in Marcus theory.

Tolerance and coupling error

The λ convergence tolerance (lambda_tol in Python, cdft_lambda_tol in Rust, default 1e-5 electrons) affects the coupling: a loosely converged \( \lambda \) produces diabatic states that miss their targets, and \( H_{ab} \) inherits that error. Tighten it before trusting a coupling. The exception is a constraint whose population barely responds to \( \lambda \), such as He₂⁺ at the localized-hole plateau with a UHF reference: on He₂⁺/def2-SVP UHF at 2.5–3.5 Å the outer loop does not converge in 40 iterations at 1e-5, while at 1e-2 (what the tests use) it stops after 3 with λ ≈ 0.91 and residuals of 1.6e-3 to 7.8e-3 electrons, so check r.max_constraint_error before trusting the state. With a UKS/PBE reference the He₂⁺ diabats do converge from λ = 0 at 1e-10 (validation_cdft_et.rs).

Cite

Wu & Van Voorhis 2006 (electron-transfer coupling from cDFT); Becke 1988 (the fragment weight partition). Full entries in References.