Numerical Integration of Slater Basis Functions Over Prolate Spheroidal Grids

TL;DR

This paper introduces a new numerical integration method using prolate spheroidal grids to improve accuracy in electronic structure simulations with large basis sets.

Contribution

The novel use of prolate spheroidal grids enables accurate integration for larger basis sets in polyatomic systems.

Findings

The new grid reduces integration errors by ~3 orders of magnitude compared to Becke partitioning.

The method is tested on polyatomic systems and shows reliability for correlated electronic structure computations.

Abstract

Slater basis functions have desirable properties that can improve electronic structure simulations, but improved numerical integration methods are needed. This work builds upon the SlaterGPU library for the evaluation of Hamiltonian matrix elements in the resolution‐of‐the‐identity approximation. In particular, a prolate spheroidal grid will provide sufficient integral accuracy to employ larger basis sets (quadruple‐zeta and greater) in practical computations involving polyatomics. To integrate 3‐center Coulomb and nuclear attraction terms, an improved grid representation around the third center is introduced. The RMSEs of the integral quantities are evaluated and compared to the previous numerical integration method used in SlaterGPU (Becke partitioning), resulting in a ~3 order of magnitude reduction in the error for 2‐center integral quantities. The procedure is generally applicable to polyatomic systems, GPU accelerated for high performance computing, and tested on self‐consistent field and full configuration interaction wavefunctions. Results for a number of 3‐atom models as well as propanediyl (C3H6) demonstrate the reliability of the new integration scheme.

Atom‐centered Slater functions are integrated over a prolate spheroidal grid, allowing the use of larger basis sets in correlated electronic structure computations.

Figures

The prolate spheroidal grid on a cartesian grid with spheroids at constant μ (red), hyperboloids at constant ν (blue), and planes at constant ϕ (green).

The process of reorienting the grid around the third center, (A) is the initial spacing of grid with the third center being located in one of the volume elements, (B) moves the μ and ν grid lines closest to the third center so that the third center coordinates align with the grid lines, and (C) splits the grid further around the third center for greater accuracy.

Overlap, kinetic energy, and 2‐center Coulomb repulsion RMSE (Ha) analyzed for ClF and OF2 using the PS method in a TZ basis.

Electron–nuclear attraction and 3‐center Coulomb repulsion RMSE (Ha) analyzed for ClF and OF2 using the PS method in a TZ basis. The grid selection and legend are the same as in Figure 3.

Comparison of the RMSE (Ha) of the electron‐nuclear attraction elements for different divisions of the third center within the PS integration grid. A single split (NSP=2) gives significant advantages over no splitting, and NSP>2 gives marginal improvements.

Six 2‐center Coulomb integrals plotted along 16 vector directions [16] as the two centers are separated. The legend indicates unit vectors for all 16 directions tested, all basis functions used have the exponent ζ=1, and all basis functions were chosen to have m=0.

HONO and LUNO orbitals and occupation numbers for the singlet and triplet states of 1,3‐propanediyl as well as the singlet–triplet gap for the vertical excitation. iFCI uses the STO QZ basis while CCSD(T) uses a GTO basis, cc‐pVQZ.

Tables

TABLE 1. RMSE (Ha) of integrals for SO2 using BP and PS integration with TZ, QZ, and 5Z basis sets. Basis sets TZ, QZ, 5Z are described in the computational details.

	TZ		QZ		5Z
	BP	PS	BP	PS	BP	PS
Overlap	1.3 × 10−08	8.0 × 10−12	7.6 × 10−09	5.4 × 10−12	1.3 × 10−08	8.9 × 10−12
Kinetic	4.6 × 10−07	9.6 × 10−11	3.1 × 10−07	3.0 × 10−10	3.0 × 10−07	2.8 × 10−10
Electron nuclear	7.7 × 10−07	3.0 × 10−07	5.8 × 10−07	3.1 × 10−07	6.0 × 10−07	3.2 × 10−07
2‐Center Coulomb	3.2 × 10−08	4.9 × 10−11	2.0 × 10−08	6.8 × 10−11	3.1 × 10−08	8.1 × 10−11
3‐Center Coulomb	8.7 × 10−09	6.6 × 10−09	5.6 × 10−09	1.3 × 10−08	7.3 × 10−09	1.6 × 10−08

