Bernoulli Expansion based Unitary Coupled-Cluster

.gitignore

@@ -71,4 +71,4 @@ _codeql_detected_source_root
 .clangd
 .vscode
 out
-run
+run

CMakeLists.txt

@@ -399,6 +399,8 @@ set(SeQuant_symb_src
 set(SeQuant_mbpt_src
         SeQuant/domain/mbpt/antisymmetrizer.cpp
         SeQuant/domain/mbpt/antisymmetrizer.hpp
+        SeQuant/domain/mbpt/bernoulli.cpp
+        SeQuant/domain/mbpt/bernoulli.hpp
         SeQuant/domain/mbpt/biorthogonalization.cpp
         SeQuant/domain/mbpt/biorthogonalization.hpp
         SeQuant/domain/mbpt/context.cpp

SeQuant/domain/mbpt/bernoulli.hpp

@@ -0,0 +1,78 @@
+#ifndef SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP
+#define SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP
+
+#include <SeQuant/core/expr.hpp>
+#include <cstddef>
+
+namespace sequant::mbpt::bernoulli {
+
+/// Tensor-level H̄ = Σ_{k=0..rank} H̄^k in the Bernoulli expansion, for
+/// σ = T−T† of rank N.
+///
+/// The Bernoulli expansion rewrites the non-terminating UCC
+/// similarity-transform series so that Bernoulli numbers appear as the
+/// expansion coefficients; the rank-by-rank operators H̄⁰..H̄⁴ are Eqs. (46)-(50)
+/// of 10.1063/1.5030344.
+///
+/// @warning Single-reference only. The N/R split expands general indices over
+/// the hole and particle spaces alone (see detail::expand_to_blocks), dropping
+/// any other base space the registry defines. That is harmless only because the
+/// single-reference projection manifolds annihilate the dropped terms. Under a
+/// multireference registry they contribute, and both the N and the R part come
+/// out wrong. Nothing checks for this.
+///
+/// The result is a tensor-level expression: coefficient tensors times
+/// normal-ordered operators, not `mbpt::op` operators. Nothing is screened out
+/// of it, so the caller projects every term.
+///
+/// @pre an HF reference: F is taken to have no occupied-virtual block, which
+/// is what keeps F out of H̄² and higher. The f_ov terms of H̄⁰ and H̄¹ are
+/// carried symbolically and vanish only on substitution.
+///
+/// @param N cluster/excitation rank (also the N/R rank cutoff)
+/// @param rank highest Bernoulli order H̄^k to include (0..4)
+/// @param skip1 exclude singles from T
+/// @throw Exception if @p rank > 4
+ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1);
+
+namespace detail {
+
+/// Applies Wick's theorem to @p expr retaining PARTIAL contractions,
+/// reducing a product of normal-ordered operators to a sum of normal-ordered
+/// operators (each = coefficient tensor × at most one residual NormalOperator;
+/// fully-contracted terms carry none). Unlike the expectation-value path it
+/// keeps operators rather than collapsing to a scalar VEV.
+/// @note @p expr is left untouched; the reduction runs on a clone.
+ExprPtr wick_reduce(const ExprPtr& expr);
+
+/// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A). NOT the bare
+/// algebraic commutator: the operator product is Wick-reduced, so contractions
+/// between A and B generate the lower-rank terms the Bernoulli expansion relies
+/// on. B's summed indices are reindexed to fresh temporaries first, making them
+/// disjoint from A's.
+ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B);
+
+/// Rewrites every general (non-base) index of the residual NormalOperator as
+/// the sum over the hole/particle base spaces it spans (occupied/virtual). The
+/// registry's other base spaces are dropped, which changes no projected
+/// quantity: the single-reference manifolds annihilate the dropped terms (see
+/// the @warning on hbar). After expansion every residual index is definite so
+/// the N/R classifier can act on it. Idempotent on block-resolved input.
+/// @pre the registry specifies both a hole and a particle space
+ExprPtr expand_to_blocks(const ExprPtr& expr);
+
+/// Block-resolved N part (O_N of 10.1063/1.5030344): the terms whose single
+/// residual NormalOperator is a pure excitation or pure de-excitation of rank ≤
+/// @p cutoff. Applies expand_to_blocks first.
+ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff);
+
+/// R part (O_R of 10.1063/1.5030344: expr minus its N part). Unlike N_part
+/// the result is NOT block-resolved; it stays in compact general-index form,
+/// which is much cheaper for the nested commutators that consume R.
+ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff);
+
+}  // namespace detail
+
+}  // namespace sequant::mbpt::bernoulli
+
+#endif  // SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP

