Accelerating seminumerical Fock-exchange calculations using mixed single- and double-precision arithmethic

The evaluation of the exact (Fock)-exchange matrix usually represents the major computational bottleneck (typically >80% of the computation time) within Hartree–Fock or hybrid-density functional theory (DFT) calculations. Traditionally, this requires the computation of the four-center-two-electron (4c2e) integral-tensor leading to a formal $O(N4)$ scaling, which is particularly problematic when combined with larger atomic-orbital (AO) basis sets, even if efficient screening techniques^1–10 are employed.

Therefore, seminumerical integration techniques, e.g., the pseudospectral method of Friesner et al.,^11–16 the chain-of-spheres-exchange (COS-X) method of Neese et al.,¹⁷ the seminumerical methods of Plessow and Weigend,¹⁸ Liu et al.,^19–21 and Kaupp and co-workers (the latter focusing more on local-hybrid functionals),^22–26 the semi-JK algorithm of Holzer,²⁷ and our sn-LinK method^28,29 are more efficient for large basis sets due to their superior $O(M3)$ formal scaling and their $O(Nbas2)$ scaling with the basis set size. Moreover, as shown in Refs. 27 and 29, seminumerical integration is perfectly suited for modern, highly parallel hardware, where performance relies heavily on the use of single-instruction-multiple-data (SIMD) instructions, which applies to both modern central processing units (CPUs) and even to a greater extent to graphic processing units (GPUs).

By default, most quantum chemistry programs execute the necessary floating point operations with double numeric precision (fp64) due to its reliable accuracy (about 10⁻¹⁶ relative precision).³⁰ However, since most “gaming grade” GPUs typically provide significantly less computational performance for fp64 operations compared to single-precision floating point (fp32) operations (e.g., the fp32:fp64 ratio for the GTX 1080Ti is 32:1), the possibility of executing as many of the necessary computation using fp32 operations needs to be explored, despite its lower numerical precision of about 10⁻⁷. In addition, the possibility of twofold speedups on CPUs also justifies the need for such a study even for pure CPU code. Investigations of pure fp32 or mixed fp32/fp64 execution have indeed already been carried out for the traditional 4c2e integral based Fock-exchange evaluation^31–36 and post-HF correlation methods,^32,37–39 but, to our knowledge, the applicability of reduced numerical precision within seminumerical integration has not yet been studied.

Therefore, we provide such a study in this work, which is organized as follows: First, we briefly review our seminumerical integration method sn-LinK from Ref. 29 in Sec. II. After reporting on the computational details in Sec. III, we explore the applicability of fp32 operations regarding performance and numerical stability in Sec. IV. This exploration is partitioned into three parts: In part 1 (Sec. IV A), we explore mixed fp32/fp64 evaluation of the three-center-one-electron (3c1e) integral tensor. In part 2 (Sec. IV B), we then explore the possibility of pure fp32 evaluation in all steps instead of only the 3c1e integral evaluation. Finally, in part 3 (Sec. IV C), we propose a specific self-consistent-field (SCF) method, denoted as i-sn-LinK, that employs incremental exchange-builds to reduce the numerical error from pure fp32 execution. Afterward, we illustrate the performance and accuracy of the so developed mixed-precision methods for a wide variety of molecules and basis sets in Sec. V and finally summarize our results in Sec. VI.

Seminumerical integration decomposes the 4c2e integral tensor

symmetrically as

employing numeric integration grids with grid points r_g and corresponding weights w_g.

Inserting this decomposition into the atomic orbital (AO) representation of the exact-exchange matrix leads to

which is evaluated in three consecutive steps:

The so-obtained exchange-matrix is finally symmetrized to account for the transpose in Eq. (3).

If only the exchange energy

is of interest (e.g., to compute the final energy after converging the SCF), the evaluation of Eq. (7) can be avoided and E_X can instead be obtained as

Equations (5) and (7) are evaluated with dense matrix–matrix multiplications employing batch-local matrices of asymptotically constant size (see Sec. 2.4 of Ref. 29 for details) utilizing highly optimized BLAS-3 libraries.

