Periodic KS-DFT

Switch from HF to DFT by using PeriodicKSOptions instead of PeriodicSCFOptions and calling the KS dispatcher run_rks_periodic_scf (the recommended user-facing entry point — mirrors run_rhf_periodic_scf from tutorial 4):

import numpy as np
from vibeqc import (
    Atom, BasisSet, KPoints, PeriodicSystem, PeriodicKSOptions,
    run_rks_periodic_scf,
)

unit_cell = [Atom(1, [0, 0, 0]), Atom(1, [0, 0, 1.4])]
sysp = PeriodicSystem(
    dim=1,
    lattice=np.diag([4.0, 30.0, 30.0]),
    unit_cell=unit_cell,
)
basis = BasisSet(sysp.unit_cell_molecule(), "pob-tzvp")
kpoints = KPoints.monkhorst_pack(sysp, [6, 1, 1])

opts = PeriodicKSOptions()
opts.functional = "PBE"
opts.lattice_opts.cutoff_bohr = 12.0
opts.lattice_opts.nuclear_cutoff_bohr = 40.0
opts.conv_tol_energy = 1e-10
# Default: opts.lattice_opts.coulomb_method = CoulombMethod.DIRECT_TRUNCATED.
# Flip to CoulombMethod.EWALD_3D for any tight 3D bulk — see "Coulomb
# method" below.

result = run_rks_periodic_scf(sysp, basis, kpoints, opts)

print(f"E_KS per cell  = {result.energy:.10f}")
print(f"  E_Coulomb    = {result.e_coulomb:.10f}")
print(f"  E_xc         = {result.e_xc:.10f}")
print(f"  E_HF_exch    = {result.e_hf_exchange:.10f}  (0 for pure DFT)")

For Γ-only calculations there is the matching run_rks_periodic_gamma_scf(sysp, basis, opts) (no kmesh argument) — internally it just dispatches the multi-k driver at a [1,1,1] mesh, since vibe-qc’s C++ stack doesn’t ship a separate Γ-only KS entry point.

Functional support

All 500+ libxc functionals work in principle. Today’s validation covers:

  • LDA (Slater + VWN5)

  • PBE — pure GGA

  • BLYP — pure GGA

Hybrid functionals (B3LYP, PBE0, …) work in the code path — the exchange_scale parameter on build_fock_2e_real_space routes the HF-exchange fraction through the lattice sum — but have not yet been validated end-to-end against CRYSTAL reference numbers.

Coulomb method

The KS dispatcher (run_rks_periodic_scf / run_rks_periodic_gamma_scf) inspects opts.lattice_opts.coulomb_method exactly like the RHF dispatcher:

  • CoulombMethod.DIRECT_TRUNCATED (default) — routes to run_rks_periodic (the original KS C++ driver, conditionally convergent in 3D). Fine for 1D / 2D and for 3D molecular-limit cells with lattice \(\gtrsim 15\) bohr.

  • CoulombMethod.EWALD_3D — full Gaussian-charge Ewald, unconditionally convergent. Shipped end-to-end:

    • Γ-only via run_rks_periodic_gamma_scf (Phase 15c-1)

    • multi-k via run_rks_periodic_scf (Phase 15c-2)

    • open-shell counterparts run_uks_periodic_gamma_scf / run_uks_periodic_scf (Phase 15c-3) — see tutorial 3 for the molecular open-shell story; the periodic UKS dispatcher follows the same pattern.

For any quantitative 3D bulk DFT (lattice \(\lesssim 12\) bohr, metal oxides, dense ionic crystals), use EWALD_3D. The periodic Becke partition takes care of the other tight-cell pathology — DFT integration weights — and the two flags compose: turn both on for tight crystals.

opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D
opts.use_periodic_becke = True       # the integration-weight fix
opts.becke_image_radius_bohr = 10.0

See the Ewald user guide for the ω-invariance story and tunable parameters; tutorial 23 for the periodic Becke partition.

Grid quality

Reuse the molecular-DFT grid controls:

opts.grid.n_radial = 99
opts.grid.n_theta  = 29
opts.grid.n_phi    = 58

The grid is built over the unit-cell atoms with the molecular Becke partition scheme. For tight unit cells where image-atom Voronoi cells intrude on the reference cell, expect small (<1e-4 Ha) boundary effects until the periodic Becke partition lands in Phase 12f.

Tip

SCF won’t converge? examples/debug/ ships nine focused diagnostic scripts. For DFT-specific cases the most relevant are scf_easy_periodic.py (Ne / Ar / diamond C must converge — known-good baseline), scf_sad_vs_hcore.py (does the v0.6.1 SAD initial guess rescue HCORE divergences?), and scf_iteration_recorder.py (4-panel forensics plot of any divergent trajectory). Tutorial 24 walks through level shift / DIIS / damping / smearing aids in detail.

