Tutorial 30: H-chain BvK equivalence — vibe-qc’s CCM ancestor¶
You’ll learn: how to set up the model system from the project author’s 2014 CCM-HF paper (Peintinger & Bredow, J. Comput. Chem. 35, 839 — free access, cover article) in vibe-qc, run periodic SCF on it, and test the Born-von-Kármán equivalence between supercell and k-mesh representations across a sweep of cell sizes — turning up a clean two-regime characterisation of where v0.7’s gauge fix is sound vs where the open multi-cell density bug still bites.
Why: the H-chain is the proof-of-concept system the 2014 paper used to demonstrate that the cyclic cluster model (finite cluster + periodic boundary conditions imposed via image-cell folding) reproduces the conventional supercell treatment to FFT precision. That paper is the methodological ancestor of the v2.x roadmap: when v2.0 lands the CCM at HF level inside vibe-qc, the paper’s headline numerical claim becomes the regression-test target. This tutorial sets up the conventional supercell side of that comparison — and the results below double as the regression-test fixture for the v0.7.x maintenance close.
Where we stand on this in v0.7.x
Conventional supercell periodic SCF in the molecular-isolated limit (large per-pair spacing) — ✓ works to ~0.2 µHa/pair in v0.7.3, exactly as Born-von-Kármán theory says it should. The v0.7 gauge fix delivers FFT-grid precision in this regime.
Conventional supercell periodic SCF with real inter-cell density overlap (small / chemically-realistic per-pair spacings) — ✗ still off, by mHa to Ha per pair, scaling roughly as 1/L. This is the tight-ionic multi-cell density bug — the v0.7.x maintenance window will close it. Until then,
numbers in this regime are not yet meaningful in absolute terms.Cyclic cluster model (the paper’s actual method) — ✗ ships with v2.0. The image-folding lattice machinery + Wigner-Seitz boundary treatment for the multi-centre integrals is the work the v2.x roadmap covers.
The system¶
A 1D linear chain of hydrogen atoms arranged as periodic H₂ pairs:
H–H··H–H··H–H··H–H ... ← H₂ pairs in a 1D chain
↑ ↑
R_S R_L R_S = short bond, R_L = long bond
Per-pair spacing is R_pair = R_S + R_L. We hold R_S = 1.4 bohr
(H₂ equilibrium) and sweep R_L from 2.5 → 18.0 bohr to walk the
chain from “tightly coupled” (interesting chemistry, gapped Peierls
ground state) to “essentially isolated H₂ molecules in vacuum”
(the regime where BvK should be perfectly trivial).
What Born-von-Kármán equivalence says¶
Periodic boundary conditions plus translation invariance imply: the same physical system gives the same per-cell energy regardless of which unit-cell representation you choose. For our H-chain:
H₂ unit cell sampled with N×1×1 k-mesh
H₈ supercell (4 H₂ pairs) sampled at Γ
H₁₆ supercell sampled at Γ
H₃₂ supercell sampled at Γ
— all describe the same infinite chain. Per-pair energies must match to numerical precision. This is the principle the cyclic cluster model formalises: fold the periodic Bloch sum onto a finite cluster, and a Γ-point cluster calculation reproduces multi-k Bloch.
The sweep — five spacings × five representations¶
A single Python file sweeps both axes and prints the table. Source:
examples/periodic/input-h-chain-bvk-sweep.py.
Run via the
vq
queue tool on planetx (the laptop OOMs on the larger cells); 11
seconds wall-clock total on a 32-CPU box.
from itertools import product
import numpy as np
import vibeqc as vq
R_S = 1.4
R_LS = [2.5, 4.0, 6.0, 10.0, 18.0] # → R_pair = 3.9, 5.4, 7.4, 11.4, 19.4
PAD = 30.0
REPRESENTATIONS = [
(1, [1, 1, 1]), # H2 cell, Γ
(1, [8, 1, 1]), # H2 cell, 8x1x1 k-mesh
(4, [1, 1, 1]), # H8 supercell, Γ
(8, [1, 1, 1]), # H16 supercell, Γ
(16, [1, 1, 1]), # H32 supercell, Γ
]
opts = vq.PeriodicSCFOptions()
opts.lattice_opts.cutoff_bohr = 15.0
opts.lattice_opts.nuclear_cutoff_bohr = 15.0
opts.conv_tol_energy = 1e-9
opts.max_iter = 80
for R_L, (n_pairs, kmesh) in product(R_LS, REPRESENTATIONS):
R_pair = R_S + R_L
cell = []
for i in range(n_pairs):
cell.append(vq.Atom(1, [i*R_pair + 0.0, 0.0, 0.0]))
cell.append(vq.Atom(1, [i*R_pair + R_S, 0.0, 0.0]))
sysp = vq.PeriodicSystem(dim=1, lattice=np.diag([n_pairs*R_pair, PAD, PAD]),
unit_cell=cell)
basis = vq.BasisSet(sysp.unit_cell_molecule(), "sto-3g")
kp = (vq.KPoints.gamma(sysp) if kmesh == [1, 1, 1]
else vq.KPoints.monkhorst_pack(sysp, kmesh))
r = vq.run_rhf_periodic_scf(sysp, basis, kp, opts)
print(f"R_pair={R_pair:5.2f} N{n_pairs:2d} k={tuple(kmesh)} "
f"E/pair={r.energy/n_pairs:14.8f} conv={r.converged}")
Results (vibe-qc 0.7.3, planetx, 2026-05-09)¶
For each per-pair spacing, the per-pair energy across the five representations:
R_pair | representation | E/pair (Ha) | conv | iter
----------|-----------------|----------------|-------|------
3.900 | H2 Γ | -2.23604272 | ✓ | 2
3.900 | H2 8x1x1 | -1.09376293 | ✓ | 15
3.900 | H8 Γ | -1.10308388 | ✓ | 7
3.900 | H16 Γ | -1.10080687 | ✓ | 8
3.900 | H32 Γ | -1.09966810 | ✓ | 9
→ max drift = 1.14 Ha/pair ✗ (1,142,280 µHa)
----------|-----------------|----------------|-------|------
5.400 | H2 Γ | -1.73083481 | ✓ | 2
5.400 | H2 8x1x1 | -1.11519942 | ✓ | 9
5.400 | H8 Γ | -1.11581740 | ✓ | 6
5.400 | H16 Γ | -1.11566741 | ✓ | 6
5.400 | H32 Γ | -1.11559240 | ✓ | 9
→ max drift = 0.62 Ha/pair ✗ (615,635 µHa)
----------|-----------------|----------------|-------|------
7.400 | H2 Γ | -1.53331149 | ✓ | 2
7.400 | H2 8x1x1 | -1.11667402 | ✓ | 6
7.400 | H8 Γ | -1.11668439 | ✓ | 5
7.400 | H16 Γ | -1.11667915 | ✓ | 5
7.400 | H32 Γ | -1.11667653 | ✓ | 5
→ max drift = 0.42 Ha/pair ✗ (416,637 µHa)
----------|-----------------|----------------|-------|------
11.400 | H2 Γ | -1.29311649 | ✓ | 2
11.400 | H2 8x1x1 | -1.11671121 | ✓ | 2
11.400 | H8 Γ | -1.11671194 | ✓ | 4
11.400 | H16 Γ | -1.11671152 | ✓ | 4
11.400 | H32 Γ | -1.11671131 | ✓ | 4
→ max drift = 0.18 Ha/pair ✗ (176,405 µHa)
----------|-----------------|----------------|-------|------
19.400 | H2 Γ | -1.11671433 | ✓ | 2
19.400 | H2 8x1x1 | -1.11671433 | ✓ | 2
19.400 | H8 Γ | -1.11671416 | ✓ | 4
19.400 | H16 Γ | -1.11671413 | ✓ | 4
19.400 | H32 Γ | -1.11671412 | ✓ | 4
→ max drift = 2.1e-7 Ha/pair ✓ µHa precision (0.21 µHa)
What just happened¶
Look at the right-hand column of the bottom block (R_pair = 19.4
bohr, R_L = 18 bohr). All five representations agree to the
eighth decimal place. Born-von-Kármán equivalence holds to
~0.2 µHa per pair — that’s the FFT-grid noise floor, exactly
what theory predicts. The v0.7 periodic-SCF gauge fix
(Löwdin's Compass) delivered.
But look at any of the other blocks. At R_pair = 11.4 bohr the drift across representations is 176 mHa/pair. At R_pair = 3.9 bohr (closest to a real chemically-bonded chain) it’s 1.14 Ha/pair — the H₂-cell-Γ energy of -2.236 Ha vs the H₃₂-cell-Γ energy of -1.100 Ha. These are not small differences. They are the fingerprint of the open multi-cell density bug that the homepage admonition flags and that the v0.7.x maintenance window will close.
The drift scales monotonically with the inverse cell length. The H₂ cell at Γ (the smallest cell) is consistently the worst-affected representation; the H₃₂ supercell (the largest single-cell representation) is consistently closer to the BvK-correct value that the multi-k H₂ representation also approximates. Larger cell → less of the bug. This is consistent with a 1/L electrostatic-self-image character.
For the H-chain physics: the per-pair energy in the chemically-bonded regime (R_pair ≤ 5–6 bohr) is what matters for modelling a real 1D crystal. Until the maintenance fix lands, the absolute number is not trustworthy, and any cohesive-energy or band-energy difference you compute is corrupted by an unknown amount of the L-dependent shift. The relative shape of the dimerisation curve (tutorial 17) is approximately correct because it’s evaluated at fixed unit-cell length, but the per-cell energies that go into a structure prediction are biased.
Regression-test fixture¶
The 19.4-bohr row is a clean regression-test target as is: every representation should give -1.116714 Ha to ~µHa precision. The remaining four rows are the regression-test fixture for the v0.7.x maintenance close: when the multi-cell-density bug fix lands, the per-pair drift at R_pair = 3.9, 5.4, 7.4, 11.4 must drop from the current 176 mHa–1.14 Ha levels to at most a few µHa. The numbers in this tutorial table are pinned to v0.7.3 — the next major release will measurably improve them, and the diff will be visible in this table when the tutorial gets re-run.
Why the paper bothered¶
The conventional supercell + dense k-mesh approach (the rows above) is mathematically equivalent to the cyclic cluster model in the limit, but the bookkeeping is fundamentally different:
Approach |
What’s evaluated |
What’s tricky |
|---|---|---|
Conventional supercell + Bloch sum |
\(\hat{F}(\mathbf{k}) = \sum_\mathbf{g} \hat{F}^{(\mathbf{g})} e^{i\mathbf{k}\cdot\mathbf{g}}\) — Fock matrix Fourier-transformed over lattice vectors |
Long-range Coulomb / exchange tails; multipole convergence; Ewald-style screening for 3D — exactly where the open v0.7.x bug lives |
Cyclic cluster model (Γ-point) |
Single Fock matrix on a finite cluster with PBC enforced via image-cell folding |
Multi-centre integrals at the Wigner-Seitz boundary — the contributing image cells must be picked carefully so each pair (μν|λσ) is counted exactly once |
The 2014 paper showed the second route can be made to give identical numerical results to the first, with cleaner extension to post-HF correlation methods. That last point is the v2.x payoff: once the CCM machinery is in place, MP2 → CCSD → CCSD(T) on solids becomes tractable because the correlation problem reduces to a molecular post-HF problem on a finite cluster.
Wider context¶
The author background and the longer arc from the 2014 paper to vibe-qc’s v2.x CCM is on the support page.
The cyclic-cluster-model roadmap entry is at
v2.0.0— HF-CCM.The Peierls-dimerisation companion (tutorial 17) uses the same H-chain framework to demonstrate the metal-insulator transition.
If you sponsor vibe-qc via GitHub Sponsors or via Ko-fi, the v2.x CCM work + the v0.7.x maintenance close that makes the small-cell rows of the table above honest are what you’re funding most directly. See the support page for the full pitch.