Implementation of Molecular Gradients for Local Hybrid Density

Article pubs.acs.org/JCTC

Implementation of Molecular Gradients for Local Hybrid Density Functionals Using Seminumerical Integration Techniques Sascha Klawohn, Hilke Bahmann,* and Martin Kaupp* Institut für Chemie, Theoretische Chemie/Quantenchemie Sekretariat C7, Technische Universität Berlin, Straße des 17, Juni 135, 10623 Berlin, Germany ABSTRACT: We present the first implementation of the derivative of the local hybrid exchange-correlation energy with respect to the displacement of nuclei in a Gaussian-type atomic basis set. This extends a recent efficient implementation of local hybrid functionals for self-consistent Kohn−Sham and linear-response TDDFT calculations into the TURBOMOLE program package. In contrast to seminumerical schemes for global exact-exchange admixtures and to the related SCF and TDDFT implementations of local hybrid functionals, additional analytical integrals have to be evaluated at each grid point in the case of molecular gradients. The overall efficiency of the present scheme is improved through prescreening with the density matrix (P-junctions), as well as with spherical overlap estimates (S-junctions). Comparative timings for structure optimizations with local vs global hybrid functionals are discussed while gauging the accuracy for S- and P-junctions using varying thresholds. Local hybrids are furthermore assessed for structure optimization and harmonic vibrational frequency calculations (using numerical second derivatives) of a selection of test systems, comparing with experimental data and some widely used density functionals.

1. INTRODUCTION The tremendous progress of Kohn−Sham density functional theory (DFT) in a large variety of research fields during the past three decades has been fueled to no small extent by important progress in developing approximations for the crucial exchangecorrelation energy.1 The introduction of exact (Hartree−Focktype) exchange (EXX) into DFT by Becke has led to the particularly popular class of hybrid functionals. These represent an attempt toward finding a compromise between the compensation of the self-interaction error by the EXX energy and simulation of some (left−right) static correlation by typical semilocal exchange functionals. Local exchange and correlation functionals depend only on the density, while semilocal functionals may also depend on the density gradient or kinetic energy density. Exact exchange on the other hand is considered fully nonlocal since it is directly derived from the underlying pair density. The initially suggested half-and-half functional2 and many of the most popular hybrid functionals, such as the well-established B3LYP3−6 or PBE07,8 functionals, may be called global hybrids, as they introduce a constant EXX admixture. This corresponds to a “one-size-fits-all philosophy”, hoping that an optimal EXX admixture for each case is not too far from the chosen mean. The first extension of the idea of hybrid functionals beyond global hybrids are range-separated hybrid functionals, where the EXX admixture depends on the interelectronic distance. The longrange corrected functionals of that type have become popular, in particular to solve a variety of problems in linear-response timedependent density functional theory (TDDFT) calculations.9−12 © XXXX American Chemical Society

Double hybrid functionals, on the other hand, are global hybrids with an MP2-type correlation contribution.13 Local hybrid functionals (local hybrids)14,15 represent a particularly flexible form of hybrid functionals. Here the EXX admixture depends on position in real-space, governed by a local mixing function (LMF) that may depend on various possible inhomogeneity parameters. The additional flexibility introduced by position-dependent EXX admixture may be used to develop increasingly accurate functionals, where the balance between a minimization of self-interaction errors and simulation of nondynamical correlation is achieved locally, at each point in space. While the local hybrids reported so far are still relatively simple and may be considered just a starting point for further development (e.g., regarding the “gauge problem” of exchangeenergy densities16−18), results for a wide variety of properties are extremely promising already at the present stage. Examples are a simultaneous accurate description of thermochemistry and reaction barriers but also of various electric or magnetic properties.15,17,19−21 While our own efforts so far have centered on t-LMFs,14 based on the ratio τW (with the von-Weizsäcker τW and the Kohn− τ

Sham local kinetic energy densities τ), or on s-LMFs,22 based on the reduced density gradient, a variety of other LMFs, that may include additional inhomogeneity parameters, have been suggested.23−31 Most recently, we have found that local hybrids offer substantial potential for obtaining accurate vertical Received: May 12, 2016

A

DOI: 10.1021/acs.jctc.6b00486 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation excitation energies in TDDFT calculations,32 in particular for Rydberg and core excitations, as well as for triplet valence excitations. While the flexibility of local hybrids makes them promising for further development, until recently efficient implementations into fast computer programs have been lacking. This has been due to the need to compute nonstandard two-electron integrals, which cannot be obtained analytically. Initial implementations used a resolution-of-the-identity (RI) in the atomic orbital basis to resolve these integrals.33,34 This RI required large, uncontracted basis sets, rendering these earlier calculations computationally inefficient. Recently, we have shown that local hybrids can be implemented efficiently for self-consistent field (SCF)35 and TDDFT32,36 calculations, with timings virtually equal to those of global hybrids, when using seminumerical integration techniques for the nonstandard EXX-based integrals. For the general application of a given exchange-correlation functional to chemical questions, the availability of analytical molecular energy gradients, i.e. the derivative of the energy w.r.t. the position of the nuclei, is a necessity. Not only are gradients needed for structure optimization but also as a basis for vibrational frequency calculations. In the latter case analytical gradients save computational effort compared to fully numerical second derivatives, and they may give more accurate results. In the next section the underlying equations for the implementation of analytical gradients of the local hybrid exchange-correlation energy using seminumerical integration schemes are derived. In this context, the employed prescreening procedure using P-junctions, modified compared to the chain-ofspheres (COSX) scheme for global hybrids,37 is explained. We further provide a preliminary assessment of two local hybrid functionals for structure optimizations and vibrational frequency calculations and discuss timings with the new implementation as well as the balance between computation time and numerical error when using varying thresholds for the S- and P-junctions.