TABLE 2. Convergence of Hartree–Fock SCF runs with BP integration. Each number indicates the eigenvalue threshold needed to control numerical stability of the overlap matrix inversion. If the value is 0 (green) all threshold values tested allow for SCF convergence. Higher values indicate instability at this step, due to lower accuracy of the overlap matrix elements. Oscillatory cases (orange) oscillated around the energy obtained through the PS method with a difference less than 0.1 Ha but did not converge. Energies which varied by more than 0.1 Ha are marked in red. In comparison, PS integration converged for all cases with the smallest eigenvalue threshold.

Mol	TZ	QZ	5Z	Mol	TZ	QZ	5Z	Mol	TZ	QZ	5Z
HF	0	0	10−5	CS	0	10−6	10−6	N3	0	0	10−6
C2	0	0	0	SiO	0	10−6	0	N2O	0	0	10−7
N2	0	0	0	ClF	0	0	0	NO2	0	0	0
CO	0	0	0	Cl2	0	0	0	SFH	0	10−7	10−7
NO	0	0	0	H2O	0	0	0	OF2	0	0	0
O2	0	0	0	HCN	0	0	0	SO2	0	10−7	10−7
HCl	0	0	0	OFH	0	0	0	CS2	0	10−7	10−7
	Success				Oscillatory				Far from Convergent

TABLE 3. Heat‐bath configuration interaction (HBCI) correlation energies (Ha) for various molecules in TZ, QZ, and 5Z basis sets.

Mol	TZ	QZ	5Z
HF	−0.2558	−0.2821	−0.3005
C2	−0.3629	−0.3861	−0.3895
N2	−0.3646	−0.3960	−0.4091
CO	−0.3477	−0.3829	−0.3935
NO	−0.6962	−0.7309	−0.7487
O2	−0.4641	−0.5138	−0.5299
HCl	−0.1910	−0.2373	−0.2417
CS	−0.3041	−0.3485	−0.3526
SiO	−0.3247	−0.3646	−0.3738
ClF	−0.4191	−0.4750	−0.4948
Cl2	−0.3490	−0.4399	−0.4482
H2O	−0.2517	−0.2821	−0.2912
HCN	−0.3402	−0.3730	−0.3811
OFH	−0.4834	−0.5374	−0.5632
N3	−0.9796	−1.0304	−1.0537
N2O	−0.6043	−0.6643	−0.6875
NO2	−1.1706	−1.2314	−1.2607
SFH	−0.4117	−0.4731	−0.4940
OF2	−0.7200	−0.7955	−0.8395
SO2	−0.6140	−0.6965	−0.7196
CS2	−0.4646	−0.5419	−0.5503

TABLE 4. Two states of CH2 (1A1 and 3B1), at various correlated levels of theory. Geometries for 1A1 and 3B1 CH2 were obtained from Ref. [81]. The bases used for the GTO calculations are the cc‐pVXZ basis sets, while the STO bases are the same as referenced in Table 2. The energy reported under CBS was obtained in the complete basis set limit. Units are in kcal/mol.

Method	TZ	QZ	5Z	CBS f	Expt.	Basis
HBCI‐STO	11.1	9.93	9.72			Slater
HBCI‐GTO	10.5	9.71	9.51			Gaussian
Expt. a					9.0
Expt. b					9.36
SF‐CIS c	20.4					Gaussian
SF‐CIS(D) c	14.1					Gaussian
EOM‐SF‐CCSD(dT) d		9.7				Gaussian
VMC e	9.92					Slater
DMC e	9.36					Slater
CMRCI + Q f				8.97		Gaussian
FCI g	11.1					Gaussian

Full text

Introduction

1

Electronic structure theory can provide insight into a countless number of chemical systems. Practical methods for molecules most often rely on wavefunctions built upon atom‐centered, single‐electron basis functions. Certain asymptotic properties of exact wavefunctions can be captured in the single‐particle basis set [1], such as cusps at the nuclei that are consistent with the Kato conditions [2, 3]. At long distances, wavefunctions decay as simple exponentials. These two ideas lead naturally to Slater type orbitals (STOs) [1], which have the form

$$[eqn]$$

N is a normalization factor, ξ is the exponent defining each basis function, n,l,m are the principal quantum numbers, and Ylm are spherical harmonics [2, 4]. STOs are known to provide accurate descriptions of polarizability, intermolecular interactions, and nuclear shielding, as these properties are sensitive to cusp or decay of the wavefunction [5, 6, 7]. The widespread use of Slater functions as electronic structure basis sets, however, has been significantly hindered due to the requirement for numerical integration.

