The standard approach for treating pair correlations self-consistently are
coupled-cluster methods where the cluster operator contains all single and
double substitutions,^{Purvis:1982} abbreviated as CCSD. CCSD yields
results that are only slightly superior to MP2 for structures and frequencies
of stable closed-shell molecules. However, it is far superior for reactive
species, such as transition structures and radicals, for which the performance
of MP2 is quite erratic.

A full textbook presentation of CCSD is beyond the scope of this manual, and several comprehensive references are available. However, it may be useful to briefly summarize the main equations. The CCSD wave function is:

$$|{\mathrm{\Psi}}_{\mathrm{CCSD}}\u27e9=\mathrm{exp}\left({\widehat{T}}_{1}+{\widehat{T}}_{2}\right)|{\mathrm{\Phi}}_{0}\u27e9$$ | (6.25) |

where the single and double excitation operators may be defined by their actions on the reference single determinant (which is normally taken as the Hartree-Fock determinant in CCSD):

$${\widehat{T}}_{1}|{\mathrm{\Phi}}_{0}\u27e9=\sum _{i}^{\mathrm{occ}}\sum _{a}^{\mathrm{virt}}{t}_{i}^{a}|{\mathrm{\Phi}}_{i}^{a}\u27e9$$ | (6.26) |

$${\widehat{T}}_{2}|{\mathrm{\Phi}}_{0}\u27e9=\frac{1}{4}\sum _{ij}^{\mathrm{occ}}\sum _{ab}^{\mathrm{virt}}{t}_{ij}^{ab}|{\mathrm{\Phi}}_{ij}^{ab}\u27e9$$ | (6.27) |

It is not feasible to determine the CCSD energy by variational minimization of ${\u27e8E\u27e9}_{\mathrm{CCSD}}$ with respect to the singles and doubles amplitudes because the expressions terminate at the same level of complexity as full configuration interaction (!). So, instead, the Schrödinger equation is satisfied in the subspace spanned by the reference determinant, all single substitutions, and all double substitutions. Projection with these functions and integration over all space provides sufficient equations to determine the energy, the singles and doubles amplitudes as the solutions of sets of nonlinear equations. These equations may be symbolically written as follows:

${E}_{\mathrm{CCSD}}$ | $=$ | $\u27e8{\mathrm{\Phi}}_{0}|\widehat{H}|{\mathrm{\Psi}}_{\mathrm{CCSD}}\u27e9$ | (6.28) | ||

$=$ | ${\u27e8{\mathrm{\Phi}}_{0}\left|\widehat{H}\right|\left(1+{\widehat{T}}_{1}+{\displaystyle \frac{1}{2}}{\widehat{T}}_{1}^{2}+{\widehat{T}}_{2}\right){\mathrm{\Phi}}_{0}\u27e9}_{C}$ | ||||

$0$ | $=$ | $\u27e8{\mathrm{\Phi}}_{i}^{a}\left|\widehat{H}-{E}_{\mathrm{CCSD}}\right|{\mathrm{\Psi}}_{\mathrm{CCSD}}\u27e9$ | (6.29) | ||

$=$ | ${\u27e8{\mathrm{\Phi}}_{i}^{a}\left|\widehat{H}\right|\left(1+{\widehat{T}}_{1}+{\displaystyle \frac{1}{2}}{\widehat{T}}_{1}^{2}+{\widehat{T}}_{2}+{\widehat{T}}_{1}{\widehat{T}}_{2}+{\displaystyle \frac{1}{3!}}{\widehat{T}}_{1}^{3}\right){\mathrm{\Phi}}_{0}\u27e9}_{C}$ | ||||

$0$ | $=$ | $\u27e8{\mathrm{\Phi}}_{ij}^{ab}\left|\widehat{H}-{E}_{\mathrm{CCSD}}\right|{\mathrm{\Psi}}_{\mathrm{CCSD}}\u27e9$ | (6.30) | ||

$=$ | $\u27e8{\mathrm{\Phi}}_{ij}^{ab}\left|\widehat{H}\right|(1+{\widehat{T}}_{1}+{\displaystyle \frac{1}{2}}{\widehat{T}}_{1}^{2}+{\widehat{T}}_{2}+{\widehat{T}}_{1}{\widehat{T}}_{2}+{\displaystyle \frac{1}{3!}}{\widehat{T}}_{1}^{3}$ | ||||

$\mathrm{\hspace{1em}}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+{\displaystyle \frac{1}{2}}{\widehat{T}}_{2}^{2}+{\displaystyle \frac{1}{2}}{\widehat{T}}_{1}^{2}{\widehat{T}}_{2}+{\displaystyle \frac{1}{4!}}{\widehat{T}}_{1}^{4}){\mathrm{\Phi}}_{0}\u27e9{}_{C}$ |

The result is a set of equations which yield an energy that is not necessarily
variational (*i.e.*, may not be above the true energy), although it is strictly
size-consistent. The equations are also exact for a pair of electrons, and, to
the extent that molecules are a collection of interacting electron pairs, this
is the basis for expecting that CCSD results will be of useful accuracy.

The computational effort necessary to solve the CCSD equations can be shown to scale with the 6th power of the molecular size, for fixed choice of basis set. Disk storage scales with the 4th power of molecular size, and involves a number of sets of doubles amplitudes, as well as two-electron integrals in the molecular orbital basis. Therefore the improved accuracy relative to MP2 theory comes at a steep computational cost. Given these scalings it is relatively straightforward to estimate the feasibility (or non feasibility) of a CCSD calculation on a larger molecule (or with a larger basis set) given that a smaller trial calculation is first performed. Q-Chem supports both energies and analytic gradients for CCSD for RHF and UHF references (including frozen-core). For ROHF, only energies and unrelaxed properties are available.