In contrast, evaluation of Eq. (6) requires the computation of the 3c-1e integral tensor

which usually represents the most expensive step. Effective integral screening techniques are therefore essential for an efficient implementation.

In order to assess the significance of a 3c1e-integral A_νλg, we consider its contribution to the total exchange energy

and to the final exchange matrix

An integral is then labeled as significant if either of the contributions is larger than a given threshold, i.e.,

For optimal performance, this selection is performed for whole batches of grid points at once, employing distance-independent (and therefore grid-point independent) integral bounds for A_νλg,¹⁰ as detailed in Secs. 2.2 and 3 of Ref. 29.

Unless stated otherwise, all calculations are performed with our FermiONs++ quantum chemistry program,^7,8 employing the “Karlsruhe” (“def2-”) basis sets⁴⁰ (the prefix “def2-” is omitted for simplicity) and the “gm4” multi-grid⁴¹ (∼5000 points per atom within the SCF and ∼15 000 points per atom for the final energy calculation). For yttrium, the respective effective core potential (ECP)⁴⁰ was employed to replace core electrons. For the evaluation of the Coulomb-interaction, the resolution-of-the-identity (RI-J) approximation^42–44 in combination with the RI-J-optimized basis set of Ref. 45 was used throughout this work. These settings are chosen to represent typical applications.

Since the computation of the exact-exchange contributions is usually the most expensive step within hybrid-DFT calculations, the discussion for Hartree–Fock calculations within this work is equally meaningful for hybrid-DFT calculations. Indeed, since only a fraction of exact Fock-exchange is employed in hybrid-DFT methods, the numerical errors from single-precision execution are proportionally lower in this case. We provide analogous results to Tables IX and XI using the PBE0 hybrid functional⁴⁶ instead of the Hartree–Fock method within the supplementary material.

For optimal performance, the CPU code was compiled with the Intel C++ compiler (ICPC) version 19.1.0⁴⁷ with all compiler-optimization enabled (“-Ofast,” “-march = native”), which is necessary to fully utilize SIMD-vector-instructions within the 3c1e integral kernels. The GPU code was compiled with NVCC-10.1 (CUDA-10.1),⁴⁸ also employing all possible compiler-optimizations (“-O3” and “-use fast math”). To provide sufficient parallel workload for each device, grid-batches of 512 grid points on CPUs, 20 480 points on the NVIDIA GV100 GPU, and 10 240 points on the 1080Ti GPU are employed.

SCF convergence is always measured by the root mean square of the DIIS error matrices (SPF–FPS).^49,50 Unless stated otherwise, all timings are given for one full exchange-matrix build averaged over all but the very first SCF cycles, since the very first Fock-build is considerably faster due to the sparsity of the superposition-of-atomic-densities (SAD) guess density matrix. For rigorous comparisons, errors in the converged energy (denoted ΔE), the root mean square deviation (RMSD) of the converged density matrix (denoted as ΔP), and the root mean square deviation of the converged nuclear forces (denoted as ΔForces) are always referenced to pure fp64-execution with all other settings being identical.

The test geometries^51,52 in this work are chosen to represent typical applications. In particular, A–T DNA-fragments [(DNA)_x] exemplify biochemistry applications, spherical water balls [(H₂O)₆₈], explicit solvent environments, LiF cutouts [(LiF)₃₆], materials science applications, and Y₂C₅H₂₀N₂₀O₁₃ inorganic complex chemistry. All geometries employed in this work are provided in the supplementary material.

In this section, the possibility of employing single-precision algebra is systematically studied, focusing on the speedups and accuracy. The total execution times corresponding to the here presented speedups are given in Sec. 2 of the supplementary material.

The evaluation and contraction of the 3c1e-integral tensor [Eq. (6)] usually represents the most time-consuming step within sn-LinK, particularly on CPUs and for smaller systems. Therefore, we first investigate the applicability of fp32 operations for this step.

