Bernoulli Expansion based Unitary Coupled-Cluster

[ajay-mk]

Adds the Bernoulli expansion of the UCC similarity-transformed Hamiltonian $\bar{H} = e^{-\sigma} H e^{\sigma}$, $\sigma = T - T^\dagger$, as an alternative to the BCH expansion in mbpt::CC. Ranks $\bar H^0$ – $\bar H^4$, Eqs. (45)–(50) of Liu, Asthana, Cheng & Mukherjee, JCP 148, 244110 (2018)

Companion PR: MPQC4 #792.

Changes

• Adds the Bernoulli-number expansion of the UCC similarity-transformed Hamiltonian $\bar H = e^{-\sigma}He^{\sigma}$ as an alternative to BCH in mbpt::CC, plus per-block $\bar H$ truncation for block-truncated EOM (qUCCSD and its IP/EA analogues).
• CC::Options::hbar_expansion: Selects BCH or Bernoulli; the former is the default. Bernoulli requires a
  unitary ansatz and hbar_comm_rank.
• New Bernoulli expansion related logic:
  • bernoulli::hbar(N, rank, skip1) assembles $\bar H^0$ – $\bar H^4$ (Eqs. (45)–(50) of Ref. 1) rank by rank. $H$ splits as $F+V$, and every operator splits into its pure excitation/de-excitation part ($O_N$) and the rest ($O_R$).
  • detail::wick_reduce — Wick-reduces while keeping partial contractions
    (needed so nesting produces operators, not scalars). Must run with use_topology(false): its weight bookkeeping is only valid on the fully-contracted path, and it defaults to on.
  • detail::N_part / detail::R_part — the $O_N$ / $O_R$ split. Only $N$ is block-resolved; $R$ stays in compact general-index form, which keeps the nested commutators that consume it cheap.
  • bernoulli::hbar() returns a tensor-level expression (BCH's hbar() returns an operator-level one). Callers have to take care of this.
• CC::eom_r: per-block H̄ truncation (block_ranks):
  • New optional argument: a row-major $K\times K$ matrix, one commutator truncation rank per block of the EOM secular matrix instead of one uniform $\bar H$.
  • Manifolds are ordered by ascending rank, so the same matrix serves EE, IP and EA. qUCCSD is {2,1,1,0}. See Liu & Cheng, JCP 155, 174102 (2021), Sec. II C; Zhang & Liu, JCTC 22, 3341 (2026), Table 1, for the IP/EA analogues.
  • Each block is the sandwich ⟨i|H¯|j⟩−δijE, not the commutator form. Under Bernoulli, each block's H¯ has its N part (the ground-state amplitude residual) removed — it only vanishes at the block's own rank when that equals the amplitude rank, so keeping it would leak a spurious residual into the off-diagonal blocks.

Tests

• tests/unit/test_mbpt_cc.cpp: Wick reduction, N/R split, $\bar H$ structure, config validation, block-truncated EOM term counts and the three rejection cases.
• tests/integration/ucc.cpp: BCH/Bernoulli term-count pins at ranks 2–3.
• Implementation is validated numerically against literature, see MPQC4 #792 for details.

qUCCSD Example

How to derive ground and excited state equations

#include <SeQuant/core/context.hpp>
#include <SeQuant/core/expr.hpp>
#include <SeQuant/core/op.hpp>
#include <SeQuant/core/tensor_canonicalizer.hpp>
#include <SeQuant/domain/mbpt/context.hpp>
#include <SeQuant/domain/mbpt/convention.hpp>
#include <SeQuant/domain/mbpt/models/cc.hpp>
#include <SeQuant/domain/mbpt/op.hpp>

#include <iostream>

using namespace sequant;
using namespace sequant::mbpt;

int main() {
  set_locale();
  set_default_context({.index_space_registry_shared_ptr = make_sr_spaces(),
                       .vacuum = Vacuum::SingleProduct,
                       .metric = IndexSpaceMetric::Unit,
                       .spbasis = SPBasis::Spinor});
  TensorCanonicalizer::set_cardinal_tensor_labels(cardinal_tensor_labels());
  set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()});

  const CC cc(2, {.ansatz = CC::Ansatz::U,
                  .hbar_comm_rank = 2,
                  .hbar_expansion = CC::HbarExpansion::Bernoulli});


  // qUCCSD: 3 commutators in Energy, 2 in ground state amplitudes; [2,1,1,0] for excited states

  // energy
  const auto energy = cc.energy(3);
  std::wcout << "E        : " << energy->size() << " terms\n";

  // t amplitudes
  const auto t = cc.t();
  for (std::size_t p = 1; p < t.size(); ++p)
    std::wcout << "residual" << p << ": " << t[p]->size() << " terms\n";

  // EE
  const auto sigma = cc.eom_r(nₚ(2), nₕ(2), {2, 1, 1, 0});
  for (std::size_t p = 1; p < sigma.size(); ++p)
    std::wcout << "sigma" << p << "   : " << sigma[p]->size() << " terms\n";

  // IP/EA use the same block ranks, with R of unequal particle/hole rank
  std::wcout << "IP sigma1: " << cc.eom_r(nₚ(1), nₕ(2), {2, 1, 1, 0})[0]->size()
             << " terms\n";
}

Important References:

• Liu, Asthana, Cheng & Mukherjee, JCP 148, 244110 (2018): Bernoulli expansion and qUCCSD equations
• J. Liu, D. A. Matthews, and L. Cheng, JCTC 18, 2281 (2022): Efficient qUCCSD implementation and higher-order methods
• Liu & Cheng, JCP 155, 174102 (2021): Reference EE results
• Zhang & Liu, JCTC 22, 3341 (2026): Reference IP/EA results