∇A E Xlh =

(3)

Since most LMFs that have been proposed so far exhibit ingredients related to meta-GGA functionals (various derivatives of the density or kinetic-energy densities, see above), the contribution to the gradient from the third and fourth integrand in eq 3 is treated by the usual numerical integration on a grid. The first two terms contain, respectively, the EXX energy density weighted with the gradient of the LMF and the gradient of the EXX energy density multiplied with the LMF. Their implementation within the seminumerical scheme is best explained deriving all equations in the atomic-orbital (AO) basis set. Turning to the first integrand of eq 3, we will assume an LMF that depends on the density ρ, its gradient squared, γ = (∇ρ)2, 1 and the kinetic-energy density, τ = 2 ∑i |∇φi|2 . The molecular orbitals are expanded in the AO basis set φi = ∑μCμiχμ, with coefficients Cμi, that are independent of r, and with contracted basis functions χμ. Hence, the derivative is

∫

∫ a(r)εXex(r) dr + ∫ [1 − a(r)]εXsl(r) dr

1 =− 2

∑∫ ij

φi(r)φi(r′)φj(r)φj(r′) |r − r′|

dr = Ca +

∫

⎡ ∑ Dμν⎢ ∂a 2∇A χμ χν ⎢⎣ ∂ρ μν

⎤ ∂a ∂a 4∇ρ∇A (∇χμ χν ) + ∇A ∇χμ ∇χν ⎥εXex dr ⎥⎦ ∂γ ∂τ

(4)

where the density matrix in the AO basis is defined as Dμν = ∑iCμiCνi, and all terms with derivatives of the orbital coefficients are gathered in Ca. The latter can be transformed to derivatives of the overlap matrix,38 which are already considered in the analytical routines of the program and therefore are not discussed here any further.39 In eq 4 the EXX energy density is multiplied with the first derivatives of the LMF. These are usually available through the preceding SCF implementation, and extension of our algorithm is straightforward. The computation of the EXX energy density and its gradients are of primary interest when discussing the efficiency of the gradient implementation in the following. For Gaussian-type orbitals (GTOs), the derivatives of the basis functions w.r.t. the nuclear positions can be replaced by their derivatives w.r.t. the spatial coordinates by changing the sign, i.e. ∇A χ = −∇χ. The computational bottleneck for local hybrid gradients arises from the contribution to the gradient that contains the derivative of the EXX energy weighted by the LMF (second term in eq 3)

(1)

where r is a real-space coordinate, and a(r) is the LMF, a continuous function controlling the contribution of EXX. Putting the gauge problem aside, LMFs are commonly constructed to be bound between 0 and 1. The EXX energy density is defined in terms of molecular orbitals {φi} as εXex (r)

∇A aεXex

+

2. THEORY The exchange-correlation energy EXC of Kohn−Sham theory is commonly divided into exchange and correlation contributions. In local hybrid functionals the exchange part is a combination of the EXX energy density εex X and a (semi)local exchange-energy density εslX E Xlh =

∫ [∇A aεXex + a∇A εXex − ∇A aεXsl + (1 − a)∇A εXsl] dr

∫ a(r)∇A εXex(r) dr = CXex − ∑ Dμκ Dνλ μνκλ

dr′ (2)

The correlation functional is usually some typical semilocal approximation. Its gradient is calculated as usual and will not concern us in the following. The explicit dependence on r will be suppressed from now on, except where needed, and the equations are derived for the closed-shell case only, keeping in mind that the extension to open shells is straightforward for exchange. Differentiating the exchange energy (eq 1) w.r.t. the nuclear positions (denoted by ∇A ) gives

⎡ ⎢ ⎢⎣

∫ a(r)∇A χμ (r)χν (r) ∫

+

∫ a(r)χμ (r)χν (r) ∫

χκ (r′)χλ (r′) |r − r′|

∇A χκ (r′)χλ (r′) |r − r′|

dr′ dr ⎤ dr′ dr⎥ ⎥⎦

(5)

where the prefactor 1 in eq 2 has been compensated with the 2 equivalence of derivatives of χμ and χν, and of χκ and χλ, reducing the initial sum over four terms to two terms. For global hybrids, i.e. for constant a(r), the two integrals on the r.h.s. of eq 5 are B


Article

Journal of Chemical Theory and Computation

routine is called with an azimuthal quantum number of one for both shells (li = 1, lj = 1). This returns three arrays (in x, y, z) consisting of preliminary integrals for every root of the Gaussian quadrature, in this example for the combinations ss, sp, ps, and pp. The number of roots needed depends on the azimuthal quantum numbers, nroots = (li + lj)/2, rounded up. Before taking the sum over all roots the preliminary integrals are multiplied to yield the actual integrals, i.e. the elements of A. For example, the product of pp for x and ss for both y and z yield the result for the (pxpx) integral of the pp shell pair. This can be written (with r just as an index for the roots) as

