Four-center Integrals Research Articles

Benchmarking the accuracy of the separable resolution of the identity approach for correlated methods in the numeric atom-centered orbitals framework.

Four-center two-electron Coulomb integrals routinely appear in electronic structure algorithms. The resolution-of-the-identity (RI) is a popular technique to reduce the computational cost for the numerical evaluation of these integrals in localized basis-sets codes. Recently, Duchemin and Blase proposed a separable RI scheme [J. Chem. Phys. 150, 174120 (2019)], which preserves the accuracy of the standard global RI method with the Coulomb metric and permits the formulation of cubic-scaling random phase approximation (RPA) and GW approaches. Here, we present the implementation of a separable RI scheme within an all-electron numeric atom-centered orbital framework. We present comprehensive benchmark results using the Thiel and the GW100 test set. Our benchmarks include atomization energies from Hartree-Fock, second-order Møller-Plesset (MP2), coupled-cluster singles and doubles, RPA, and renormalized second-order perturbation theory, as well as quasiparticle energies from GW. We found that the separable RI approach reproduces RI-free HF calculations within 9meV and MP2 calculations within 1meV. We have confirmed that the separable RI error is independent of the system size by including disordered carbon clusters up to 116 atoms in our benchmarks.

High-Performance Evaluation of High Angular Momentum 4-Center Gaussian Integrals on Modern Accelerated Processors.

We present a high-performance evaluation method for 4-center 2-particle integrals over Gaussian atomic orbitals with high angular momenta (l ≥ 4) and arbitrary contraction degrees on graphical processing units (GPUs) and other accelerators. The implementation uses the matrix form of McMurchie-Davidson recurrences. Evaluation of the four-center integrals over four l = 6 (i) Gaussian AOs in double precision (FP64) on an NVIDIA V100 GPU outperforms the reference implementation of the Obara-Saika recurrences (Libint) running on a single Intel Xeon core by more than a factor of 1000, easily exceeding the 73:1 ratio of the respective hardware peak FLOP rates while reaching almost 50% of the V100 peak. The approach can be extended to support AOs with even higher angular momenta; for lower angular momenta (l ≤ 3), additional improvements will be reported elsewhere. The implementation is part of an open-source LibintX library freely available at github.com:ValeevGroup/LibintX.

Expansion of multicenter Coulomb integrals in terms of two-center integrals

Three- and four-center Coulomb integrals in the solid spherical harmonic Gaussian basis are solved by expansion in terms of two-center integrals. The two-electron Gaussian product rule, coupled with the addition theorem for solid spherical harmonics, reduces four-center Coulomb integrals into a linear combination of two-center Coulomb integrals and one-center overlap integrals. With this approach, three- and four-center Coulomb integrals can be reduced to the same form of two-center integrals. Resulting two-center Coulomb integrals can be further simplified into a simpler form, which can be related to the Boys function. Multi-center Coulomb integrals are solved hierarchically: simple two-center Coulomb integrals are used for calculation of more complicated two-center Coulomb integrals, which are used in the calculation of multicenter integrals.

Tight distance-dependent estimators for screening two-center and three-center short-range Coulomb integrals over Gaussian basis functions.

We derive distance-dependent estimators for two-center and three-center electron repulsion integrals over a short-range Coulomb potential, erfc(ωr12)/r12. These estimators are much tighter than the ones based on the Schwarz inequality and can be viewed as a complement to the distance-dependent estimators for four-center short-range Coulomb integrals and for two-center and three-center full Coulomb integrals previously reported. Because the short-range Coulomb potential is commonly used in solid-state calculations, including those with the Heyd-Scuseria-Ernzerhof functional and with our recently introduced range-separated periodic Gaussian density fitting, we test our estimators on a diverse set of periodic systems using a wide range of the range-separation parameter ω. These tests demonstrate the robust tightness of our estimators, which are then used with integral screening to calculate periodic three-center short-range Coulomb integrals with linear scaling in system size.

Optimizing Molecular Geometries in Strong Magnetic Fields.