While STOs correctly describe physical properties, their features—such as steepness near the nucleus and long tails—make them difficult to numerically integrate. While some Slater integrals are available analytically [8, 9, 10], the electronic Hamiltonian for polyatomics requires Slater integrals which have no known analytical form. Gaussian Type Orbitals (GTOs), in contrast, are generally analytically integrable [11, 12, 13, 14, 15] but cannot precisely capture nuclear cusps or long‐range decays. Improved integration protocols could start to bridge the practicality gap between STOs and GTOs, making STOs more readily usable in a wide range of electronic structure theories. This work therefore builds upon our recent efforts [16] to provide highly accurate Slater integrals, especially for their use in wavefunction simulations.

Integration of Slater functions has often been performed using atom‐centered integration grids, for instance as done in the Amsterdam density functional (ADF) program [17, 18, 19, 20, 21], and more recently through the SlaterGPU library [16]. These integrals use products of radial and angular grids on each atom, where the grid weights are adjusted to avoid overcounting in regions where the grids are overlapping. While various partitioning methods are possible [22], Becke's method of fuzzy Voronoi cells [23] is probably the most well‐known, due to its widespread use in the integration of quantities related to density functional theory (DFT) [24, 25]. Becke partitioning (BP) in ADF has produced integrals required for the GW approximation [26, 27, 28, 29, 30]. In SlaterGPU, BP has been used in RKS calculations to produce Kohn–Sham potentials [31].

Two of the authors recently introduced SlaterGPU [16], an algorithm for accelerating STO integrations on graphics processing units (GPUs). This GPU library uses parallel numerical integration with efficient vector operators based on mixed‐precision arithmetic to keep execution costs low. Equally important to the parallelized code was the use of the resolution‐of‐the‐identity (RI) [32, 33, 34] approximation to reduce the complexity of Coulomb integrals. Coulomb integration with RI can be performed in three dimensions in an STO basis, saving a great amount of computational time compared to six‐dimensional integration. Altogether, the original SlaterGPU implementation provided a practical means to perform electronic structure computations at the Hartree–Fock (HF), Complete Active Space Self‐Consistent Field (CASSCF), and full configuration interaction (FCI) levels for double‐ and triple‐zeta STO basis sets. The atomic Becke/Voronoi grids of SlaterGPU, however, are insufficiently precise to treat larger basis sets (e.g., for pentuple‐zeta, polarized basis sets).

The quality of numerical integration scheme is tied closely to the coordinate system and quadrature that together form the integration grid [23]. Becke showed that an orthogonal curvilinear coordinate system known as prolate spheroidal (PS) coordinates can be effective for integrals involving two centers [35]. PS coordinates have an origin between two foci and are composed of a radial part μ and two angular parts ν and ϕ, where μ is composed of spheroids that encompass these foci, ν is composed of hyperboloids, and ϕ is a plane on the focal axis. An interesting property of PS coordinates is that the μ,ν degrees of freedom are closely related to r1, r2, the distances to the two nuclei that define the coordinate system. Close to the nucleus, r1 and r2 are quadratic in μ,ν, meaning that Slater functions become Gaussians (exp−ξr→exp−ξaμ2+ν2). This coordinate transformation reduces the steepness of the STOs near the nuclei in the PS coordinate space, without compromising the Slater shape expected by the Hamiltonian.

While PS coordinates are a natural generalization of spherical coordinates from atoms to diatomics, there is no general procedure to handle cases with more than two atoms within this coordinate system. Fortunately, under the RI approximation, only 3‐center integrals are required, so a complete generalization to many‐center integration is not required. A path forward to use PS coordinates in practical STO integrations for molecules is therefore conceivable, as long as a careful choice of grid discretization around the third center is made.

The SlaterGPU library and its extension to PS coordinates are designed to generate all of the terms necessary for a nonrelativistic electronic Hamiltonian expressed in atom‐centered Slater orbitals, for example, for HF, DFT, and post‐self‐consistent‐field (SCF) methods. The goal of this work is to show that polyatomic systems can be handled within the PS coordinate system, with sufficient accuracy and efficiency to allow for the evaluation of larger basis sets. Not only will this allow increased accuracy in post‐SCF correlated methods, it also will be useful in deriving Kohn–Sham orbitals and exchange‐correlation potentials to high accuracy [31, 36].

Theory

2

PS Coordinates

2.1

PS coordinates are an orthogonal curvilinear coordinate system in three dimensions (μ, ν, ϕ). The coordinates are defined with respect to two foci, which in this case will be positions of two basis functions. If the two foci are defined in Cartesian coordinates at the points 0,0,a and 0,0,−a then the relation between Cartesian and PS coordinates is defined by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Any two‐center integral can be performed after rotation and translation of this grid, assuming distance 2a between the centers (integration with more than two atoms is discussed later on). Radial distances from the foci are

