Coupled-Cluster Class

Coupled-cluster (CC) theory is one of the most accurate and widely used quantum chemistry methods for describing electron correlation in molecular systems. It represents the many-electron wavefunction using an exponential ansatz:

\[|\Psi_{\text{CC}}\rangle = e^{\hat{T}}|\Phi_0\rangle\]

where \(|\Phi_0\rangle\) is a reference determinant (typically Hartree-Fock), and \(\hat{T}\) is a cluster operator that generates excited determinants. The cluster operator is typically expanded as:

\[\hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3 + \ldots\]

where \(\hat{T}_n\) generates \(n\)-fold excited determinants. For computational tractability, the cluster operator is usually truncated. For example, CCSD includes only single and double excitations (\(\hat{T} = \hat{T}_1 + \hat{T}_2\)).

The CC class provides a convenient interface for setting up and processing CC equations using SeQuant’s symbolic algebra engine. It supports various CC formulations, including traditional, unitary, and orbital-optimized ansätze.

Overview

The CC class can be used to derive:

  • Ground state amplitude equations

  • λ (de-excitation) amplitude equations — Lagrange multipliers conjugate to the ground-state amplitudes, needed for properties, analytic gradients and perturbative corrections

  • Equation-of-motion (EOM) CC equations for excited states

  • Response equations for properties and perturbations

Expressions are generated in spin-orbital basis and can be post-processed using SeQuant’s spin-tracing capabilities. See Spin Tracing of Expressions for more details.

Multireference contexts are not fully supported yet: only hbar(), energy() and t() are available, and only with the BCH expansion. Because the amplitude operators do not commute with each other there, the expansion does not terminate and grows exponentially with the commutator rank, so hbar_comm_rank must be set explicitly, as for the unitary ansatz; the constructor throws otherwise.

Ansatz Options

The CC class supports several CC ansätze through the CC::Ansatz enum:

  • Ansatz::T: Traditional CC ansatz, where the wavefunction is represented as \(|\Psi_{\text{CC}}\rangle = e^{\hat{T}}|\Phi_0\rangle\). This is the standard approach used in most implementations.

  • Ansatz::oT: Orbital-optimized traditional ansatz. Singles amplitudes (\(\hat{T}_1\)) are excluded from the cluster operator, with orbital optimization performed instead.

  • Ansatz::U: Unitary CC ansatz, where the wavefunction is represented as \(|\Psi_{\text{UCC}}\rangle = e^{\hat{T} - \hat{T}^\dagger}|\Phi_0\rangle\). Particularly useful for quantum computing applications.

  • Ansatz::oU: Orbital-optimized unitary ansatz, combines both unitary and orbital optimized ansätze.

Key Methods

Ground State Amplitudes

std::vector<ExprPtr> t(size_t commutator_rank = 4,
                       size_t pmax = std::numeric_limits<size_t>::max(),
                       size_t pmin = 0);

Derives the equations for the \(t\) amplitudes (\(\langle \Phi_P|\bar{H}|\Phi_0 \rangle = 0\)) up to specified excitation levels.

std::vector<ExprPtr> λ(size_t commutator_rank = 4);

Derives the equations for the \(\lambda\) de-excitation amplitudes.

Coupled-Cluster Response

std::vector<ExprPtr> tʼ(size_t rank = 1, size_t order = 1);
std::vector<ExprPtr> λʼ(size_t rank = 1, size_t order = 1);

Derives perturbed amplitude equations for response theory calculations.

Equation-of-Motion Coupled-Cluster

std::vector<ExprPtr> eom_r(
    nₚ np, nₕ nh,
    const std::vector<std::size_t>& block_ranks = {});
std::vector<ExprPtr> eom_l(nₚ np, nₕ nh);

Derives equation-of-motion coupled-cluster (EOM-CC) equations for excited states. The eom_r method generates equations for the right eigenvectors, while eom_l generates equations for the left eigenvectors. Traditional CC always uses the connected H̄R product. UCC always uses Hamiltonian-matrix assembly, subtracting the scalar part of each diagonal block’s Hamiltonian.

Ground-state cluster amplitudes, EOM vectors, and projectors follow the active sequant::mbpt::NormalizationConvention throughout CC and UCC equation construction, including H̄ in the BCH and Bernoulli paths. The caller selects the convention before deriving equations and supplies amplitudes in that convention. UCC EOM requires Symmetric and rejects Default: its blocks then form \(\mathbf{D}\mathbf{H}\mathbf{D}\) with \(D_{II}=1/\sqrt{m_I}\) (\(m_I\) as in Operator Interface for Symbolic Many-Body Algebra), which is Hermitian when block_ranks is symmetric. Under Default they would form the non-Hermitian \(\mathbf{H}\mathbf{D}^2\).

For UCC, the optional block_ranks argument gives the per-block truncation orders as a row-major matrix over the EOM manifolds: nested-commutator order for BCH and \(\bar{H}^{k}\) order for Bernoulli. When omitted, every block uses the configured H̄ rank. A positive hbar_singles_comm_rank applies the additional singles transform even to rank-0 BCH blocks. Traditional CC does not support block ranks.

Examples

The following examples demonstrate how to use the CC class to derive CC equations for various ansätze and excitation levels.

From this point onward, assume the following namespaces are imported, and sequant::Context and sequant::mbpt::Context are set up as shown.

using namespace sequant;
using namespace sequant::mbpt;
set_default_context({.index_basis_registry_shared_ptr = make_min_sr_spaces(),
                     .vacuum = Vacuum::SingleProduct});
set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()});

Ground State CC Amplitude Equations

// Traditional CCSD
auto t_eqs = CC{2}.t();
std::wcout << "T1: " << to_latex(t_eqs[1]) << "\n"
           << "T2: " << to_latex(t_eqs[2]) << "\n";

