An efficient implementation is presented for analytic gradients of the coupled-cluster singles and doubles (CCSD) method with the density-fitting approximation, denoted DF-CCSD. Frozen core terms are also included. When applied to a set of alkanes, the DF-CCSD analytic gradients are significantly accelerated compared to conventional CCSD for larger molecules. The efficiency of our DF-CCSD algorithm arises from the acceleration of several different terms, which are designated as the “gradient terms”: computation of particle density matrices (PDMs), generalized Fock-matrix (GFM), solution of the Z-vector equation, formation of the relaxed PDMs and GFM, back-transformation of PDMs and GFM to the atomic orbital (AO) basis, and evaluation of gradients in the AO basis. For the largest member of the alkane set (C_{10}H_{22}), the computational times for the gradient terms (with the cc-pVTZ basis set) are 2582.6 (CCSD) and 310.7 (DF-CCSD) min, respectively, a speed up of more than 8-folds. For gradient related terms, the DF approach avoids the usage of four-index electron repulsion integrals. Based on our previous study [U. Bozkaya, J. Chem. Phys. **141**, 124108 (2014)], our formalism completely avoids construction or storage of the 4-index two-particle density matrix (TPDM), using instead 2- and 3-index TPDMs. The DF approach introduces negligible errors for equilibrium bond lengths and harmonic vibrational frequencies.

## I. INTRODUCTION

Tensor decompositions of the four-index electron repulsion integrals (ERIs) have been of great interest in modern computational chemistry.^{1–21} One of the popular approximations is density fitting (DF), which is also known as the resolution of the identity (RI).^{1–10,15–18} In the DF approach the four-index ERIs are expressed in terms of three-index tensors. Another integral factorization method is Cholesky decomposition (CD) of the four-index integral tensor.^{11–16} The DF and CD approximations are very successful to reduce the computational cost of integral transformations and the storage requirements of four-index two-electron integrals (TEIs). In the DF approach, TEIs of the primary basis set are written in terms of a pre-defined auxiliary basis set, while in the CD approximation the three-index tensors are generated directly from the primary basis set.

In the last four decades, coupled-cluster methods^{22–26} have proven to be highly reliable for the computation of molecular properties. Analytic gradient methods greatly facilitate the location of stationary points on potential energy surfaces and enable the computation of one-electron properties.^{27–30} In 1984, Adamowicz *et al.*^{31} reported the initial algebraic expressions for analytic gradients of the coupled-cluster singles and doubles (CCSD) method,^{32} following the Z-vector method of Handy and Schaefer.^{33} In a 1985 study, Fitzgerald *et al.*^{34} implemented analytic gradients for the coupled-cluster doubles (CCD) method. The first implementation of analytic gradients for the closed-shell CCSD method was reported by Scheiner *et al.*^{35} Later, several CCSD analytic gradient papers were also reported.^{36–40}

Analytic gradients for electron correlation methods with the DF approximation have been reported for second-order Møller–Plesset perturbation theory (MP2) with conventional Hartree-Fock (HF) reference (RI-MP2),^{41–43} coupled cluster model CC2,^{44,45} local MP2 (DF-LMP2),^{46} complete active space second-order perturbation theory (DF-CASPT2),^{47} MP2 with DF-HF reference (DF-MP2),^{17} orbital-optimized MP2 (DF-OMP2),^{18} and local time-dependent coupled cluster response theory (DF-TD-LCC2).^{48} However, analytic gradients cannot be computed with the CD approximation unless one employs an approach like that of Aquilante *et al.*^{49}

Energies for the CCSD method with the DF and CD approximations have been reported in previous studies.^{15,50–52} In this research, analytic gradient expressions are presented for the density-fitted CCSD and CCD methods (DF-CCSD and DF-CCD), where the density-fitting has been applied to both reference and correlation energies. The equations presented (including frozen core terms) have been implemented in a new computer code and added to the Dfocc module^{16–18,53,54} of the Psi4 package.^{55} The computational time of the DF-CCSD gradients is compared with that of the conventional CCSD method (from Molpro2012.1).^{56–58} The DF-CCSD and DF-CCD methods are applied to equilibrium geometries and vibrational frequencies.

## II. DF-CCSD VS CCSD FOR ANALYTIC GRADIENTS

Before reporting detailed equations for DF-CCSD analytic gradients, we would like to highlight the differences between DF-CCSD and CCSD for computation of analytic gradients. Theoretically, DF-CCSD analytic gradients can be obtained by replacing the conventional ERIs (including derivative integrals) with DF tensors. At the implementation level, the gradient formalism can be more efficient by using the DF intermediates. Further, we need to take derivatives of the 2-index and 3-index integrals instead of the 4-index ERIs. Moreover, with the DF approach the 4-index two-particle density matrix (TPDM) can be replaced with 2-index and 3-index TPDMs, which greatly simplifies the computation of analytic gradients. Further, with the new definition of the 3-index TPDMs, the cost of the constructing the generalized Fock matrix (GFM) is substantially reduced, from *O*(*N*^{5}) to *O*(*N*^{4}). Moreover, the computational time for the solution of the Z-vector equation^{33} is significantly reduced with the DF intermediates.^{17} Finally, we would like to note that the orbital relaxation terms are computed more efficiently with the DF-HF reference.^{17}