identical as the integration variables can be interchanged. With local hybrids this simplification is not possible. Due to the nonlinear dependence of the LMF a(r) on the density and related meta-GGA quantities, the integrals in eq 5 are not solvable analytically. Employing numerical integration, the first integration variable is transformed to a set of grid points (r → rg). No special grids are needed, as standard DFT grids available in any quantum chemistry package are sufficient.35,37 For each grid point the inner integrals are calculated analytically, thus defining the matrix elements A κλg =

′g = A κλ

∫ ∫

χκ (r′)χλ (r′) |rg − r′|

dr′

∇A χκ (r′)χλ (r′) |rg − r′|

nroots

(6)

(px px ) =

r

dr′

=

G2ex =

∇A Gu(x , Ax ) = ∇A [(x − Ax )u exp( −α(x − Ax )2 )]

∑ ag wg ∑ Dμκ ∇A χμg A κλg Dνλχνg g

μνκλ

= 2αGu + 1 − uGu − 1

(8)

μνκλ

(11)

where α is the exponential factor of the Gaussian function, x is a spatial coordinate, and Ax is the position of the nucleus on which the basis function is centered and whose displacement is considered. That is, for gradients the number of integrals per pair doubles. Additionally, the azimuthal quantum number of each shell is artificially increased by one, thus increasing the computational demand. This routine is the most time-consuming step of the calculation. To ease the computational burden, we employ a scheme that reduces the number of calls to the integral routine for the matrix elements. This is achieved by building the elements of A and A′ simultaneously using the Rys scheme. While for the symmetrical A only half of the off-diagonal elements need to be be calculated (see above), there are two paths to choose from for the asymmetrical A′: a) loop over all shell pairs, invoking the underlying integral routine for the elements Aκλg ′ that constitute a given shell pair ij. That is, in one call only integrals that contain derivatives of basis functions belonging to the i-shell are calculated; b) loop only over shell pairs ij with j ≤ i (as for the calculation of A) and calculate off-diagonal elements Aκλg ′ and Aλκg ′ together. The integral routine is invoked only once for all integrals that contain derivatives of basis functions belonging to both the i- and j-shell. In both cases, Aκλg can be acquired as a byproduct. Option b), which we chose, increases the virtual azimuthal quantum number of a given shell pair by one as compared to option a), thereby raising the computational demand for the integral routine. This is, however, overcompensated because option b) results in only half the calls to the integral routine as compared to option a). Revisiting the above example for a pp shell pair, we call the routine as if it was a dd shell pair (li = 2, lj = 2). This gives rise to the preliminary integrals ss, sp, sd, ps, pp, pd, ds, dp, and dd. The elements of A can be calculated from these as described

∑ ag wg ∑ Dμκ χμg A κλ′ g Dνλχνg g

(10)

Accordingly, ss, ps, and sp are multiplied to give the (py pz) integral and so on. Thus, all integral elements for the given shell pair can be calculated, and A is constructed by looping over all shell pairs. As the matrix is symmetrical, only the diagonal and lower-triangle elements of A are calculated, and the latter are multiplied by two to include the upper-triangle values. Turning to the computation of A′, we remember that the derivative of a one-dimensional GTO primitive Gu with a Cartesian quantum number u located at a nucleus A spawns two primitives, one with a lower and one with a higher quantum number

(7)

Note that, while A is symmetrical, A′ is not because the derivative only affects one of the basis functions. Also, each element A′κλg is in fact a three-dimensional vector. As compared to the energy gradient of a global hybrid functional, we calculate therefore more analytical integrals on each grid point. Consequently only the prefactor for the computation time increases, while the overall scaling of the seminumerical scheme is not affected. Inserting the definitions of eq 6 and eq 7 into the r.h.s. of eq 5 and applying the aforementioned numerical integration, we arrive at two terms corresponding to the two double integrals G1ex

∑ (pp)rx ·(ss)ry ·(ss)rz

(9)

Here, the weights wg replace the integration variable dr. For efficient memory management the numerical integration is subdivided into batches of about 100 grid points which are evaluated completely before the next batch is loaded. The elements of A and A′ from eq 6 and eq 7 are calculated using the Rys scheme.40 For example, as illustrated in Figure 1, to get only the A elements of two nonidentical p shells i and j, the

Figure 1. Calculating elements of A arising from a pair of two different p shells. Left: preliminary integrals in x, y, and z for multiple roots (only two are shown); top right: two example calculations; bottom right: overview of all results for the pp example. The preliminary integrals are multiplied in specific patterns to give the actual integrals, e.g. pp, ss, and ss yield the (pxpx) integral (yellow), ss, ps, and sp give (pypz) (blue), and so on. The products have to be calculated for each root and then summed up over those roots (see eq 10). C


Article

Journal of Chemical Theory and Computation above. For the A′ elements we first collect intermediate integrals (center matrices in Figure 2), which are then combined to form

Figure 2. Calculating elements of A′ arising from a pair of two different p shells. Left: preliminary integrals in x, y and z (only x is shown) for multiple roots (only one is shown); center: intermediate integrals; top right: example calculations; bottom right: overview of all results for the pp example (∂ summarizes ∂x, ∂y, and ∂z). The intermediate integrals sp, pp, and dp are crafted from the preliminary integrals as in Figure 1. The intermediate integrals of pp (gray) represent the elements of A. Those from sp and dp are combined according to eq 11 to yield the gradient of pp in all three directions, including the exponential factor α and the Cartesian quantum number as a prefactor, here 1. For example, each of (dxxpz), (dxy pz), and (dxz pz) is combined with (spz) to give the gradients (∂xpxpz), (∂y pxpz), and (∂z pxpz), which can be denoted as a vector (∂pxpz). This example only depicts the calculation of elements A′κλg, but ′ elements (e.g., (px∂pz)) with the remaining preliminary integrals the Aλκg can also be created.