We can assign each 3c1e-integral (or batch of integrals) an upper bound to its contribution to the final exchange-energy [Eq. (11)] and exchange-potential [Eq. (12)]. These upper bounds not only enable the linear-scaling evaluation of the 3c1e integrals via computing only the significant subset of the tensor but also enable a finer-grained partitioning of the 3c1e tensor into three instead of only two categories:

• Category one contains the most significant 3c1e integrals (⁠$ενλgE/K≥ϑE/Kfp64$⁠). These are computed with fp64 operations.

• Category two contains all 3c1e integrals that are too significant to be completely neglected, but not so significant that they require fp64 evaluation (⁠$ϑE/Kfp64>ενλgE/K≥ϑE/Kfp32$⁠). These are computed and accumulated with fp32 operations.

• Category three contains all the insignificant integrals (⁠$ενλgE/K<ϑE/Kfp32$⁠). These are not computed at all, even with pure fp64-execution.

Compared to the original sn-LinK from Ref. 29, we just split the set of significant integrals into two subsets (category one and two), resulting in only a slight adjustment of the pseudocode for the 3c1e-integral evaluation, which is given in Fig. 1.

The complete discarding of all category three integrals (line 3 of Fig. 1) is just standard integral screening, so we can employ the same “parent” thresholds as in the original sn-LinK, i.e., $ϑKfp32=ϑK=10−8$⁠, $ϑEfp32=ϑE=10−11$ within the SCF, and $ϑEfp32=ϑE=10−12$ in the final energy calculation. The screening-errors from these thresholds are usually well below 1 µE_h. This integral screening is not altered within this work, and instead, the speedups presented here are solely caused by a more efficient handling of the category two integrals by using a lower precision arithmetic.

We decided to choose our second set of thresholds, which distinguish between fp32- and fp64-execution (l.5 and 8 of Fig. 1) $ϑE/Kfp64$ dynamically, based on the original thresholds ϑ_E/K. Indeed, since the relative numerical precision of fp32 numbers is about 10⁻⁷, one would naturally choose $ϑE/Kfp64=107ϑE/K$⁠. In Table I, we test a variety of threshold-multipliers and find that we can employ even tighter fp64-thresholds, since still only a very small fraction of 3c1e-integrals need to be computed with fp64 operations. Therefore, we decided to use a fp64 threshold multiplier of 10⁵ for all further calculations, providing virtually perfect accuracy regarding both the converged energy and the converged density matrix, while over 97% of all integrals are computed using fp32 instructions. As expected, this results in nearly 2× speedup for the evaluation of the 3c1e integrals.

We further evaluate the robustness and possible speedups of mixed fp32/fp64 3c1e-integral evaluation for different molecular systems and basis sets in Table II. Mixed precision evaluation provides essentially error-free results for both the converged energy and the density matrix and converges within exactly the same amount of SCF cycles. That is, mixed-precision evaluation of the 3c1e integrals is essentially error-free while providing close to 2× speedups (1.73–1.97×) for the 3c1e integral evaluation, in line with the expected 2× performance gains for CPUs.

As presented in Sec. 1 of the supplementary material, this speedup is only achieved if vector-instructions (AVX2 in this work) are utilized. This matches our expectations, since one 256-bit vector instruction processes four fp64 but eight fp32 values, whereas one scalar instruction always processes 1 number, regardless of the precision.

However, the 2×-speedup for the 3c1e integral evaluation does not perfectly translate to the full K-build, since the other two steps [Eqs. (5) and (7)] have not been accelerated. As depicted in Table III, this effect is even more pronounced with GPU evaluation because the 3c1e integral evaluation contributes less to the total runtime on GPUs (compare also with the discussion in Sec. 5.2.1 of Ref. 29). For this reason, only minor speedups from employing the mixed precision evaluation of the 3c1e integrals are obtained for the GV100 GPU (0.96–1.79×) and the speedups on the GTX 1080Ti (1.6× to 2.5×) are far away from the theoretical value of 32×.

To summarize part 1, the 3c1e integral tensor is partitioned so that only the most significant contributions are evaluated with fp64 operations. The introduced errors are negligible while providing nearly 2× speedups in the integral-evaluation part. However, when employing GPUs with low fp64-performance, the BLAS-3 steps [Eqs. (5) and (7)] also need to be performed with fp32 operations to achieve optimal performance.