## III. THEORETICAL APPROACH

### A. Density-fitting

In the DF approximation the atomic-orbital (AO) basis TEIs are written as

the DF factor $b\mu \nu Q$ is defined as

where

and

where {*χ*_{μ}(**r**)} and {*φ _{P}*(

**r**)} are members of the primary and auxiliary basis sets, respectively.

The molecular-orbital (MO) basis integrals can be written similarly,

where $bpqQ$ is a MO basis DF factor.

### B. DF-CCSD energy and amplitude equations

For the orbital indexing, the conventional notation is used: *i*, *j*, *m*, *n* for occupied orbitals; *a*, *b*, *e*, *f* for virtual orbitals; and *p*, *q*, *r*, *s*, *t*, *u*, *v*, *w* for general spin orbitals. The DF-CCSD correlation energy can be written as follows:

where $H\u02c6N$ is the normal-ordered Hamiltonian operator,^{25,59} |0〉 is the reference determinant, and $T\u02c6$ is the cluster excitation operator. For the CCSD wave function,

where $T\u02c61$ and $T\u02c62$ are single- and double-excitation operators, respectively.

where $tia$ and $tijab$ are the single and double excitation amplitudes, and $a\u02c6\u2020$ and $i\u02c6$ are the creation and annihilation operators, respectively.

The DF-CCSD correlation energy can be written explicitly as follows:

where 〈*ij*| |*ab*〉_{DF} is the antisymmetrized TEI within the DF approximation and *f _{pq}* is the Fock matrix.

The amplitude equations can be written as

where $\u3008\Phi ia|$ and $\u3008\Phi ijab|$ are singly and doubly excited Slater determinants, respectively.

The explicit form of the *T*_{1} amplitude equation with the DF approximation can be written as

and the *T*_{2} amplitude equation with the DF approximation can be expressed as

Intermediates involved in Eqs. (14) and (15) are reported in Appendix A. Further, in all equations *P*_{±}(*pq*) is defined by

with $P(pq)$ acting to permute the indices *p* and *q*. *δ _{pq}* denotes Kronecker delta and the denominator tensors are defined by

### C. DF-CCSD-Λ energy functional (unrelaxed Lagrangian)

In order to obtain a variational energy functional ($L\u02dc$) it is convenient^{60,61} to introduce a Lagrangian (DF-CCSD-Λ functional) as follows:

where $H\u02c6$ is the Hamiltonian operator and $\Lambda \u02c6$ is the coupled-cluster de-excitation operator. For the CCSD wave function,

where $\Lambda \u02c61$ and $\Lambda \u02c62$ are the single and double de-excitation operators, respectively.

where $\lambda ai$ and $\lambda abij$ are the single and double de-excitation amplitudes, respectively.

The CCSD *T*-amplitude equations are obtained by requiring that $L\u02dc$ be stationary with respect to Λ-amplitudes, while the stationary requirement with respect to *T*-amplitudes leads to the Λ-amplitude equations^{36,38,39,62,63}

where *E* is the CCSD total energy.

Starting from the formulation of Gauss and Stanton^{62} for the Λ-amplitudes with the conventional integrals, we obtain our equations for the DF-CCSD wave function. We would like to note that our expressions are somewhat different than those of Gauss and Stanton. To avoid the formation or the storage of 3-virtual index intermediates we significantly reformulate the Λ-amplitude equations. In our formulation the $Wabei$ and $Wamef$ intermediates of Gauss and Stanton^{62} are not present. The explicit form of the Λ_{1}-amplitude equation with the DF approximation is given by

and the explicit form of Λ_{2}-amplitude equation with the DF approximation is given by

Intermediates involved in Eqs. (25) and (26) are reported in Appendix B.

### D. DF-CCSD analytic gradients

Since DF-CCSD energy is not stationary with respect to orbital rotation parameters and coupled-cluster amplitudes, one needs to consider their derivatives in the gradient equation. For derivatives of the orbital rotation parameters, one should have solved the perturbation dependent coupled-perturbed Hartree-Fock (CPHF) equations. However, as it is first shown by Handy and Schaefer for configuration interaction (CI) gradients,^{33} one may solve a perturbation independent equation (Z-vector equation) instead. Later, a more systematic approach, called as the Lagrangian approach, was developed by Helgaker and Jørgensen for analytic gradient methods.^{60,61,64,65} With the Lagrangian approach, it is also possible to replace the perturbation dependent amplitude derivative equations with the *λ*-amplitude equations, as in the case of orbital rotation parameters. The details of the Lagrangian approach are out of scope for our paper. Readers who are interested may refer to the references provided. Hence, for the DF-CCSD wave function, the following Lagrangian can be employed:

