Periodic HF on a 1D H chain¶

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

Run the baseline¶

Use a working engine installation and an empty run directory outside the checkout; see Running vibe-qc. Coordinates and lattice lengths below are in bohr. The calculation uses a three-dimensional periodic array with transverse vacuum, not a true one-dimensional Coulomb boundary condition. The chosen vacuum and mesh are starting settings to converge, not certified values for an isolated chain.

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

# Unit cell: H2 oriented along the chain, 4.0 bohr between unit cells.
unit_cell = [Atom(1, [0, 0, 0]), Atom(1, [1.4, 0, 0])]

# Chain axis on x. A wire that needs k sampling is run as a dim=3 cell
# with vacuum on y and z (15 bohr here); converge the energy with it.
lattice = np.diag([4.0, 15.0, 15.0])

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

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

result = run_periodic_job(
    sysp,
    basis,
    method="RHF",
    jk_method="gdf",
    kpoints=kpoints,
    conv_tol_energy=1e-10,
    output="h-chain-rhf-k4-vac15",
)
if not result.converged:
    raise RuntimeError("Periodic RHF did not converge; inspect the output and SCF trace")
print(f"E_RHF per cell = {result.energy:.10f}  ({result.n_iter} iters)")

run_periodic_job is the recommended user-facing dispatcher. The chain is modelled as a dim=3 cell with transverse vacuum because the native GDF route (Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644)) has no consistent Coulomb gauge for a true dim=1 cell: PeriodicSystem(dim=1, ...) with jk_method="gdf" raises NotImplementedError at Gamma and with a k mesh, and jk_method="auto" refuses a dim=1 system before any setup. (The Gamma-only form used to run, but its total energy depended on the lattice-sum cutoff.) The cost of the vacuum-box form grows with the box volume, so start with a modest vacuum and check how the energy per cell changes as it grows. The historical run_rhf_periodic_scf direct-truncated branch is retained only for vacuum-padded molecular-limit diagnostics, not quantitative 1D solids. For a Gamma-only calculation, pass KPoints.gamma(sysp).

Turn the baseline into a convergence study¶

Keep the atomic geometry and chain repeat length fixed while studying one numerical choice at a time. For each point, rebuild the PeriodicSystem, basis and k-point object together if the lattice changes, and use a distinct output stem. A k-point object belongs to the cell from which it was constructed.

Study	Change	Hold fixed
Chain sampling	mesh along x	transverse mesh components at 1, vacuum, basis and route
Transverse separation	y and z lattice lengths	chain geometry, x mesh and numerical controls
Lattice/fit support	the selected route’s documented cutoffs	cell, basis and k mesh

Record energy per simulation cell, increments, convergence and resource use. Changing vacuum or k sampling can change cost substantially. Do not interpret an SCF stopping criterion as convergence of the vacuum, basis or lattice sums. See periodic method routes and lattice-cutoff convergence.

Other ways to build a k-mesh, `vibeqc.KPoints`¶

The KPoints builder fronts every common k-point input mode in a single API. The following are alternative constructor examples, not a sequence of calculations to append to the chain run. Three-dimensional meshes and auto-density meshes can add unwanted transverse sampling to a vacuum cell; keep [N, 1, 1] for the chain study. Symmetry and band-path helpers also need their optional dependencies. See the k-point guide.

The five most useful for everyday work:

# doc-audit: skip - alternative constructors using sysp and imports from the baseline
# 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, n_kpts_per_atom=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 conventions behind these modes are the Monkhorst-Pack mesh (Monkhorst, Hendrik J. and Pack, James D. (1976). *Special points for Brillouin-zone integrations.* Physical Review B, **13**(12), 5188–5192. doi:[10.1103/PhysRevB.13.5188](https://doi.org/10.1103/PhysRevB.13.5188)), spglib symmetry detection (Togo, Atsushi, Shinohara, Kohei, and Tanaka, Isao (2024). *Spglib: a software library for crystal symmetry search.* Science and Technology of Advanced Materials: Methods, **4**(1), 2384822. doi:[10.1080/27660400.2024.2384822](https://doi.org/10.1080/27660400.2024.2384822)), the AFLOW KPPRA density (Curtarolo, Stefano, Setyawan, Wahyu, et al. (2012). *AFLOW: An automatic framework for high-throughput materials discovery.* Computational Materials Science, **58**, 218–226. doi:[10.1016/j.commatsci.2012.02.005](https://doi.org/10.1016/j.commatsci.2012.02.005)) and the HPKOT band path (Hinuma, Yoyo, Pizzi, Giovanni, et al. (2017). *Band structure diagram paths based on crystallography.* Computational Materials Science, **128**, 140–184. doi:[10.1016/j.commatsci.2016.10.015](https://doi.org/10.1016/j.commatsci.2016.10.015)).

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 that images don’t overlap. The native GDF route raises NotImplementedError for a dim=1 cell, at Gamma and with a k mesh, so the example above pads a dim=3 cell with vacuum instead. dim=2 treats cols 0 and 1 as in-plane lattice vectors; col 2 is synthesized bookkeeping, not a vacuum gap, and the slab energy is invariant to it (build these with vq.slab_2d). 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).
lattice_opts.gamma_only_0, off by default. Set it only when you mean the molecular limit: one unit cell, no periodic image coupling. Image cells are selected by |g| <= cutoff_bohr, so if your shortest lattice vector is longer than the cutoff, no image survives and the calculation is a free-boundary cluster wearing a periodic label: it converges, it returns the isolated molecule’s energy, and it is exactly k-independent, so refining the mesh only confirms a stable wrong answer. The driver therefore refuses that cell list unless you set this flag. If the refusal surprises you, raise cutoff_bohr or use a primitive cell with a matched k mesh; if you really did mean a molecule in a vacuum box (the big-box case below), set the flag. It declares intent and changes no sum.