Theory

The Kohn-Sham equations at a single \(\mathbf{k}\)-point look identical to the molecular KS equations of tutorial 02:

\[ \mathbf{F}^\text{KS}(\mathbf{k}) \, \mathbf{C}(\mathbf{k}) = \mathbf{S}(\mathbf{k}) \, \mathbf{C}(\mathbf{k}) \, \boldsymbol{\varepsilon}(\mathbf{k}), \]

with \(\mathbf{F}^\text{KS}(\mathbf{k})\) assembled from its real-space lattice sum via the Bloch transform of tutorial 04. The core ideas for both machines — the Bloch summation, the Monkhorst-Pack k-integration, the direct-truncation caveat for the 3D Coulomb sum — are identical to periodic HF and documented there.

What’s new in periodic KS-DFT is the exchange-correlation contribution. The density at a point is now Bloch-summed:

\[ \rho(\mathbf{r}) = \frac{1}{N_\mathbf{k}} \sum_\mathbf{k} \sum_n^\text{occ}(\mathbf{k}) |\psi_{n\mathbf{k}}(\mathbf{r})|^2, \]

where the sum is over \(N_\mathbf{k}\) sampled k-points and over the occupied bands at each. The numerical integration of \(V^\text{xc}_{\mu\nu}(\mathbf{T})\) on the DFT grid then needs to account for periodic images of the basis functions. vibe-qc’s current implementation builds the grid from the unit-cell atomic centers using the standard molecular Becke partition. This is accurate for well-separated periodic cells (the molecular-limit regime) but produces small partition-boundary artefacts for tight cells — the periodic Becke partition (Phase 12f on the roadmap) is the fix.

Hybrid DFT in periodic systems

Hybrid functionals (B3LYP, PBE0, …) mix a fraction \(\alpha\) of HF exchange into the XC energy. In a periodic system this translates to an \(\mathcal{O}(N_\text{bf}^4)\) lattice sum of the exchange-type two-electron integrals

\[ K_{\mu\nu}(\mathbf{T}) = \sum_{\lambda\sigma, \mathbf{T}', \mathbf{T}''} P_{\lambda\sigma}(\mathbf{T}'' - \mathbf{T}') \, (\mu \mathbf{0} \, \sigma \mathbf{T}' \, | \, \nu \mathbf{T} \, \lambda \mathbf{T}''), \]

instead of the \(\mathcal{O}(N_\text{bf}^3)\) sum of local XC matrix elements for pure DFT. That’s 10–100× more expensive — the reason solid-state DFT defaults to pure GGAs (PBE, BLYP). vibe-qc’s build_fock_2e_real_space(..., exchange_scale=α) routes the hybrid exchange through the same lattice-summed ERI machinery used by HF.

Basis-set considerations for solids

Two basis-set issues specific to periodic calculations, summarized from tutorial 13:

  1. Linear dependencies. Diffuse molecular basis functions (aug-cc-pVXZ, 6-311++G**) overlap strongly across neighboring periodic cells, pushing the overlap matrix toward singularity. vibe-qc’s canonical-orthogonalisation path handles mild cases; stronger contamination needs a redesigned basis.

  2. Solid-state basis families. pob-TZVP (Peintinger-Vilela Oliveira-Bredow 2013) is the standard starting point — it’s a systematic re-derivation of TZVP with the smallest-exponent primitives either removed or replaced, and ships with vibe-qc.

Resources

~250 MB peak RAM, ~2 s on one core (Apple M2 baseline) for the H₂-chain PBE / 6-k example. Add ~1 s per k-point for the XC grid build inside each iteration (75-radial × 302-angular Lebedev shells per atom by default). Use the periodic Becke partition on tight cells to avoid the over-counting of integration weights that hurts molecular-Becke DFT.

References

  • Periodic KS-DFT. The theoretical framework is Kohn-Sham 1965 (see tutorial 02 for full citation) applied to crystalline orbitals per Pisani-Dovesi 1980 (see tutorial 04).

  • Periodic Becke partition. [CITATION MISSING: no standardised published derivation of the Becke atomic partitioning scheme extended to periodic systems — the CRYSTAL implementation is documented in the CRYSTAL manual but lacks a clean primary reference.]

  • Solid-state basis sets (pob-TZVP). M. F. Peintinger, D. V. Oliveira, and T. Bredow, “Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations,” J. Comput. Chem. 34, 451 (2013).

  • CRYSTAL reference. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith- Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Brandenburg, and R. Dovesi, “CRYSTAL23: a program for computational solid state physics and chemistry,” J. Chem. Theory Comput. 19, 6891 (2023).