where the operator $Z\u02c6$ is defined as

and let us recall that the HF MO gradient is

Differentiating the Lagrangian in Eq. (27) with respect to orbital rotation parameters we can obtain the following linear equation, which defines the Z-vector:^{33}

where **A** is the HF MO Hessian and **w** is the MO gradient of the post-HF method considered. The details of the orbital response equation (30) are given in our recent paper.^{17}

The gradient of energy can be written as follows:^{60,61,64–70}

when we insert the explicit expressions for integral derivatives into Eq. (31), the first derivative equation can be cast into the following form:^{17}

where $\gamma pqeff$, $\Gamma PQeff$, $\Gamma \u02dcpqQ(eff)$, and $Fpqeff$ are the relaxed one-particle density matrix (OPDM), 2- and 3-index two-particle density matrix (TPDM), and GFM. The unrelaxed OPDM,^{71,72}

the unrelaxed GFM,^{16–18}

and, unrelaxed 2- and 3-index TPDMs are defined by^{17}

The usage of 3-index TPDMs was also considered in other previous studies in the context of DF-MP2, CD-HF, and CD-MP2 gradients.^{41–43,73,74} However, in these studies the employed 3-index TPDM is equivalent to our $\Gamma \u02dcpqQ$. In a 2014 study,^{17} it has been demonstrated that with the definition of $\Gamma pqQ$, 3-index TPDM can be readily evaluated from the conventional 4-index TPDM and the [**J**^{−1/2}]_{PQ} factor can be excluded until the final stage of the gradient computations. Therefore, it is more convenient for higher-order methods, such as DF-CCSD.

The gradient is evaluated by back-transforming the relaxed PDMs and the relaxed GFM into the AO basis and contracting against the appropriate AO derivative integrals, as usual^{75–77}

where *F*_{μν}, *γ*_{μν}, and $\Gamma \mu \nu Q$ are the AO basis GFM, OPDM, and 3-index TPDM, respectively. The 2-index TPDM can be written as

hence, the final gradient equation in the AO basis can be written as follows:^{17}

All explicit equations for the GFM and the orbital response contributions to GFM and PDMs, as well as the detailed discussion of frozen-core gradients, were presented in our recent paper.^{17} Hence, we will not repeat those detailed discussions.

### E. Unrelaxed particle density matrices

The OPDM can be decomposed as follows:

where $\gamma pqref$ and $\gamma pqcorr$ are the reference and correlation contributions to OPDM. The reference OPDM is given by

where *δ _{ij}* denotes the Kronecker delta.

Similar to the OPDM case, the TPDM can be decomposed as follows:

where $\Gamma pqQ(ref)$ and $\Gamma pqQ(corr)$ are the reference and correlation contributions to TPDM, and $\Gamma pqQ(sep)$ is the separable part of TPDM. The reference TPDM is given by

where $bijQ$ is evaluated using the reference auxiliary (DF-REF) basis, for example, the JK-FIT basis set.^{6} The unique non-zero blocks of the separable TPDMs can be written as

where $bpqQ$ is again evaluated in the DF-REF basis, and intermediates are given by

We would like to note that for the reference and separable parts of PDMs the integrals of the DF-REF basis (e.g., the JK-FIT basis set) should be employed.

#### 1. Correlation part of unrelaxed particle density matrices

The unique blocks of the correlation OPDM can be written as

Similarly, the unique blocks of the correlation TPDM can be expressed as

Intermediates involved in Eqs. (56)–(61) are reported in Appendix C. Finally, we would like to note that for the correlation parts of PDMs, the integrals of the DF-CC basis (e.g., the cc-pVXZ-RI basis sets) should be employed.

### F. Gradient algorithm steps

The algorithm steps to compute the analytic gradients are illustrated as follows:

Perform the self-consistent field (SCF) iterations to obtain the DF-HF energy and the MO coefficients.

Compute DF-CC integrals in the AO basis and transform them to the MO basis.

Perform DF-CCSD

*t*-amplitude iterations and compute the DF-CCSD energy.Perform DF-CCSD

*λ*-amplitude iterations to obtain the converged*λ*-amplitudes.Compute the unrelaxed PDMs and GFM.

Compute MO gradient and form the effective MO gradient if frozen-core orbitals are present.

^{17}Solve Eq. (30) to obtain the Z-vector (DF-SCF integrals should be used).

^{17}Compute the effective OPDM, GFM, and 3-index TPDM ($\Gamma pqQ(sep)$ should be considered).

