Is coupled cluster variational for two electrons?

Question

I know that coupled cluster (CC) is not variational for the general case. However, if we only have two electrons with one nucleus, CCSD should be exact for this system like full configuration interaction (CI). Since full CI is variational, can we argue that CCSD is also variational for this special case?

Answer 1

No. The reason for this is not to be found in the excitations, but in the evaluation of the method, i.e. the working equations. $$%Introducing some shortcuts \require{cancel} \newcommand{\op}[1]{\mathrm{#1}} %\op{H} \newcommand{\bracket}[2]{\left\langle#1\middle|#2\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\order}[1]{^{(#1)}} %E_n\order{1} \newcommand{\overlap}[3]{\mathcal{#1}_{#2}\order{#3}} %\overlap{S}{m}{1} \newcommand{\integral}[3]{\mathcal{#1}_{#2,#3}} %\integral{V}{i}{j} $$

The variational principle

In short: a trial wave function always has a higher energy than the exact wave function.

Assume we know the exact solutions to the Schrödinger equation, which form a complete basis and there is one wave function $\Psi_0$ which leads to the lowest energy $E_0$. \begin{align} \op{H}\Psi_i &= E_i\Psi_i & \{i \in \mathbb{N}, i\geq 0\}\tag{1}\label{schrodinger} \end{align}

We choose the solutions to be orthonormal (boundary conditions). $$\bracket{\Psi_i}{\Psi_j} = \delta_{ij} = \begin{cases}1, & i=j \\ 0, & i\neq j \end{cases} \tag2\label{orthonorm}$$

The trial wave function $\Phi$ can be expressed as a linear combination of the complete set of wave functions. $$\Phi = \sum_i^\infty a_i\Psi_i\tag3\label{expansion}$$

The expectation energy of the trial wave function is given through $$W = \frac{\bra{\Phi}\op{H}\ket{\Phi}}{\bracket{\Phi}{\Phi}}. \tag4\label{expen}$$

With $\eqref{expansion}$, $\eqref{orthonorm}$, and $\eqref{schrodinger}$ we can expand and simplify. \begin{align} W &= \frac{ \sum_i^\infty\sum_j^\infty a_i a_j \bra{\Psi_i}\op{H}\ket{\Psi_j} }{ \sum_i^\infty\sum_j^\infty a_i a_j \bracket{\Psi_i}{\Psi_j} }\\ &= \frac{ \sum_i^\infty\sum_j^\infty a_i a_j E_j\bracket{\Psi_i}{\Psi_j} }{ \sum_i^\infty\sum_j^\infty a_i a_j \bracket{\Psi_i}{\Psi_j} }\\ &= \frac{ \sum_i^\infty a_i^2 E_i\bracket{\Psi_i}{\Psi_j} + \sum_{j\neq i}^\infty a_i a_j E_j\bracket{\Psi_i}{\Psi_j} }{ \sum_i^\infty a_i^2 \bracket{\Psi_i}{\Psi_j} + \sum_{j\neq i}^\infty a_i a_j \bracket{\Psi_i}{\Psi_j} }\\ &= \frac{ \sum_i^\infty a_i^2 E_i }{ \sum_i^\infty a_i^2 } & | -E_0\\ W - E_0 &= \frac{ \sum_i^\infty a_i^2 (E_i - E_0) }{ \sum_i^\infty a_i^2 } \\ W - E_0 &= \frac{ \sum_{i>0}^\infty a_i^2 (E_i - E_0) }{ \sum_{i>0}^\infty a_i^2 } + \frac{a_0^2 (E_0 - E_0)}{a_0^2}\tag{4'}\label{notnegative}\\ W - E_0 &\geq 0\tag5\label{varprin} \end{align}

In $\eqref{notnegative}$ the last term is zero. Since $E_i > E_0$ (by definition) and $a_i^2$ is always greater or equal zero, the expectation value of the trial wave function must always have a higher energy than the exact solution of the Schrödinger equation, i.e. $\eqref{varprin}$.

Configuration Interaction

The CI wave function is set up as a linear combination of determinants from a reference calculation. This reference calculation can in principle be any other method the only requirement is that it also forms a complete (orthonormal) basis.
When Hartree-Fock is chosen, we know that we obtain a complete set and that this method itself is variational. The CI trial wave function can then be expressed as a linear combination of solutions of the HF wave function. $$\Psi^\mathrm{CI} = \sum_i^\infty a_i\Phi_i = a_0\Phi^\mathrm{HF} + \sum_{i>0}^\infty a_i\Phi_i$$ An alternative expression of the above is using an "excitation operator" $\op{T}$ to generate the higher energy determinants from the reference. $$\Psi^\mathrm{CI} = (1 + \op{T})\Phi^\mathrm{HF}$$

The energy of the CI trial wave function can be evaluated with the Lagrange multiplier method, which requires, that the wave function is normalised (adding zero). $$L = \bra{\Psi^\mathrm{CI}}\op{H}\ket{\Psi^\mathrm{CI}} - \lambda\left[\bracket{\Psi^\mathrm{CI}}{\Psi^\mathrm{CI}}-1\right]$$

Since the boundary conditions did not change, the variational principle still holds.

Coupled Cluster

Because of the inherent complexity of the CC ansatz, one key element to evaluate to CC energy is missing: the requirement of the normalised CC trial wave function.
We generate the trial wave function with an exponential approach (rather than a linear) expanded into a Taylor series. \begin{align} \Psi^\mathrm{CC} &= \mathrm{e}^{\op{T}}\Phi^\mathrm{HF}\\ &= \left[\op{1} + \op{T}_1 + \left(\op{T}_2 + \frac12\op{T}_1^2\right) + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right) + \cdots \right] \Phi^\mathrm{HF}\\ &= \Phi^\mathrm{HF} + \left[\op{T}_1 + \left(\op{T}_2 + \frac12\op{T}_1^2\right) + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right) + \cdots \right] \Phi^\mathrm{HF}\\ \end{align}

Or to shorten the journey $$\Psi^\mathrm{CC} = \Phi^\mathrm{HF} + \sum_i^\infty t_i \Phi_i.$$

One can easily see what kind of mess it would become requiring the wave function to be normalised $\bracket{\Psi^\mathrm{CC}}{\Psi^\mathrm{CC}}=1$, so the Lagrange ansatz won't work. Instead we choose the wave function to be orthonormal to the reference, which in the Hartree-Fock (complete orthonormal basis) takes away a lot of work. \begin{align} \bracket{\Phi^\mathrm{HF}}{\Psi^\mathrm{CC}}&= 1& \implies\quad \bracket{\Phi^\mathrm{HF}}{\Phi_i} &= 0 \end{align}

We can write the CC-Schrödinger equation as $$\op{H}\mathrm{e}^{\op{T}}\Phi^\mathrm{HF} = E^\mathrm{CC}\mathrm{e}^{\op{T}}\Phi^\mathrm{HF}$$ and evaluate the energy as \begin{align} \bra{\Phi^\mathrm{HF}}\op{H}\mathrm{e}^{\op{T}}\ket{\Phi^\mathrm{HF}} &= E^\mathrm{CC} \bracket{\Phi^\mathrm{HF}}{\mathrm{e}^{\op{T}}\Phi^\mathrm{HF}}\\ E^\mathrm{CC} &= \bra{\Phi^\mathrm{HF}}\op{H}\mathrm{e}^{\op{T}}\ket{\Phi^\mathrm{HF}}\\ E^\mathrm{CC} &= E^\mathrm{HF} + \bra{\Phi^\mathrm{HF}}\op{H} \ket{\left[\op{T}_1 + \left(\op{T}_2 + \frac12\op{T}_1^2\right) + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right) + \cdots \right]\Phi^\mathrm{HF}}.\\ \end{align}