An efficient implementation of geometrical derivatives at the Hartree–Fock (HF) and current-density functional theory (CDFT) levels is presented for the study of molecular structure in strong magnetic fields. The required integral derivatives are constructed using a hybrid McMurchie–Davidson and Rys quadrature approach, which combines the amenability of the former to the evaluation of derivative integrals with the efficiency of the latter for basis sets with high angular momentum. In addition to its application to evaluating derivatives of four-center integrals, this approach is also applied to gradients using the resolution-of-the-identity approximation, enabling efficient optimization of molecular structure for many-electron systems under a strong magnetic field. The CDFT contributions have been implemented for a wide range of density functionals up to and including the meta-GGA level with current-density dependent contributions and (range-separated) hybrids for the first time. Illustrative applications are presented to the OH and benzene molecules, revealing the rich and complex chemistry induced by the presence of an external magnetic field. Challenges for geometry optimization in strong fields are highlighted, along with the requirement for careful analysis of the resulting electronic structure at each stationary point. The importance of correlation effects is examined by comparison of results at the HF and CDFT levels. The present implementation of molecular gradients at the CDFT level provides a cost-effective approach to the study of molecular structure under strong magnetic fields, opening up many new possibilities for the study of chemistry in this regime.

Multicomponent Coupled Cluster Singles and Doubles with Density Fitting: Protonated Water Tetramers with Quantized Protons.

Nuclear quantum effects such as zero-point energy are important for describing a wide range of chemical properties. The nuclear-electronic orbital (NEO) approach incorporates such effects into quantum chemistry calculations by treating specified nuclei, typically protons, quantum mechanically on the same level as electrons. Herein, both the traditional and t1-transformed NEO coupled cluster with singles and doubles (NEO-CCSD) methods are implemented with a density fitting (DF) scheme for approximating the four-center two-particle integrals. The enhanced computational efficiency enables calculations on larger molecules with multiple quantum protons. The NEO-DF-CCSD method predicts proton affinities within chemical accuracy. Its application to protonated water tetramers with all nine protons treated quantum mechanically produces the qualitatively correct ordering of the isomer energies, which are strongly influenced by the zero-point energy contributions inherently included in NEO energy calculations. This work showcases the capabilities of the NEO-DF-CCSD method and provides the foundation for future developments and applications.

Efficient evaluation of exact exchange for periodic systems via concentric atomic density fitting.

The evaluation of the exact [Hartree-Fock (HF)] exchange operator is a crucial ingredient for the accurate description of the electronic structure in periodic systems through ab initio and hybrid density functional approaches. An efficient formulation of periodic HF exchange in a linear combination of atomic orbitals representation presented here is based on the concentric atomic density fitting approximation, a domain-free local density fitting approach in which the product of two atomic orbitals is approximated using a linear combination of fitting basis functions centered at the same nuclei as the AOs in that product. A significant reduction in the computational cost of exact exchange is demonstrated relative to the conventional approach due to avoiding the need to evaluate four-center two-electron integrals, with sub-millihartree/atom errors in absolute HF energies and good cancellation of fitting errors in relative energies. The novel aspects of the evaluation of the Coulomb contribution to the Fock operator, such as the use of real two-center multipole expansions and spheropole-compensated unit cell densities, are also described.

Computing Four-Center Two-Electron Coulomb Integrals Using Exponential Transformations and Trapezoidal Rule

The numerical evaluations of the four-center two-electron Coulomb integrals are among the most time-consuming computations involved in molecular electronic structure calculations. In the present paper we extend the double exponential (DE) transform method, previously developed for the numerical evaluation of threecenter one-electron molecular integrals [J. Lovrod, H. Safouhi, Molecular Physics (2019) DOI:10.1030/0026867.2019.1619854], to four-center two-electron integrals. The fast convergence properties analyzed in the numerical section illustrate the advantages of the new approach.

Evaluating two-electron-repulsion integrals over arbitrary orbitals using zero variance Monte Carlo: Application to full configuration interaction calculations with Slater-type orbitals