^{17}Compute $\Gamma \u02dcpqQ$ from Eq. (36) for $\Gamma pqQ(ref)$, $\Gamma pqQ(sep)$, and $\Gamma pqQ(corr)$ separately. For $\Gamma pqQ(ref)$ and $\Gamma pqQ(sep)$ the DF-SCF basis [

**J**^{−1/2}]_{PQ}should be used, while the DF-CC basis [**J**^{−1/2}]_{PQ}tensor should be used for $\Gamma pqQ(corr)$.Combine $\Gamma pqQ(ref)$ and $\Gamma pqQ(sep)$ to form $\Gamma pqQ(RefSep)$.

Backtransform the effective OPDM, GFM, $\Gamma pqQ(RefSep)$, and $\Gamma pqQ(corr)$ to the AO basis, and compute the 2-index TPDM in the AO basis [Eqs. (39)–(42)].

Compute the gradients in the AO basis using Eq. (44). The one-electron and overlap derivatives are computed once, the first two terms of Eq. (44). However, two-electron and metric derivatives are computed in the DF-SCF and DF-CC basis, separately, using $\Gamma \mu \nu Q(RefSep)/\Gamma PQ(RefSep)$ and $\Gamma \mu \nu Q(corr)/\Gamma PQ(corr)$, respectively.

### G. Specific notes on the implementation

At first, we would like to note that in our implementation we do not store any 3- and 4-virtual four-index tensors either in the core memory or on the disk. For example, we need to build *W _{abef}* intermediate for

*T*

_{2}and Λ

_{2}amplitude equations and the

*V*intermediate for the 3-index TPDM. In such cases we build these tensors on the fly, and add their contributions directly to the corresponding amplitudes or TPDM. In our direct algorithms we assume that there is enough memory to store $Wbefa$ or $Vbefa$ type tensors for a fixed virtual index

_{abef}*a*, requiring dynamic arrays with the size of

*V*

^{3}(where

*V*denotes virtual orbitals). Further, whenever 4-index integrals including two or fewer virtual indices (such as

*OOOO*,

*OOOV*,

*OOVV*, and

*OVOV*type integrals, where

*O*denotes occupied orbitals) are necessary, they are built on the fly from the DF integral tensors and they are stored in the core memory.

Another significant point is the evaluation of the most expensive term of *T*_{2} and Λ_{2} amplitude equations. For example, for the closed-shell case the most expensive term can be written as

hereafter, this term will be called as *particle-particle ladder* (PPL). Following the previous studies of Saebø and Pulay,^{78} and Scuseria, Janssen, and Schaefer^{79} we employ the following algorithm for the evaluation of PPL:

where *S* is the symmetric component, while *A* is the anti-symmetric component. Now let us define

where *S* and *A* have the following symmetry properties:

hence, we can always keep *i* ≥ *j* and *a* ≥ *b*.

The pseudo code for our PPL algorithm is given as

Further, we would like to note that our present DF-CCSD gradient code assumes that there is enough memory for amplitudes, DF integrals, and CC contractions. With a future implementation of frozen-natural orbitals (FNO)^{15,80–82} we may reduce the memory requirements of the present code and the computational cost (especially of the PPL term, as discussed in Ref. 15). Finally, we note that our DF-CCSD gradient code makes use of shared-memory parallelism through threaded BLAS calls as well as OpenMP parallelization of all tensor manipulations.

### H. Verification of DF-CCSD gradients code

Our brand new DF-CCSD energy, DF-CCSD-Λ, and DF-CCSD-PDM codes are separately tested against the conventional CCSD (in the ccenergy module), CCSD-Λ (in the cclambda module), and CCSD-PDM (in the ccdensity module) codes of Psi4 program.^{55} To test the DF-CCSD energy code with respect to the conventional CCSD code, we use large auxiliary basis sets for small primary basis sets (for example, cc-pV5Z-RI for cc-pVDZ) at first, and we match energies up to 4-5 decimal places.

Further, we also implemented Cholesky decomposed CCSD codes: CD-CCSD energy, CD-CCSD-Λ, and CD-CCSD-PDM codes, but not the gradients. Except for the 3-index integrals used, the DF-CCSD and CD-CCSD codes are exactly the same. With very tight CD tolerances such as 10^{−14} and tight energy convergence criteria, such as 10^{−10}, we match CD-CCSD and CCSD energies up to 9-10 decimal places. The same strategy is also applied to the CD-CCSD-Λ code, where pseudo-energies are compared, and CD-CCSD-PDM code, where energies computed via PDMs are compared. Further, we would like to note that PDM codes for DF/CD and conventional CCSD are very different since in the former one we use 3-index TPDM whereas in the latter one 4-index TPDMs are used. Hence, the verification of DF/CD-CCSD-PDM code requires a great effort. For this purpose we “hacked” Crawford’s ccdensity code and compared 4-index TPDM terms one by one with their contributions in the 3-index TPDM. Most of the 4-index TPDM terms contribute to the different blocks of 3-index TPDM reported here.