Direct native Gamma operator¶

run_rhf_periodic_gamma_scf resolves physical support for its 3D native RHF calculation. A single Gamma quadrature point gives the same real-space density matrix at every lattice displacement. Both density indices therefore participate in the J/K image sum. Treating the inner density as zero away from the home cell while folding every outer image does not give the derivative of the corresponding Gamma energy.

The finite operator uses the following support, with the requested radii:

Overlap, kinetic and nuclear attraction retain AO pairs whose centers are within cutoff_bohr.
Nuclear charges and ECP projector centers lie within nuclear_cutoff_bohr of each AO-pair midpoint.
Each two-electron product has center separation at most cutoff_bohr; the two product midpoints are within eri_interaction_cutoff_bohr, or cutoff_bohr when that option is zero. The enumeration includes the translated counterparts of all retained terms and their density blocks.
Direct ionic repulsion uses the same physical nuclear radius.

These domains are closed under the integral permutations and common translations required for a Hermitian operator and reciprocal density response. The builder returns raw J/K sums; it does not repair them by averaging matrix triangles. The public build_jk_gamma_physical function uses this two-electron domain without changing the supplied options.

This is a formulation correction to the previous independent translation balls and home-cell-only inner density, so finite-cutoff energies can change. It is still a finite direct approximation: cutoff and k-point convergence, and the choice of electrostatic treatment, require separate validation. The explicitly named build_jk_gamma_molecular_limit helper and lower-dimensional native Gamma driver retain their existing formulation. GDF, BIPOLE and the native multi-k driver are separate routes.

The density convention follows Pisani and Dovesi (1980), Equations (1), (5)-(8), and Causa et al. (1988), Equations (1)-(2). The latter also explains why finite Gamma sampling does not establish the correct asymptotic exchange density of a bulk crystal.

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 with lattice_opts.gamma_only_0 = True produces the same energy as molecular RHF, machine precision across dim ∈ {1, 2, 3} and multiple k-mesh sizes. Without the declaration the driver refuses the box, by design; see the option above.
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-metric Ewald backend. The supported route uses exact erfc-screened real-space ERIs plus reciprocal long-range Coulomb; it does not activate the fail-closed quartet multipole prototype. See the Ewald user guide for details on the full-metric FFT long-range solver and the FFTDF/GDF parity track.

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). Periodic SCF convergence: damping, DIIS, and level shifts 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¶