A Monte Carlo method for evaluating multicenter two-electron-repulsion integrals over any types of orbitals (Slater, Sturmian, finite-range, numerical, etc.) is presented. The approach is based on a simple and universal (orbital-independent) gaussian sampling of the two-electron configuration space and on the use of efficient zero-variance Monte Carlo estimators. Quite remarkably, it is shown that the high level of accuracy required on two-electron integrals to make Hartree-Fock (HF) and configuration interaction (CI) calculations feasible can be achieved. A first zero-variance estimator is built by introducing a gaussian approximation of the orbitals and by evaluating the two-electron integrals using a correlated sampling scheme for the difference between exact and approximate orbitals. A second one is based on the introduction of a general coordinate transformation. The price to pay for this simple and general Monte Carlo scheme is the high computational cost required. However, we argue that the great simplicity of the algorithm, its embarrassingly parallel nature, its ideal adaptation to modern computational platforms and, most importantly, the possibility of using more compact and physically meaningful basis sets nevertheless make the method attractive. HF and near full CI (FCI) calculations using Slater-type orbitals (STOs) are reported for Be, CH4, and [H2N(CH)NH2]+ (a simple model of cyanine). To the best of our knowledge, our largest FCI calculation involving 18 active electrons distributed among 90 orbitals for the cyanine molecule is the most extensive molecular calculation performed so far using pure STOs (no gaussian approximation, even for the challenging four-center two-electron integrals).

Range-Separated Hybrid Functionals with Variational Fitted Exact Exchange

This work presents a variationally fitted long-range exact exchange algorithm that can be used for the computation of range-separated hybrid density functionals in the linear combination of Gaussian type orbital (LCGTO) approximation. The obtained LCGTO energy and gradient expressions are free of four-center integrals and employ modified three-center integral recurrence relations to obtain optimal computational performance. The accuracy and performance of selected range-separated hybrid functionals with variational fitted long-range exact exchange are analyzed and discussed. A parallel implementation is presented and benchmarked.

Robust Periodic Fock Exchange with Atom-Centered Gaussian Basis Sets.

In this work, we present a robust implementation of the periodic Fock exchange for atom-centered Gaussian-type orbitals (GTOs). We discuss the divergence, appearing in the formulation of the periodic Fock exchange in the case of a finite number of k-points, and compare two schemes that remove it. These are the minimum image convention (MIC) and the truncated Coulomb interaction (TCI) that we use here in combination with k-meshes. We observe artifacts in four-center integrals of GTOs, when evaluated in the TCI scheme. They carry over to the exchange and density matrices of Hartree-Fock calculations for TCI but are absent in MIC. At semiconducting and insulating systems, we show that both MIC and TCI yield the same energies for a sufficiently large supercell or k-mesh, but the self-consistent field algorithm is more stable for MIC. We therefore conclude that the MIC is superior to TCI and validate our implementation by comparing not only to other GTO-based calculations but also by demonstrating excellent agreement with results of plane-wave codes for sufficiently large Gaussian basis sets.

Distance-including rigorous upper bounds and tight estimates for two-electron integrals over long- and short-range operators.

We introduce both rigorous and non-rigorous distance-dependent integral estimates for four-center two-electron integrals derived from a distance-including Schwarz-type inequality. The estimates are even easier to implement than our so far most efficient distance-dependent estimates [S. A. Maurer et al., J. Chem. Phys. 136, 144107 (2012)] and, in addition, do not require well-separated charge-distributions. They are also applicable to a wide range of two-electron operators such as those found in explicitly correlated theories and in short-range hybrid density functionals. For two such operators with exponential distance decay [e-r12 and erfc(0.11⋅r12)/r12], the rigorous bound is shown to be much tighter than the standard Schwarz estimate with virtually no error penalty. The non-rigorous estimate gives results very close to an exact screening for these operators and for the long-range 1/r12 operator, with errors that are completely controllable through the integral screening threshold. In addition, we present an alternative form of our non-rigorous bound that is particularly well-suited for improving the PreLinK method [J. Kussmann and C. Ochsenfeld, J. Chem. Phys. 138, 134114 (2013)] in the context of short-range exchange calculations.

Dispersion-Corrected Spin-Component-Scaled Double-Hybrid Density Functional Theory: Implementation and Performance for Non-covalent Interactions.