In contrast to the integral evaluation, computing the BLAS-3 steps [Eqs. (5) and (7)] with mixed precision instruction is not straightforward because they are best evaluated with a single call of a highly optimized dense linear algebra routine for optimal performance. Although algorithms for mixed precision linear-algebra have been proposed in the literature,⁵³ we prefer to, if possible, avoid them for the sake of simplicity.

Therefore, the possibility of pure fp32 evaluation, i.e., evaluation of Eqs. (5)–(7) with only fp32 operations and accumulation of the batch-local fp32-K-matrices into a global fp64-K-matrix, needs to be explored. Note that in the final energy calculation, Eqs. (5) and (6) are evaluated with single-precision, but the accumulation in Eq. (9) is performed with double-precision.⁵⁴ 

As depicted in Table IV, pure fp32 evaluation results in up to 7.4× speedups compared to pure fp64 execution. This speedup is still significantly less than the theoretically possible 32× speedups because the code is not purely limited by the floating point throughput alone and the demand for memory bandwidth and local storage (cache) is only reduced by a factor of two.

More important, however, is the fact that pure fp32 execution leads to a significant deterioration of the accuracy, with errors up to 90 µEh. In contrast, the exchange matrix from pure fp32 evaluation is much more accurate than the errors in the energy indicate. That is, even though converging the SCF using fp32 only can lead to quite significant errors in the converged density matrix of up to 4 × 10⁻⁵ a.u., computing the final energy from this inaccurate matrix employing mixed precision integral evaluation instead of pure fp32-execution (denoted as “fp32^*”) leads to significantly smaller errors (1.8 µEh at most) while providing the performance of full fp32 execution (i.e., up to 7× speedups) within the SCF.

The situation is analogous to adaptive integration grids, i.e., employing three times coarser grids within the SCF than for the final energy calculation provides virtually the same result as performing the whole calculation with the large grid.⁴¹ That is, the exchange energy is generally more sensitive to numeric errors than the exchange matrix.

However, as presented in Table V, executing all K-builds within the SCF with pure fp32 arithmetic deteriorates the SCF convergence behavior measurably, meaning that very tight convergence, e.g., to 10⁻⁸, is impossible due to the increase in numerical fluctuations. Although the level of precision is sufficient for many applications (e.g., ab initio molecular dynamics), a truly error-free method similar to part 1 (i.e., no measurable change in the SCF convergence behavior and as little change in the converged energy and density matrix as possible) is our goal in the following.

The idea of incremental Fock/exchange-builds, i.e., computing K[ΔP] and then incrementally updating K instead of always recomputing the full exchange matrix, was initially proposed by Almloef et al.⁵⁵ and later improved by Haeser and Ahlrichs⁵⁶ to allow for tighter density-matrix-based screening of the 4c2e-integrals within the SCF. However, because only small increments to K are computed, incremental Fock-builds have also been shown to reduce the numerical error introduced by fp32-evaluation, since the absolute error is proportional to the magnitude of the contribution.³¹ In part 2 (Sec. IV B), we showed that pure fp32 exchange-builds are already quite accurate and allow for unproblematic convergence to about 10⁻⁶. Therefore, the possibility to perform the last exchange-builds incrementally should be explored.

We thus propose i-sn-LinK, a special SCF scheme given in Fig. 2. First, the SCF is converged to a relative loose threshold with non-incremental fp32 K-builds. Since the difference density is comparatively large for these early SCF steps, there are no significant performance gains from incremental updates in these steps. In this work, we choose a convergence criterion of 10⁻⁵ to ensure that convergence to this point is always reached.

