11.1 Characteristics of the Implementation and Computational Demands

In CCSD the ground–state energy is (as for CC2) evaluated as \[\begin{align} E_{\mathrm{CC}} & = \langle\mathrm{HF}| H |\mathrm{CC}\rangle = \langle\mathrm{HF}| H \exp(T) |\mathrm{HF}\rangle ~~, \end{align}\](11.1) where the cluster operator \(T = T_1 + T_2\) consist of linear combination of single and double excitations: \[\begin{align} T_1 & = \sum_{ai} t_{ai} \tau_{ai} ~, \\ T_2 & = \frac{1}{2}\sum_{aibj} t^{ij}_{ab} \tau_{aibj} ~. \end{align}\](11.2–11.3) In difference to CC2, the cluster amplitudes \(t_{ai}\) and \(t_{aibj}\) are determined from equations which contain no further approximations apart from the restriction of \(T\) to single and double excitations: \[\begin{align} \Omega_{\mu_1} & = \langle\mu_{1} | \hat{\tilde{H}} + [\hat{\tilde{H}},T_2] |\mathrm{HF}\rangle = 0 ~~, \\ \Omega_{\mu_2} & = \langle\mu_{2} | \hat{\tilde{H}} + [\hat{\tilde{H}},T_2] + [[\hat{H},T_2],T_2]|\mathrm{HF}\rangle = 0 ~~, \end{align}\](11.4–11.5) where again \[\hat{\tilde{H}} = \exp(-T_1) \hat{H} \exp(T_1),\] and \(\mu_{1}\) and \(\mu_{2}\) are, respectively, the sets of all singly and doubly excited determinants. For MP3 the energy is computed from the first-order amplitudes (\(t^{(1)}_{\mu}\)) as \[\begin{align} E_{\mathrm{MP3},tot} = & E_{\mathrm{HF}}+E_{\mathrm{MP2}} + E_{\mathrm{MP3}} \\ = & \langle \mathrm{HF}| \hat{H} + [\hat{H},T_2^{(1)}] | \mathrm{HF} \rangle + \sum_{\mu_2} t_{\mu_2}^{(1)} \langle \mu_2 | [\hat{W},T^{(1)}_2] | \mathrm{HF} \rangle \end{align}\](11.6–11.7) with \(\hat{W} = \hat{H}-\hat{F}\). To evaluate the fourth-order energy one needs in addition to the first-order also the second-order amplitudes, which are obtained from the solution of the equations \[\begin{align} \langle \mu_{1}|[\hat{F},T_1^{(2)}] + [\hat{W},T^{(1)}_2]|\mathrm{HF}\rangle = & \, 0 \\ \langle \mu_{2}|[\hat{F},T_2^{(2)}] + [\hat{W},T^{(1)}_2]|\mathrm{HF}\rangle = & \, 0 \\ \langle \mu_{3}|[\hat{F},T_3^{(2)}] + [\hat{W},T^{(1)}_2]|\mathrm{HF}\rangle = & \, 0 \end{align}\](11.8–11.10) From these the fourth-order energy correction is computed as: \[\begin{align} E_{\mathrm{MP4}} = & \sum_{\mu_2} t_{\mu_2}^{(1)} \langle \mu_2 | [\hat{W},T^{(2)}_1 + T^{(2)}_2 + T^{(2)}_3] + [[\hat{W},T^{(1)}_2],T^{(1)}_2] | \mathrm{HF} \rangle ~. \end{align}\](11.11) Eqs. (11.5) and (11.7) – (11.11) are computational much more complex and demanding than the corresponding doubles equations for the CC2 model. If \({\cal N}\) is a measure for the system size (e.g. the number of atoms), the computational costs (in terms of floating point operations) for CCSD calculations scale as \({\cal O}({\cal N}^6)\). If for the same molecule the number of one-electron basis functions \(N\) is increased the costs scale with \({\cal O}({\cal N}^4)\). (For RI-MP2 and RI-CC2 the costs scale with the system size as \({\cal O}({\cal N}^5)\) and with the number of basis functions as \({\cal O}({\cal N}^3)\).) The computational costs for an MP3 calculations are about the same as for one CCSD iteration. For MP4 the computational costs are comparable to those for two CCSD iteration plus the costs for the perturbation triples correction (see below).

Explicitly-correlated CCSD(F12) methods:

In explicitly-correlated CCSD calculations the double excitations into products of virtual orbitals, described by the operator \(T_2 = \frac{1}{2}\sum_{aibj} t^{ij}_{ab} \tau_{aibj}\), are augmented with double excitations into the explicitly-correlated pairfunctions (geminals) which are described in Sec. 9.5: \[\begin{align} T & = T_1 + T_2 + T_{2'} \\ T_{2'} & = \frac{1}{2} \sum_{ijkl} c^{kl}_{ij} \tau_{kilj} \end{align}\](11.12–11.13) where \(\tau_{kilj} | ij\rangle = \hat{Q}_{12} f_{12} | kl \rangle\) (for the definition \(\hat{Q}_{12}\) and \(f_{12}\) see Sec. 9.5). This enhances dramatically the basis set convergence of CCSD calculations [21]. Without any further approximations than those needed for evaluating the necessary matrix elements, this extension of the cluster operator \(T\) leads to the CCSD-F12 method. CCSD(F12) is an approximation [290, 21] to CCSD-F12 which neglects certain computationally demanding higher-order contributions of \(\hat{T}_{2'}\). This reduces the computational costs dramatically, while the accuracy of CCSD(F12) is essentially identical to that of CCSD-F12 [291, 292]. In the CCSD(F12) approximation the amplitudes are determined from the equations: \[\begin{align} \Omega_{\mu_1} & = \langle\mu_{1} | \hat{\tilde{H}} + [\hat{\tilde{H}},T_2+T_{2'}] |\mathrm{HF}\rangle = 0 ~~, \\ \Omega_{\mu_2} & = \langle\mu_{2} | \hat{\tilde{H}} + [\hat{\tilde{H}},T_2+T_{2'}] + [[\hat{H},T_2+2T_{2'}],T_2]|\mathrm{HF}\rangle = 0 ~~, \\ \Omega_{\mu_{2'}} & = \langle\mu_{2'} | [\hat{F},T_{2'}] + \hat{\tilde{H}} + [\hat{\tilde{H}},T_2] |\mathrm{HF}\rangle = 0 ~~. \end{align}\](11.14–11.16) Similar as for MP2-F12, also for CCSD(F12) the coefficients for the doubles excitations into the geminals, \(c^{kl}_{ij}\) can be determined from the electronic cusp conditions using the rational generator (also known as SP or fixed amplitude) approach. In this case Eq. (11.16) is not solved. To account for this, the energy is then computed from a Lagrange function as: \[\begin{align} E_{\mathrm{CCSD(F12)-SP}} & = L_{\mathrm{CCSD(F12)}} = \langle\mathrm{HF}| H |\mathrm{CC}\rangle + \sum_{\mu_{2'}} c_{\mu_{2'}} \Omega_{\mu_{2'}} \end{align}\](11.17) This is the recommended approach which is used by default if not any other approach has been chosen with the examp option in $rir12 (see Sec. 9.5 for further details on the options for F12 calculations; note that the examp noinv option should not be combined with CCSD calculations). CCSD(F12)-SP calculations are computationally somewhat less expensive that CCSD(F12) calculations which solve Eq. (11.16), while both approaches are approximately similar accurate for energy differences.

The SP approach becomes in particular very efficient if combined with the neglect of certain higher-order explicitly-correlated contributions which have a negligible effect on the energies but increase the costs during the CC iterations. The most accurate and recommended variant is the CCSD(F12*) approximation [22], which gives essentially identical energies as CCSD(F12)-SP. Also available are the CCSD[F12] (Ref. [22]), CCSD-F12a (Ref. [293]) and CCSD-F12b (Ref. [294]) approximations as well as the perturbative corrections CCSD(2)\(_{\overline{\mathrm{F12}}}\) and CCSD(2*)\(_{\overline{\mathrm{F12}}}\) (see Refs. [295, 296, 22]) and CCSD\(_{\mathrm{(F12*)}}\) (Ref. [247]). These approximations should only be used with ansatz 2 and the SP approach (i.e. fixed geminal amplitudes). Note that the CCSD-F12c method which is available in the molpro package is the same as CCSD(F12*). Since CCSD(F12*) is the original name that was introduced when the method was proposed for the first time in Ref. [22] users are asked to use in publications this abbreviation and cite the original reference [22].

For MP3 the approximations (F12*), and (F12) to a full F12 implementation become identical: they include all contributions linear in the coefficients \(c^{kl}_{ij}\). The explicitly-correlated MP4 method MP4(F12*) is defined as fourth-order approximation to CCSD(F12*)(T). Note that MP4(F12*) has to be used with the SP or fixed amplitude approach for the geminal coefficients \(c^{kl}_{ij}\). MP3(F12*) and MP4(F12*) are currently only available for closed-shell or unrestricted Hartree-Fock reference wavefunctions.