The implementation of 300 combinations of generalized gradient approximation/local density approximation exchange-correlation dispersion-corrected spin-component-scaled double-hybrid (DSD) density functional theory (DFT) methods has been carried out and the performance assessed against several DFT and post-Hartree-Fock methods, enabling further advancements toward the long-standing challenge of accurate prediction of interaction energies and associated properties. The resulting framework is flexible and has been further extended to include the resolution of identity (RI) approximation for solving the critical four-center two-electron repulsion integrals in the basis of the Kohn-Sham orbitals for cost effectiveness. To evaluate the performance of this set of new cost-effective methods, denoted as RI-DSD-DFTs, seven validation data sets were designed to cover a broad range of non-covalent interactions with characteristic stabilizing contributions. Inclusion of the perturbative treatment of correlation effects is shown to significantly improve the description of weak interactions. The set of DSD-DFTs provide interaction energies with root-mean-square deviations and mean absolute errors within 0.5 kcal/mol. The cost-effective RI-DSD-DFT counterparts deviate by less than 0.18 kcal/mol on average with only 2% of the computational cost.

Hybrid and Constrained Resolution-of-Identity Techniques for Coulomb Integrals.

The introduction of auxiliary bases to approximate molecular orbital products has paved the way to significant savings in the evaluation of four-center two-electron Coulomb integrals. We present a generalized dual space strategy that sheds a new light on variants over the standard density and Coulomb-fitting schemes, including the possibility of introducing minimization constraints. We improve in particular the charge- or multipole-preserving strategies introduced respectively by Baerends and Van Alsenoy that we compare to a simple scheme where the Coulomb metric is used for lowest angular momentum auxiliary orbitals only. We explore the merits of these approaches on the basis of extensive Hartree-Fock and MP2 calculations over a standard set of medium size molecules.

Accurate localized resolution of identity approach for linear-scaling hybrid density functionals and for many-body perturbation theory

A key component in calculations of exchange and correlation energies is the Coulomb operator, which requires the evaluation of two-electron integrals. For localized basis sets, these four-center integrals are most efficiently evaluated with the resolution of identity (RI) technique, which expands basis-function products in an auxiliary basis. In this work we show the practical applicability of a localized RI-variant (‘RI-LVL’), which expands products of basis functions only in the subset of those auxiliary basis functions which are located at the same atoms as the basis functions. We demonstrate the accuracy of RI-LVL for Hartree–Fock calculations, for the PBE0 hybrid density functional, as well as for RPA and MP2 perturbation theory. Molecular test sets used include the S22 set of weakly interacting molecules, the G3 test set, as well as the G2–1 and BH76 test sets, and heavy elements including titanium dioxide, copper and gold clusters. Our RI-LVL implementation paves the way for linear-scaling RI-based hybrid functional calculations for large systems and for all-electron many-body perturbation theory with significantly reduced computational and memory cost.

Hermitian “chemical” Hamiltonian: an alternative ab initio method

Some previous results of the present author are combined in order to develop a Hermitian version of the “Chemical Hamiltonian Approach.” In this framework the second quantized Born–Oppenheimer Hamiltonian is decomposed into one- and two-center components, if some finite basis corrections are omitted. (No changes are introduced into the one- and two-center integrals, while projective expansions are used for the three- and four-center ones, which become exact only in the limit of complete basis sets.) The total molecular energy calculated with this Hamiltonian can then presented as a sum of the intraatomic and diatomic energy terms which were introduced in our previous “chemical energy component analysis” scheme. The corresponding modified Hartree–Fock–Roothaan equations are also derived; they do not contain any three- and four-center integrals, while the non-empirical character of the theory is conserved. This scheme may be useful also as a “layer” in approaches like ONIOM.

The cyclic cluster model at hartree–fock level