// Lambda equations
auto l_eqs = CC{2}.λ();
std::wcout << "λ1: " << to_latex(l_eqs[1]) << "\n"
           << "λ2: " << to_latex(l_eqs[2]) << "\n";

// Unitary CCSD
auto Ut_eqs = CC(2, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 4})
                  .t();  // Use 4th-order commutator expansion
std::wcout << "T1 (UCC): " << to_latex(Ut_eqs[1]) << "\n"
           << "T2 (UCC): " << to_latex(Ut_eqs[2]) << "\n";

// Orbital-optimized CCSD
auto oT_eqs = CC(2, {.ansatz = CC::Ansatz::oT}).t();
std::wcout << "T2 (oT): " << to_latex(oT_eqs[2]) << "\n";

EOM-CC Equations

// EE-EOM-CCSD (excitation energy)
auto r_eqs = CC{2}.eom_r(nₚ(2), nₕ(2));
std::wcout << "R1: " << to_latex(r_eqs[1]) << "\n"
           << "R2: " << to_latex(r_eqs[2]) << "\n";

//  EE-EOM-CCSD Left eigenvectors
auto ee_l_eqs = CC{2}.eom_l(nₚ(2), nₕ(2));
std::wcout << "L1: " << to_latex(ee_l_eqs[1]) << "\n"
           << "L2: " << to_latex(ee_l_eqs[2]) << "\n";

// IP-EOM-CCSD (ionization potential)
auto ip_eqs = CC{2}.eom_r(nₚ(0), nₕ(1));
std::wcout << "IP-R1: " << to_latex(ip_eqs[0]) << "\n";

// EA-EOM-CCSD (electron attachment)
auto ea_eqs = CC{2}.eom_r(nₚ(1), nₕ(0));
std::wcout << "EA-R1: " << to_latex(ea_eqs[0]) << "\n";

Response and Perturbation Equations

// First-order perturbed amplitude equations
auto t_pt = CC{2}.tʼ(1, 1);
std::wcout << "t1 perturbed: " << to_latex(t_pt[1]) << "\n"
           << "t2 perturbed: " << to_latex(t_pt[2]) << "\n";

// First-order perturbed Lambda amplitude equations
auto l_pt = CC{2}.λʼ(1, 1);
std::wcout << "λ1 perturbed: " << to_latex(l_pt[1]) << "\n"
           << "λ2 perturbed: " << to_latex(l_pt[2]) << "\n";

Advanced Usage

Truncating the Commutator Expansion

The similarity-transformed Hamiltonian is built by mbpt::lst, see Lie Similarity Transformation. For traditional CC with two-body Hamiltonians and the reference equal to the Wick vacuum, the commutator expansion is truncated at 4th order. For unitary CC or a reference that differs from the Wick vacuum the expansion does not terminate, so the commutator rank must be set explicitly; for other Hamiltonians you may need to set it too:

// Unitary CCSDT with 6th-order commutator expansion
auto Ut_ccsdt = CC(3, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 6}).t();

For unitary BCH expansions, Options::hbar_singles_comm_rank applies an additional singles-only similarity transform to H̄ generated by \(\sigma_1 = T_1 - T_1^\dagger\) after the primary expansion. A zero rank disables this transform; a positive rank requires a unitary ansatz with singles amplitudes enabled and is not supported with Bernoulli expansions. Perturbed amplitude equations and RDMs are not supported with this option.

Using \(\bar{H}\) outside the CC class

CC::hbar() is public, but the form of the expression it returns depends on the ansatz and on the reference. For a non-unitary ansatz with the reference equal to the Wick vacuum, each commutator is written as a connected product (see Lie Similarity Transformation), which equals the commutator only once the operators are connected when taking the expectation value; evaluated with empty connectivity it retains disconnected terms. For a unitary ansatz, or when the reference differs from the Wick vacuum (under any vacuum, since then no connectivity makes a connected product reproduce the commutator), H̄ is built from explicit commutators and is self-contained; imposing connectivity on it would drop terms that must survive.

CC::hbar_connections() returns the connectivity that matches the form H̄ was built with, default_op_connections() or empty, and the class uses it (or a superset) for its own equations. Pass it whenever you evaluate CC::hbar() yourself:

// Evaluate H̄ outside the class with the connectivity that matches its form
CC cc{2};
auto hbar = cc.hbar();
auto r2_eq =
    op::ref_av(op::P(nₚ(2)) * hbar, {.connect = cc.hbar_connections()});

The public operator-level and tensor-level ref_av and vac_av functions default to empty connectivity constraints; omitting the options argument is equivalent to passing {}. Alternatively, build H̄ with explicit commutators using mbpt::lst with its default options, which needs no connectivity.

Spin Tracing of Expressions

Equations generated by the CC class are in spin-orbital basis. SeQuant also provides capability to transform these equations into spin-traced forms; see Spin Tracing for the underlying concepts (the Spin quantum number and the general spintrace() function).

Make sure to include the <SeQuant/domain/mbpt/spin.hpp> header to access spin-tracing functions.

// Generate CCSD R2 equation
auto ccsd_eqs = CC{2}.t();
auto t2_eq = t_eqs[2];

// Convert to the closed-shell spin-traced form
auto t2_cs = closed_shell_CC_spintrace(t2_eq);
std::wcout << "Closed-shell spin-traced CCSD-R2: " << to_latex(t2_cs) << "\n";

// Convert to the open-shell spin-traced form
auto t2_os = open_shell_CC_spintrace(
    t2_eq);  // vector of expressions: {ɑɑ, ɑβ, ββ} blocks
std::wcout << "Open-shell spin-traced CCSD-R2:\n";
for (const auto& eq : t2_os) {
  std::wcout << to_latex(eq) << "\n";
}