$$[eqn]$$ $$[eqn]$$

As pointed out by Becke [35] for small r1 or r2, a Taylor expansion shows that each of these distances is quadratic in μ,ν. Slater functions therefore can be efficiently integrated due to the avoidance of the cusp near the nucleus, while maintaining the correct physical cusp shape. The ν and ϕ coordinates span 0≤ν≤π and 0≤ϕ≤2π. The ν coordinate moves between the two foci, and the ϕ coordinate rotates around the central axis with cylindrical symmetry. The μ coordinate spans 0≤μ<∞, and therefore is mapped onto a finite range in practice (Figure 1).

The prolate spheroidal grid on a cartesian grid with spheroids at constant μ (red), hyperboloids at constant ν (blue), and planes at constant ϕ (green).

The PS grid must be discretized to perform numerical integration, as Slater functions will be evaluated at discrete points. Uniform grids in the ν and ϕ coordinates are plausible, though the μ coordinate deserves more consideration. For example, the μ grid points should be more concentrated near the nuclei, where the most rapid changes in basis functions occur. Therefore for 0<t1<1 we have

$$[eqn]$$

where C1 and the maximum value of t1 together fix the maximum value of μ. The t1 coordinates, when uniformly divided, lead to more μ points near the nuclei. Here, the ν and ϕ grid points are also spaced evenly. Related discretization methods for other coordinate systems can be found in Refs. [16, 35, 37, 38].

Having introduced a method to divide PS coordinate space into discrete volume elements, quadrature within each volume element completes the integration scheme. Due to the way the grid is constructed—where grid discretization places nuclei only at the edges of each volume element—quadrature points will never be evaluated at a nucleus. This can be done using Gauss‐Legendre quadrature for all dimensions, which approximates an integral over a volume element without placing points on boundaries [39]. Gauss–Legendre quadrature is exact for polynomials of 2n−1 order, allowing rapid convergence with size of quadrature grid.

A generic integral involving a pair of Slater functions has the form

$$[eqn]$$

Applying discretization and quadrature to this generic integral yields

$$[eqn]$$

where functions will be evaluated at points xij, O^ is the operator of interest, wxij is the weight function associated with the quadrature, Q is the number of quadrature points within a volume element, and M is the number of volume elements. The weights are

$$[eqn]$$

where Gμij,Gνij,Gϕij are the weights from the Gaussian quadrature. The volume element is determined by the spatial extent of the discretized cells.

$$[eqn]$$

This two‐center integration scheme will be shown below to be highly effective.

Treatment of a Third Center

2.2

Two‐center integration using PS coordinates is efficient since volume elements can be naturally distributed based on the positions of the two centers. Quadrature is also facilitated around these centers due to the grid lines at μ=0 and ν=0 or π. The third center, however, will sit somewhere in an arbitrarily sized volume element, with no particular location relative to the PS grid lines (Figure 2A). To achieve accurate quadrature, the grid lines should be placed to intersect the third nucleus, and the volume elements subdivided. Our grid lines are therefore shifted to accommodate the third center, specifically by moving the nearest volume element borders in μ and ν (the position of ϕ is trivial, as the three centers will be placed within the same plane before generating the grid). After the grid lines are moved, the 8 volume elements surrounding the third center are then further divided (Figure 2B →C). Figure 2 shows a single division surrounding the third center, where 4 cells are divided in half along μ and ν. Including the ϕ degree of freedom, the 8 neighboring cells become 64 cells. This division can be increased to create more cells as necessary to achieve higher precision. The parameter N_SP_ defines how many cells are used to discretize around the third center.

The process of reorienting the grid around the third center, (A) is the initial spacing of grid with the third center being located in one of the volume elements, (B) moves the μ and ν grid lines closest to the third center so that the third center coordinates align with the grid lines, and (C) splits the grid further around the third center for greater accuracy.

Implementation

2.3

The integration grid is set up by enumerating over the discrete volume elements and their weights. First, the angular grid lines are defined for ν and ϕ by dividing π and 2π by Nν and Nϕ, respectively. The radial component that determines μ (cf. Equation 12) is divided in the t1 transform uniformly through

$$[eqn]$$