Figure 3. Algorithm to calculate elements of A and A′ for a shell pair ij (without consideration of P-junctions, see Figure 4).

We turn now to the P-junctions. As shown in eq 8 and eq 9, both A and A′ are multiplied with the basis functions, the density matrix, and the weight at each grid point. We calculate two intermediates that combine these values Fκg =

wg

∑ Dμκ χκg μ

Fκ′g =

wg

∑ Dμκ ∇A χκg μ

(12)

(13)

With these the EXX contributions to the gradient in eq 8 and eq 9 become

the actual A′ elements. For example, to compute the derivatives of the first shell in a pp pair, we construct a dp pseudopair for the higher and an sp pseudopair for the lower quantum number. Then we subtract specific l-decreased elements from the lincreased ones, multiplying with the exponential factors and quantum numbers according to eq 11, to arrive at the lowertriangle A′ element. We repeat this process in a “mirrored” way with pd and ps to get the upper-triangle elements respresenting the derivatives of the second shell in the pp pair from the same routine call. The algorithm is roughly depicted in Figure 3. For both A and A′, the diagonal elements are treated in a special way, taking advantage of their identical basis-function parameters and the lack of a “mirrored” element, to save computation time. Prescreening with S- and P-Junctions. To speed up the gradient calculations, we use a similar prescreening as in the SCF implementation of local hybrids,35 probing the overlap of basisfunction shells (S-junctions) and a term related to the density matrix (P-junctions) to skip the evaluation of certain A and A′ elements. The S-junctions follow the COSX algorithm37 for exact exchange in global hybrids. Starting at an initial distance from the nucleus, estimated from the most diffuse primitive, we probe each basis function by increasing the distance and evaluating its value (not its derivative) until it is smaller than a given threshold. Constructing spheres with these distances as radii, only shell pairs whose spheres overlap are evaluated. Accordingly, a lower threshold results in more shell pairs to be evaluated. The threshold is adjusted in negative powers of ten, e.g. 10−5.

G1ex =

∑ ag ∑ Fκ′gA κλg Fλg g

G2ex =

κλ

(14)

∑ ag ∑ FκgA κλ′ g Fλg g

κλ

(15)

For each shell, the values of F or F′ are compared to a given threshold, and if any elements are above that threshold for the current batch of grid points, the shell is marked as mandatory. This results in two lists of junctions that indicate non-negligible elements of A and A′ upon multiplication with F or F′. The Pjunction threshold is distinct but analogous to that for the Sjunctions. For each shell pair to be calculated, the combination of primary (i) and secondary (j) shell is evaluated to see if the pair can be skipped, or at least if the virtual azimuthal quantum number of the integral routine call can be lowered for this pair. The procedure is illustrated in Figure 4.

3. COMPUTATIONAL DETAILS We implemented our method into a development version of the TURBOMOLE program package.41,42 Unless noted otherwise, def2-TZVP43,44 basis sets were employed throughout this work. To examine the effect of basis-set size on computation time, def2SVP44,45 and def2-QZVP46 basis sets were evaluated additionally. All calculations were done with the RI approximation for the Coulomb energy39 using default auxiliary basis sets.47 Out of a larger available number of LMFs, for the purpose of evaluating the efficient implementation of the molecular D


Article


NumForce module, to obtain the second derivatives. For consistency, this was also done for the reference calculations with other functionals. Additional calculations with fully analytical second derivatives, computed using the aoforce module within the RI approximation,59,61−64 were performed to gauge the accuracy of numerical differentiation. The impact on mean errors is, however, marginal (below 1 cm−1 for any GGA and global hybrid functional), and only results with numerical derivatives will be compared below. S- and P-junctions were not used for these calculations. To evaluate computational efficiency aspects, we measured the application of S- and P-junctions by timing gradient calculations for unoptimized linear alkanes (CnH2n+2 with n ∈ {1, ..., 20}). All timings were done using a single CPU core (Intel i3-4130 CPU @ 3.40 GHz). We used the timing mechanism available in the TURBOMOLE programs. The initial structures were created with C−C distances of 145.0 pm, H−C distances of 108.9 pm, and angles of 109.471 degrees. A single SCF was run on each structure. Afterward the gradient was calculated with thresholds for S-junctions and P-junctions varying from 10−4 to 10−8 or without any prescreening. Grid size 1 and def2-TZVP basis sets were used for these calculations. Subsequently, the timing measurements were extended to adamantane (C10H16) as a more compact case. The initial structure parameters are the same as for the linear alkanes above. Here we also investigated the time for a complete structure optimization to an energy threshold of 10−6 au. For the local hybrid functional we tested different P-junction thresholds between 10−4 and 10−6 and also distinct ones for SCF and gradient calculations. In one case both S- and P-junction thresholds were set to 10−5. For comparison we also measured the CPU time for gradients of the seminumerical exact-exchange senex65 algorithm for global hybrids in TURBOMOLE, with grid size 1 and its default S- and P-junction thresholds (which are not directly comparable to the present ones). To estimate the influence of basis-set size we compared def2-SVP, def2-TZVP, and def2-QZVP basis sets, using the same computer as above and grid size 1. For all calculations with def2-QZVP basis sets, a grid point batch size of 70 was used, otherwise it was 100.