Finally, we verified our overall analytic gradients code with respect to the numerical gradients (with the 5-point formulas).

## IV. RESULTS AND DISCUSSION

Results from the CCSD and DF-CCSD methods were obtained for a set of alkanes for comparison of the computational cost for single-point analytic gradient computations. Geometries for the considered alkanes were reported in our recent study.^{16} For the alkanes Dunning’s correlation-consistent polarized valence triple-*ζ* basis set (cc-pVTZ) was employed with the frozen core approximation.^{83,84} For the cc-pVTZ primary basis sets, cc-pVTZ-JKFIT^{6} and cc-pVTZ-RI^{85} auxiliary basis sets were employed for reference and correlation energies, respectively.

Additionally, the CCD, CCSD, DF-CCD, and DF-CCSD methods were applied to a set of small molecules^{86} for comparison of geometries and vibrational frequencies. For computations of geometries and frequencies, Dunning’s correlation-consistent polarized valence quadruple-*ζ* basis set (cc-pVQZ) was employed with the frozen core approximation.^{83,84} For the cc-pVQZ primary basis sets, cc-pVQZ-JKFIT^{6} and cc-pVQZ-RI^{85} auxiliary basis sets were employed for reference and correlation energies, respectively. The DF-CCD and DF-CCSD computations were performed with our present codes, which are available in the Psi4 program package,^{55} while the conventional CCSD computations were carried out with the Molpro2012.1 program.^{56–58}

### A. Efficiency of DF-CCSD analytic gradients

To assess the efficiency of the DF-CCSD analytic gradients, we consider a set of alkanes.^{16} The total computational (wall) times for single-point frozen-core analytic gradient computations using the DF-CCSD and CCSD methods are presented graphically in Figure 1. The DF-CCSD method noticeably reduces the total computational cost compared to CCSD; there is a ∼20% reduction compared to CCSD (from Molpro2012.1) for the largest member of the alkanes set.

Next, we analyze where this cost saving is coming from. Figure 1 also shows the total computational time broken down into three major components: (1) the total time to converge the CCSD energy and amplitudes, (2) the total time to converge the Λ equations, and (3) the cost of the remaining gradient terms (defined as above to include the formation of the PDMs and GFM, the solution of the Z-vector equation, the formation of the relaxed PDMs and GFM, back-transformation of the relaxed PDMs and GFM to the AO basis, and evaluation of the gradients in the AO basis). We see that 46% of the time is spent in the amplitude equations, 47% is spent in the Λ equations, and 7% is spent on the remaining gradient terms. Clearly, DF-CCSD demonstrates a very large speedup in the gradient terms, while it is actually somewhat slower for the other terms.

The reason Molpro is faster than the current DF-CCSD code for the solution of the amplitude equations is that it uses an alternative, integral-direct AO-basis (IDAOB) algorithm for the PPL term. This algorithm, like our DF-CCSD algorithm, also avoids costly I/O of 〈*ab*|*cd*〉 integrals. It has the disadvantage that the AO-basis integrals are computed each iteration, but it has the compensating advantage that efficient prescreening can substantially reduce the number of computed integrals due to the sparsity of the AO basis. In previous studies, it was argued that the integral-direct AO-basis algorithm is expected to approach linear scaling at large molecular limits.^{87} All scaling analyses are based on a large ratio of number of virtual vs. number of occupied orbitals, such that the following Equations (74)–(76) will represent the leading computational cost. For the closed-shell case, the PPL term can be evaluated in the AO basis as follows:^{87}

where

and

Similar to the MO basis algorithm, one can adapt the symmetric-antisymmetric tensor decompositions for the AO basis algorithm. With this algorithm, the cost of PPL may be estimated as $14O2N4+2O2VN2+2O2V2N$. Hence, for the PPL term, the theoretical cost of the IDAOB algorithm is roughly similar to the cost of our MO-basis algorithm with DF integrals outlined in Section III G. The main difference is that the DF algorithm involves the costly *V*^{4}*N _{aux}* construction of the 〈

*ab*|

*cd*〉 integrals while the IDAOB algorithm replaces this with the direct computation of the AO integrals. The locality of the AO basis can result in a lower cost for the IDAOB algorithm.

The linear alkanes considered for our timing tests are ideal for sparsity in the AO-basis electron repulsion integrals. For the largest member of the alkanes set, C_{10}H_{22}, the number of active occupied orbitals, virtual orbitals, and auxiliary basis functions are *O* = 31, *V* = 567, and *N _{aux}* = 1470. For the largest test case, C

_{10}H

