Periodic systems I. Density fitting

These notes are the preliminary calculations I ran to build familiarity with periodic systems, lightly cleaned up for publishing.

The four-index electron-repulsion integral is a main ingredient for HF, MP2, RPA, hybrid functionals, and most of quantum chemistry. It is written as

$$(pq|rs) = \int \int d\mathbf{r}\, d\mathbf{r}' \, \frac{\phi_p(\mathbf{r})\phi_q(\mathbf{r})\phi_r(\mathbf{r}')\phi_s(\mathbf{r}')}{r_{12}}$$

For periodic systems these integrals are computationally expensive, so we approximate them with density fitting. We can read the integral as the repulsion between two generalized electron densities $\rho_{pq} = \phi_p \phi_q$, such that

$$(pq|rs) = \int \int d\mathbf{r}\, d\mathbf{r}' \, \rho_{pq}(\mathbf{r}) \frac{1}{r_{12}} \rho_{rs}(\mathbf{r}')$$

The approximation is to expand that density in an auxiliary basis,

$$\tilde{\rho}_{pq}(\mathbf{r}) = \sum_P^{N_\text{fit}} d_P^{pq} \, \chi_P(\mathbf{r})$$

where the $d_P^{pq}$ are fitting coefficients, obtained by minimizing some functional of the fitting error (for the Coulomb metric, the error in the electric field). One standard form is

$$d_Q^{pq} = \sum_P (pq|P)\, [\mathbf{J}^{-1}]_{PQ}$$

with the three-center integral

$$(pq|P) = \int \int d\mathbf{r}\, d\mathbf{r}' \, \phi_p(\mathbf{r})\phi_q(\mathbf{r}) \frac{1}{r_{12}} \chi_P(\mathbf{r}')$$

and $\mathbf{J}$ the Coulomb repulsion matrix of the auxiliary functions $\chi$. The four-index integral then factorizes as

$$\begin{gathered} (pq|rs) = \sum_Q d_Q^{pq} (Q|rs) = \sum_{PQ} (pq|P)\,[\mathbf{J}^{-1}]_{PQ}\,(Q|rs) \end{gathered}$$

which resembles a resolution of identity in quantum mechanics, so density fitting (DF) is also called the resolution-of-identity (RI) approximation.

DF-MP2

The bottleneck of an MP2 calculation is the $\mathcal{O}(N^5)$ integral transformation from AO to MO basis, $(ia|jb) = \sum_{\mu\nu\rho\sigma} C_\mu^i C_\nu^a C_\rho^j C_\sigma^b (\mu\nu|\rho\sigma)$. In the DF scheme this becomes $\mathcal{O}(N_\text{occ}^2 N_\text{vir}^2 N_\text{aux})$, much less than the conventional $\mathcal{O}(N_\text{occ} N_\text{AO}^4)$ where $N_\text{AO}$ dominates.

DF in PySCF

There are four flavors of DF implemented in PySCF. Gaussian density fitting (GDF) uses Gaussian-type orbitals as the auxiliary functions. Fast-Fourier-transform density fitting (FFTDF) uses plane waves, and is also known as the Gaussian-and-plane-wave (GPW) approach. Mixed density fitting (MDF) uses both. Range-separated density fitting (RSDF) is an efficient implementation of GDF that splits the Coulomb metric into short-range and long-range parts with an error function. The short-range part decays fast with distance, so its lattice sum converges quickly when evaluated directly in real space; the long-range part is smooth, so its Fourier components die off quickly and the sum converges in few terms in reciprocal space. Splitting the metric lets each piece be evaluated where it is cheap, the same logic as Ewald summation.

By default PySCF uses FFTDF in a periodic calculation. This can be changed, for example to GDF, by altering the object as

mf = scf.KRHF(cell, kpts).density_fit()

Only GDF and MDF can be used for all-electron calculations (that is, without freezing the core), because FFTDF would need an impractically large plane-wave basis to describe the sharp core orbitals accurately.

Diamond

The first benchmark is an energy calculation on the diamond primitive cell. Use the primitive cell when possible, since the result is the same as the conventional cell while the conventional cell costs $(n^c_\text{atom}/n^p_\text{atom})^\kappa$ times more, where $\kappa$ is the scaling exponent of the method. The calculation is carried out with the following parameters.