Figure 4. Scheme for P-junctions. F and F′ for primary (i) and secondary (j) shells are compared to the threshold to possibly skip all for a primary shell (skip i), all for a secondary shell (skip j), or to lower the virtual azimuthal quantum number for the integral routine by one (li + 1 or lj + 1 instead li + 1,lj + 1). Horizontal arrows (orange) denote that all values of a shell are below the threshold, vertical ones (blue) that at least one is not.

gradient, two simple one-parameter LMFs are used as representative examples: the t-LMF15 (see section 1) defined as a(r) = c

τW |∇ρ|2 =c 8ρτ τ

(16)

22

and the s-LMF given by ⎞ ⎛ |∇ρ| ⎟ a(r) = erf (ds) = erf ⎜d 2 1/3 4/3 ⎝ 2(3π ) ρ ⎠

(17)

To be consistent with two local hybrids previously optimized for thermochemistry and kinetics,15,22 we use Slater exchange,48,49 VWN correlation50 and set the constant prefactors c = 0.48 and d = 0.22. Henceforth we will refer to these specific local hybrids as “t-lh” and “s-lh”. For comparison we used the following functionals: BP8650−52 and PBE53,54 as examples for generalized-gradient approximations (GGA), as well as TPSSh,53,55,56 B3LYP,3−6 PBE0,7,8 and BHLYP2,4,51 as examples for global hybrids. This is the first implementation of molecular gradients for local hybrid functionals; therefore, we also present a preliminary assessment of local hybrids (with the above-mentioned functionals as representatives) for structure optimizations using a set of small molecules of main group elements by Zhao and Truhlar57 (MGBL19 test set) and for a 3d transition-metal test set by Bühl and Kabrede.58 A large grid size59 of 5 was applied for this purpose. As an even more critical test, we computed vibrational frequencies for a set of small molecules (the F2 subset by Scott and Radom60). The structures were optimized, and the frequencies were calculated with a grid size m5 (i.e., a medium grid size 3 during the SCF but a large grid size 5 for the last iteration and the gradient). Furthermore, the SCF convergence criterium was set to 10−9, and the gradient threshold was set to 10−5 during the structure optimization. S- and P-junctions were not used for these calculations. In some cases, frequencies from different irreducible representations are very close, and the order may thus differ from one functional to another. We have therefore compared the calculated to experimental frequencies in numerical order without attempting to match representations. This avoids favoring a given method that is used for the initial assignment. Since analytical second derivatives so far are not available for local hybrid functionals, we used the numerical differentiation of analytical gradients, that is provided by TURBOMOLE’s

4. RESULTS As we evaluate only two simple local hybrid functionals in the present work, the aim of the following section is not a general survey of the performance of this class of functionals for structures and vibrational frequencies. We focus on computational efficiency and numerical aspects. 4.1. Main-Group and Transition-Metal Structure Test Sets. Figure 5 shows mean and mean absolute errors (ME, MAE) of computed bond lengths with the local hybrids and with some other functionals (see section 3), compared to the experimental values of the MGBL19 set. The results of t-lh are comparable to the B3LYP ones, with MAEs of 0.58 and 0.59 pm and maximum errors (MAX) of −2.68 pm (F2) and −2.65 pm (Cl2), respectively. The s-lh results are slightly worse with MAE 0.64 pm and MAX 3.05 pm (F2). TPSSh performs somewhat better, and PBE0 performs slightly worse. BHLYP, chosen as its large exact-exchange admixture of 50% is close to the maximum of 48% in the selected t-LMF, exhibits the largest errors (with generally negative ME), whereas the GGA functionals perform moderately well. Figure 6 shows the results for the set of 3d transition-metal complexes. Here TPSSh has the lowest MAE. The MAEs of t-lh and B3LYP are similar (1.73 pm vs 1.68 pm). The GGA E


Article


hybrids. Only BHLYP requires notably more scaling, while after scaling performance is comparable to the other functionals. These preliminary results suggest that both local hybrids perform similarly for main-group vibrational frequencies as established global hybrids. The IR intensities (not shown) of t-lh and B3LYP are also similar. Studies of vibrational spectra simulations with a larger range of local hybrids and for more difficult cases, needed for a systematic evaluation of performance, are currently done in our lab. 4.3. Timings for Linear Alkanes with S- and PJunctions. Figure 7 and Figure 8 provide timings for Figure 5. Mean error (ME) and mean absolute error (MAE) for bond lengths (in pm) of main-group structure test set MGBL19, comparing two local hybrids (t-lh, s-lh) and a few other functionals.

Figure 7. Relative CPU time of a local hybrid gradient calculation for nalkanes as a function of chain length, and with different thresholds for Sjunctions (in negative powers of ten), compared to results without Sjunctions (no P-junctions have been used here). The kinks of graphs 5 and 6 for C4H10 are artifacts caused by rounding the timings to seconds for times longer than a minute.

Figure 6. Mean error (ME) and mean absolute error (MAE) for bond lengths (in pm) of the set of 3d transition-metal complexes, comparing two local hybrids (t-lh, s-lh) and a few other functionals.