From this we can see that the evaluation of the expectation value of the energy is significantly different from what is required for the variational principle $\eqref{expen}$.

Without going further into details, the dropping of the requirement of a normalised trial wave function leads to not obeying the variational principle.

Additional considerations

It is often said that Full-CC is equivalent to Full-CI, which is a convenient half truth. It is extremely important that this gets more true with a better description of the reference wave function, i.e. with larger basis sets. Full-CI is exact at the complete basis set limit (for the electronic Schrödinger equation). This is not necessarily true for Full-CC.

Since we are in computational chemistry hardly ever concerned about absolute energies, the non-variational nature of the method is by far outperformed by the "infinite" approach generating the "excited" determinants. Relative energies don't suffer significantly (compared with the size consistency and extensivity) under the non-variational approach.

Answer 2

My previous answer was downvoted (-3!), but I maintain that in the particular case where there are only two electrons, CCSD is equivalent to CISD and Full-CI. If the discussion was about asking if for a one-electron system Coupled Cluster gives an upper bound to the energy, everyody would agree that it is true. I agree with all that was said in the top answers for the general case, but with two electrons things are a bit special since there are no triples and quadruples.

1) The CCSD equations can be rewritten in such a way to be solved with an intermediate Hamiltonian in the space of singly and doubly excited determinants.[1] In the special case where there are 2 electrons the space is the FCI space, and the equations of FCI, CISD and CCSD are strictly equivalent : the dressing matrix $\Delta$ is zero so the effective Hamiltonian is equal to the Hamiltonian.

