  In single-reference coupled-cluster (CC) methods, one has to solve a set of
non-linear polynomial equations in order to determine the so-called amplitudes
which are then used to compute the energy and other properties. Although it is
of common practice to converge to the (lowest-energy) ground-state solution, it
is also possible, thanks to tailored algorithms, to access higher-energy roots
of these equations which may or may not correspond to genuine excited states.
Here, we explore the structure of the energy landscape of variational CC (VCC)
and we compare it with its (projected) traditional version (TCC) in the case
where the excitation operator is restricted to paired double excitations
(pCCD). By investigating two model systems (the symmetric stretching of the
linear \ce{H4} molecule and the continuous deformation of the square \ce{H4}
molecule into a rectangular arrangement) in the presence of weak and strong
correlations, the performance of VpCCD and TpCCD are gauged against their
configuration interaction (CI) equivalent, known as doubly-occupied CI (DOCI),
for reference Slater determinants made of ground- or excited-state Hartree-Fock
orbitals or state-specific orbitals optimized directly at the VpCCD level. The
influence of spatial symmetry breaking is also investigated.