1. GTH-HF-rev pseudopotential
2. GTH-cc-pVDZ basis (not in PySCF’s preloaded pbc basis sets; parse it with gto.mol.basis.parse)
3. k-point meshes of size $n^3$, with $n = 1\!-\!4$ for FFTDF, $1\!-\!5$ for GDF, and $1\!-\!6$ for RSDF

When running MP2 at large k-meshes, don’t store the $T_2$ amplitudes, which would eat memory at large meshes. This is achieved by

mmp = mp.KMP2(mf)
mmp.kernel(with_t2=False)

Energies

HF energy (Eh)	GDF	RSDF	FFTDF
$n=1$	−10.2128685	−10.2128685	−10.21286713
$n=2$	−10.95686485	−10.95686485	−10.95687283
$n=3$	−11.02220654	−11.02220654	−11.02221347
$n=4$	−11.02867244	−11.02867244	−11.02868016
$n=5$	−11.02830944	−11.02830944
$n=6$		−11.02734337

Taking FFTDF as numerically exact for this AO basis, the tables below list the errors $|E_\mathrm{GDF/RSDF} - E_\mathrm{FFTDF}|$.

HF energy error (Eh)	GDF	RSDF
$n=1$	0.0000013715	0.0000013726
$n=2$	0.0000079775	0.0000079782
$n=3$	0.0000069266	0.0000069244
$n=4$	0.0000077172	0.0000077172
MAE	0.0000059982	0.0000059981
RMSE	0.0000065775	0.0000065772

MP2 corr. energy error (Eh)	GDF	RSDF
$n=1$	0.0000270229	0.0000270241
$n=2$	0.0000070320	0.0000070301
$n=3$	0.0000043544	0.0000043547
$n=4$	0.0000040885	0.0000040886
MAE	0.0000106245	0.0000106244
RMSE	0.0000142773	0.0000142777

The errors are all on the order of $10^{-6}$ hartree, the same magnitude as the HF energy errors reported in the RSDF paper (Ye and Berkelbach, J. Chem. Phys. 154, 131104 (2021)). Even for MP2 correlation energies, GDF and RSDF stay on the same order. There is no significant accuracy difference between GDF and RSDF as a function of mesh size when compared against FFTDF.

CPU time

HF CPU time (s)	$N_k$	FFTDF	GDF	RSDF
$n=1$	1	10.3	2390.46	231.92
$n=2$	8	143.69	3255.22	339.53
$n=3$	27	1529.01	5150.75	699.31
$n=4$	64	8294.77	8953.04	1914.79
$n=5$	125		40943.18	29613.32
$n=6$	216			89105.05

While FFTDF is fast at small $N_k$, its cost catches up to GDF’s by $N_k = 64$, consistent with the known $\mathcal{O}(N_k^2)$ scaling of FFTDF. GDF and RSDF grow more slowly over this range, with RSDF keeping the absolute time almost 10-fold below GDF from $N_k = 1$ to $27$. (GDF and RSDF timings include a DF warm-up, since their integrals are not computed on the fly; that warm-up is also why these timings shouldn’t be compared too literally against the near-linear small-mesh scaling reported in the RSDF paper, which was measured on different systems and hardware.)

MP2 CPU time (s)	$N_k$	FFTDF	GDF	RSDF
$n=1$	1	0.49	0.13	0.23
$n=2$	8	261.59	2.23	2.45
$n=3$	27	9198.81	52.48	45.42
$n=4$	64	120043.76	496.46	525.25
$n=5$	125		15230.4	15535.94
$n=6$	216			66361.81

For MP2 the gap is dramatic. The bottleneck is the integral transformation, and GDF/RSDF work with three-center rather than four-center integrals, so their times sit far below FFTDF’s. The strong $N_k$ dependence itself comes from conservation of crystal momentum: every orbital index in $(i\mathbf{k}_i\, a\mathbf{k}_a | j\mathbf{k}_j\, b\mathbf{k}_b)$ carries a k-point label, and momentum conservation fixes only one of the four, leaving a triple sum over k-points, an $\mathcal{O}(N_k^3)$ factor on top of the per-cell orbital work. Effective exponents fitted over a small range of meshes come out steeper than the formal ones (my fit suggested roughly $N_k^5$-like growth for GDF/RSDF here) because the range is short and the warm-up and transformation costs haven’t separated cleanly.