Analytic energy gradients for the coupled-cluster singles and doubles method with the density-fitting approximation

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.³¹ reported the initial algebraic expressions for analytic gradients of the coupled-cluster singles and doubles (CCSD) method,³² following the Z-vector method of Handy and Schaefer.³³ In a 1985 study, Fitzgerald et al.³⁴ 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.³⁵ 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),⁴⁶ complete active space second-order perturbation theory (DF-CASPT2),⁴⁷ MP2 with DF-HF reference (DF-MP2),¹⁷ orbital-optimized MP2 (DF-OMP2),¹⁸ and local time-dependent coupled cluster response theory (DF-TD-LCC2).⁴⁸ However, analytic gradients cannot be computed with the CD approximation unless one employs an approach like that of Aquilante et al.⁴⁹ 

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.⁵⁵ 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.

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⁵) to O(N⁴). Moreover, the computational time for the solution of the Z-vector equation³³ is significantly reduced with the DF intermediates.¹⁷ Finally, we would like to note that the orbital relaxation terms are computed more efficiently with the DF-HF reference.¹⁷ 

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

the DF factor $bμν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.

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ˆN$ is the normal-ordered Hamiltonian operator,^25,59 |0〉 is the reference determinant, and $Tˆ$ is the cluster excitation operator. For the CCSD wave function,

where $Tˆ1$ and $Tˆ2$ are single- and double-excitation operators, respectively.

where $tia$ and $tijab$ are the single and double excitation amplitudes, and $aˆ†$ and $iˆ$ 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 $〈Φia|$ and $〈Φijab|$ are singly and doubly excited Slater determinants, respectively.

The explicit form of the T₁ amplitude equation with the DF approximation can be written as

and the T₂ 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

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

where $Hˆ$ is the Hamiltonian operator and $Λˆ$ is the coupled-cluster de-excitation operator. For the CCSD wave function,

where $Λˆ1$ and $Λˆ2$ are the single and double de-excitation operators, respectively.

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

The CCSD T-amplitude equations are obtained by requiring that $L˜$ 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⁶² 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⁶² are not present. The explicit form of the Λ₁-amplitude equation with the DF approximation is given by

and the explicit form of Λ₂-amplitude equation with the DF approximation is given by

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

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,³³ 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ˆ$ 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:³³ 

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.¹⁷ 

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:¹⁷ 

where $γpqeff$⁠, $ΓPQeff$⁠, $Γ˜pqQ(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¹⁷ 

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 $Γ˜pqQ$⁠. In a 2014 study,¹⁷ it has been demonstrated that with the definition of $Γ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 $Γμν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:¹⁷ 

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.¹⁷ Hence, we will not repeat those detailed discussions.

The OPDM can be decomposed as follows:

where $γpqref$ and $γ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 $ΓpqQ(ref)$ and $ΓpqQ(corr)$ are the reference and correlation contributions to TPDM, and $Γ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.⁶ 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.

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.

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

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

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

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

4. Perform DF-CCSD λ-amplitude iterations to obtain the converged λ-amplitudes.

5. Compute the unrelaxed PDMs and GFM.

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

7. Solve Eq. (30) to obtain the Z-vector (DF-SCF integrals should be used).¹⁷ 

8. Compute the effective OPDM, GFM, and 3-index TPDM (⁠$ΓpqQ(sep)$ should be considered).¹⁷ 

9. Compute $Γ˜pqQ$ from Eq. (36) for $ΓpqQ(ref)$⁠, $ΓpqQ(sep)$⁠, and $ΓpqQ(corr)$ separately. For $ΓpqQ(ref)$ and $Γ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 $ΓpqQ(corr)$⁠.

10. Combine $ΓpqQ(ref)$ and $ΓpqQ(sep)$ to form $ΓpqQ(RefSep)$⁠.

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

12. 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 $ΓμνQ(RefSep)/ΓPQ(RefSep)$ and $ΓμνQ(corr)/ΓPQ(corr)$⁠, respectively.

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₂ and Λ₂ amplitude equations and the V_abef 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 a, requiring dynamic arrays with the size of V³ (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₂ and Λ₂ 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,⁷⁸ and Scuseria, Janssen, and Schaefer⁷⁹ 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