functionals perform well for this test set, as had been noted before,58 whereas PBE0 and in particular BHLYP are somewhat inferior. The MAE of s-lh (1.93 pm) lies between t-lh and PBE0. In summary, the selected local hybrids have a similar accuracy for molecular structures as other commonly used functionals like B3LYP for the chosen test sets, while they have been shown to be more accurate for a larger range of properties.15,17,19,66 4.2. Main-Group Vibrational Frequencies. Let us turn to the typically more critical case of vibrational frequencies, using the F2 set. As is commonly done, for each functional we determined a scaling factor λ = ∑νthνexpt/∑ν2th that minimizes the root-mean-square deviation between calculated and experimental frequencies. The scaling factor compensates for a general overestimation in calculated frequencies, which is only in part caused by a given functional and to a larger extent by the harmonic approximation.60 Table 1 lists the scaling factors and the statistics without and with scaling. Overall, scaling factors and errors before and after scaling are very similar for t-lh, s-lh, and most of the global

Figure 8. Relative CPU time of a local hybrid gradient calculation for nalkanes as a function of chain length, and with different thresholds for Pjunctions (in negative powers of ten), compared to results without Pjunctions (no S-junctions have been used here). The kink of graph 6 for C4H10 is an artifact caused by rounding the timings to seconds for times longer than a minute.

Table 1. Scaling Factor λ and Errors (in cm−1) for the F2 Vibrational Frequency Test Seta

a

functional

λ

ME

MAE

MAX

SME

SMAE

SMAX

BHLYP PBE0 t-lh s-lh B3LYP TPSSh

0.934339 0.960138 0.960182 0.967318 0.967821 0.968549

111.3 60.0 60.9 48.0 47.5 45.4

113.2 63.5 64.6 52.4 52.1 51.6

319.8 216.6 215.8 206.5 210.4 210.7

−0.8 −6.0 −5.1 −5.7 −5.4 −6.2

20.4 24.9 22.1 23.1 20.3 20.6

179.3 169.5 170.6 169.7 174.1 175.2

ME: mean error; MAE: mean absolute error; MAX: absolute maximum error; SME, SMAE, SMAX: the corresponding scaled results. F


Article


th

S

P

for functionals of the same family are almost the same). The timings are given relative to those of the corresponding calculation with the global hybrid (TPSSh) using the efficient analytical gradient of TURBOMOLE’s rdgrad module. We additionaly provide data for TPSSh obtained with the senex module, which also uses a seminumerical treatment of the EXX integrals (see section 3). As expected, the GGA gradient calculation is much faster than that with the global hybrid, which in turn is faster than the current implementation for local hybrids. However, the seminumerical implementation for local and, as pointed out earlier, for global hybrids scales better with basis-set size than the analytical implementation for global hybrids. Thus, while the local hybrid gradient takes 7.5 times longer than a standard TPSSh calculation with the small def2-SVP basis sets, the factor decreases to 3.3 with def2-TZVP and to 2.6 with def2-QZVP. We also confirm that the effect of using P-junctions (shaded area on the bar) becomes more notable with increasing basis-set size. The senex algorithm for global hybrids (provided with standard settings for S- and P-junctions) performs even better, with factors of 1.5 (def2-SVP), 0.64 (def2-TZVP), and 0.62 (def2-QZVP). Due to the additional integrals needed for A′ with local hybrids (see section 2), the corresponding factor decreases less quickly for this system size. Table 3 lists the computation times and cycles of a complete structure optimization of adamantane using t-lh with varying P-

10−4 10−5 10−6 10−7 10−8

6 ·10−05 9 ·10−09 2 ·10−09 4 ·10−12 2 ·10−14

5 ·10−05 4 ·10−07 1 ·10−09 2 ·10−12 1 ·10−15

Table 3. Absolute CPU Time for a Full Structure Optimization (Excluding the Initial SCF) of Adamantane with a Local Hybrid Using def2-QZVP Basis Set and Different PJunction Thresholds (no S-Junctions)a

computation of the local hybrid gradient as a function of alkane chain length and with different thresholds for S- and P-junctions, respectively. While the overall appearance of the two graphs is similar, the magnitude of the time savings due to prescreening by S- and P-junctions is notably different. In both cases, the percentage saving increases with chain length and thus with system size. However, S-junctions (Figure 7) are less efficient for prescreening in this case than P-junctions (Figure 8). Taking reasonably conservative and accurate (see below) thresholds of 10−5 for both cases, S-junction savings converge to about 7% for longer chains, whereas the reduction in computation time due to P-junction does not seem to level off much even at 20 carbon atoms, where it amounts already to almost 40%. With tighter thresholds, the savings are less, and they start at larger chain lengths. Turning to the effects of S- and P-junctions on numerical accuracy, we note that the errors of the gradients with S- and Pjunctions relative to calculations without prescreening remain approximately constant with chain length. Table 2 provides Table 2. Mean Absolute Errors (MAE) of Local Hybrid Alkane Molecular Gradients (Averaged over the Chain Length) for Different Thresholds (th) of S- and P-Junctionsa MAE (∇AE)/au

a

The reference values are gradients without prescreening.

P threshold SCF

grad  10−4

