Periodic HF on a 1D H chain

The simplest non-trivial periodic system: a one-dimensional chain of H₂ dimers.

import numpy as np
import vibeqc as vq
from vibeqc import (
    Atom, BasisSet, KPoints, PeriodicSystem, PeriodicSCFOptions,
    run_rhf_periodic_scf,
)

# Unit cell: H2 oriented along the chain, 4.0 bohr between unit cells.
unit_cell = [Atom(1, [0, 0, 0]), Atom(1, [0, 0, 1.4])]
lattice = np.diag([4.0, 30.0, 30.0])   # 1D axis on x; y, z are vacuum

sysp = PeriodicSystem(
    dim=1,
    lattice=lattice,
    unit_cell=unit_cell,
)

basis = BasisSet(sysp.unit_cell_molecule(), "pob-tzvp")
kpoints = KPoints.monkhorst_pack(sysp, [6, 1, 1])

opts = PeriodicSCFOptions()
opts.conv_tol_energy = 1e-10
opts.lattice_opts.cutoff_bohr = 12.0         # μν overlap cutoff
opts.lattice_opts.nuclear_cutoff_bohr = 40.0  # longer-range nuclear sum
# Default: opts.lattice_opts.coulomb_method = CoulombMethod.DIRECT_TRUNCATED.
# For quantitative 3D bulk, set CoulombMethod.EWALD_3D and the dispatcher
# routes to the Ewald backend automatically — see the
# [Ewald user guide](../user_guide/ewald.md).

result = run_rhf_periodic_scf(sysp, basis, kpoints, opts)
print(f"E_RHF per cell = {result.energy:.10f}  ({result.n_iter} iters)")

run_rhf_periodic_scf is the Coulomb-method dispatcher — the recommended user-facing entry point. It accepts a :class:vibeqc.KPoints object (or, for back-compat, a raw BlochKMesh) and inspects opts.lattice_opts.coulomb_method to route to the appropriate backend (run_rhf_periodic for direct truncation, the multi-k Ewald backend for EWALD_3D). For a Γ-only run, use run_rhf_periodic_gamma_scf(sysp, basis, opts) (no kmesh argument) or KPoints.gamma(sysp).

Other ways to build a k-mesh — vibeqc.KPoints

The KPoints builder fronts every common k-point input mode in a single API. The five most useful for everyday work:

# 1. Monkhorst–Pack (classical convention — auto-shifts even meshes by ½)
kp = KPoints.monkhorst_pack(sysp, [6, 1, 1])

# 2. Γ-centred — required for hex/trigonal cells (the classical MP shift
#    breaks 3-fold symmetry); the MP constructor will refuse non-zero
#    shifts on those spacegroups with an actionable error.
kp = KPoints.gamma_centred(sysp, [6, 6, 6])

# 3. Symmetry-reduced to the irreducible BZ (via spglib).
#    Requires `vq.attach_symmetry(sysp)` first.
vq.attach_symmetry(sysp)
kp = KPoints.monkhorst_pack(sysp, [4, 4, 4], symmetry=True)
# Equivalent two-step form: kp.symmetry_reduce()

# 4. Density-based auto-mesh (KPPRA / kspacing / VASP `Auto`).
#    Picks the mesh from a target k-point density rather than from
#    a fixed [N, N, N] shape.
kp = KPoints.from_kppra(sysp, kppra=1000)        # AFLOW / Curtarolo
kp = KPoints.from_kspacing(sysp, kspacing=0.04)   # Materials Project / ASE
kp = KPoints.auto(sysp, length=40)                # VASP Auto

# 5. Band path — auto-detects Bravais lattice via spglib, returns
#    the canonical Hinuma 2017 (HPKOT) k-path via seekpath.
band = KPoints.band_path(sysp)
print(' → '.join(label for _, label in band.labels))

The legacy free-function vq.monkhorst_pack(sysp, [N, M, L]) is still supported and returns a raw BlochKMesh — the periodic SCF entry points accept both. For new code, prefer KPoints: it carries the symmetry-reduction state, supports band paths, and gives an actionable error on the hex/trigonal MP-shift trap.

End-to-end demo across every mode (Γ-only, MP scan, IBZ reduction, KPPRA scan, HPKOT band path) on cubic Mg: examples/periodic/input-k-mesh-convergence.py. The full reference is user_guide/k_points.md.

