Periodic systems II. The thermodynamic limit

Simulating an extended system means approximating a system of infinite size, also known as the thermodynamic limit (TDL).

Supercell vs k-point sampling

In practice we cannot compute an infinite system, so we approach the limit by extrapolating finite systems of increasing size. There are two ways to grow the system, by either

1. increase from a unit cell to larger and larger supercells, or
2. sample the first Brillouin zone with an increasingly fine k-point mesh.

These two are mathematically equivalent as long as

$$N_k \times n_\text{atm} = N_\text{atmSc}$$

where $n_\text{atm}$ is the number of atoms per unit cell, $N_k$ the total number of k-points in a k-sampled calculation, and $N_\text{atmSc}$ the total atoms in the supercell. They are equivalent because enlarging the supercell shrinks its Brillouin zone by the same factor, so a finer mesh on the unit cell’s Brillouin zone ends up sampling the same set of points.

As a check, a $(2\times 2\times 2)$ diamond supercell gives −87.6549188020541 hartree; dividing by the supercell size reproduces the $(2\times 2\times 2)$ k-point sampled value of −10.95686485 (the k-sampled energy being reported per cell, $E = \frac{1}{N_k}\sum_{k} E(k)$).

Computational cost

The cost of a periodic calculation scales in general as

$$\text{Cost} \sim N_k^{\,n}\, n_\text{atm}^{\,m}$$

where the exponents depend on the method. HF has $n=2$, $m=3$, MP2 has $n=3$, $m=5$, and CCSD has $n=3$, $m=6$. Comparing the two routes to the same system size (doubling each cell dimension, a supercell at the $\Gamma$ point with $N_k=1$ versus a $2^3$ mesh on the unit cell), the supercell approach costs a factor of $\mathcal{N}^2$ more for HF, where $\mathcal{N}$ is the supercell atom count, and $\mathcal{N}^2$ and $\mathcal{N}^3$ more for MP2 and CCSD. Translational symmetry is doing real work; k-point sampling wins, and the rest of these notes uses it.

Extrapolating in the mesh size

To reach the TDL we extrapolate to $N_k = \infty$, assuming the energies are converging. The HF and MP2 correlation energies of diamond, plotted against the reciprocal mesh size below.

The extrapolation variable follows from the decay of the leading finite-size error. With the Madelung-constant correction applied (next section), the remaining HF and correlation errors fall off like the inverse volume of the sampled cell, $N_k^{-1}$. Without it, the exchange error decays like an inverse length, $N_k^{-1/3}$, and that slower decay is the one that matters at large $N_k$. Since the correction is on here, I extrapolate linearly in $N_k^{-1}$ using the two largest meshes ($N_k = 5^3$ and $6^3$), executed as

inf_scf = interpolate.interp1d((Nk**(-1))[:2],
    scf_e[-2:], fill_value='extrapolate')

with interpolate from scipy. The extrapolated energies are shown below.

These TDL values get used later for the cohesive energy of diamond.

Ewald correction

Because we approximate an infinite system, there are finite-size errors from the kinetic, Coulomb, exchange, and correlation terms, each decaying at its own rate with $N_k$. The dominant one at the HF level comes from an integrable divergence at $\mathbf{G} = \mathbf{0}$ in the nonlocal exact-exchange energy. The standard fix is the “Ewald” scheme (Broqvist, Alkauskas, and Pasquarello, Phys. Rev. B 80, 085114 (2009)), which shifts the occupied-orbital energies by the Madelung constant of the system; PySCF applies it by default.

With the correction turned off, both HF and MP2 converge visibly more slowly toward the TDL values obtained above. Treating finite-size errors explicitly matters for both accuracy and efficiency, since a corrected calculation at a modest mesh is worth an uncorrected one at a much larger mesh.