2) I don't agree with the argument that the fact that the energy is computed by projection implies that the energy is non-variational. What makes the result usually approximate is the fact that the wave function is not an eigenfunction of the Hamiltonian. But when it is, computing the energy by projection is exact. In the special case where we have 2 electrons, the CISD is the Full-CI. So if one projects the equations on the singles and doubles, the Full space is considered and the FCI wave function is obtained (up to a normalization factor). In the FCI, the wave function is an eigenfunction of the Hamiltonian, so computing the energy by projection gives exactly the FCI energy (in a less numerically stable manner in practice).

3) $\Psi_{\rm CCSD} = \exp^{T_1+T_2} \Phi_{\rm HF} = (1 + (T_1 + T_2) + \frac{1}{2}(T_1+T_2)^2 + \dots) \Phi_{\rm HF}$

With only 2 electrons,

$\Psi_{\rm CCSD} = (1 + T_1 + T_2 + \frac{1}{2}T_1^2) \Phi_{\rm HF}$

and $\Psi_{\rm FCI} = (1+T_1+T_2)$

Both $T_1^2\Phi_{\rm HF}$ and $T_2 \Phi_{\rm HF}$ span the same determinant space, so finding the $\Psi_{\rm CCSD}$ and the $\Psi_{\rm FCI}$ which minimize the energy will give the same solution in the basis of determinants. All the additional constraints of Coupled Cluster are on the amplitudes of the Triples and Quadruples, which are absent here.

---

[1] Eigenvalue problem formulation of coupled-cluster expansions through intermediate Hamiltonians, I Nebot-Gil, J Sanchez-Marin, J.L Heully, J.P Malrieu, D Maynau, Chemical Physics Letters, Volume 234, Issues 1–3, 3 March 1995, Pages 45-49, https://doi.org/10.1016/0009-2614(95)00026-Z

Note: In the paper the energies they obtain are slightly different from standard CCSD, but this is due to an approximation on the $T_1\times T_2$ terms, which don't exist with only two electrons.

Answer 3

In the most general case, CC can be understood simply as a prescription for a trial wave function (ansatz) that uses the exponential excitation operators. This ansatz can be then optimized variationally, and this is variational CC.^[1] It is much more expensive than the common CC that uses another approximation (besides the ansatz), which linearizes the CC amplitude equations in a certain way. In most cases, this is a good approximation.

"Variational" means that a method has some ansatz for the wave function which is optimized (varied) to get the minimal energy. This is then guaranteed to be greater or equal to the true energy by the variational principle. In this sense, the common CCSD is not variational even for two-electron systems: there is no energy that would be optimized in the procedure.

Perhaps you could understand "variational" as a method that is in some way guaranteed to always give an energy greater or equal to the true energy. In general, the non-variational CC is not such a method. But for two-electron systems, the CCSD energy is exact, so technically, it is variational in this other derived sense. But as it is exact, which is a stronger statement, I see little sense in stressing the "variational" aspect. Also, I know of no method that would be variational in this second sense, but not in the first sense.