MAEs for all alkanes studied. These depend appreciably on the thresholds used. Considering an accuracy of 10−6 au for the gradient as reasonable for most purposes, we see that thresholds of 10−5 for both S- and P-junctions provide sufficient accuracy but still allow for favorable timings (see above). If we want to be even more conservative, thresholds of 10−6 may be used. 4.4. Timings (with S- and P-Junctions) for Adamantane. Figure 9 provides timings for a single gradient calculation of adamantane. In addition to a local hybrid with the present implementation (t-lh), we have chosen TPSSh as an example of a global hybrid and PBE as an example GGA functional (timings

−4

−5

10 10−5

10 10−4 10−5

−5

10 10−6

10−6 10−5

10−6 S,P th. = 10−5

cycles SCF

struc

t/h

ΔE, kJ/mol

50 116 84 72 50 50 47 50 50

10 15 11 17 10 10 9 10 10

13.7 13.8 12.3 12.6 10.4 12.3 10.1 12.8 9.8

10+0 10+0 10−1 10−2 10−2 10−3 10−4 10−2

t/t0 +1% −10% −8% −24% −10% −26% −6% −28%

a

The energy difference and timing ratio refer to the optimization without P-junctions. For comparison, timings with both the recommended value of S- and P-junctions (10−5) are given as well.

junction thresholds (no S-junctions), plus one result with both Sand P-junction thresholds set to 10−5. While the computation time of a single SCF cycle and gradient calculation decreases with looser thresholds, for thresholds of 10−4 the overall time for a structure optimization increases due to inaccuracies in the intermediate gradients. The total-energy error after structure convergence decreases to 0.01 kJ/mol upon using P-junction thresholds of 10−5 or lower, and the computation time decreases by 24%. If S-junction thresholds are additionally set to 10−5, the error remains the same, but the time is lowered further by about 4%, resulting in an absolute CPU time of 9.8 h compared to 13.7 h without prescreening. This also confirms previous findings that P-junctions are more important than S-junctions. Additional variations with tighter thresholds than 10−5 for either the SCF or the gradient are likely insignificant within our measurement accuracy for a full optimization process. Finally our value of 10−5

Figure 9. Relative CPU times for a gradient calculation of adamantane with def2-SVP, def2-TZVP, and def2-QZVP basis sets, respectively. The time of TPSSh is set to 1 with each basis set. The shaded area of t-lh depicts the time savings obtained with S- and P-junction thresholds of 10−5. sxTPSSh stands for a TPSSh calculation using the senex module with standard prescreenings. G


Article

Journal of Chemical Theory and Computation Notes

for both S- and P-junctions suggested above is confirmed as a good compromise between accuracy and efficiency for most applications.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Toni M. Maier for helpful discussions and the TURBOMOLE developers consortium for access to the development version of the program package.

5. CONCLUSIONS This work provides the first implementation of molecular gradients for local hybrid functionals in a Gaussian-type atomic basis set, in a development version of the TURBOMOLE program package. The gradients are derived representatively with meta-GGA ingredients as needed for a t-LMF or s-LMF. Extension to other LMFs is, however, straightforward on the basis of an SCF implementation. This opens the way for structure optimization and vibrational frequency calculation (by numerical differentiation of analytical gradients) for this promising class of density functionals. As a first assessment, we have applied two simple one-parameter local hybrid functionals to two test sets for molecular structures (including both main-group and 3d elements) and to a vibrational-frequency test set of small maingroup systems. Both local hybrids have been optimized previously for thermochemistry, and they perform comparably well as some standard global hybrids. More detailed evaluations of a larger variety of local hybrid functionals for more challenging cases of structures and vibrational spectra should reveal areas where local hybrids provide distinct advantages over more common functionals, as they have been found to do for other ground-state properties and for electronic excitation energies. As for our previous SCF and TDDFT implementations of local hybrids, we have used seminumerical techniques to solve the nonstandard two-electron integrals arising with local hybrids. While our recent SCF implementation exhibits virtually the same timings for local and global hybrid functionals (using the seminumerical implementation of exact exchange), we have shown that the energy gradient requires additional two-center electron repulsion integrals that may be omitted for global hybrids. Within the seminumerical scheme the computational effort of a gradient calcuation with local hybrids is thus always larger than with a global hybrid although this effect is expected to diminish for large systems and basis sets. For global hybrid functionals the scaling of seminumerical exact exchange is favorable over fully analytical integration schemes, which has been shown numerically as well. We note that this pertains also to our gradient implementation for local hybrid functionals. To nevertheless obtain an efficient implementation also for small and medium-sized systems, we have made extensive use of prescreenings with S- and P-junctions, as introduced in Neese’s COSX algorithm for global hybrids. The effect of such prescreenings on computation time and numerical accuracy has been evaluated for a series of linear alkanes and for adamantane. A threshold of 10−5 for both S- and P-junctions provides good numerical accuracy but still allows a significant reduction of computation time (by about 30% for a complete structure optimization of adamantane with def2-QZVP basis sets). Since the seminumerical treatment becomes more competitive to analytical integration schemes for global hybrids, and the prescreening more effective with increasing system size and basis sets, structure optimizations for large molecules with local hybrid functionals have become possible.

■

■

REFERENCES

(1) Becke, A. D. J. Chem. Phys. 2014, 140, 18A301. (2) Becke, A. D. J. Chem. Phys. 1993, 98, 1372−1377. (3) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (4) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (5) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200−206. (6) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (7) Burke, K.; Ernzerhof, M.; Perdew, J. P. Chem. Phys. Lett. 1997, 265, 115−120. (8) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6170. (9) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. J. Chem. Phys. 2004, 120, 8425−8433. (10) Jacquemin, D.; Perpéte, E. A.; Scuseria, G. E.; Ciofini, I.; Adamo, C. J. Chem. Theory Comput. 2008, 4, 123−135. (11) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. J. Chem. Phys. 2008, 128, 044118. (12) Isegawa, M.; Peverati, R.; Truhlar, D. G. J. Chem. Phys. 2012, 137, 244104. (13) Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2011, 7, 291− 309. (14) Jaramillo, J.; Scuseria, G. E.; Ernzerhof, M. J. Chem. Phys. 2003, 118, 1068−1073. (15) Bahmann, H.; Rodenberg, A.; Arbuznikov, A. V.; Kaupp, M. J. Chem. Phys. 2007, 126, 011103. (16) Arbuznikov, A. V.; Kaupp, M. J. Chem. Phys. 2014, 141, 204101. (17) Theilacker, K.; Arbuznikov, A. V.; Kaupp, M. Mol. Phys. 2016, 114, 1118−1127. (18) Maier, T. M.; Haasler, M.; Arbuznikov, A. V.; Kaupp, M. Phys. Chem. Chem. Phys. 2016, DOI: 10.1039/C6CP00990E. (19) Kaupp, M.; Bahmann, H.; Arbuznikov, A. V. J. Chem. Phys. 2007, 127, 194102. (20) Arbuznikov, A. V.; Kaupp, M. Int. J. Quantum Chem. 2011, 111, 2625−2638. (21) Theilacker, K.; Arbuznikov, A. V.; Bahmann, H.; Kaupp, M. J. Phys. Chem. A 2011, 115, 8990−8996. (22) Arbuznikov, A. V.; Kaupp, M. Chem. Phys. Lett. 2007, 440, 160− 168. (23) Johnson, E. R. J. Chem. Phys. 2014, 141, 124120. (24) de Silva, P.; Corminboeuf, C. J. Chem. Phys. 2015, 142, 074112. (25) Schmidt, T.; Kraisler, E.; Makmal, A.; Kronik, L.; Kümmel, S. J. Chem. Phys. 2014, 140, 18A510. (26) Haunschild, R.; Odashima, M. M.; Scuseria, G. E.; Perdew, J. P.; Capelle, K. J. Chem. Phys. 2012, 136, 184102. (27) Perdew, J. P.; Staroverov, V. N.; Tao, J.; Scuseria, G. E. Phys. Rev. A: At., Mol., Opt. Phys. 2008, 78, 052513. (28) Arbuznikov, A. V.; Kaupp, M. J. Chem. Phys. 2008, 128, 214107. (29) Janesko, B. G.; Scuseria, G. E. J. Chem. Phys. 2008, 128, 084111. (30) Arbuznikov, A. V.; Bahmann, H.; Kaupp, M. J. Phys. Chem. A 2009, 113, 11898−11906. (31) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Phys. 2014, 141, 144104. (32) Maier, T. M.; Bahmann, H.; Arbuznikov, A. V.; Kaupp, M. J. Chem. Phys. 2016, 144, 074106. (33) Arbuznikov, A. V.; Kaupp, M.; Bahmann, H. J. Chem. Phys. 2006, 124, 204102. (34) Janesko, B. G.; Krukau, A. V.; Scuseria, G. E. J. Chem. Phys. 2008, 129, 124110.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. H


Article

Journal of Chemical Theory and Computation (35) Bahmann, H.; Kaupp, M. J. Chem. Theory Comput. 2015, 11, 1540−1548. (36) Maier, T. M.; Bahmann, H.; Kaupp, M. J. Chem. Theory Comput. 2015, 11, 4226−4237. (37) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Chem. Phys. 2009, 356, 98−109. (38) Pople, J. A.; Gill, P. M. W.; Johnson, B. G. Chem. Phys. Lett. 1992, 199, 557−560. (39) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 242, 652−660. (40) Dupuis, M.; Rys, J.; King, H. F. J. Chem. Phys. 1976, 65, 111−116. (41) TURBOMOLE V6.6 2014, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007. Available from http://www. turbomole.com (accessed July 1, 2016). (42) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165−169. (43) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152. (44) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297− 3305. (45) Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571− 2577. (46) Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003, 119, 12753−12762. (47) Weigend, F. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (48) Dirac, P. A. M. Proc. R. Soc. London, Ser. A 1929, 123, 714−733. (49) Slater, J. C. Phys. Rev. 1951, 81, 385−390. (50) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (51) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (52) Perdew, J. P. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (53) Perdew, J. P.; Wang, Y. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 13244−13249. (54) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (55) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (56) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. J. Chem. Phys. 2003, 119, 12129−12137. (57) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101. (58) Bühl, M.; Kabrede, H. J. Chem. Theory Comput. 2006, 2, 1282− 1290. (59) Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346−354. (60) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502−16513. (61) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119−124. (62) Deglmann, P.; Furche, F. J. Chem. Phys. 2002, 117, 9535−9538. (63) Deglmann, P.; Furche, F.; Ahlrichs, R. Chem. Phys. Lett. 2002, 362, 511−518. (64) Deglmann, P.; May, K.; Furche, F.; Ahlrichs, R. Chem. Phys. Lett. 2004, 384, 103−107. (65) Plessow, P.; Weigend, F. J. Comput. Chem. 2012, 33, 810−816. (66) Kaupp, M.; Arbuznikov, A. V.; Bahmann, H. Z. Phys. Chem. 2010, 224, 545−567.

I


Implementation of Molecular Gradients for Local Hybrid Density

Recommend Documents