Then, one full K-build is performed with mixed fp32/fp64 precision (i.e., employing part 1) and the SCF is subsequently converged further with incremental fp32 K-builds. These incremental K-builds are always built from the difference density matrix ΔP referenced to the last full K-build, not to the previous density matrix. This removes the possibility of incremental error accumulation, further improving the numerical stability in this way.⁵⁷ In order to reduce errors from the density-matrix including integral screening when operating on incremental matrices, we employ two orders of magnitude tighter screening thresholds ϑ_E/K for these incremental builds, which, due to the elements of ΔP being much smaller than P, leads to a similar performance as full K-builds with looser thresholds. The final energy (also optionally nuclear forces) is then computed using mixed fp32/fp64 precision, analogous to part 2.

In contrast to the works of Almloef et al.⁵⁵ and Haeser and Ahlrichs,⁵⁶ the aim of our incremental scheme i-sn-LinK is only to improve the convergence behavior while using as many fp32 K-builds as possible, not to improve the tightness of the integral-selection.

The performance, accuracy, and SCF convergence of i-sn-LinK are presented in Table VI. Overall, we obtain virtually the same accuracy and numerical stability as for pure fp64 execution while obtaining up to 5.2× speedups (averaged over the whole calculation). As expected, the speedups are higher if more SCF cycles need to be performed because the cost of the two remaining mixed-precision K-builds, which comprise 36%–47% of the computation time, is comparatively less impactful in this situation.

However, the most significant result of Table VI is the fact that i-sn-LinK nearly always converges within the same amount of SCF cycles as pure fp64 execution, proving the numerical stability of the method.

In order to compare the different approaches discussed in Secs. IV A–IV C, a brief summary of them is given in Table VII. In this context, the mixed-precision 3c1e integral method (“mp”) of part 1 represents the simplest introduction of reduced numerical precision; that is, only the 3c1e integral evaluation part [Eq. (6)] is adjusted. This approach has virtually no impact on the SCF convergence behavior or the accuracy of the final result. Since the 3c1e integral evaluation is by far the most significant bottleneck for CPU execution, accelerating the other steps cannot provide significant speedups. Hence, we recommend to only employ the “mp” approach on CPUs.

On GPUs, however, the BLAS-3 steps with fp64 execution can become significant. Therefore, we investigated the possibility of pure fp32-execution in Sec. IV B (part 2) and found that the so-converged SCF can yield surprisingly accurate energies, if the final energy was evaluated with higher precision. However, the convergence behavior and the quality of the converged density matrix were measurably deteriorated.

These shortcomings were then addressed in Sec. IV C (i-sn-LinK; part 3) by introducing incremental K-builds, which update a once-computed high-precision K-matrix. This recovers the accuracy and stable convergence behavior of the “mp” approach, but only one non-fp32 K-build has to be performed within the SCF. The final energy is, of course, never computed with pure single-precision, since this leads to unacceptably large errors (cf. “fp32” in Table IV). Because the i-sn-LinK scheme represents the best trade-off between accuracy and performance on GPUs, i.e., close to pure-fp32 performance and essentially pure-fp64 accuracy, we recommend the i-sn-LinK method for GPUs.

Since deterioration in SCF convergence stability was the main reason to develop the i-sn-LinK method of part 3 (Sec. IV C) and since it was also shown to be problematic in other works on single-precision execution,^31–36 we provide a more detailed investigation on difficult molecules regarding SCF stability in Table VIII. In order to exemplify the worst-case scenario, we selected the ten molecules with the worst convergence behavior of the ASCDB benchmark set⁵⁸ and employ the def-TZVPPD basis set, which, due to the additional diffuse basis-functions, requires even higher numerical precision.

The results verify that both the mixed precision 3c1e evaluation and i-sn-LinK have no impact on the stability of the SCF convergence, even when considering very difficult electronic structures. Furthermore, the largest error of 0.83 µE_h occurs for the most difficult structure (cisHO₃) and is still smaller than the SCF convergence criterion of 10⁻⁶.

To illustrate the practical applicability of the mixed-precision methods developed in this work, a range of systems are tested with up to quadruple-ζ basis sets in Table IX.

First, the reliability and accuracy of the two methods is proven again, i.e., nearly all molecules converge within the same amount of SCF steps (10–12 steps) regardless of the method used for the Fock-exchange. Furthermore, the converged energies match the fp64 results to <0.1 µE_h even in the case of heavy elements (e.g., yttrium-complexes in Table XI), the nuclear forces are accurate within 0.7 µE_h a₀⁻¹, and the converged density matrices are accurate within 10⁻⁶ a.u. except for C₆₀.