The sections below build up the periodic-HF machinery the code runs: Bloch’s theorem and the Bloch-summed AO basis, the k-diagonal generalised eigenvalue problem it produces, Monkhorst-Pack sampling of the Brillouin zone, the Γ-only and molecular-limit checks, and the conditionally convergent 3D Coulomb sum that motivates the Ewald backend.

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; Pisani, Cesare and Dovesi, Roberto (1980). *Exact-exchange Hartree-Fock calculations for periodic systems. I. Illustration of the method.* International Journal of Quantum Chemistry, **17**(3), 501–516. doi:[10.1002/qua.560170311](https://doi.org/10.1002/qua.560170311), Pisani, Cesare, Dovesi, Roberto, and Roetti, Carla (1988). *Hartree-Fock Ab Initio Treatment of Crystalline Systems.* Lecture Notes in Chemistry, Vol. 48, Springer, **48**. doi:[10.1007/978-3-642-93385-1](https://doi.org/10.1007/978-3-642-93385-1)).

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 (Pisani, Cesare, Dovesi, Roberto, and Roetti, Carla (1988). *Hartree-Fock Ab Initio Treatment of Crystalline Systems.* Lecture Notes in Chemistry, Vol. 48, Springer, **48**. doi:[10.1007/978-3-642-93385-1](https://doi.org/10.1007/978-3-642-93385-1), Dovesi, R., Orlando, R., et al. (2014). *CRYSTAL14: A program for the ab initio investigation of crystalline solids.* International Journal of Quantum Chemistry, **114**, 1287–1317. doi:[10.1002/qua.24658](https://doi.org/10.1002/qua.24658)).

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 (Monkhorst, Hendrik J. and Pack, James D. (1976). *Special points for Brillouin-zone integrations.* Physical Review B, **13**(12), 5188–5192. doi:[10.1103/PhysRevB.13.5188](https://doi.org/10.1103/PhysRevB.13.5188)) constructs a uniform grid of k-points offset from the BZ origin, optionally reduced to the irreducible wedge by the space-group operations (Dovesi, Roberto (1986). *On the role of symmetry in the ab initio Hartree-Fock linear-combination-of-atomic-orbitals treatment of periodic systems.* International Journal of Quantum Chemistry, **29**, 1755–1774. doi:[10.1002/qua.560290608](https://doi.org/10.1002/qua.560290608)). 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 the parameter-free Gilat-Raubenheimer net, a tetrahedron-family integrator (metallic BZ integration), 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 (Ewald, Paul Peter (1921). *Die Berechnung optischer und elektrostatischer Gitterpotentiale.* Annalen der Physik, **369**(3), 253-287. doi:[10.1002/andp.19213690304](https://doi.org/10.1002/andp.19213690304), Dovesi, R., Pisani, C., et al. (1983). *Treatment of Coulomb interactions in Hartree-Fock calculations of periodic systems.* Physical Review B, **28**, 5781–5792. doi:[10.1103/PhysRevB.28.5781](https://doi.org/10.1103/PhysRevB.28.5781); see Madelung constants via Ewald summation for the full story). Mixed-boundary 1D and 2D systems also require dimension-specific Coulomb kernels. No public route has one for 1D yet (pad a dim=3 cell with vacuum); 2D uses SLAB_EWALD_2D (Parry, D. E. (1975). *The electrostatic potential in the surface region of an ionic crystal.* Surface Science, **49**, 433–440. doi:[10.1016/0039-6028(75)90362-3](https://doi.org/10.1016/0039-6028(75)90362-3), de Leeuw, S. W. and Perram, J. W. (1979). *Electrostatic lattice sums for semi-infinite lattices.* Molecular Physics, **37**(4), 1313–1322. doi:[10.1080/00268977900100951](https://doi.org/10.1080/00268977900100951)) or explicit slab GDF. DIRECT_TRUNCATED is only a vacuum-padded molecular-limit diagnostic. For 3D, set opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D; the dispatcher routes to vibe-qc’s full-metric Ewald backend, with exact real-space ERIs rather than the fail-closed quartet multipole prototype (validated to ω-invariance to 1e-8 Ha across H₂ / LiH / MgO / Ne in the bulk benchmark suite). The native FFT long-range solver supports skew-cell metrics; FFTDF/GDF parity is the active production path for tight 3D ionic crystals.

For the full menu of periodic Coulomb/exchange routes that tame this sum (BIPOLE: Saunders, V. R., Freyria-Fava, C., et al. (1992). *On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions.* Molecular Physics, **77**, 629–665. doi:[10.1080/00268979200102671](https://doi.org/10.1080/00268979200102671); GDF: Sun, Qiming, Berkelbach, Timothy C., et al. (2017). *Gaussian and plane-wave mixed density fitting for periodic systems.* Journal of Chemical Physics, **147**(16), 164119. doi:[10.1063/1.4998644](https://doi.org/10.1063/1.4998644); GPW / GAPW: Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1997). *A hybrid Gaussian and plane wave density functional scheme.* Molecular Physics, **92**(3), 477-487. doi:[10.1080/00268979709482119](https://doi.org/10.1080/00268979709482119), Lippert, Gerald, Hutter, Jürg, and Parrinello, Michele (1999). *The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations.* Theoretical Chemistry Accounts, **103**(2), 124-140. doi:[10.1007/s002140050523](https://doi.org/10.1007/s002140050523)) and which external program validates each, see periodic-SCF methods and the periodic JK routes & parity policy.

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.
Gaussian-density 3D electrostatics. 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). This supplies electrostatic-potential, penetration, and spheropole context, not the implemented quartet dispatcher.
Periodic bipolar quartet expansion. C. Pisani, R. Dovesi, and C. Roetti, Hartree-Fock Ab Initio Treatment of Crystalline Systems, Lecture Notes in Chemistry 48, Springer (1988), Chapter II.4c. DOI 10.1007/978-3-642-93385-1.