What the options mean

  • PeriodicSystem.dim — number of periodic axes. dim=1 makes only lattice.col(0) a lattice vector; the other two columns are interpreted as vacuum with enough separation (30 bohr here) that images don’t overlap. dim=2 treats cols 0 and 1 as in-plane lattice vectors. dim=3 is full bulk.

  • lattice_opts.cutoff_bohr — the real-space distance beyond which (μν) pair integrals are dropped from the lattice sum. Bump it up until your energy stops changing.

  • lattice_opts.nuclear_cutoff_bohr — separate, usually larger, cutoff for the nuclear-attraction term (1/r tail decays more slowly than AO overlap).

Choosing a k-mesh

For the SCF to represent an insulator correctly, the mesh needs to sample enough of the Brillouin zone that the occupied manifold stays separated from the virtual. Start with [N, 1, 1] for 1D, [N, N, 1] for 2D, [N, N, N] for 3D, and increase N until the energy per cell stops changing. For small H-chain tests a 4×1×1 mesh is usually plenty; for covalent 3D crystals in TZ basis, 4×4×4 is a reasonable starting point.

What’s validated today

  • Molecular limit: a big-box (50 bohr) unit cell produces the same energy as molecular RHF — machine precision across dim ∈ {1, 2, 3} and multiple k-mesh sizes.

  • Bloch-sum machinery: kinetic-only folding equivalence between 1-cell × K k-points and K-cell × Γ, machine precision.

  • Non-trivial 1D SCF convergence at tight cell spacings.

3D bulk calculations with the direct-truncated Coulomb sum are cutoff-dependent (the infinite Coulomb series is conditionally convergent in 3D). For quantitative bulk results set opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D — the dispatcher then routes to vibe-qc’s full Saunders–Dovesi multipolar Ewald backend. See the Ewald user guide for details and the orthorhombic-cell constraint.

Tip

SCF won’t converge? examples/debug/ ships a kit of focused diagnostic scripts (Ne / Ar / diamond C known-good baselines, trajectory recorder, convergence-aid sweep, PySCF cross-check). Tutorial 24 walks through the level-shift / DIIS / SAD-guess aids the SCF exposes. Every dispatcher accepts progress=True (v0.5.1) for live per-iter logging.

Theory

Bloch’s theorem

For a Hamiltonian that commutes with every translation by a lattice vector \(\mathbf{T}\), the eigenfunctions can be labeled by a crystal momentum \(\mathbf{k}\) in the first Brillouin zone and a band index \(n\):

\[ \psi_{n\mathbf{k}}(\mathbf{r}) = e^{i \mathbf{k} \cdot \mathbf{r}} \, u_{n\mathbf{k}}(\mathbf{r}), \qquad u_{n\mathbf{k}}(\mathbf{r} + \mathbf{T}) = u_{n\mathbf{k}}(\mathbf{r}). \]

The lattice-periodic function \(u_{n\mathbf{k}}\) is expanded in a Bloch-summed basis: take each atom-centered Gaussian \(\chi_\mu\) in the reference unit cell and form the infinite symmetry-adapted combination

\[ \chi_{\mu\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{T}} e^{i \mathbf{k} \cdot \mathbf{T}} \, \chi_\mu(\mathbf{r} - \mathbf{T}). \]

These Bloch-summed AOs are the natural basis for crystalline orbitals in the Roothaan-style framework (CRYSTAL, vibe-qc’s periodic driver, QE’s hybrid-DFT implementation all do this).

The generalised eigenvalue problem at every k

Plug the Bloch-summed basis into the HF variational equation and the resulting matrices are \(\mathbf{k}\)-diagonal: for each \(\mathbf{k}\) in the sampled mesh, you get an independent generalised eigenvalue problem

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

where each k-space matrix is a Fourier transform of the real- space lattice-summed matrix:

\[ \mathbf{F}(\mathbf{k}) = \sum_{\mathbf{T}} e^{i \mathbf{k} \cdot \mathbf{T}} \, \mathbf{F}(\mathbf{T}), \qquad \mathbf{F}(\mathbf{T})_{\mu\nu} = \int \chi_\mu(\mathbf{r}) \, \hat{F} \, \chi_\nu(\mathbf{r} - \mathbf{T}) \, \mathrm{d}^3 r. \]