_{22}, the total time to converge the CCSD amplitude equations is 2125.82 min for DF-CCSD, and 1258.64 min for conventional CCSD in Molpro2012.1 using the IDAOB algorithm, a difference of 41% [the total number of CCSD iterations is 11 and 12 for CCSD and DF-CCSD, respectively].

The use of frozen natural orbitals (FNOs) greatly reduces the computational cost of the PPL term in DF-CCSD, simply by reducing the number of retained virtual orbitals, *V*. On the other hand, the use of FNOs will have less of an effect on the cost of the IDAOB algorithm, which still has to compute all the non-negligible AO-basis ERIs. Hence, the use of FNOs may help to make DF-CCSD faster relative IDAOB-based conventional CCSD for the solution of the amplitude equations. As a point of comparison, computing the CCSD energy for C_{10}H_{22} using the FNO-DF-CCSD code of DePrince and Sherrill^{15} requires 1740.77 min on the same hardware, with a conservative FNO cutoff of 10^{−5}, while it requires only 306.17 min with modest FNO cutoff of 10^{−4}.

Next, we consider the solution of the CCSD Λ-amplitudes. The conventional CCSD code of Molpro is again faster than our present DF-CCSD implementation because the Λ equations contain a PPL term analogous to the one discussed above for the amplitude equations. For the largest member of the alkanes set, C_{10}H_{22}, the total computational time to converge the Λ equations is 2198.12 (DF-CCSD) and 1733.79 (CCSD) min, a difference of 21.1% [the total number of CCSD Λ-amplitude iterations is 12 and 11 for CCSD and DF-CCSD, respectively]. Once again, the FNO approach will substantially speed up the PPL term of the Λ-amplitudes in the DF-CCSD algorithm.

Finally we consider the evaluation of the remaining gradient related terms, which is the most time consuming part for the conventional CCSD code of Molpro. By gradient related terms we mean the computations of PDMs, GFM, the solution of the Z-vector equation, formation of the relaxed PDMs and GFM, back-transformation of relaxed PDMs and GFM to the AO basis, and evaluation of gradients in the AO basis. As shown in Figure 1, the DF-CCSD method substantially reduces the computational cost compared to CCSD. For the largest member of the alkanes set, the computational time for the gradient terms is 2582.59 (CCSD) and 310.67 (DF-CCSD) min, respectively. Hence, our DF-CCSD code is more than 8-fold faster than the conventional CCSD code of Molpro for the gradient terms for C_{10}H_{22}. This speed increase is due to several reasons. As it is obvious from Eq. (32), we only need to consider the 3-index integrals and 3-index TPDM instead of the 4-index one for the evaluation of gradients. Hence, the scaling of the gradient equation (32) is reduced to *O*(*N*^{3}) from *O*(*N*^{4}); costs of back transformation for TPDM and the contraction of TPDM with the 4-index AO-basis ERI derivative are significantly reduced due to a reduction in the I/O costs. Further, with our definition of the 3-index TPDM, the GFM can be readily evaluated.^{17} In case of the conventional integrals one needs to consider four-index integrals including four virtual orbitals (〈*ab*| |*cd*〉 type integrals), while in case of DF formalism one just needs to three-index DF factors.^{17} More specifically, with the conventional integrals the most expensive term for GFM is^{71,86,88,89}

which requires the simultaneous handling of 4-virtuals integrals and TPDM. However, with the DF approach this problematic term is replaced by

hence, for the most expensive term the cost of GFM can be reduced to *V*^{3}∗*N _{aux}* from

*V*

^{5}, and instead of out-of-core storage of 4-virtuals integrals and TPDM one can use in-core memory storage for the DF factors and three-index TPDM. For the entire GFM, with the DF approach the computational cost is reduced to the order of

*O*(

*N*

^{4}) instead of

*O*(

*N*

^{5}). Further, as it is demonstrated in our previous study,

^{17}the computational cost for the solution of the Z-vector equation can be dramatically reduced with the DF approach. Finally, we illustrate the scaling with respect to the number of computer cores for the DF-CCSD gradients for the alkane test set in Figure 2. As shown in Figure 2, our algorithm scales fairly well to six cores for the large molecules.

### B. Accuracy of DF-CCSD analytic gradients