We present the implementation of the cyclic cluster model (CCM) formalism at the Hartree-Fock (HF) level. In contrast to other periodic models, the CCM is a Γ-point approach. Integration is carried out in real space within a finite interaction area determined by the size and the shape of the cluster that corresponds to a supercell of the solid, surface, or polymer. Particular care has to be taken for the proper treatment of three- and four-center integrals that involve basis functions located at the boundaries of the Wigner-Seitz supercell, which defines the interaction region. The similarity between the CCM formalism and molecular approaches allows in principle the application of sophisticated post HF methods to solid-state problems with only moderate modification of the molecular code. We show for selected model systems, that with our approach, the electronic structure and energetics obtained by the conventional supercell model is fully reproduced.

Efficient two-component self-consistent field procedures and gradients: implementation in TURBOMOLE and application to

The implementation of gradients of two-component density functional (DFT) and hybrid-DFT energies is presented as well as several improvements with respect to our previous version in the self-consistent field (SCF) procedure: Calculation of the Hartree–Fock (HF) exchange with four-centre integrals, use of difference density formalism and partial exploitation of symmetry. Efficiency is demonstrated by comparison to recently reported data and by calculation of energy and gradient of Au147. As an example of use, structure parameters of Au− 20 were optimised and subsequently vibration frequencies were calculated numerically from the analytical gradients at two-component level. For this example, spin–orbit coupling leads to higher degeneracy of vibrational levels by lifting degeneracy of electronic levels, as the latter prevents Jahn–Teller distortion.

Attractive electron–electron interactions within robust local fitting approximations

An analysis of Dunlap's robust fitting approach reveals that the resulting two-electron integral matrix is not manifestly positive semidefinite when local fitting domains or non-Coulomb fitting metrics are used. We present a highly local approximate method for evaluating four-center two-electron integrals based on the resolution-of-the-identity (RI) approximation and apply it to the construction of the Coulomb and exchange contributions to the Fock matrix. In this pair-atomic resolution-of-the-identity (PARI) approach, atomic-orbital (AO) products are expanded in auxiliary functions centered on the two atoms associated with each product. Numerical tests indicate that in 1% or less of all Hartree-Fock and Kohn-Sham calculations, the indefinite integral matrix causes nonconvergence in the self-consistent-field iterations. In these cases, the two-electron contribution to the total energy becomes negative, meaning that the electronic interaction is effectively attractive, and the total energy is dramatically lower than that obtained with exact integrals. In the vast majority of our test cases, however, the indefiniteness does not interfere with convergence. The total energy accuracy is comparable to that of the standard Coulomb-metric RI method. The speed-up compared with conventional algorithms is similar to the RI method for Coulomb contributions; exchange contributions are accelerated by a factor of up to eight with a triple-zeta quality basis set. A positive semidefinite integral matrix is recovered within PARI by introducing local auxiliary basis functions spanning the full AO product space, as may be achieved by using Cholesky-decomposition techniques. Local completion, however, slows down the algorithm to a level comparable with or below conventional calculations.

Distance-dependent Schwarz-based integral estimates for two-electron integrals: Reliable tightness vs. rigorous upper bounds

A new integral estimate for four-center two-electron integrals is introduced that accounts for distance information between the bra- and ket-charge distributions describing the two electrons. The screening is denoted as QQR and combines the most important features of the conventional Schwarz screening by Häser and Ahlrichs published in 1989 [J. Comput. Chem. 10, 104 (1989)] and our multipole-based integral estimates (MBIE) introduced in 2005 [D. S. Lambrecht and C. Ochsenfeld, J. Chem. Phys. 123, 184101 (2005)]. At the same time the estimates are not only tighter but also much easier to implement, so that we recommend them instead of our MBIE bounds introduced first for accounting for charge-distance information. The inclusion of distance dependence between charge distributions is not only useful at the SCF level but is particularly important for describing electron-correlation effects, e.g., within AO-MP2 theory, where the decay behavior is at least 1/R(4) or even 1/R(6). In our present work, we focus on studying the efficiency of our QQR estimates within SCF theory and demonstrate the performance for a benchmark set of 44 medium to large molecules, where savings of up to a factor of 2 for exchange integrals are observed for larger systems. Based on the results of the benchmark set we show that reliable tightness of integral estimates is more important for the screening performance than rigorous upper bound properties.