The CPU time for a CCSD(F12) calculation is approximately the sum of the CPU time for an MP2-F12 calculation with the same basis sets plus that of a conventional CCSD calculation multiplied by \((1+N_{CABS}/N)\), where \(N\) is the number of basis and \(N_{CABS}\) the number of complementary auxiliary basis (CABS) functions (typically \(N_{CABS} \approx 2-3 N\)). If the geminal coefficients are determined by solving Eq. (11.16) instead of using fixed amplitudes, the costs per CCSD(F12) iteration increase to \(\approx (1+2N_{CABS}/N)\) times the costs for a conventional CCSD iteration. Irrespective how the geminal coefficients are determined, the disc space for CCSD(F12) calculations are approximately a factor of \(\approx (1+2N_{CABS}/N)\) larger than the disc space required for a conventional CCSD calculation. Note that this increase in the computational costs is by far outweighted by the enhanced basis set convergence.

In combination with the CCSD(F12*) approximation (and also CCSD[F12], CCSD-F12a, CCSD-F12b, CCSD(2)\(_{\overline{\mathrm{F12}}}\), CCSD(2*)\(_{\overline{\mathrm{F12}}}\) and CCSD\(_{\mathrm{(F12*)}}\)) the CPU time for the SP approach is only about 20% or less longer than for a conventional CCSD calculation within the same basis set.

CC calculations with restricted open-shell (ROHF) references:

The MP2 and all CC calculations for ROHF reference wavefunctions are done by first transforming to a semi-canonical orbital basis which are defined by the eigenvectors of the occupied/occupied and virtual/virtual blocks of the Fock matrices of alpha and beta spin. No spin restrictions are applied in the cluster equations. This approach is sometimes also denoted as ROHF-UCCSD.

Note that if a frozen-core approximation is used, the semicanonical orbitals depend on whether the block-diagonalization of the Fock matrices is done in space of all orbitals or only in the space of the correlated valence orbitals. The two approaches lead thus to slightly different energies, but neither one is more valid or more accurate than the other. Both schemes are available through the specification of core res/can. The default is core can where the full occupied block is diagonalised. The same scheme is used e.g. in the CFOUR program suite, but other codes as e.g. the implementation in MOLPRO use a block-diagonalization restricted to the active valence space. The semi-canonical orbitals are written to the $sc_alph and $sc_bet data groups.

Brueckner Calculations and KS-CCSD(T):

Brueckner orbital optimisation in place of single excitations is available, which is useful when the HF reference relaxes significantly in the presence of correlation (for example for calculations on heavy elements). The default way a calculation proceeds is by first building the Fock matrix for the input starting orbitals, semi-canonicalising these orbitals and then optimising the (non-frozen) orbitals at the same time as solving for the doubles amplitudes, such that the singles residual vanishes at convergence. The default is therefore UBCCD for open-shell references. The final orbitals are written in the $bcc_mos or $bcc_alph, $bcc_bet data groups. It is possible through options in the $ricc2 data group to skip the initial semi-canonicalisation and to control the nature of the frozen orbitals. It is also possible to request restricted open-shell BCCD where the Brueckner optimisation applies only to \(t_1^\alpha + t_1^\beta\) such that the alpha and beta doubly occupied orbitals of a starting ROHF reference remain identical to each other during the optimisation, and \(t_1^\alpha - t_1^\beta\) is no longer zero at convergence. Since it is necessary to semi-canonicalise the converged orbitals in order to perform a (T) calculation, the final orbitals are always in semi-canonical form if (T) is requested. Since the F12 integrals depend on the orbitals, the iterative BCCD(F12*) approximation is very expensive. Instead, it is recommended to use the perturbative equivalent BCCD(T)\(_{\mathrm{(F12*)}}\).

As an alternative to BCCD theory, it is possible (but not recommended) to used Kohn-Sham orbitals to define the reference. When using orbitals other than converged HF orbitals, it is necessary to add the $non-canonical MOs data group to the control file. When performing an F12 calculation, the option r12orb arb should be set in the $rir12 data group.

Perturbative triples corrections:

To achieve ground state energies with a high accuracy that systematically surpasses the accuracy MP2 and DFT calculations for reaction and binding energies, the CCSD model has to be combined with a perturbative correction for connected triples. The recommended approach for the correction is the CCSD(T) model \[\begin{equation} E_{\mathrm{CCSD(T)}} = E_{\mathrm{CCSD}} + E^{(4)}_{DT} + E^{(5)}_{ST} \end{equation}\](11.18) which includes the following two terms: \[\begin{align} E^{(4)}_{DT} & = \sum_{\mu_2} t_{\mu_2}^{CCSD} \langle \mu_2 |[H,T_3^{(2)}]|\mathrm{HF}\rangle \\ E^{(5)}_{ST} & = \sum_{\mu_1} t_{\mu_1}^{CCSD} \langle \mu_2 |[H,T_3^{(2)}]|\mathrm{HF}\rangle \end{align}\](11.19–11.20) where the approximate triples amplitudes are evaluated as: \[\begin{align} t^{(2)}_{aibjck} = - \frac{\langle {{ijk}\atop{abc}}|[\hat{H},T_2]|\mathrm{HF}\rangle }{ \epsilon_a-\epsilon_i + \epsilon_b-\epsilon_j + \epsilon_c-\epsilon_k } \end{align}\](11.21) In the literature one also finds sometimes the approximate triples model CCSD[T] (also denoted as CCSD+T(CCSD)), which is obtained by adding only \(E^{(4)}_{DT}\) to the CCSD energy. Usually CCSD(T) is slightly more accurate than CCSD[T], although for closed-shell or spin-unrestricted open-shell reference wavefunctions the energies of both models, CCSD(T) and CCSD[T] model, are correct through 4.th order perturbation theory. For a ROHF reference, however, \(E^{(5)}_{ST}\) contributes already in 4.th order and only the CCSD(T) model is correct through 4.th order perturbation theory.

Integral-direct implementation and resolution-of-the-identity approximation:

The computationally most demanding (in terms floating point operations) steps of a CCSD calculation are related to two kinds of terms. One of the most costly steps is the contraction \[\begin{equation} \Omega^B_{aibj} = \sum_{cd} t^{ij}_{cd} (ac|bd) \end{equation}\](11.22) where \(a\), \(b\), \(c\), and \(d\) are virtual orbitals. For small molecules with large basis sets or basis sets with diffuse functions, where integral screening is not effective, it is the time-determing step and can most efficiently be evaluated with a minimal operation count of \(\tfrac{1}{4} O^2 V^2\) (where \(O\) and \(V\) are number of, respectively occupied and virtual orbitals), if the 4-index integrals \((ac|bd)\) in the MO are precalculated and stored on file before the iterative solution of the coupled-cluster equation, 11.4 and 11.5, and the full permutational symmetries of \(t^{ij}_{cd}\), \((ac|bd)\), and \(\Omega_{aibj}\) are exploited. For larger systems, however, the storage and I/O of the integrals \((ac|bd)\) leads to bottlenecks. Alternatively, this contribution can be evaluated in an integral-direct way as \[\begin{equation} t^{ij}_{\kappa \lambda } = \sum_{cd} t^{ij}_{cd} C_{\kappa c} C_{\lambda d} , \qquad \Omega^B_{\mu i \nu j} = \sum_{\kappa\lambda} t^{ij}_{\kappa \lambda} (\mu\kappa|\nu\lambda) , \qquad \Omega^B_{a i b j} = \sum_{\mu\nu} \Omega^B_{\mu i \nu j} C_{\mu a} C_{\nu b} \end{equation}\](11.23) which, depending on the implementation and system, has formally a 2–3 times larger operation count, but allows to avoid the storage and I/O bottlenecks by processing the 4-index integrals on-the-fly without storing them. Furthermore, integral screening techniques can be applied to reduce the operation count for large systems to an asymptotic scaling with \({\cal O}({\cal N}^4)\).

In TURBOMOLE only the latter algorithm is presently implemented. (For small systems other codes will therefore be faster.)

The other class of expensive contributions are so-called ring terms (in some publications denoted as C and D terms) which involve contractions of the doubles amplitudes \(t_{aibj}\) with several 4-index MO integrals with two occupied and two virtual indices, partially evaluated with \(T_1\)-dependent MO coefficients. For these terms the implementation in TURBOMOLE employs the resolution-of-the-identity (or density-fitting) approximation (with the cbas auxiliary basis set) to reduce the overhead from integral transformation steps. Due this approximation CCSD energies obtained with TURBOMOLE will deviate from those obtained with other coupled-cluster programs by a small RI error. This error is usually in the same order or smaller than the RI error for a RI-MP2 calculation for the same system and basis sets.

The RI approximation is also used to evaluate the 4-index integrals in the MO basis needed for the perturbative triples corrections.

Disc space requirements:

In difference to CC2 and MP2, the CCSD model does no longer allow to avoid the storage of double excitation amplitudes (\(t^{ij}_{ab}\)) and intermediates of with a similar size. Thus, also the disc space requirements for CCSD calculations are larger than for RI-MP2 and RI-CC2 calculations for the same system. For a (closed-shell) CCSD ground state energy calculation the amount of disc space needed can be estimated roughly as \[\begin{equation} N_{disc} \approx \Big( 4 N^3 + (4 +m_{DIIS}) O^2 N^2 \Big)/(128\times 1000) ~\text{MBytes} ~, \end{equation}\](11.24) where \(N\) is the number of basis functions, \(O\) the number of occupied orbitals and \(m_{DIIS}\) the number of vectors used in the DIIS procedure (by default 10, see Sec. 25.2.23 for details).

For (closed-shell) CCSD(T) calculations the required disc space is with \[\begin{equation} N_{disc} \approx \Big( 5 N^3 + 5 O^2 N^2 + ON^3\Big)/(128\times 1000) ~ \text{MBytes} ~, \end{equation}\](11.25) somewhat larger.

For calculations with an open-shell (UHF or ROHF) reference wavefunction the above estimates should be multiplied by factor of 4.

Memory requirements:

The CCSD and CCSD(T) implementation in TURBOMOLE uses multi-pass algorithms to avoid strictly the need to store arrays with a size of \({\cal O}(N^3)\) or \({\cal O}(O^2N^2)\) or larger as complete array in main memory. Therefore, the minimum memory requirements are relatively low—although it is difficult to give accurate estimates for them.

On should, however, be aware that, if the amount of memory provided to the program in the data group $maxcor becomes too small compared to \(O^2N^2/(128\times 1000)\) MBytes, loops will be broken in many small batches at the cost of increased I/O and a decrease in performance. As mentioned above, it is recommended to set $maxcor to 66–75% of the physical core memory available for the calculation.

Important options:

The options to define the orbital and the auxiliary basis sets, the maximum amount of allocatable core memory ($maxcor), and the frozen-core approximation ($maxcor) have been mentioned above and described in the chapters on MP2 and CC2 calculations. Apart from this, CCSD and CCSD(T) calculations require very little additional input.

Relevant are in particular some options in the $ricc2 data group:

$ricc2
  # ccsd
  ccsd(t)
  conv=7
  oconv=6
  mxdiis=10
  maxiter=25

The options ccsd and ccsd(t) request, respectively, CCSD and CCSD(T) calculations. Since CCSD(T) requires the cluster amplitudes from a converged CCSD calculation, the option ccsd(t) includes a CCSD calculation or a restart of the latter if the option ccsd is also set in the input¹. bccd(t) requests the Brueckner coupled cluster variant BCCD(T).

The number given for mxdiis defines the maximum number of vectors included in the DIIS procedure for the solution of the cluster equations. As mentioned above, it has some impact on the amount of disc space used by a CCSD calculation. Unless disc space becomes a bottleneck, it is not recommended to change the default value.

With maxiter one defines the maximum number of iterations for the solution of the cluster equations. If convergence is not reached within this limit, the calculation is stopped. Usually 25 iterations should be sufficient for convergence. Only in difficult cases with strong correlation effects more iterations are needed. It is recommended to increase this limit only if the reason for the strong correlation effects is known. (Since one reason could also be an input error as e.g. unreasonable geometries or orbital occupations or a wrong basis set assignment. With diffuse basis functions it is sometimes neccessary to tighten the integral screening threshold in $scftol.) If oscillatory convergence is observed, adding a level shift can help convergence.

The two parameters conv and oconv define the convergence thresholds for the iterative solution of the cluster equations. Convergence is assumed if the change in the energy (with respect to the previous iteration) is smaller than \(10^{-\mathrm{conv}}\) and the euclidian norm of the residual (the so-called vector function) smaller than \(10^{-\mathrm{oconv}}\). If conv is not given in the data group $ricc2 the threshold for changes in the energy is set to the value given in $denconv (by default \(10^{-6}\)). If oconv is not given in the data group $ricc2 the threshold for the residual norm is by default set to 10 \(\times\) conv. With the default settings for these thresholds, the energy will thus be converged until changes drop below \(10^{-7}\) Hartree, which typically ensures an accuracy of about \(1~\mu\mathrm{H}\). These setting are thus rather tight and conservative even for the calculation of highly accurate reaction energies. If for your application larger uncertainties for the energy are tolerable, it is recommended to use less tight thresholds, e.g. conv=6 or conv=5 for an accuracy of, respectively, at least 0.01 mH (0.03 kJ/mol) or 0.1 mH (0.3 kJ/mol). The settings for conv (and oconv) have not only an impact on the number of iterations for the solution of the cluster equations, but as they determine the thresholds for integral screening also on the costs for the individual iterations.