In this section we assess the accuracy of the DF-CCSD method. For this purpose we consider equilibrium geometries and harmonic vibrational frequencies. Geometry optimizations and vibrational frequency computations are carried out with a 10^{−8} energy convergence tolerance. In order to assess the performance of DF-CCD and DF-CCSD for bond lengths, we consider a set of small molecules employed in previous studies.^{86,90} Table I reports the bond lengths of the molecules considered in increasing order. Errors for DF-CCSD with respect to conventional CCSD are depicted in Figure 3. The mean absolute error (Δ_{mae}), the standard deviation of errors (Δ_{std}), and the maximum absolute error (Δ_{max}) are ≤0.0001, 0.0001, and 0.0001 Å, respectively. Further, for the bond angles of molecules considered, Δ_{mae}, Δ_{std}, and Δ_{max} values are 0.059°, 0.083°, and 0.206°, respectively. The performance of the CCSD and DF-CCSD methods is essentially the same, as in the case of MP2/DF-MP2.^{17} Hence, the error introduced by the DF approach is quite negligible. A similar situation also occurs for the CCD/DF-CCD case, errors for DF-CCD with respect to conventional CCD are depicted in the supplementary material.^{91} Further, we consider transition-states for three reactions^{92,93} to assess the accuracy of DF-CCSD for geometries (reported in the supplementary material^{91}). As in the case of minima, the Δ_{mae} value of DF-CCSD is ≤10^{−4} Å for bond lengths when compared to conventional CCSD.

. | Molecule . | Bond . |
---|---|---|

1 | HF | R_{HF} |

2 | H_{2}O | R_{OH} |

3 | H_{2}O_{2} | R_{OH} |

4 | HOF | R_{OH} |

5 | HNC | R_{NH} |

6 | Acetamide | R_{NH1} |

7 | Acetamide | R_{NH2} |

8 | NH_{3} | R_{NH} |

9 | N_{2}H_{2} | R_{NH} |

10 | HNO | R_{NH} |

11 | C_{2}H_{2} | R_{CH} |

12 | HCN | R_{CH} |

13 | Propene | R_{CH1} |

14 | Benzene | R_{CH} |

15 | C_{2}H_{4} | R_{CH} |

16 | Propene | R_{CH2} |

17 | Propene | R_{CH3} |

18 | Acetamide | R_{CH1} |

19 | CH_{4} | R_{CH} |

20 | Acetamide | R_{CH2} |

21 | Propene | R_{CH4} |

22 | Acetamide | R_{CH3} |

23 | Propene | R_{CH5} |

24 | Propene | R_{CH6} |

25 | N_{2} | R_{NN} |

26 | CH_{2}O | R_{CH} |

27 | CO | R_{CO} |

28 | HCN | R_{CN} |

29 | CO_{2} | R_{CO} |

30 | HNC | R_{CN} |

31 | HNO | R_{NO} |

32 | CH_{2}O | R_{CO} |

33 | C_{2}H_{2} | R_{CC} |

34 | Acetamide | R_{CO} |

35 | BH | R_{BH} |

36 | N_{2}H_{2} | R_{NN} |

37 | C_{2}H_{4} | R_{CC} |

38 | Propene | R_{CC1} |

39 | Acetamide | R_{CN} |

40 | Benzene | R_{CC} |

41 | HOF | R_{OF} |

42 | H_{2}O_{2} | R_{OO} |

43 | Propene | R_{CC2} |

44 | Acetamide | R_{CC} |

. | Molecule . | Bond . |
---|---|---|

1 | HF | R_{HF} |

2 | H_{2}O | R_{OH} |

3 | H_{2}O_{2} | R_{OH} |

4 | HOF | R_{OH} |

5 | HNC | R_{NH} |

6 | Acetamide | R_{NH1} |

7 | Acetamide | R_{NH2} |

8 | NH_{3} | R_{NH} |

9 | N_{2}H_{2} | R_{NH} |

10 | HNO | R_{NH} |

11 | C_{2}H_{2} | R_{CH} |

12 | HCN | R_{CH} |

13 | Propene | R_{CH1} |

14 | Benzene | R_{CH} |

15 | C_{2}H_{4} | R_{CH} |

16 | Propene | R_{CH2} |

17 | Propene | R_{CH3} |

18 | Acetamide | R_{CH1} |

19 | CH_{4} | R_{CH} |

20 | Acetamide | R_{CH2} |

21 | Propene | R_{CH4} |

22 | Acetamide | R_{CH3} |

23 | Propene | R_{CH5} |

24 | Propene | R_{CH6} |

25 | N_{2} | R_{NN} |

26 | CH_{2}O | R_{CH} |

27 | CO | R_{CO} |

28 | HCN | R_{CN} |

29 | CO_{2} | R_{CO} |

30 | HNC | R_{CN} |

31 | HNO | R_{NO} |

32 | CH_{2}O | R_{CO} |

33 | C_{2}H_{2} | R_{CC} |

34 | Acetamide | R_{CO} |

35 | BH | R_{BH} |

36 | N_{2}H_{2} | R_{NN} |

37 | C_{2}H_{4} | R_{CC} |

38 | Propene | R_{CC1} |

39 | Acetamide | R_{CN} |

40 | Benzene | R_{CC} |

41 | HOF | R_{OF} |

42 | H_{2}O_{2} | R_{OO} |

43 | Propene | R_{CC2} |

44 | Acetamide | R_{CC} |