The SCF loop alternates between two representations: build the Fock contributions in real space (lattice-summed one- and two-electron integrals), transform to \(\mathbf{k}\)-space, diagonalize per \(\mathbf{k}\), inverse-transform the density back to real space. This “real-space build, \(\mathbf{k}\)-space diagonalize” pattern is what CRYSTAL pioneered for Gaussian crystalline orbitals.

Monkhorst-Pack k-sampling

Integrals over the Brillouin zone — including the density

\[ \rho(\mathbf{r}) = \sum_n^{\text{occ}} \int_\text{BZ} |\psi_{n\mathbf{k}}(\mathbf{r})|^2 \, \frac{\mathrm{d}^3 k}{(2\pi)^3} \]

— are replaced by weighted sums over a discrete mesh. The Monkhorst-Pack scheme constructs a uniform grid of k-points offset from the BZ origin, optionally reduced to the irreducible wedge by the space-group operations. monkhorst_pack(sysp, [N, M, L]) builds the grid with \(N \times M \times L\) points along each reciprocal-lattice direction.

For insulators, convergence with \(N\) is exponential — the error shrinks as \(e^{-\alpha N}\) because the occupied bands are well separated from the virtual ones. For metals, convergence degrades to a power law because the Fermi surface cuts bands; smearing or tetrahedron methods (not yet in vibe-qc) accelerate convergence there.

The Γ-only limit and molecular-limit validation

Setting the mesh to [1, 1, 1] samples only \(\mathbf{k} = \mathbf{0}\) (the Γ point), equivalent to using purely real Bloch sums. This is the cheapest periodic calculation and the right default for a large unit cell where a denser mesh doesn’t change anything.

The molecular-limit test — make the cell so big that periodic images don’t interact, then the per-cell energy must equal the isolated-molecule HF energy — is vibe-qc’s cleanest correctness check. Run it at every k-mesh size: the per-cell energy must be independent of the mesh (because in the molecular limit, no \(\mathbf{k}\)-dependent structure exists to resolve). vibe-qc’s regression suite asserts this at machine precision across dim ∈ {1, 2, 3}.

The Coulomb-sum problem in 3D

The last term of the Fock matrix — the Coulomb contribution from the infinite lattice of nuclei and electrons — is where periodic HF gets hard in 3D. The per-unit-cell sum

\[ \sum_{\mathbf{T}} \frac{1}{|\mathbf{r} - \mathbf{T}|} \]

is conditionally convergent: direct truncation gives shape- dependent answers (see tutorial 06 for the full story). 1D and 2D are absolutely convergent and DIRECT_TRUNCATED is fine for them. For 3D, set opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D; the dispatcher routes to vibe-qc’s full Saunders–Dovesi multipolar Ewald backend (validated to ω-invariance to 1e-8 Ha across H₂ / LiH / MgO / Ne in the bulk benchmark suite). The Ewald path currently requires orthorhombic cells; triclinic support is on the roadmap.

Resources

~150 MB peak RAM, ~2 s on one core (Apple M2 baseline) for the H₂-chain pob-TZVP / 6-k example. Cost scales linearly in the k-mesh size and as \(\mathcal{O}(N_\text{bf}^4)\) in the per-cell basis. EWALD_3D adds roughly 10–30% wall time for the FFT-based long-range piece on top of the direct-truncated baseline.

References

  • Bloch’s theorem. F. Bloch, “Über die Quantenmechanik der Elektronen in Kristallgittern,” Z. Phys. 52, 555 (1929).

  • CRYSTAL-style Gaussian crystalline orbitals. C. Pisani and R. Dovesi, “Exact-exchange Hartree-Fock calculations for periodic systems. I. Illustration of the method,” Int. J. Quantum Chem. 17, 501 (1980); C. Pisani, R. Dovesi, and C. Roetti, Hartree-Fock Ab Initio Treatment of Crystalline Systems, Lecture Notes in Chemistry 48, Springer (1988).

  • Monkhorst-Pack k-mesh. H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations,” Phys. Rev. B 13, 5188 (1976).

  • Textbook, crystalline electronic structure. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Cengage Learning reprint (1976), chapter 8 for Bloch’s theorem, chapter 10 for the tight-binding formulation that most closely matches Gaussian crystalline orbitals.

  • Saunders-Dovesi multipolar 3D Coulomb. V. R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco, and C. Roetti, “On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions,” Mol. Phys. 77, 629 (1992).