For the C₆₀-Fullerene, larger deviations (∼10⁻⁵ a.u.) are, however, not directly due to fp32-execution but are instead caused by the slightly different integral screening in i-sn-link, where ΔP is employed instead of P within the incremental updates. This is supported by Table X, where tighter thresholds are employed for both the fp64-reference and the i-sn-LinK calculations in order to remove the effect of the integral screening. The remaining deviations are then only caused by the reduced numerical precision in i-sn-LinK. These remaining errors are significantly smaller for the most challenging molecules such as C₆₀ or (DNA)₄/QZVP, proving that those larger deviations in Table IX were indeed caused by the altered integral screening and not by the reduced numerical precision.

Summarizing the accuracy comparisons, the two tested mixed-precision methods “mp” and i-sn-LinK lead to essentially negligible errors (ΔE < 0.1 µEh, ΔP ≤ 10⁻⁶, ΔForces <1 µEh a₀⁻¹), which are multiple orders of magnitude smaller than “chemical accuracy”(∼1000 µEh).

Meanwhile, the methods lead to considerable speedups compared to pure fp64-execution, as presented in Table XI: 1.4–1.8× speedups are obtained from the mixed-precision method on CPUs, around 2× speedups are obtained from the mixed-precision method on GPUs, and up to 6.9× speedups are obtained with the incremental i-sn-LinK method on GPUs. The speedups of i-sn-LinK are typically larger for more expensive computations (larger molecules and larger basis sets) because the BLAS-3 steps [Eqs. (5) and (7)] are comparatively more expensive in these situations.

Comparing the CPU and GPU performance, we note that our incremental i-sn-LinK method is essential for an effective GPU acceleration with gaming GPUs. That is, pure fp64 execution is about as fast on the GPU as on CPUs,⁵⁹ but i-sn-LinK is up to 4.4× faster on the GPU [333 vs 1475 s for (S₈)₂₀] than on CPUs.

This work presents a systematic study of the applicability of single-precision instructions to evaluate the Fock-exchange matrix within seminumerical integration schemes. First, we demonstrated that only a very small fraction of the 3c1e integrals needs to be evaluated with fp64 instructions to still provide virtually the same result while providing nearly 2× speedups for this step on CPUs. We then demonstrated that pure fp32 execution still leads to surprisingly accurate results, as long as the final energy is computed at higher accuracy, although the SCF convergence behavior is somewhat deteriorated. Finally, we proposed the i-sn-LinK method, where incremental exchange-builds are employed for the last SCF steps, which removes all instabilities and inaccuracies from pure fp32 execution while requiring only a single high precision K-matrix build. This method provides up to 7× speedups for the whole SCF on typical “gaming” GPUs (1080Ti) without any significant impact on the numerical stability or accuracy and is therefore essential for an effective GPU acceleration using more affordable “gaming” GPUs.

Finally, we note that the present study on single-precision execution within seminumerical integration is also relevant for the evaluation of local-hybrid functionals, where the most time consuming step, i.e., the evaluation of the local exchange contributions, is identical to the seminumerical expression for the exchange matrix. Thus, we expect the performance benefits presented in this work to also translate to this novel and exciting class of functionals.

See the supplementary material for the effect of vector-instructions on the timings in Table II, absolute timings associated with the speedups of Tables II–IV and VI, analogous results to Tables IX and XI for the PBE0-functional, and all xyz-structures employed in this work.

The authors acknowledge financial support by the “Deutsche Forschungsgemeinschaft” (DFG) via the Grant No. SFB 1309-32587107 and the Cluster of Excellence “Munich Center for Quantum Science and Technology” (MCQST, Grant No. EXC2111-390814868). C.O. further acknowledges financial support as a Max-Planck-Fellow at the MPI-FKF, Stuttgart.

The data that support the findings of this study are available within the article and its supplementary material.