Reference

1. Troy Van Voorhis and Martin Head-Gordon. J. Chem. Phys. 2000, 113, 8873.

Answer 4

Yes, in the particular case of 2 electrons, CCSD is strictly equivalent to CISD since there can't be any excitation higher that singles and doubles. CISD is the Full-CI of a 2-electron system, so CCSD is variational for 2 electrons.

Answer 5

Here is a partial bibliography of papers that refer to CCSD as exact for two electrons, meaning equivalent to FCI, or CISD, for the same basis set.

• https://doi.org/10.1016/0009-2614(93)85124-7: "The QCISD and CCSD transition energies are size consistent, and exact for two-electron systems"
• https://doi.org/10.1063/1.5128795: "CCSD is exact for 2-electron systems."
• https://doi.org/10.1063/1.5094788: "An important feature of these methods is that they are exact for systems with a number of electrons smaller or equal to the highest excitation level (i.e., CCSD is exact for two-electron systems)."
• https://doi.org/10.1021/acs.jctc.9b00674: "In contrast, it is clear that CCSD (despite being exact for an ensemble of isolated two electron pairs)..."
• https://doi.org/10.1016/0009-2614(93)E1398-Z: "This CPMET-C based excited state method is size-consistent, but unlike CCSD and QCISD, is not exact for separated pairs of electrons."
• https://doi.org/10.3389/fmats.2019.00123: "For two-electron systems, 2-fold excited Slater determinants can be generated at most and therefore CCSD theory becomes exact. Higher orders of coupled cluster theory become exact for systems with the corresponding number of electrons."

It is left as an exercise to the reader if a method that is equivalent to FCI is variational.

Answer 6

Proposition. Consider a quantum chemical system with $n$ electrons in $N$ spin orbitals. For any FCI wave function

$$| \Psi^{\mathrm{CI}} \rangle = \left( 1 + \sum_{k=1}^n T_k \right) |\Phi^\mathrm{HF}\rangle $$ there exists an excitation operator $\tilde{T}=\sum_k^n \tilde{T}_k$ such that $$| \Psi^{\mathrm{CI}} \rangle = e^\tilde{T} |\Phi^\mathrm{HF}\rangle. $$

Therefore, for an $n$-electron system, the range of configurations accessible to CC-n is the entire space of $n$-electron Slater determinants.

Proof.

The operator $T=\sum_{k=1}^n T_k$ is nilpotent of degree $\le n$ and the matrix logarithm of $1 + T$ therefore exists. It is given by the Mercator series [1]

$$ \log(1+T) = T - \frac{T^2}{2} +\frac{T^3}{3}-\frac{T^4}{4}+\cdots+(-1)^{n+1}\frac{T^n}{n}.$$

Define $ \tilde{T} \equiv \log(1+T) $.

Because $\tilde{T}$ is a sum of excitation operators, it is itself an excitation operator. Then verify

\begin{align} e^\tilde{T} |\Phi^\mathrm{HF}\rangle &= e^{\log(1+T)} |\Phi^\mathrm{HF}\rangle \\ &= (1+T) |\Phi^\mathrm{HF}\rangle \\ &= |\Psi^\mathrm{CI}\rangle. \end{align}

Because the $k$th term in the series for $\tilde{T}$ is composed of $k$-excitations and higher, this produces the following formulae for the $\tilde{T}_k$:

\begin{align} \tilde{T}_1 &= T_1 \\ \tilde{T}_2 &= T_2 - \frac{1}{2}T_1^2 \\ \tilde{T}_3 &= T_3 - \frac{1}{2}(2 T_1 T_2) + \frac{1}{3} T_1^3 \\ \tilde{T}_4 &= T_4 - \frac{1}{2}(T_2^2 + 2 T_1 T_3) + \frac{1}{3}(3T_1^2 T_2) - \frac{1}{4} T_1^4 \\ &\mathrm{etc.} \end{align}

[1] https://math.stackexchange.com/questions/3555422/logarithm-of-i-n-n-for-n-nilpotent