where Nμ is the number of μ grid separations. Once the initial 2‐center grid is constructed, the third center is considered in the grid. The third center is placed at ϕ=0 for alignment within the ϕ grid line. Next, the μ and ν grid lines closest to the third center are moved to intersect the third center, as shown in Figure 2. Finally, each volume element defined by the above grid lines is divided up via three‐dimensional quadrature to give the full integration grid (Algorithms 1 and 2).Algorithm 1Coordinate grid discretization: Nμ, Nν, and Nϕ, are the number of divisions over a coordinate. μ, ν, and ϕ represent a point in the center of a volume element and dμ, dν, and dϕ are the distance from one edge of the volume element to the other edge of the volume element.

Algorithm 2Quadrature grid initialization: Qμ, Qν, and Qϕ, are the quadrature points of a respective coordinate and wμ, wν, and wϕ are the associated weights.

Computational Details

3

Most results in this work utilize an all‐electron triple‐zeta basis set with polarization functions (denoted TZ), while others utilize quadruple‐ and pentuple‐zeta basis sets (denoted QZ and 5Z, respectively). Basis sets were constructed to be even tempered (with exponents ζ = αβ ^ n ^, n = 0, −1, · · ·, −N) where N depends on the angular momentum and row on the periodic table. Auxiliary functions were generated in a combinatorial fashion from the original basis by adding together the angular momentum ℓ as well as the exponents of all pairs of functions on each atom. For each channel ℓ, the minimum and maximum summed exponents (ζ max and ζ min) were selected to define the range for the auxiliary basis. The auxiliary basis is then generated through an even‐tempered procedure with ζ max and ζ min as its limits [40, 41, 42, 43]. Finally, all m degrees of freedom for each ℓ were enumerated. The basis sets are provided in Table S1.

Figures 3, 4, 5, 6 and Table 1 all analyze integral matrices directly. Figures 3, 4, 5 and Table 1 do so through analysis of the root mean square error (RMSE) which is defined by

$$[eqn]$$

gi are matrix elements in the integral matrices, gref is a matrix produced from an accurate calculation with a large grid, and g is a matrix produced from some less accurate calculation.

Overlap, kinetic energy, and 2‐center Coulomb repulsion RMSE (Ha) analyzed for ClF and OF2 using the PS method in a TZ basis.

Electron–nuclear attraction and 3‐center Coulomb repulsion RMSE (Ha) analyzed for ClF and OF2 using the PS method in a TZ basis. The grid selection and legend are the same as in Figure 3.

Comparison of the RMSE (Ha) of the electron‐nuclear attraction elements for different divisions of the third center within the PS integration grid. A single split (NSP=2) gives significant advantages over no splitting, and NSP>2 gives marginal improvements.

Six 2‐center Coulomb integrals plotted along 16 vector directions [16] as the two centers are separated. The legend indicates unit vectors for all 16 directions tested, all basis functions used have the exponent ζ=1, and all basis functions were chosen to have m=0.

Coulomb integrals take the form

$$[eqn]$$

Classical and nonclassical (exchange) Coulomb integrals are computed under the RI approximation, where 4‐index integrals are determined as follows

$$[eqn]$$

These two 3‐index terms can be expressed as

$$[eqn]$$

and simplified to

$$[eqn]$$

More information on the form of this potential is discussed in the original SlaterGPU paper [16], while alternatives for Coulomb resolution are in principle also possible [44, 45]. The analytical 1‐center Slater potential has been derived in a previous work [46], this allows for direct use of the expression VCP without fitting the density to the Coulomb operator.