SeQuant/domain/mbpt/models/cc.cpp

@@ -1,8 +1,10 @@
+#include <SeQuant/core/container.hpp>
 #include <SeQuant/core/expr.hpp>
 #include <SeQuant/core/rational.hpp>
 #include <SeQuant/core/reserved.hpp>
 #include <SeQuant/core/runtime.hpp>
 #include <SeQuant/core/utility/macros.hpp>
+#include <SeQuant/domain/mbpt/bernoulli.hpp>
 #include <SeQuant/domain/mbpt/context.hpp>
 #include <SeQuant/domain/mbpt/convention.hpp>
 #include <SeQuant/domain/mbpt/models/cc.hpp>
@@ -16,6 +18,7 @@
 #include <new>
 #include <stdexcept>
 #include <utility>
+#include <vector>
 
 namespace {
 // alias reserved labels for readability
@@ -42,14 +45,21 @@ CC::CC(size_t n, const Options& opts)
       screen_(opts.screen),
       use_topology_(opts.use_topology),
       hbar_comm_rank_(opts.hbar_comm_rank),
-      pertbar_comm_rank_(opts.pertbar_comm_rank) {
+      pertbar_comm_rank_(opts.pertbar_comm_rank),
+      hbar_expansion_(opts.hbar_expansion) {
   if (unitary())
     SEQUANT_ASSERT(hbar_comm_rank_,
                    "CC: hbar_comm_rank is required for unitary ansatz");
   if (ansatz_ == Ansatz::oT || ansatz_ == Ansatz::oU)
     SEQUANT_ASSERT(skip_singles_,
                    "CC: skip_singles must be true for orbital-optimized "
                    "ansatz");
+  if (hbar_expansion_ == HbarExpansion::Bernoulli) {
+    SEQUANT_ASSERT(ansatz_ == Ansatz::U,
+                   "CC: Bernoulli expansion requires the U ansatz");
+    SEQUANT_ASSERT(hbar_comm_rank_,
+                   "CC: Bernoulli expansion requires hbar_comm_rank");
+  }
 }
 
 CC::Ansatz CC::ansatz() const { return ansatz_; }
@@ -60,6 +70,8 @@ bool CC::unitary() const {
 
 std::optional<size_t> CC::hbar_comm_rank() const { return hbar_comm_rank_; }
 
+CC::HbarExpansion CC::hbar_expansion() const { return hbar_expansion_; }
+
 bool CC::skip_singles() const { return skip_singles_; }
 
 bool CC::screen() const { return screen_; }
@@ -69,13 +81,22 @@ bool CC::use_topology() const { return use_topology_; }
 ExprPtr CC::hbar(std::optional<size_t> truncation_rank) const {
   const auto truncation = truncation_rank.value_or(hbar_comm_rank_.value_or(4));
 
+  if (hbar_expansion_ == HbarExpansion::Bernoulli)
+    return bernoulli::hbar(N, truncation, skip_singles());
+
   // for a non-unitary ansatz this is the cheaper connected-product form, which
   // is only equivalent to the commutator once the caller supplies operator
   // connectivity to ref_av (see lst_options() and the @warning on hbar())
   return mbpt::lst(H(), T(N, skip_singles()), truncation, lst_options());
 }
 
 ExprPtr CC::energy(std::optional<size_t> comm_rank) const {
+  // Bernoulli: the hbar expansion is at tensor level, call the tensor level
+  // ref_av directly. No connectivity or screening.
+  if (hbar_expansion_ == HbarExpansion::Bernoulli) {
+    const auto erank = comm_rank.value_or(hbar_comm_rank_.value());
+    return op::tensor::ref_av(this->hbar(erank));
+  }
   // <0|H̄|0>: reference expectation value of H̄ at the requested commutator
   // truncation. No projector ⇒ this is the energy. ref_av applies the
   // connectivity (empty for unitary, default otherwise).
@@ -86,7 +107,19 @@ ExprPtr CC::energy(std::optional<size_t> comm_rank) const {
 
 std::vector<ExprPtr> CC::t(size_t pmax, size_t pmin) const {
   pmax = (pmax == std::numeric_limits<size_t>::max() ? N : pmax);
-  SEQUANT_ASSERT(pmax >= pmin && "pmax should be >= pmin");
+  SEQUANT_ASSERT(pmax >= pmin, "pmax should be >= pmin");
+
+  // Bernoulli: the hbar expansion is at tensor level, project and call the
+  // tensor level ref_av directly.
+  if (hbar_expansion_ == HbarExpansion::Bernoulli) {
+    const auto hbar = this->hbar();
+    std::vector<ExprPtr> result(pmax + 1);
+    for (std::int64_t p = pmax; p >= static_cast<std::int64_t>(pmin); --p) {
+      const auto projected = (p != 0) ? op::tensor::P(nₚ(p)) * hbar : hbar;
+      result.at(p) = op::tensor::ref_av(projected);
+    }
+    return result;
+  }
 
   // 1. construct hbar(op) in canonical form
   auto hbar = this->hbar();
@@ -137,11 +170,11 @@ std::vector<ExprPtr> CC::t(size_t pmax, size_t pmin) const {
 }
 
 std::vector<ExprPtr> CC::λ() const {
-  SEQUANT_ASSERT(!unitary() && "there is no need for CC::λ for unitary ansatz");
+  SEQUANT_ASSERT(!unitary(), "there is no need for CC::λ for unitary ansatz");
 
   // construct hbar
   const auto commutator_rank = hbar_comm_rank_.value_or(4);
-  SEQUANT_ASSERT(commutator_rank >= 1 && "CC::λ: hbar_comm_rank must be >= 1");
+  SEQUANT_ASSERT(commutator_rank >= 1, "CC::λ: hbar_comm_rank must be >= 1");
   auto hbar = this->hbar(commutator_rank -
                          1);  // -1 because of the connection with the projector
 
@@ -205,6 +238,9 @@ std::vector<ExprPtr> CC::λ() const {
 }
 
 ExprPtr CC::rdm(size_t rank, std::optional<size_t> comm_rank) const {
+  SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli,
+                 "CC::rdm: the Bernoulli expansion is not supported yet");
+
   // 1. replacement operator {ã^{p_1..p_r}_{p_{r+1}..p_{2r}}} (see op::ã); its
   // indices are free, so they become the free indices of γ.
   auto replacer = op::ã(rank);
@@ -237,15 +273,17 @@ ExprPtr CC::rdm(size_t rank, std::optional<size_t> comm_rank) const {
 
 std::vector<ExprPtr> CC::tʼ(size_t rank, size_t order,
                             std::optional<size_t> nbatch) const {
-  SEQUANT_ASSERT(order == 1 &&
+  SEQUANT_ASSERT(order == 1,
                  "sequant::mbpt::CC::tʼ(): only first-order perturbation is "
                  "supported now");
-  SEQUANT_ASSERT(rank == 1 &&
+  SEQUANT_ASSERT(rank == 1,
                  "sequant::mbpt::CC::tʼ(): only one-body perturbation "
                  "operator is supported now");
   if (unitary())
     SEQUANT_ASSERT(pertbar_comm_rank_,
                    "pertbar_comm_rank must be specified for unitary ansatz");
+  SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli,
+                 "CC::tʼ: the Bernoulli expansion is not supported yet");
 
   // construct h1_bar
   // truncate h1_bar at rank 2 for one-body perturbation operator and at rank 4
@@ -297,15 +335,14 @@ std::vector<ExprPtr> CC::tʼ(size_t rank, size_t order,
 
 std::vector<ExprPtr> CC::λʼ(size_t rank, size_t order,
                             std::optional<size_t> nbatch) const {
-  SEQUANT_ASSERT(order == 1 &&
+  SEQUANT_ASSERT(order == 1,
                  "sequant::mbpt::CC::λʼ(): only first-order perturbation is "
                  "supported now");
-  SEQUANT_ASSERT(rank == 1 &&
+  SEQUANT_ASSERT(rank == 1,
                  "sequant::mbpt::CC::λʼ(): only one-body perturbation "
                  "operator is supported now");
-  SEQUANT_ASSERT(!unitary() &&
-                 "there is no need for CC::λʼ for unitary ansatz");
-  SEQUANT_ASSERT(ansatz_ == Ansatz::T &&
+  SEQUANT_ASSERT(!unitary(), "there is no need for CC::λʼ for unitary ansatz");
+  SEQUANT_ASSERT(ansatz_ == Ansatz::T,
                  "CC::λʼ: only traditional ansatz is supported");
 
   // construct hbar
@@ -364,13 +401,104 @@ namespace {
 constexpr Normalization eom_norm = Normalization::SquareRoot;
 }  // namespace
 
-std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
-  SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order");
+// Per-block-truncated EOM sigma equations. For the qUCCSD ranks see
+// 10.1063/5.0062090 Sec. II C, Eqs. (29)-(48); for the IP/EA analogues,
+// 10.1021/acs.jctc.5c01991 Fig. 1 (Table 1 there maps out which H̄ components
+// enter each block, not the commutator ranks they are truncated at).
+//
+// Each block is the sandwich <i|H̄|j> of Eq. (7). Eq. (10) writes H̄ as
+// E_gr + a normal-ordered remainder and builds the blocks from the remainder
+// alone, so here the diagonal carries an explicit -<0|H̄|0> instead.
+std::vector<ExprPtr> CC::eom_r_blocked(
+    nₚ np, nₕ nh, const std::vector<size_t>& block_ranks) const {
+  SEQUANT_ASSERT(unitary(), "eom_r_blocked requires a unitary ansatz");
+
+  std::vector<std::pair<std::int64_t, std::int64_t>> manifolds;
+  for (std::int64_t rp = np, rh = nh; rp >= 0 && rh >= 0; --rp, --rh) {
+    if (rp == 0 && rh == 0) break;
+    manifolds.emplace_back(rp, rh);
+    if (rp == 0 || rh == 0) break;
+  }
+
+  std::ranges::reverse(manifolds);
+  const auto K = manifolds.size();
+  // empty means uniform truncation at hbar_comm_rank in every block
+  const std::vector<size_t> ranks =
+      block_ranks.empty() ? std::vector<size_t>(K * K, hbar_comm_rank().value())
+                          : block_ranks;
+  SEQUANT_ASSERT(ranks.size() == K * K,
+                 "CC::eom_r: block_ranks must be a K x K row-major matrix, "
+                 "K = number of projection manifolds");
+
+  // Bernoulli H̄ is tensor-level, BCH H̄ operator-level; the bra/ket/vev trio
+  // below must match it. Empty connectivity, as everywhere on the unitary path.
+  const bool tensor_level = hbar_expansion_ == HbarExpansion::Bernoulli;
+
+  // One H̄ per distinct truncation order, reduced to its R part (Bernoulli
+  // only: the operator-level BCH H̄ has no N/R split to take). The N part is
+  // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes only
+  // at the amplitude rank, so a block truncated below that rank would keep it.
+  // The diagonal is untouched either way: an N operator of rank r shifts the
+  // manifold rank by r, so it never lands on a diagonal block, and it has no
+  // reference expectation value. Off the diagonal removing it is a no-op only
+  // where the block rank equals hbar_comm_rank; below that rank the terms are
+  // off-shell, so the numbers change too. That is the point: Eqs. (41)-(47) of
+  // 10.1063/5.0062090 carry no such intermediate.
+  container::map<size_t, ExprPtr> hbars;
+  for (const auto k : ranks) {
+    auto [it, fresh] = hbars.try_emplace(k);
+    if (!fresh) continue;  // deriving H̄ twice for one rank is not cheap
+    it->second = hbar(k);
+    if (tensor_level) it->second = bernoulli::detail::R_part(it->second, N);
+  }
+  auto bra_of = [tensor_level](std::int64_t p, std::int64_t h) {
+    return tensor_level ? op::tensor::δl(nₚ(p), nₕ(h)) : op::δl(nₚ(p), nₕ(h));
+  };
+  auto ket_of = [tensor_level](std::int64_t p, std::int64_t h) {
+    return tensor_level ? op::tensor::r(nₚ(p), nₕ(h), eom_norm)
+                        : op::r(nₚ(p), nₕ(h), eom_norm);
+  };
+  auto vev = [tensor_level, this](const ExprPtr& e) {
+    return tensor_level ? op::tensor::ref_av(e)
+                        : op::ref_av(e, {.connect = {},
+                                         .screen = screen_,
+                                         .use_topology = use_topology_});
+  };
+
+  using std::min;
+  std::vector<ExprPtr> result(min(np, nh) + 1);
+  for (size_t i = 0; i < K; ++i) {
+    const auto [bp, bh] = manifolds[i];
+    const auto bra = bra_of(bp, bh);
+    auto acc = std::make_shared<Sum>();
+    for (size_t j = 0; j < K; ++j) {
+      const auto [kp, kh] = manifolds[j];
+      const auto& hbar_ij = hbars.at(ranks.at(i * K + j));
+      const auto ket = ket_of(kp, kh);
+      acc->append(vev(bra * hbar_ij * ket));
+      // -<0|H̄^(k_ii)|0>, written as <i|r_i H̄|0> so Wick keeps E's summed
+      // indices disjoint from the block's external ones.
+      if (i == j) acc->append(ex<Constant>(-1) * vev(bra * ket * hbar_ij));
+    }
+    result.at(static_cast<size_t>(min(bp, bh))) = simplify(ExprPtr{acc});
+  }
+  return result;
+}
+
+std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh,
+                               const std::vector<size_t>& block_ranks) const {
+  SEQUANT_ASSERT(np > 0 || nh > 0, "Unsupported excitation order");
   if (np != nh)
     SEQUANT_ASSERT(
-        get_default_context().spbasis() != SPBasis::Spinfree &&
+        get_default_context().spbasis() != SPBasis::Spinfree,
         "spin-free basis does not yet support non particle-conserving cases");
 
+  // Bernoulli always takes the blocked path: the uniform one below commutes H̄
+  // with an operator-level R, which a tensor-level H̄ cannot take part in. An
+  // empty matrix there means uniform truncation at hbar_comm_rank.
+  if (!block_ranks.empty() || hbar_expansion_ == HbarExpansion::Bernoulli)
+    return eom_r_blocked(np, nh, block_ranks);
+
   // construct hbar
   const auto hbar = this->hbar();
 
@@ -412,9 +540,9 @@ std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
 }
 
 std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) const {
-  SEQUANT_ASSERT(!unitary() &&
+  SEQUANT_ASSERT(!unitary(),
                  "there is no need for CC::eom_l for unitary ansatz");
-  SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order");
+  SEQUANT_ASSERT(np > 0 || nh > 0, "Unsupported excitation order");
 
   if (np != nh)
     SEQUANT_ASSERT(