Next, we consider harmonic vibrational frequencies to assess the accuracy of DF-CCSD. For this purpose we consider the molecules shown in Table I. The Δ_{mae}, Δ_{std}, and Δ_{max} values are 0.5, 0.9, and 2.5 cm^{−1}, respectively. Thus, results for CCSD and DF-CCSD are essentially the same, as in the case of bond lengths, and the error introduced by the DF approach is again negligible. A similar situation also holds for the CCD/DF-CCD case, as shown in the supplementary material.^{91} For the transition states considered, errors of the DF-CCSD method are again quite negligible (Δ_{mae} = 0.2 cm^{−1}, Δ_{std} = 0.3 cm^{−1}, and Δ_{max} = 0.9 cm^{−1}) compared to conventional CCSD.

## V. CONCLUSIONS

In this study, an efficient implementation of frozen core analytic gradients for the CCSD method with the DF approximation, which is denoted as DF-CCSD, has been reported. The computational time of the DF-CCSD single point analytic gradients has been compared with that of the conventional CCSD (from Molpro2012.1 package^{56–58}). The DF-CCSD method noticeably reduces the computational cost compared to CCSD, with a ∼20% reduction for the largest member of alkane test set considered, C_{10}H_{22} in a cc-pVTZ basis set. This cost savings result from a dramatic reduction in the cost of the “gradient terms,” i.e., the computation of the PDMs, GFM, solution of the Z-vector equation, formation of the relaxed PDMs and GFM, back-transformation of the relaxed PDMs and GFM to the AO basis, and evaluation of gradients in the AO basis. The DF-CCSD method is substantially more efficient than CCSD for these terms because four-index integrals and TPDMs can be avoided, leading to reduced scaling of certain terms and reduced I/O and storage requirements. For C_{10}H_{22}, the computational times for the gradient terms is 2583 (CCSD) and 311 (DF-CCSD) min, i.e., more than an 8-fold speedup for DF-CCSD for these terms.

The computational cost of the amplitude and Λ equations can also be reduced by density fitting, especially on machines with slow or network-attached disks, because then one can avoid costly I/O operations on the 〈*ab*|*cd*〉 integrals that enter into the particle-particle ladder (PPL) terms. On the other hand, density fitting can be more costly than conventional CCSD for the PPL terms in the amplitude and Λ equations when the conventional code uses an IDAOB algorithm. In particular, re-computing the electron repulsion integrals in the IDAOB CCSD algorithm can be faster than the construction of the MO-basis 〈*ab*|*cd*〉 integrals from the density-fitting intermediates in the DF-CCSD algorithm. This is especially true for larger and/or more extended molecules like linear alkanes, which exhibit significant sparsity in the AO integrals. For C_{10}H_{22}, the amplitude equations took 1258.64 (CCSD) and 2125.82 (DF-CCSD) min to converge, a difference of 40.8%, while the Λ equations took 2198.12 (CCSD) and 1733.79 (DF-CCSD) min, a difference of 21.1%. On the other hand, using FNOs one can dramatically speed up the cost of the PPL term in DF-CCSD.^{15} We intend to investigate DF-CCSD gradients with FNOs in future work.

Finally, the errors introduced by the density-fitting approximation have been assessed by applying DF-CCSD to compute the equilibrium geometries and harmonic vibrational frequencies of a variety of molecules. Our results show that the DF approach introduces negligible errors compared to CCSD: equilibrium bond lengths are within 10^{−4} Å, and harmonic vibrational frequencies exhibit mean absolute errors of 0.5 cm^{−1} for the cases considered.

## Acknowledgments

This research was supported by the Scientific and Technological Research Council of Turkey (Grant No. TÜBİTAK-114Z786), the European Cooperation in Science and Technology (Grant No. CM1405), and the U.S. National Science Foundation (Grant Nos. ACI-1147843 and CHE-1300497). U.B. also acknowledges the support from Turkish Academy of Sciences, Outstanding Young Scientist Award (TÜBA-GEBİP 2015). We thank Professor Eugene DePrince of Florida State University for many helpful discussions. U.B. thanks the Department of Chemistry, Middle East Technical University for hosting a visit during the summer season.

### APPENDIX A: SPIN-ORBITAL EQUATIONS FOR DF-CCSD INTERMEDIATES

### APPENDIX B: SPIN-ORBITAL EQUATIONS FOR DF-CCSD-Λ INTERMEDIATES

### APPENDIX C: SPIN-ORBITAL EQUATIONS FOR CORRELATION PDM INTERMEDIATES

## REFERENCES

*ab initio*programs, 2012, see http://www.molpro.net.

By gradient terms we mean the computations of PDMs, GFM, solution of the Z-vector equation, formation of the relaxed PDMs and GFM, back-transformation of relaxed PDMs and GFM to the AO basis, and evaluation of gradients in the AO basis. All these tasks are performed by our general DF gradient code,^{17} which is called as Dfgrad.