All CI computations in Table 3 were run in a neutral state (singlet spin for all, except doublet spin for NO, N_3_, and NO_2_). The geometries for the diatomic systems in the table were obtained from the Computational Chemistry Comparison and Benchmark DataBase (https://cccbdb.nist.gov/). Triatomic geometries were optimized using Q‐Chem version 5.2 [47], geometries for specific species and information about methods used to obtain these geometries can be found in Table S2.

Heat‐bath configuration interaction (HBCI), used in Tables 3 and 4, is a select‐CI approach that returns a close approximation to the full CI limit [48, 49, 50, 51]. The energy thresholds ε1 and ε2 control the extent of recovery of correlation through variational (ε1) and perturbative (ε2) steps. These thresholds control convergence by screening electronic configurations based on estimating how much an external configuration may contribute to the total energy, and adding in components with the highest importance. In Table 3 HBCI was performed using ε1=1×10−4 Ha and ε2=1×10−7 Ha except in cases which possess greater than 15 valence electrons (N_2_O, NO_2_, OF_2_, SO_2_, CS_2_), which used ε1=5×10−4 Ha and ε2=5×10−7 Ha. The iFCI method, used in Figure 7, truncates the search for configurations in the Hilbert space further by defining localized molecular orbitals as base units for correlation. A many‐body expansion combines these units to systematically recover correlation from a reference state (a valence bond, perfect‐pairing wave function), ensuring convergence to full CI as the expansion level, n, is increased [52, 53, 54]. This allows for polynomial scaling of the iFCI method while maintaining a similar accuracy to HBCI. Here, n=3 recovers the majority of the correlation energy and simplifies the computation of the 24 electron in 225 orbital all electron (core + valence) space. See Refs. [53, 55, 56] for further details of this approach. The HBCI solver in the iFCI calculation uses energy thresholds of ε1,doubles=5×10−4 Ha, ε1,singles=2.5×10−4, and ε2=1×10−7 Ha. The frozen core approximation has been used in all calculations, besides in Table 4 and Figure 7, where all electrons were correlated. The CCSD(T) method [57, 58] was also performed in Figure 7, applied to the same 1,3‐propanediyl system, core correlation was included and a UHF reference was utilized.

HONO and LUNO orbitals and occupation numbers for the singlet and triplet states of 1,3‐propanediyl as well as the singlet–triplet gap for the vertical excitation. iFCI uses the STO QZ basis while CCSD(T) uses a GTO basis, cc‐pVQZ.

For comparisons to GTOs for the methylene and 1,3‐propanediyl systems, the cc‐pVXZ family was used, the cc‐pVXZ‐RIFIT auxiliary basis was utilized in calculations on methylene, no auxiliary basis was used for 1,3‐propanediyl. Here X is T, Q, and 5 for polarized triple‐, quadruple‐, and pentuple‐zeta basis sets, respectively. HBCI [51] was used as a representative electronic structure method, which closely approximates the full CI energy. The Nvidia HPC SDK 25.5 compiler suite with OpenACC was used to compile SlaterGPU and HBCI. The calculations using GTOs on the CH_2_ system in Table 4 were run using Q‐Chem version 5.2 [47]. The CCSD(T) calculations used in Figure 7 were run using ORCA 6.0.1 [59, 60, 61, 62].

When discretizing the initial PS grid, the scalar C1 transforms the overall grid size, depending on the spacing between the nuclei. This is set to

$$[eqn]$$

The angular components are left untransformed.

There are five variables (Nμ, Nν, Nϕ, quadrature points (Q), and third‐center split (N SP)) which need to be chosen to specify the grid for PS integration. To demonstrate convergence with respect to grid size, we vary the radial grid size as well as the angular grid size but fix Q=4 and NSP=3 unless otherwise mentioned. Table S3 shows the choice of grid sizes for the radial and angular grids used in Figures 3, 4, 5. For comparisons in Table 1 between the BP integration (with a grid composed of 5810 angular points and 120 radial points per atom) and the similarly sized PS grid (Nμ:26Nν:32Nϕ:14Q:4NSP:3), these results are compared to a large PS grid (Nμ:80Nν:70Nϕ:50Q:4NSP:3). Table 2 uses the same BP grid as Table 1. Figure 6 as well as Tables 2, 3, 4 use the same large PS grid as Table 1. The iFCI calculations on the 1,3‐propanediyl system in Figure 7 use a PS grid with the following parameters Nμ:80Nν:60Nϕ:40Q:4NSP:4.

Results and Discussion

4

As a starting point, two molecules were selected to demonstrate and benchmark the new PS integration method: ClF and OF_2_. The 2‐ and 3‐center integrals needed for the electronic Hamiltonian under the RI approximation were computed using a range of grid sizes. By varying the number of radial and angular discretization points, the convergence and numerical accuracy of PS integration will be discussed and then compared to integrals from the original SlaterGPU method [16]. The Slater basis for the initial tests is of triple‐ζ quality, including s, p, and d angular momentum functions and exponents ranging from 0.75 to 22.0. The auxiliary basis contains up to g functions.

Figures 3 and 4 show that the PS integration technique produces low errors and smooth convergence with respect to radial and angular discretization for all Hamiltonian elements for ClF. The most challenging cases are the electron‐nuclear attraction integrals, due to the singularity at each nucleus. Regardless, the smallest integration grid for 2‐center Coulomb (Figure 3) reaches an RMSE of order 10−7, and increased grid discretization lowers errors to order 10−12 RMSE. These errors are lower than those of the original SlaterGPU method (a detailed comparison is given later on). The remaining integrals—the overlap, kinetic, electron nuclear and 3‐center Coulomb integrals—all show excellent convergence. Obtaining RMSE below 10−10 does not require especially large grids. The largest grids investigated here show RMSEs around 10−13−10−14, which is close to what is possible with double precision arithmetic. These results confirm two‐center integrals are well matched to numerical PS integration [63, 64].

For the three atom OF_2_ system, the Coulomb and the electron‐nuclear attraction integrals are somewhat more challenging for PS integration. This can be seen in comparison to the 2‐center integrals for ClF in Figure 4, revealing decreased accuracy for the 3‐center integrals. Increasing the grid sizes can systematically lower the RMSE, bringing the errors from the smallest to the largest grid from 10−6 to 10−8 for electron‐nuclear attraction, and 10−7 to 10−9 for the 3‐center Coulomb. As will be shown later on, these errors are sufficiently low to allow high‐quality wavefunction simulations to be performed. Since the nonrelativistic electronic Hamiltonian under the RI approximation involves only integrals with up to three centers, these RMSE values are expected to also apply to polyatomic systems. That is, even for a polyatomic system, the grid does not need to be extended beyond three centers.

As discussed in Section 2.2 and motivated by Figure 4, additional discretization in the PS grid around the third center may be helpful for numerical accuracy. The NSP parameter controls this discretization, so the errors in the electron‐nuclear attraction integrals were analyzed for 1≤NSP≤4. Figure 5 shows that NSP=2 yields noticeable improvements, but NSP>2 has little utility, at least for the equilibrium geometry of the OF_2_ molecule. Since the spacing between μ grid lines grows with μ, the volume elements around the third center grow at large distances. Therefore we expect that higher NSP might have more utility for centers which are farther from each other. To test this hypothesis, a nonequilibrium geometry for OF_2_ was created by moving one fluorine atom to ~30 times its equilibrium bond distance. NSP=3 significantly improved accuracy over NSP=2, but NSP=4 provided little additional utility. For a molecule with long distances between atoms, the NSP discretization scheme may therefore be useful up to about NSP=3. Quantities supporting these results are given in Figure S1.

To test the applicability of PS integration on variation in the molecular geometry, a set of paths involving changing nuclear positions were considered. As the geometry changes, the integral values should be smooth and lacking any artifacts from the numerical integration scheme. To do this, two centers are defined: one at the origin and the other displaced by a unit vector in one of the 16 directions. The second center was moved in one of these vector directions until the distance between the centers was 6 Bohr. Figure 6 illustrates the value of a 2‐center Coulomb integral along this path. All integral profiles are smooth, even for difficult cases with high angular momentum, for instance a pair with f and h functions (ℓ=3 and ℓ=5) on the bottom right corner of Figure 6. A larger list of basis combinations can be found in Figure S2 that show related results with similar accuracy.

The analysis so far suggests that the PS integration scheme is able to accurately integrate 2‐ and 3‐center quantities of the electronic Hamiltonian. Now, we compare the PS method with the original BP method of the SlaterGPU code. This comparison was performed with relatively similar grid sizes of ~700,000 points. The RMSE obtained for the electron‐nuclear and 3‐center Coulomb integrals when comparing BP and PS methods is within the same order of magnitude (Table 1). The overlap, kinetic, and 2‐center Coulomb integrals have improvements in the error of up to three orders of magnitude for PS over BP integration. The improved accuracy of the PS method for 2‐center integrals will now be shown to be important for SCF calculations.

The ability to precisely compute integrals with the PS integration method will be a significant advantage to practical electronic structure simulations. Therefore, a variety of polyatomic 2‐ and 3‐atom molecules were examined to analyze the ability of the PS and BP methods to handle larger basis sets, including basis sets with relatively high angular momentum functions. The TZ, QZ, and 5Z basis sets include two d functions (TZ), three d and one f (QZ), and four d and two f (5Z), respectively. As integral accuracy is lower for BP compared to PS integration, there may be noticeable differences in practical electronic structure computations. As a simple test, self‐consistent field (SCF) computations at the Hartree–Fock (HF) level were performed. The PS method allows for SCF convergence of all systems in all basis sets explored in Table 2. The BP method, however, performs well for the TZ basis set but is less reliable for the larger basis sets.

To understand why the larger basis sets are much more challenging, we considered the smallest eigenvalues of the overlap matrix for the SO_2_ system. For the TZ, QZ, and 5Z basis sets, these values are 1.3 × 10^−5^, 9.8 × 10^−8^, and 7.1 × 10^−8^, respectively. The latter two values are close to the limit of the BP integration accuracy, meaning one or more degrees of freedom in the basis may be poorly behaved. The mixed precision used in the evaluation of BP integrals is believed to contribute to the degradation in performance [16]. PS integration, having accuracies of around 10^−12^ RMSE, experiences no difficulty through 5Z basis sets.

Further tests of the utility of the PS integrals were done for the same molecules in Table 2, this time using a correlated post‐HF method. High‐quality wavefunctions were computed using TZ, QZ, and 5Z basis sets at the HBCI level of theory [48, 49, 50, 51]. HBCI provides a close approximation to full CI, recovering all static and dynamic correlation available to the basis. Given that all degrees of orbital freedom are accessed by HBCI, errors in the integrals can cause serious problems in the reliability of the CI procedure. The HBCI calculations in Table 3 illustrate the accuracy of the integrals and reference orbitals generated using the PS method. In particular, from TZ to QZ, additional correlation energies of −4.22 mHa/electron on average were found, and from QZ to 5Z, −1.29 mHa/electron, suggesting the relevance of high‐quality basis sets for converging these wavefunctions. PS integrals, therefore, were instrumental in achieving accurate results from large Slater basis sets at a correlated level.

Having found that large Slater basis sets can be used in correlated wavefunction computations, we test the PS integration method further on two statically correlated systems. The first case is a well‐studied benchmark [65, 66, 67, 68, 69, 70, 71] for electronic structure methods: the singlet‐triplet spin gap of methylene (CH_2_). The CH_2_ transition between ^1^A_1_ and ^3^B_1_ states were obtained via extrapolation of a series of HBCI computations to the FCI limit with the following HBCI convergence parameters ε1=2×10−4,1×10−4,0.5×10−4 Ha and ε2=5×10−7 Ha. Further details are available in Figures S3 and S4. The HBCI‐STO results obtained are close to Diffusion Monte Carlo (DMC) and internally contracted multireference configuration interaction (CMRCI + Q). Both the HBCI‐STO and HBCI‐GTO calculations behave well and come within 0.36–0.72 and 0.15–0.51 kcal/mol of the experimental spin gap (Table 4).

PS integrals via SlaterGPU are applicable to molecular systems with more than three atoms, as the electronic Hamiltonian can be constructed by computing integrals over subsets of three centers at a time. As an example, incremental FCI (iFCI) [53] computations were run on 1,3‐propanediyl to obtain the vertical singlet–triplet gap. iFCI is able to recover impressive amounts of correlation energy at polynomial scaling: all core and valence electrons were correlated within a 225 orbital active space (full CI). When comparing the singlet–triplet gap energies of the iFCI method with CCSD(T), we observe that the predicted energies for the vertical excitation differ by more than 1 kcal/mol. The predicted energy difference can be explained by the use of an unrestricted reference for CCSD(T), as spin contamination can be a problem for open‐shell systems [72, 73], in this case likely causing a cancelation of errors. Based on the natural orbitals obtained from the iFCI calculation, the 1,3‐propanediyl system has two radical centers located on carbons 1 and 3. Simultaneously, the natural orbital occupation numbers in the ground singlet state indicate the system has significant biradicaloid character (Figure 7).

Conclusion

5

Leveraging the advancements made in the original SlaterGPU paper [16], this work reports an efficient numerical integration scheme to allow for high‐accuracy construction of electronic structure matrix elements. The two‐center PS integrals are around three orders of magnitude more accurate than the prior integration scheme (see Table 1), with the 3‐center integrals being of similar accuracy. The increased accuracy of the overlap integrals in particular allows for the use of larger basis sets in practical electronic structure simulations, since the integrals generated via BP are inadequate for the convergence of QZ and 5Z basis SCF simulations (cf. Table 2). Examples of high precision CI computations demonstrate the practical utility of the new scheme (Tables 3 and 4, Figure 7).

The method presented excels in calculations which require an accurate description of the cusp and exponential wavefunction tails. This is notably beneficial in inverse DFT applications where the cusp conditions can greatly impact the exchange‐correlation potential [31, 36]. The PS scheme should allow for less numerical artifacts and improve the output obtained from the inverse Kohn–Sham problem [20, 74, 75, 76, 77, 78, 79, 80].

Funding

This work was supported by the U.S. Department of Energy (DE‐SC0022241).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1: The information contains basis set zeta values for all centers used in this paper, geometries for all systems, grid sizes used to generate Figures 3, 4, 5. Plot of N SP convergence for the nonequilibrium OF_2_ geometry. Plots of all basis combinations for 2‐center Coulomb integrals varying with displacement along 16 vectors. Plots that show extrapolation to FCI for Table 4 in GTO and STO basis sets.