From 802711e469d786c63024f4c9c8cdaac64fdbc4b6 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 00:57:34 -0400 Subject: [PATCH 01/37] feat(mbpt): operator-valued Wick reduction for the Bernoulli expansion Adds the bottom layer of the Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian: a Wick reduction that retains partial contractions, so a product of normal-ordered operators reduces to a sum of normal-ordered operators rather than collapsing to a scalar vacuum average, and the normal-ordered commutator built on it. WickTheorem::use_topology is disabled explicitly rather than left alone: it defaults to ON (wick.hpp), and its one-representative-times-multiplicity bookkeeping is only exercised by the fully-contracted path. On this partial-contraction path it rescales terms whose amplitude pairs are symmetric, which leaves vacuum averages correct while corrupting projections onto excited manifolds. wick_commutator reindexes B's summed indices to fresh temporaries before forming A*B, since A and B are independently constructed and may otherwise share labels, which would fuse two independent summations. --- CMakeLists.txt | 2 + SeQuant/domain/mbpt/bernoulli.cpp | 68 +++++++++++++++++++++++++++++++ SeQuant/domain/mbpt/bernoulli.hpp | 29 +++++++++++++ 3 files changed, 99 insertions(+) create mode 100644 SeQuant/domain/mbpt/bernoulli.cpp create mode 100644 SeQuant/domain/mbpt/bernoulli.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..4a6a6213a8 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp new file mode 100644 index 0000000000..619f94d5f8 --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -0,0 +1,68 @@ +#include + +#include +#include +#include +#include + +// Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian +// H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). Because σ mixes excitation +// and de-excitation the plain BCH series does not terminate; the Bernoulli +// expansion rewrites it so that Bernoulli numbers are the expansion +// coefficients, leaving the final truncation at a chosen commutator rank as the +// only approximation. +// +// This file builds that expansion in three layers: an operator-valued Wick +// reduction (here), the N/R operator split, and the rank-by-rank assembly of +// H̄. All equation references are to 10.1063/1.5030344 (Sec. III B). + +namespace sequant::mbpt::bernoulli { + +namespace detail { + +/// Operator-valued Wick reduction (see header): reduces a product of +/// normal-ordered operators to a sum of normal-ordered operators, retaining +/// partial contractions so the result is an operator, not a scalar VEV. +ExprPtr wick_reduce(ExprPtr expr) { + simplify(expr); + // full_contractions(false) is the whole point: it yields the normal-ordered + // operator form rather than the scalar VEV. Otherwise mirrors + // mbpt::tensor::expectation_value_impl. See core/wick.hpp. + FWickTheorem wick{expr}; + // use_topology MUST be disabled explicitly -- it defaults to ON + // (wick.hpp: `bool use_topology_ = true`), so merely not asking for it is not + // enough. It counts one representative per symmetry-equivalent contraction + // class times a multiplicity; the weight bookkeeping is exercised by the + // fully-contracted (vacuum-average) path, not by this partial-contraction + // one. With it on, <0|H̄³|0> stays correct but 16 of the 332 coefficients of + // <μ|H̄³|0> are rescaled by 2, 1/2, 3, 8/3 or 2/3 -- every one of them a term + // with a symmetric amplitude pair. With it off, all 332 signatures match + // pdaggerq exactly. Cost of giving up the optimisation: the rank-3-amplitude + // derivation goes 6.1 s -> 7.6 s wall. + wick.use_topology(false).full_contractions(false); + auto result = wick.compute(/*count_only=*/false, + /*skip_input_canonicalization=*/true); + simplify(result); + return result; +} + +/// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A) (see header). +ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { + // Disjoin B's (bound) indices from A's before forming the product: A and B + // are independently constructed operator expressions whose summed indices are + // local to each. If they happen to share labels (e.g. a block-resolved R/N + // part, which carries definite a/i/o/g indices, commuted with sigma, which + // also uses a/i), the naive product A*B would identify two independent + // summations, corrupting the contraction. Reindexing B to globally-fresh + // temporaries makes the two index sets disjoint; canonicalization restores + // tidy labels afterward. + container::map repl; + for (const auto& idx : get_used_indices(B)) + repl.emplace(idx, Index::make_tmp_index(idx.space())); + const auto Bd = repl.empty() ? B : transform_expr(B, repl); + return wick_reduce(simplify(A * Bd - Bd * A)); +} + +} // namespace detail + +} // namespace sequant::mbpt::bernoulli diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp new file mode 100644 index 0000000000..6153f3d2ef --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -0,0 +1,29 @@ +#ifndef SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP +#define SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP + +#include +#include + +namespace sequant::mbpt::bernoulli { + +namespace detail { + +/// Operator-valued Wick reduction: 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 × one residual +/// NormalOperator). Unlike the expectation-value path it keeps operators rather +/// than collapsing to a scalar VEV. +ExprPtr wick_reduce(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); + +} // namespace detail + +} // namespace sequant::mbpt::bernoulli + +#endif // SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP From 5fb0700099434d60e3fb7d81342c98cb1e047c50 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:01:29 -0400 Subject: [PATCH 02/37] feat(mbpt): N/R operator split for the Bernoulli expansion Adds the second layer: the split of an operator O into O_N, "the non-diagonal part containing all the excitation and de-excitation operators" (defined above Eq. (43) of 10.1063/1.5030344), and the rank-preserving remainder O_R = O - O_N. The Bernoulli expansion's inner commutators carry N/R subscripts, so every nesting level needs this classification. Classification needs definite index spaces, so expand_to_blocks first rewrites each general index of the residual NormalOperator as a sum over the base spaces it spans. Only the hole and particle spaces are expanded over: in the single-reference setting the remaining base spaces are empty, so restricting to those keeps the expansion 2-way per index instead of compounding across the nested commutators. That makes the routine single-reference only, which the header warns about. The rank cutoff mirrors pdaggerq (nt_bra > bernoulli_excitation_level -> R) rather than the paper's uncapped O_N, since that is the convention defining qUCCSD and the one the numbers are validated against; terms above the cutoff fall to R rather than being dropped. --- SeQuant/domain/mbpt/bernoulli.cpp | 229 +++++++++++++++++++++++++++++- SeQuant/domain/mbpt/bernoulli.hpp | 17 +++ 2 files changed, 243 insertions(+), 3 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 619f94d5f8..ff27196249 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -1,20 +1,100 @@ #include +#include #include +#include +#include #include #include +#include #include +#include +#include +#include + // Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian // H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). Because σ mixes excitation // and de-excitation the plain BCH series does not terminate; the Bernoulli // expansion rewrites it so that Bernoulli numbers are the expansion // coefficients, leaving the final truncation at a chosen commutator rank as the -// only approximation. +// only approximation. H is split as F (Fock, rank-preserving) + V +// (fluctuation potential), and every operator O is split into O_N (all +// excitation and de-excitation operators) and O_R = O − O_N. // // This file builds that expansion in three layers: an operator-valued Wick -// reduction (here), the N/R operator split, and the rank-by-rank assembly of -// H̄. All equation references are to 10.1063/1.5030344 (Sec. III B). +// reduction, the N/R operator split (here), and the rank-by-rank assembly of +// H̄. All equation references are to 10.1063/1.5030344 (Sec. III B); the N/R +// split and the UCC amplitude condition V̄_N = 0 are defined above and at +// Eq. (43). + +namespace { + +/// Returns the single residual fermionic NormalOperator carried by @p term, or +/// nullptr when it has none (a pure scalar / fully-contracted term). Every term +/// produced by wick_reduce is either a bare NormalOperator or a Product with +/// exactly one NormalOperator factor times tensor coefficients. +const sequant::NormalOperator* find_nop( + const sequant::ExprPtr& term) { + using namespace sequant; + if (term.is>()) + return &term.as>(); + if (term.is()) { + const NormalOperator* found = nullptr; + for (const auto& f : term.as().factors()) + if (f.is>()) { + // The one-residual-operator invariant is load-bearing: N/R + // classification reads this operator alone, so a second one would be + // silently ignored and misclassify the term. + SEQUANT_ASSERT(!found && + "find_nop: term carries >1 NormalOperator; wick_reduce " + "is expected to leave exactly one residual operator"); + found = &f.as>(); + } + return found; + } + return nullptr; +} + +/// Classifies one block-resolved term as N or R (Cancellation #2). A term is N +/// iff its single residual NormalOperator is a pure excitation (all creators +/// pure-unoccupied AND all annihilators pure-occupied) or a pure de-excitation +/// (the mirror), with rank ≤ @p cutoff. A term with no residual NormalOperator +/// is rank-preserving, hence R. +/// +/// Rank > @p cutoff falls to R rather than being dropped. The paper's O_N +/// (above Eq. (43): "containing all the excitation operators and de-excitation +/// operators in O") carries no rank cutoff, but this mirrors pdaggerq (`nt_bra +/// > bernoulli_excitation_level -> R`), which is the convention that defines +/// qUCCSD and the one our numbers are validated against. The choice is a real +/// degree of freedom: at cutoff = 2 the HF rank-3 energy is +2.507 mEh, versus +/// +2.428 mEh with no cutoff. +bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { + using namespace sequant; + auto isr = get_default_context().index_space_registry(); + const auto* nop = find_nop(term); + if (!nop) return false; // no residual operator => rank-preserving => R + const auto ncre = ranges::distance(nop->creators()); + const auto nann = ranges::distance(nop->annihilators()); + if (static_cast(std::max(ncre, nann)) > cutoff) return false; + auto all_unocc = [&](auto&& ops) { + return ranges::all_of(ops, [&](const auto& o) { + return isr->is_pure_unoccupied(o.index().space()); + }); + }; + auto all_occ = [&](auto&& ops) { + return ranges::all_of(ops, [&](const auto& o) { + return isr->is_pure_occupied(o.index().space()); + }); + }; + const bool pure_exc = + all_unocc(nop->creators()) && all_occ(nop->annihilators()); + const bool pure_deexc = + all_occ(nop->creators()) && all_unocc(nop->annihilators()); + return pure_exc || pure_deexc; +} + +} // namespace namespace sequant::mbpt::bernoulli { @@ -63,6 +143,149 @@ ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { return wick_reduce(simplify(A * Bd - Bd * A)); } +namespace { + +/// Core of expand_to_blocks for input already in wick_reduce'd form (a +/// simplified sum of coefficient × single-NormalOperator terms). Skipping the +/// reduction is an identity: wick_reduce is idempotent (terms with a single +/// residual NormalOperator admit no further contractions). @p expr is not +/// mutated. +ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { + auto isr = get_default_context().index_space_registry(); + const auto& bases = isr->base_spaces(); + + auto is_base_space = [&](const IndexSpace& sp) { + return ranges::any_of(bases, [&](const auto& b) { return b == sp; }); + }; + + // Only the hole and particle base spaces are physically populated in the + // single-reference qUCCSD setting; the other base spaces the complete space + // nominally spans (e.g. the SR "o"/"g" blocks) are empty, so splitting a + // general index into them only multiplies the term count without changing any + // projected quantity. Restricting to {hole, particle} keeps the expansion + // 2-way (occupied/virtual) per index -- exactly what the N/R classifier needs + // -- instead of splitting every base space, which compounds across the nested + // commutators. Falls back to all base spaces if the registry defines no + // hole/particle split. + const auto& hole_t = isr->hole_space(/*nulltype_ok=*/true); + const auto& particle_t = isr->particle_space(/*nulltype_ok=*/true); + auto physical = [&](const IndexSpace& b) { + return hole_t.includes(b.type()) || particle_t.includes(b.type()); + }; + + auto expand_term = [&](const ExprPtr& term) -> ExprPtr { + // collect the residual NormalOperator's distinct general (non-base) indices + const auto* nop = find_nop(term); + if (!nop) return term; // pure scalar/contraction: nothing to split + container::svector gens; + for (const auto& op : nop->creann()) { + if (!is_base_space(op.index().space()) && + ranges::none_of(gens, [&](const auto& g) { return g == op.index(); })) + gens.push_back(op.index()); + } + if (gens.empty()) return term; + // candidate base spaces per general index: base b is a sub-block of the + // general space iff its type bits are included and its quantum numbers + // match (stay within the same spin sector). + container::svector> choices; + for (const auto& g : gens) { + container::svector c, c_all; + for (const auto& b : bases) + if (b.qns() == g.space().qns() && g.space().type().includes(b.type())) { + c_all.push_back(b); + if (physical(b)) c.push_back(b); + } + choices.push_back(c.empty() ? c_all : c); + } + // cartesian product of assignments => sum of transformed terms; + // accumulate via Sum::append (linear) rather than operator+, which + // deep-copies the accumulated Sum on every call (quadratic) + auto sum = std::make_shared(); + container::svector idx(gens.size(), 0); + for (;;) { + container::map repl; + for (std::size_t k = 0; k < gens.size(); ++k) + // fresh ordinal (not gens[k].ordinal()): reusing the general index's + // ordinal would collide with any pre-existing definite index of the + // same base space and ordinal already in the term (e.g. an a/i index + // from an amplitude in a commutator result), producing a duplicate + // index. A globally-unique temporary is disjoint by construction; + // canonicalization restores tidy labels. + repl.emplace(gens[k], Index::make_tmp_index(choices[k][idx[k]])); + sum->append(transform_expr(term, repl)); + // increment mixed-radix counter over the assignments + std::size_t k = 0; + for (; k < gens.size(); ++k) { + if (++idx[k] < choices[k].size()) break; + idx[k] = 0; + } + if (k == gens.size()) break; + } + return ExprPtr{sum}; + }; + + ExprPtr out; + if (expr.is()) { + auto out_sum = std::make_shared(); + for (const auto& t : expr.as()) out_sum->append(expand_term(t)); + out = out_sum->empty() ? ex(0) : ExprPtr{out_sum}; + } else { + // clone: expand_term may return its argument, which simplify would mutate + out = expand_term(expr)->clone(); + } + simplify(out); + return out; +} + +/// Keeps only the N terms of block-resolved @p bx (shared tail of N_part and +/// N_part_reduced). +ExprPtr keep_N_terms(const ExprPtr& bx, std::size_t cutoff) { + if (bx.is()) { + auto out = std::make_shared(); + for (const auto& t : bx.as()) + if (is_N_term(t, cutoff)) out->append(t); + return out->empty() ? ex(0) : simplify(ExprPtr{out}); + } + return is_N_term(bx, cutoff) ? bx : ex(0); +} + +/// N_part for input already in wick_reduce'd form. +ExprPtr N_part_reduced(const ExprPtr& reduced, std::size_t cutoff) { + return keep_N_terms(expand_to_blocks_reduced(reduced), cutoff); +} + +/// R_part for input already in wick_reduce'd form. +ExprPtr R_part_reduced(const ExprPtr& reduced, std::size_t cutoff) { + return simplify(reduced - N_part_reduced(reduced, cutoff)); +} + +} // namespace + +/// Identity expansion of every general index into its base sub-blocks (see +/// header): after expansion every residual index is definite, so the N/R +/// classifier can act on it. +ExprPtr expand_to_blocks(const ExprPtr& expr_in) { + return expand_to_blocks_reduced(wick_reduce(expr_in->clone())); +} + +/// N part of @p expr at truncation @p cutoff (see header): block-resolve, then +/// keep only the pure excitation / de-excitation terms. +ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff) { + return keep_N_terms(expand_to_blocks(expr), cutoff); +} + +/// R part of @p expr at truncation @p cutoff (see header): the reduced operator +/// minus its N part. Because expand_to_blocks is an identity +/// (N ⊎ R = expr as operators), R = expr − N holds exactly while keeping expr +/// in its compact (general-index) form -- only N is block-resolved. This is +/// verified equivalent to the fully block-resolved remainder: ref_av([R,σ]) is +/// identical either way. Keeping expr compact makes the nested commutators that +/// consume R operate on far fewer terms. +ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { + auto reduced = wick_reduce(expr->clone()); + return R_part_reduced(reduced, cutoff); +} + } // namespace detail } // namespace sequant::mbpt::bernoulli diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index 6153f3d2ef..ee993adacc 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -22,6 +22,23 @@ ExprPtr wick_reduce(ExprPtr expr); /// 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), an +/// identity in the single-reference setting where the other base spaces are +/// empty. After expansion every residual index is definite so the +/// N/R classifier can act on it. Idempotent on block-resolved input. +ExprPtr expand_to_blocks(const ExprPtr& expr); + +/// Block-resolved N part (O_N of 10.1063/1.5030344, defined above Eq. (43)): +/// 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); + +/// Block-resolved R (rank-preserving remainder) part: +/// expand_to_blocks(expr) minus N_part(expr, cutoff). +ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff); + } // namespace detail } // namespace sequant::mbpt::bernoulli From eae176729f6d6cde9fa65b60ee4aa2368ac3804b Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:02:03 -0400 Subject: [PATCH 03/37] =?UTF-8?q?feat(mbpt):=20rank-by-rank=20assembly=20o?= =?UTF-8?q?f=20the=20Bernoulli=20H=CC=84,=20Eqs.=20(45)-(50)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the top layer: hbar(N, rank, skip1) sums H̄⁰..H̄^rank of 10.1063/1.5030344 Eq. (45), each order transcribed from its equation with the published coefficients and per-level N/R subscripts. Bernoulli numbers B₁=-1/2, B₂=1/12, B₃=0, B₄=-1/720 (Eq. 40) enter as those coefficients; a subscript R/N on a commutator means "form the commutator, then keep only its R/N part before the next nesting", which is what the split from the previous commit provides. Two cancellations from the paper are relied on and noted in place: F enters H̄ only at first order (stated just below Eq. (50)), and the higher orders carry only R-subscripted inner commutators. Every term is a nested commutator whose prefix is shared with other terms, within a rank and across ranks, so nest() memoizes each prefix (keyed by the base operator plus the tags applied so far). The nine rank-4 terms have only 3 distinct level-1 and 6 distinct level-2 nodes. Reusing a memoized ExprPtr is safe because expression composition deep-copies its operands. Contributions accumulate through Sum::append rather than chained operator+, which deep-copies the whole accumulated Sum on every call and is quadratic in the term count at high rank. --- SeQuant/domain/mbpt/bernoulli.cpp | 132 ++++++++++++++++++++++++++++-- SeQuant/domain/mbpt/bernoulli.hpp | 20 +++++ 2 files changed, 144 insertions(+), 8 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index ff27196249..0aae2d1978 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -4,29 +4,38 @@ #include #include #include +#include +#include #include #include #include #include +#include #include #include #include +#include +#include + // Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian // H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). Because σ mixes excitation // and de-excitation the plain BCH series does not terminate; the Bernoulli // expansion rewrites it so that Bernoulli numbers are the expansion // coefficients, leaving the final truncation at a chosen commutator rank as the -// only approximation. H is split as F (Fock, rank-preserving) + V -// (fluctuation potential), and every operator O is split into O_N (all -// excitation and de-excitation operators) and O_R = O − O_N. +// only approximation. H is split as F (Fock, rank-preserving) + V (fluctuation +// potential), and every operator O is split into O_N (all excitation and +// de-excitation operators) and O_R = O − O_N. At a converged RHF/UHF reference +// two cancellations hold: F survives only in H̄¹, and the higher orders carry +// only R-subscripted inner commutators. // -// This file builds that expansion in three layers: an operator-valued Wick -// reduction, the N/R operator split (here), and the rank-by-rank assembly of -// H̄. All equation references are to 10.1063/1.5030344 (Sec. III B); the N/R -// split and the UCC amplitude condition V̄_N = 0 are defined above and at -// Eq. (43). +// All equation references are to 10.1063/1.5030344 (Sec. III B): superoperator +// inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, B₃=0, B₄=−1/720 +// Eq. (40); the N/R split and the UCC amplitude condition V̄_N = 0 above and at +// Eq. (43); the iterative recursion for V̄ Eq. (44); the rank-by-rank operators +// H̄⁰..H̄⁴ Eqs. (45)-(50). Cancellation #1 (F enters only H̄¹) is stated just +// below Eq. (50). namespace { @@ -288,4 +297,111 @@ ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { } // namespace detail +/// Assembles H̄ order by order (see header), summing H̄⁰..H̄^rank of Eq. (45). +/// Each H̄^k below is transcribed from its equation, verified term-by-term +/// against the published coefficients and N/R subscripts. A subscript R/N on a +/// commutator means "form the commutator, then keep only its R/N part before +/// the next nesting". +ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { + if (rank > 4) + throw Exception("bernoulli::hbar: only ranks 0..4 are implemented"); + + using detail::N_part; + using detail::N_part_reduced; + using detail::R_part; + using detail::R_part_reduced; + using detail::wick_commutator; + const auto cutoff = N; + const auto F = op::tensor::F(); + const auto V = op::tensor::h(2); + const auto T = op::tensor::T(N, skip1); + const auto sigma = simplify(T - adjoint(T)); // σ = T − T† + auto c = [&](rational num, const ExprPtr& e) { + return ex(num) * e; + }; + + // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} + // with a per-level N/R/A partition tag applied after each commutator ('A' = + // no filter). nest(p0, f) evaluates such a node, memoizing every prefix + // (key = p0 + tags applied so far): the terms share prefixes both within a + // rank (the 9 rank-4 terms have only 3 distinct level-1 and 6 level-2 nodes) + // and across ranks (each H̄^k node is a prefix of H̄^{k+1} nodes), so the memo + // avoids recomputing them. Reusing a memoized ExprPtr is safe: expression + // composition deep-copies operands (Product/Sum append clone), so + // wick_commutator does not mutate its arguments. Commutator outputs are + // already wick_reduce'd, so the reduced-input N/R filters apply. + container::map memo; + auto nest = [&](char p0, const char* f) -> ExprPtr { + std::string key{p0}; + auto it = memo.find(key); + if (it == memo.end()) { + ExprPtr base = (p0 == 'N') ? N_part(V, cutoff) + : (p0 == 'R') ? R_part(V, cutoff) + : V; + it = memo.emplace(std::move(key), std::move(base)).first; + } + ExprPtr op = it->second; + for (int i = 0; f[i] != '\0'; ++i) { + key = it->first + f[i]; + it = memo.find(key); + if (it == memo.end()) { + auto cx = wick_commutator(op, sigma); + ExprPtr filtered = (f[i] == 'R') ? R_part_reduced(cx, cutoff) + : (f[i] == 'N') ? N_part_reduced(cx, cutoff) + : cx; + it = memo.emplace(std::move(key), std::move(filtered)).first; + } + op = it->second; + } + return op; + }; + + // accumulate H̄ contributions via Sum::append (linear) rather than chained + // operator+ / operator+=, each of which deep-copies the accumulated Sum -- + // quadratic in the (large) term count at high rank + auto acc = std::make_shared(); + acc->append(simplify(F + V)); // H̄⁰ = F + V [Eq. (46)] + if (rank >= 1) { + // H̄¹ = [F,σ] + ½[V,σ] + ½[V_R,σ] [Eq. (47)]. F enters H̄ ONLY here + // (Cancellation #1, stated just below Eq. (50)). + acc->append(wick_commutator(F, sigma)); + acc->append(c({1, 2}, nest('A', "A"))); + acc->append(c({1, 2}, nest('R', "A"))); + } + if (rank >= 2) { + // H̄² = 1/12[[V_N,σ],σ] + ¼[[V,σ]_R,σ] + ¼[[V_R,σ]_R,σ] [Eq. (48)] + acc->append(c({1, 12}, nest('N', "AA"))); + acc->append(c({1, 4}, nest('A', "RA"))); + acc->append(c({1, 4}, nest('R', "RA"))); + } + if (rank >= 3) { + // H̄³ = 1/24[[[V_N,σ],σ]_R,σ] + ⅛[[[V,σ]_R,σ]_R,σ] + ⅛[[[V_R,σ]_R,σ]_R,σ] + // − 1/24[[[V,σ]_R,σ],σ] − 1/24[[[V_R,σ]_R,σ],σ] [Eq. (49)] + acc->append(c({1, 24}, nest('N', "ARA"))); + acc->append(c({1, 8}, nest('A', "RRA"))); + acc->append(c({1, 8}, nest('R', "RRA"))); + acc->append(c({-1, 24}, nest('A', "RAA"))); + acc->append(c({-1, 24}, nest('R', "RAA"))); + } + if (rank >= 4) { + // H̄⁴ = Eq. (50), the nine order-4 terms produced by the recursion Eq. (44), + // V̄^{k+1} = σ̂F + X̂⁻¹(σ̂)e^{σ̂}V − Σ_{n≠0} B_n σ̂^n V̄_R^{k}. F is absent here + // (Cancellation #1). Listed in the paper's order; the outermost tag is + // always A. Verified term-by-term against Eq. (50): coefficients + // +1/16, +1/16, +1/48, −1/48, −1/48, −1/144, −1/48, −1/48, −1/720 and the + // N/R subscripts all match as transcribed. + acc->append(c({1, 16}, nest('R', "RRRA"))); + acc->append(c({1, 16}, nest('A', "RRRA"))); + acc->append(c({1, 48}, nest('N', "ARRA"))); + acc->append(c({-1, 48}, nest('A', "RARA"))); + acc->append(c({-1, 48}, nest('R', "RARA"))); + acc->append(c({-1, 144}, nest('N', "ARAA"))); + acc->append(c({-1, 48}, nest('A', "RRAA"))); + acc->append(c({-1, 48}, nest('R', "RRAA"))); + acc->append(c({-1, 720}, nest('N', "AAAA"))); + } + ExprPtr result{std::move(acc)}; + return simplify(result); +} + } // namespace sequant::mbpt::bernoulli diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index ee993adacc..e5845222dc 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -6,6 +6,26 @@ 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 relies on expanding general +/// indices over the hole and particle spaces alone (see +/// detail::expand_to_blocks), which is an identity only when the remaining base +/// spaces are empty. Under a multireference registry the R part comes out wrong +/// silently -- there is no check for this. +/// +/// @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 { /// Operator-valued Wick reduction: applies Wick's theorem to @p expr retaining From d69340bc67b0c2e34370ac6ed7241d2576208500 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:03:02 -0400 Subject: [PATCH 04/37] =?UTF-8?q?feat(mbpt):=20select=20the=20Bernoulli=20?= =?UTF-8?q?H=CC=84=20expansion=20from=20the=20CC=20engine?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds CC::Options::hbar_expansion (BCH by default, Bernoulli opt-in) and dispatches CC::hbar() to bernoulli::hbar() when it is selected. Two constructor assertions guard the combination: the Bernoulli expansion is defined for the unitary ansatz only, and it requires an explicit hbar_comm_rank, since CC::hbar() otherwise falls back to rank 4 and would silently select the most expensive and least exercised order. CC::energy() takes the plain reference expectation value under this expansion: the tensor-level H̄ is already fully expanded, so no operator connectivity remains to constrain. Its comm_rank argument defaults to the amplitude rank and is passed explicitly for the qUCCSD [2|3] split, where the energy is taken at H̄³ while the amplitudes stop at H̄². --- SeQuant/domain/mbpt/models/cc.cpp | 39 ++++++++++++++++++++++++++++++- SeQuant/domain/mbpt/models/cc.hpp | 13 +++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 36a4ab1108..4793984a08 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -41,7 +42,8 @@ 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"); @@ -49,6 +51,14 @@ CC::CC(size_t n, const Options& opts) SEQUANT_ASSERT( skip_singles_ && "CC: skip_singles must be true for orbital-optimized ansatz"); + if (hbar_expansion_ == HbarExpansion::Bernoulli) { + SEQUANT_ASSERT(unitary(), + "CC: Bernoulli expansion requires a unitary ansatz"); + // without hbar_comm_rank CC::hbar() falls back to rank 4, silently + // selecting the most expensive (and least exercised) order + SEQUANT_ASSERT(hbar_comm_rank_, + "CC: Bernoulli expansion requires hbar_comm_rank"); + } } CC::Ansatz CC::ansatz() const { return ansatz_; } @@ -59,6 +69,8 @@ bool CC::unitary() const { std::optional 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_; } @@ -68,6 +80,9 @@ bool CC::use_topology() const { return use_topology_; } ExprPtr CC::hbar(std::optional 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()) @@ -75,6 +90,14 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { } ExprPtr CC::energy(std::optional comm_rank) const { + // Bernoulli: the tensor-level H̄ is already fully expanded, so there is no + // operator connectivity left to constrain -- take the plain reference + // expectation value. The energy rank defaults to the amplitude rank; + // pass comm_rank explicitly for the qUCCSD [2|3] split (energy at H̄³). + if (hbar_expansion_ == HbarExpansion::Bernoulli) { + const auto erank = comm_rank.value_or(*hbar_comm_rank_); + 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). @@ -87,6 +110,20 @@ std::vector CC::t(size_t pmax, size_t pmin) const { pmax = (pmax == std::numeric_limits::max() ? N : pmax); SEQUANT_ASSERT(pmax >= pmin && "pmax should be >= pmin"); + // Bernoulli: project the tensor-level H̄ (built at the amplitude + // rank = hbar_comm_rank_) onto each manifold hbar(); + std::vector result(pmax + 1); + for (std::int64_t p = pmax; p >= static_cast(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(); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 278e8081c4..eafd9e755c 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -32,6 +32,13 @@ class CC { oU }; + enum class HbarExpansion { + /// standard Baker-Campbell-Hausdorff commutator expansion + BCH, + /// Bernoulli expansion, 10.1063/1.5030344 (unitary ansatz only) + Bernoulli + }; + /// Configuration options for CC class struct Options { SEQUANT_DESIGNATED_INIT_ONLY; @@ -54,6 +61,8 @@ class CC { /// perturbation operator; must be specified if unitary ansatz is used in /// perturbed amplitude derivation std::optional pertbar_comm_rank = std::nullopt; + /// choice of H̄ expansion; Bernoulli requires a unitary ansatz + HbarExpansion hbar_expansion = HbarExpansion::BCH; }; /// @brief constructs CC engine with default options (traditional ansatz, @@ -76,6 +85,9 @@ class CC { /// not set [[nodiscard]] std::optional hbar_comm_rank() const; + /// @return the choice of H̄ expansion + [[nodiscard]] HbarExpansion hbar_expansion() const; + /// @return true if singles amplitudes are excluded from \f$ \hat{T} \f$ and /// \f$ \hat{\Lambda} \f$ [[nodiscard]] bool skip_singles() const; @@ -191,6 +203,7 @@ class CC { bool use_topology_ = true; std::optional hbar_comm_rank_ = std::nullopt; std::optional pertbar_comm_rank_ = std::nullopt; + HbarExpansion hbar_expansion_ = HbarExpansion::BCH; /// @return the `LSTOptions` this engine uses for every `mbpt::lst()` call /// @note The choice of commutator representation is really a question of From 7dca5621ffb42c9eb093e47ad198dd378586ec4f Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:03:30 -0400 Subject: [PATCH 05/37] test(mbpt): tests for Bernoullie UCC equations Pins the derived equations at Bernoulli ranks 1-3 for the unitary ansatz: term counts for the energy and for the singles/doubles residuals, plus the guards on invalid configurations (Bernoulli with a non-unitary ansatz, and Bernoulli without an explicit hbar_comm_rank). The rank-3 numbers (46 energy, 32 singles, 38 doubles terms) are the ones cross-checked term-by-term against pdaggerq, so a change here means the derivation changed. --- tests/unit/test_mbpt_cc.cpp | 157 ++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index f08b8c574b..52f276f98f 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -5,7 +5,9 @@ #include #include #include +#include #include +#include #include #include "catch2_sequant.hpp" @@ -14,6 +16,16 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { using namespace sequant; using namespace sequant::mbpt; + auto has_tensor = [](const ExprPtr& e, std::wstring label) { + bool found = false; + e->visit( + [&](const ExprPtr& n) { + if (n.is() && n.as().label() == label) found = true; + }, + /*atoms_only=*/true); + return found; + }; + SECTION("sr_tcc") { SECTION("t") { // TCC R1 @@ -50,6 +62,151 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { } // SECTION("λ") } + SECTION("bernoulli_wick") { + using namespace sequant; + using namespace sequant::mbpt; + // [V, T2] is antisymmetric: [A,B] == -[B,A] after Wick reduction + const auto V = op::tensor::h(2); + const auto T2 = op::tensor::t(2); // rank-2 excitation, tensor form + const auto ab = bernoulli::detail::wick_commutator(V, T2); + const auto ba = bernoulli::detail::wick_commutator(T2, V); + REQUIRE_THAT(ab, EquivalentTo(simplify(ex(-1) * ba))); + // wick_reduce of a bare (already normal-ordered) operator is itself + REQUIRE_THAT(bernoulli::detail::wick_reduce(V), EquivalentTo(V)); + // Wick reduction adds contractions beyond the naive V*T2 - T2*V + REQUIRE(bernoulli::detail::wick_commutator(V, T2) != ex(0)); + REQUIRE_THAT(bernoulli::detail::wick_commutator(V, T2), + !EquivalentTo(simplify(V * T2 - T2 * V))); + } + + SECTION("bernoulli_expand_to_blocks") { + using namespace sequant; + using namespace sequant::mbpt; + const auto V = op::tensor::h(2); // general g + const auto Vx = bernoulli::detail::expand_to_blocks(V); + // identity on each manifold: the expansion changes no physical content + for (const auto n : {1, 2}) + REQUIRE_THAT(op::tensor::ref_av(op::tensor::P(nₚ(n)) * Vx), + EquivalentTo(op::tensor::ref_av(op::tensor::P(nₚ(n)) * V))); + REQUIRE(Vx.is()); + REQUIRE(Vx.as().size() > 1); + REQUIRE_THAT(bernoulli::detail::expand_to_blocks(Vx), + EquivalentTo(Vx)); // idempotent + // no general index survives: every residual index is occ or uocc + auto isr = get_default_context().index_space_registry(); + Vx->visit( + [&](const ExprPtr& n) { + if (!n.is>()) return; + for (const auto& o : + n.as>().creann()) { + const auto& sp = o.index().space(); + REQUIRE((isr->is_pure_occupied(sp) || isr->is_pure_unoccupied(sp))); + } + }, + /*atoms_only=*/true); + } + + SECTION("bernoulli_N_R_split") { + using namespace sequant; + using namespace sequant::mbpt; + const auto V = op::tensor::h(2); // fluctuation potential g (general) + const auto Vn = bernoulli::detail::N_part(V, 2); + const auto Vr = bernoulli::detail::R_part(V, 2); + // N ⊎ R reconstructs V. R stays in compact general-index form, so check the + // identity on the manifolds rather than symbolically. + const auto NR = simplify(Vn + Vr); + for (const auto n : {1, 2}) + REQUIRE_THAT(op::tensor::ref_av(op::tensor::P(nₚ(n)) * NR), + EquivalentTo(op::tensor::ref_av(op::tensor::P(nₚ(n)) * V))); + REQUIRE(Vn != ex(0)); + REQUIRE(Vr != ex(0)); + // N is idempotent; R has no pure-exc/deexc content + REQUIRE_THAT(bernoulli::detail::N_part(Vn, 2), EquivalentTo(Vn)); + REQUIRE_THAT(bernoulli::detail::N_part(Vr, 2), + EquivalentTo(ex(0))); + } + + SECTION("bernoulli_hbar_structure") { + using namespace sequant; + using namespace sequant::mbpt; + // Cancellation #1: F appears only in H̄¹, so rank r − rank r−1 is F-free + // for r ≥ 2. + auto h0 = bernoulli::hbar(2, 0, false); + auto h1 = bernoulli::hbar(2, 1, false); + auto h2 = bernoulli::hbar(2, 2, false); + auto h3 = bernoulli::hbar(2, 3, false); + auto has_f = [&](const ExprPtr& e) { return has_tensor(e, L"f"); }; + REQUIRE(has_f(simplify(h1 - h0))); // [F,σ] + REQUIRE_FALSE(has_f(simplify(h2 - h1))); + REQUIRE_FALSE(has_f(simplify(h3 - h2))); + REQUIRE_THAT(h0, // H̄⁰ = F + V, Eq. (46) + EquivalentTo(simplify(op::tensor::F() + op::tensor::h(2)))); + } + + SECTION("bernoulli_config_validation") { + using namespace sequant; + using namespace sequant::mbpt; + // only ranks 0..4 are implemented. The CC-level preconditions (unitary + // ansatz, hbar_comm_rank set) are SEQUANT_ASSERTs per the class convention, + // so they are not testable here -- their behavior depends on + // SEQUANT_ASSERT_BEHAVIOR. + REQUIRE_THROWS_AS(bernoulli::hbar(2, 5, false), Exception); + } + + SECTION("bernoulli_quccsd") { + using namespace sequant; + using namespace sequant::mbpt; + const CC::Options opts{.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 2, + .hbar_expansion = CC::HbarExpansion::Bernoulli}; + CC cc(2, opts); + + // energy through H̄³, amplitudes through H̄² (hbar_comm_rank) + const auto E = cc.energy(3); + REQUIRE(E); + REQUIRE_THAT(E, !EquivalentTo(ex(0))); + const auto amps = cc.t(); + REQUIRE(amps.size() == 3); + REQUIRE(amps[1]); + REQUIRE(amps[2]); + REQUIRE_THAT(amps[1], !EquivalentTo(ex(0))); + REQUIRE_THAT(amps[2], !EquivalentTo(ex(0))); + + // energy at rank 3 vs amplitudes at rank 2, so unlike BCH/UCC the + // energy()==t()[0] invariant intentionally does not hold + REQUIRE_THAT(cc.energy(3), !EquivalentTo(cc.t().at(0))); + + // Reference expectation values of H̄¹ and H̄² (Eqs. (47), (48) of + // 10.1063/1.5030344), taken as successive-rank differences. + const auto E0 = op::tensor::ref_av(bernoulli::hbar(2, 0, false)); + const auto E1 = op::tensor::ref_av(bernoulli::hbar(2, 1, false)); + const auto E2 = op::tensor::ref_av(bernoulli::hbar(2, 2, false)); + const auto E1_contrib = simplify(E1 - E0); + const auto E2_contrib = simplify(E2 - E1); + + // <0|H̄¹|0>: g-content is exactly 1/8 σ_ij^ab + h.c.; the remainder + // is the [F,σ] Brillouin terms, which vanish at RHF. + const auto E1_g_closed = deserialize( + L"1/8 t{a_1,a_2;i_1,i_2}:A-N-S * g{i_1,i_2;a_1,a_2}:A-C-S " + L"+ 1/8 t⁺{i_1,i_2;a_1,a_2}:A-N-S * g{a_1,a_2;i_1,i_2}:A-C-S"); + const auto E1_brillouin = simplify(E1_contrib - E1_g_closed); + REQUIRE_FALSE(has_tensor(E1_brillouin, L"g")); + REQUIRE(has_tensor(E1_brillouin, L"f")); + + // <0|H̄²|0> = 1/12 σ_i^a σ_j^b + h.c. + REQUIRE_THAT(E2_contrib, + EquivalentTo(L"1/12 t{a_1;i_1}:A-N-S * t{a_2;i_2}:A-N-S " + L"* g{i_1,i_2;a_1,a_2}:A-C-S " + L"+ 1/12 t⁺{i_1;a_1}:A-N-S * t⁺{i_2;a_2}:A-N-S " + L"* g{a_1,a_2;i_1,i_2}:A-C-S")); + + // characterization goldens: term counts, frozen from a numerically + // validated run to catch changes to hbar/projection + REQUIRE(size(E) == 46); + REQUIRE(size(amps[1]) == 32); + REQUIRE(size(amps[2]) == 38); + } + SECTION("energy") { // CC::energy() must equal the p==0 element of CC::t() for both ansätze. const auto N = 2; From c806ceac1cb2239afa92591a2261c827782bd70a Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 27 Jul 2026 01:41:58 -0400 Subject: [PATCH 06/37] docs(mbpt): correct the Bernoulli comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Corrections: - Eq. (45) is the assembly H̄ = Σ_k H̄^k; (46)-(50) are the rank-by-rank operators. The file header attributed (45)-(50) to the latter. - expand_to_blocks: the SR "o"/"g" base spaces are not empty, so the old justification for dropping them was wrong. They are droppable because the single-reference projection annihilates those terms. The header @warning said the same wrong thing. - R_part's result is not block-resolved; it stays in compact general-index form. The header claimed the opposite. - wick_reduce leaves at most one residual NormalOperator, not exactly one -- fully-contracted terms carry none, which find_nop already handled. - The memo shares level-1 nodes across ranks; it is not a prefix relation. - The use_topology rescaling set is {2, 1/2, 1/3, 8/3, 2/3}; the comment said 3 where it should have said 1/3. --- SeQuant/domain/mbpt/bernoulli.cpp | 73 ++++++++++++++----------------- SeQuant/domain/mbpt/bernoulli.hpp | 25 ++++++----- 2 files changed, 48 insertions(+), 50 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 0aae2d1978..521d297e59 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -33,16 +33,16 @@ // All equation references are to 10.1063/1.5030344 (Sec. III B): superoperator // inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, B₃=0, B₄=−1/720 // Eq. (40); the N/R split and the UCC amplitude condition V̄_N = 0 above and at -// Eq. (43); the iterative recursion for V̄ Eq. (44); the rank-by-rank operators -// H̄⁰..H̄⁴ Eqs. (45)-(50). Cancellation #1 (F enters only H̄¹) is stated just -// below Eq. (50). +// Eq. (43); the iterative recursion for V̄ Eq. (44); the assembly +// H̄ = Σ_k H̄^k Eq. (45); the rank-by-rank operators H̄⁰..H̄⁴ Eqs. (46)-(50). +// Cancellation #1 (F enters only H̄¹) is stated just below Eq. (50). namespace { /// Returns the single residual fermionic NormalOperator carried by @p term, or /// nullptr when it has none (a pure scalar / fully-contracted term). Every term -/// produced by wick_reduce is either a bare NormalOperator or a Product with -/// exactly one NormalOperator factor times tensor coefficients. +/// produced by wick_reduce is either a bare NormalOperator or a Product with at +/// most one NormalOperator factor times tensor coefficients. const sequant::NormalOperator* find_nop( const sequant::ExprPtr& term) { using namespace sequant; @@ -57,7 +57,7 @@ const sequant::NormalOperator* find_nop( // silently ignored and misclassify the term. SEQUANT_ASSERT(!found && "find_nop: term carries >1 NormalOperator; wick_reduce " - "is expected to leave exactly one residual operator"); + "is expected to leave at most one residual operator"); found = &f.as>(); } return found; @@ -75,9 +75,8 @@ const sequant::NormalOperator* find_nop( /// (above Eq. (43): "containing all the excitation operators and de-excitation /// operators in O") carries no rank cutoff, but this mirrors pdaggerq (`nt_bra /// > bernoulli_excitation_level -> R`), which is the convention that defines -/// qUCCSD and the one our numbers are validated against. The choice is a real -/// degree of freedom: at cutoff = 2 the HF rank-3 energy is +2.507 mEh, versus -/// +2.428 mEh with no cutoff. +/// qUCCSD. The cutoff is a real degree of freedom, not a formality: it moves +/// the correlation energy at the sub-mEh level. bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { using namespace sequant; auto isr = get_default_context().index_space_registry(); @@ -118,16 +117,15 @@ ExprPtr wick_reduce(ExprPtr expr) { // operator form rather than the scalar VEV. Otherwise mirrors // mbpt::tensor::expectation_value_impl. See core/wick.hpp. FWickTheorem wick{expr}; - // use_topology MUST be disabled explicitly -- it defaults to ON - // (wick.hpp: `bool use_topology_ = true`), so merely not asking for it is not - // enough. It counts one representative per symmetry-equivalent contraction - // class times a multiplicity; the weight bookkeeping is exercised by the - // fully-contracted (vacuum-average) path, not by this partial-contraction - // one. With it on, <0|H̄³|0> stays correct but 16 of the 332 coefficients of - // <μ|H̄³|0> are rescaled by 2, 1/2, 3, 8/3 or 2/3 -- every one of them a term - // with a symmetric amplitude pair. With it off, all 332 signatures match - // pdaggerq exactly. Cost of giving up the optimisation: the rank-3-amplitude - // derivation goes 6.1 s -> 7.6 s wall. + // use_topology must be disabled explicitly: it defaults to ON (wick.hpp: + // `bool use_topology_ = true`; the doc block at wick.hpp:154 claims otherwise + // and is stale), so not asking for it is not enough. It keeps one + // representative per symmetry-equivalent contraction class and multiplies by + // the class size, a weight bookkeeping that only holds on the + // fully-contracted path -- not on this partial-contraction one, where it + // silently rescales the terms carrying a symmetric amplitude pair. The vacuum + // expectation value stays correct either way, so the damage shows up only + // under projection. Turning it back on costs correctness, not just speed. wick.use_topology(false).full_contractions(false); auto result = wick.compute(/*count_only=*/false, /*skip_input_canonicalization=*/true); @@ -167,15 +165,14 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { return ranges::any_of(bases, [&](const auto& b) { return b == sp; }); }; - // Only the hole and particle base spaces are physically populated in the - // single-reference qUCCSD setting; the other base spaces the complete space - // nominally spans (e.g. the SR "o"/"g" blocks) are empty, so splitting a - // general index into them only multiplies the term count without changing any - // projected quantity. Restricting to {hole, particle} keeps the expansion - // 2-way (occupied/virtual) per index -- exactly what the N/R classifier needs - // -- instead of splitting every base space, which compounds across the nested - // commutators. Falls back to all base spaces if the registry defines no - // hole/particle split. + // Split each general index over the hole and particle base spaces only. A + // general index also spans the registry's other base spaces (under the SR + // convention the frozen-core "o" and inactive-virtual "g"), but terms landing + // in those are annihilated by the single-reference projection onto the + // hole/particle manifolds, so dropping them changes no projected quantity. + // This keeps the expansion 2-way per index instead of 4-way, which otherwise + // compounds across the nested commutators. Falls back to all base spaces if + // the registry defines no hole/particle split. const auto& hole_t = isr->hole_space(/*nulltype_ok=*/true); const auto& particle_t = isr->particle_space(/*nulltype_ok=*/true); auto physical = [&](const IndexSpace& b) { @@ -287,9 +284,8 @@ ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff) { /// minus its N part. Because expand_to_blocks is an identity /// (N ⊎ R = expr as operators), R = expr − N holds exactly while keeping expr /// in its compact (general-index) form -- only N is block-resolved. This is -/// verified equivalent to the fully block-resolved remainder: ref_av([R,σ]) is -/// identical either way. Keeping expr compact makes the nested commutators that -/// consume R operate on far fewer terms. +/// equivalent to the fully block-resolved remainder, and keeping expr compact +/// makes the nested commutators that consume R operate on far fewer terms. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { auto reduced = wick_reduce(expr->clone()); return R_part_reduced(reduced, cutoff); @@ -298,9 +294,8 @@ ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { } // namespace detail /// Assembles H̄ order by order (see header), summing H̄⁰..H̄^rank of Eq. (45). -/// Each H̄^k below is transcribed from its equation, verified term-by-term -/// against the published coefficients and N/R subscripts. A subscript R/N on a -/// commutator means "form the commutator, then keep only its R/N part before +/// Each H̄^k below is a direct transcription of its equation. A subscript R/N on +/// a commutator means "form the commutator, then keep only its R/N part before /// the next nesting". ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { if (rank > 4) @@ -325,9 +320,9 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { // no filter). nest(p0, f) evaluates such a node, memoizing every prefix // (key = p0 + tags applied so far): the terms share prefixes both within a // rank (the 9 rank-4 terms have only 3 distinct level-1 and 6 level-2 nodes) - // and across ranks (each H̄^k node is a prefix of H̄^{k+1} nodes), so the memo - // avoids recomputing them. Reusing a memoized ExprPtr is safe: expression - // composition deep-copies operands (Product/Sum append clone), so + // and across ranks (all four ranks share the same three level-1 nodes), so + // the memo avoids recomputing them. Reusing a memoized ExprPtr is safe: + // expression composition deep-copies operands (Product/Sum append clone), so // wick_commutator does not mutate its arguments. Commutator outputs are // already wick_reduce'd, so the reduced-input N/R filters apply. container::map memo; @@ -387,9 +382,7 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { // H̄⁴ = Eq. (50), the nine order-4 terms produced by the recursion Eq. (44), // V̄^{k+1} = σ̂F + X̂⁻¹(σ̂)e^{σ̂}V − Σ_{n≠0} B_n σ̂^n V̄_R^{k}. F is absent here // (Cancellation #1). Listed in the paper's order; the outermost tag is - // always A. Verified term-by-term against Eq. (50): coefficients - // +1/16, +1/16, +1/48, −1/48, −1/48, −1/144, −1/48, −1/48, −1/720 and the - // N/R subscripts all match as transcribed. + // always A. acc->append(c({1, 16}, nest('R', "RRRA"))); acc->append(c({1, 16}, nest('A', "RRRA"))); acc->append(c({1, 48}, nest('N', "ARRA"))); diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index e5845222dc..5670a295e0 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -14,11 +14,12 @@ namespace sequant::mbpt::bernoulli { /// 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 relies on expanding general -/// indices over the hole and particle spaces alone (see -/// detail::expand_to_blocks), which is an identity only when the remaining base -/// spaces are empty. Under a multireference registry the R part comes out wrong -/// silently -- there is no check for this. +/// @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 silently -- there is no check for this. /// /// @param N cluster/excitation rank (also the N/R rank cutoff) /// @param rank highest Bernoulli order H̄^k to include (0..4) @@ -30,9 +31,10 @@ namespace detail { /// Operator-valued Wick reduction: 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 × one residual -/// NormalOperator). Unlike the expectation-value path it keeps operators rather -/// than collapsing to a scalar VEV. +/// 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. ExprPtr wick_reduce(ExprPtr expr); /// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A). NOT the bare @@ -55,8 +57,11 @@ ExprPtr expand_to_blocks(const ExprPtr& expr); /// @p cutoff. Applies expand_to_blocks first. ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff); -/// Block-resolved R (rank-preserving remainder) part: -/// expand_to_blocks(expr) minus N_part(expr, cutoff). +/// R (rank-preserving remainder) part: wick_reduce(expr) minus +/// N_part(expr, cutoff). Unlike N_part the result is NOT block-resolved -- it +/// stays in compact general-index form, which is exact here because +/// expand_to_blocks is an identity, and much cheaper for the nested +/// commutators that consume R. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff); } // namespace detail From 3086ffeaf4298df252647b502cd279cb1cd68a49 Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 27 Jul 2026 08:03:08 -0400 Subject: [PATCH 07/37] test(mbpt): cleanup Bernoulli related unit tests --- tests/unit/test_mbpt_cc.cpp | 76 ++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 43 deletions(-) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 52f276f98f..5f732607ee 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -129,27 +129,47 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { SECTION("bernoulli_hbar_structure") { using namespace sequant; using namespace sequant::mbpt; + // Equation references are to 10.1063/1.5030344, Sec. III B. // Cancellation #1: F appears only in H̄¹, so rank r − rank r−1 is F-free // for r ≥ 2. auto h0 = bernoulli::hbar(2, 0, false); auto h1 = bernoulli::hbar(2, 1, false); auto h2 = bernoulli::hbar(2, 2, false); - auto h3 = bernoulli::hbar(2, 3, false); auto has_f = [&](const ExprPtr& e) { return has_tensor(e, L"f"); }; REQUIRE(has_f(simplify(h1 - h0))); // [F,σ] REQUIRE_FALSE(has_f(simplify(h2 - h1))); - REQUIRE_FALSE(has_f(simplify(h3 - h2))); REQUIRE_THAT(h0, // H̄⁰ = F + V, Eq. (46) EquivalentTo(simplify(op::tensor::F() + op::tensor::h(2)))); + + // Reference expectation values of H̄¹ and H̄², Eqs. (47) and (48), taken as + // successive-rank differences. + const auto E0 = op::tensor::ref_av(h0); + const auto E1 = op::tensor::ref_av(h1); + const auto E2 = op::tensor::ref_av(h2); + const auto E1_contrib = simplify(E1 - E0); + const auto E2_contrib = simplify(E2 - E1); + + // <0|H̄¹|0>: g-content is exactly 1/8 σ_ij^ab + h.c.; the remainder + // is the [F,σ] Brillouin terms, which vanish at RHF. + const auto E1_g_closed = deserialize( + L"1/8 t{a_1,a_2;i_1,i_2}:A-N-S * g{i_1,i_2;a_1,a_2}:A-C-S " + L"+ 1/8 t⁺{i_1,i_2;a_1,a_2}:A-N-S * g{a_1,a_2;i_1,i_2}:A-C-S"); + const auto E1_brillouin = simplify(E1_contrib - E1_g_closed); + REQUIRE_FALSE(has_tensor(E1_brillouin, L"g")); + REQUIRE(has_tensor(E1_brillouin, L"f")); + + // <0|H̄²|0> = 1/12 σ_i^a σ_j^b + h.c. + REQUIRE_THAT(E2_contrib, + EquivalentTo(L"1/12 t{a_1;i_1}:A-N-S * t{a_2;i_2}:A-N-S " + L"* g{i_1,i_2;a_1,a_2}:A-C-S " + L"+ 1/12 t⁺{i_1;a_1}:A-N-S * t⁺{i_2;a_2}:A-N-S " + L"* g{a_1,a_2;i_1,i_2}:A-C-S")); } SECTION("bernoulli_config_validation") { using namespace sequant; using namespace sequant::mbpt; - // only ranks 0..4 are implemented. The CC-level preconditions (unitary - // ansatz, hbar_comm_rank set) are SEQUANT_ASSERTs per the class convention, - // so they are not testable here -- their behavior depends on - // SEQUANT_ASSERT_BEHAVIOR. + // only ranks 0..4 are implemented. REQUIRE_THROWS_AS(bernoulli::hbar(2, 5, false), Exception); } @@ -161,50 +181,20 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { .hbar_expansion = CC::HbarExpansion::Bernoulli}; CC cc(2, opts); - // energy through H̄³, amplitudes through H̄² (hbar_comm_rank) - const auto E = cc.energy(3); - REQUIRE(E); - REQUIRE_THAT(E, !EquivalentTo(ex(0))); + // amplitudes through H̄² (hbar_comm_rank) const auto amps = cc.t(); REQUIRE(amps.size() == 3); - REQUIRE(amps[1]); - REQUIRE(amps[2]); REQUIRE_THAT(amps[1], !EquivalentTo(ex(0))); REQUIRE_THAT(amps[2], !EquivalentTo(ex(0))); - // energy at rank 3 vs amplitudes at rank 2, so unlike BCH/UCC the - // energy()==t()[0] invariant intentionally does not hold - REQUIRE_THAT(cc.energy(3), !EquivalentTo(cc.t().at(0))); - - // Reference expectation values of H̄¹ and H̄² (Eqs. (47), (48) of - // 10.1063/1.5030344), taken as successive-rank differences. - const auto E0 = op::tensor::ref_av(bernoulli::hbar(2, 0, false)); - const auto E1 = op::tensor::ref_av(bernoulli::hbar(2, 1, false)); - const auto E2 = op::tensor::ref_av(bernoulli::hbar(2, 2, false)); - const auto E1_contrib = simplify(E1 - E0); - const auto E2_contrib = simplify(E2 - E1); - - // <0|H̄¹|0>: g-content is exactly 1/8 σ_ij^ab + h.c.; the remainder - // is the [F,σ] Brillouin terms, which vanish at RHF. - const auto E1_g_closed = deserialize( - L"1/8 t{a_1,a_2;i_1,i_2}:A-N-S * g{i_1,i_2;a_1,a_2}:A-C-S " - L"+ 1/8 t⁺{i_1,i_2;a_1,a_2}:A-N-S * g{a_1,a_2;i_1,i_2}:A-C-S"); - const auto E1_brillouin = simplify(E1_contrib - E1_g_closed); - REQUIRE_FALSE(has_tensor(E1_brillouin, L"g")); - REQUIRE(has_tensor(E1_brillouin, L"f")); - - // <0|H̄²|0> = 1/12 σ_i^a σ_j^b + h.c. - REQUIRE_THAT(E2_contrib, - EquivalentTo(L"1/12 t{a_1;i_1}:A-N-S * t{a_2;i_2}:A-N-S " - L"* g{i_1,i_2;a_1,a_2}:A-C-S " - L"+ 1/12 t⁺{i_1;a_1}:A-N-S * t⁺{i_2;a_2}:A-N-S " - L"* g{a_1,a_2;i_1,i_2}:A-C-S")); - - // characterization goldens: term counts, frozen from a numerically - // validated run to catch changes to hbar/projection - REQUIRE(size(E) == 46); REQUIRE(size(amps[1]) == 32); REQUIRE(size(amps[2]) == 38); + +#ifndef SEQUANT_SKIP_LONG_TESTS + const auto E = cc.energy(3); + REQUIRE_THAT(E, !EquivalentTo(amps.at(0))); + REQUIRE(size(E) == 46); +#endif // !defined(SEQUANT_SKIP_LONG_TESTS) } SECTION("energy") { From 40ab0d4569432ff3e6efd65c24f63825a19bb890 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:38 -0400 Subject: [PATCH 08/37] refactor(mbpt): guard the Bernoulli partition tags, simplify the memo key MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A character outside {A,N,R} in a nest() tag string silently read as 'A' (no filter) and would have yielded the wrong H̄; the whole Eq. (46)-(50) transcription lives in these strings, so assert on every tag. Grow the memo key in place instead of deriving it from the memo iterator: container::map is a flat_map, whose insertions invalidate iterators, so the read-back was correct only by the accident of no insertion happening in between. Include for std::max and range/v3's primitives for ranges::distance rather than relying on them arriving transitively. --- SeQuant/domain/mbpt/bernoulli.cpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 521d297e59..5bd7690521 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -15,7 +15,9 @@ #include #include #include +#include +#include #include #include @@ -327,24 +329,30 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { // already wick_reduce'd, so the reduced-input N/R filters apply. container::map memo; auto nest = [&](char p0, const char* f) -> ExprPtr { + SEQUANT_ASSERT((p0 == 'A' || p0 == 'N' || p0 == 'R') && + "bernoulli::hbar: partition tag must be one of A, N, R"); + // grow `key` in place rather than deriving it from the memo iterator: + // container::map is a flat_map, whose insertions invalidate iterators std::string key{p0}; auto it = memo.find(key); if (it == memo.end()) { ExprPtr base = (p0 == 'N') ? N_part(V, cutoff) : (p0 == 'R') ? R_part(V, cutoff) : V; - it = memo.emplace(std::move(key), std::move(base)).first; + it = memo.emplace(key, std::move(base)).first; } ExprPtr op = it->second; for (int i = 0; f[i] != '\0'; ++i) { - key = it->first + f[i]; + SEQUANT_ASSERT((f[i] == 'A' || f[i] == 'N' || f[i] == 'R') && + "bernoulli::hbar: partition tag must be one of A, N, R"); + key += f[i]; it = memo.find(key); if (it == memo.end()) { auto cx = wick_commutator(op, sigma); ExprPtr filtered = (f[i] == 'R') ? R_part_reduced(cx, cutoff) : (f[i] == 'N') ? N_part_reduced(cx, cutoff) : cx; - it = memo.emplace(std::move(key), std::move(filtered)).first; + it = memo.emplace(key, std::move(filtered)).first; } op = it->second; } From 98e882fdd8aced0552868e630d899bd9b1d8f4b3 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:46 -0400 Subject: [PATCH 09/37] =?UTF-8?q?docs(mbpt):=20note=20that=20the=20Bernoul?= =?UTF-8?q?li=20H=CC=84=20bypasses=20CC::ref=5Fav?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CC::Options::screen and use_topology reach the derivation only through CC::ref_av(); the Bernoulli path calls op::tensor::ref_av() directly and so picks up that function's own defaults instead. --- SeQuant/domain/mbpt/models/cc.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index eafd9e755c..dd55202189 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -61,7 +61,11 @@ class CC { /// perturbation operator; must be specified if unitary ansatz is used in /// perturbed amplitude derivation std::optional pertbar_comm_rank = std::nullopt; - /// choice of H̄ expansion; Bernoulli requires a unitary ansatz + /// choice of H̄ expansion; Bernoulli requires a unitary ansatz. + /// @note the Bernoulli H̄ is assembled at the tensor level and does not go + /// through CC::ref_av(), which is what forwards `screen` and + /// `use_topology`; it calls `op::tensor::ref_av()` with that function's own + /// defaults instead HbarExpansion hbar_expansion = HbarExpansion::BCH; }; From f5c15a237b8679a21d0fadfe5064dbfbb5853530 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:56 -0400 Subject: [PATCH 10/37] test(mbpt): take the has_tensor label by const reference --- tests/unit/test_mbpt_cc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 5f732607ee..2be3187f0a 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -16,7 +16,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { using namespace sequant; using namespace sequant::mbpt; - auto has_tensor = [](const ExprPtr& e, std::wstring label) { + auto has_tensor = [](const ExprPtr& e, const std::wstring& label) { bool found = false; e->visit( [&](const ExprPtr& n) { From af21188a80e8fb2eabc04bd2fa478ea8782ffd70 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:56 -0400 Subject: [PATCH 11/37] test(mbpt): UCC equation-derivation integration test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The srcc.cpp analogue for the unitary ansatz, covering both H̄ expansions. CC::t() yields the whole equation set in one derivation -- element 0 the energy, element R the residual -- and the term counts are pinned so a change in either expansion fails ctest. Registered variants run in seconds; the Bernoulli H̄⁴ pins are recorded but left out of ctest, since that configuration takes ~2 minutes against sub-second times for everything else in this directory. --- tests/integration/CMakeLists.txt | 4 + tests/integration/ucc.cpp | 136 +++++++++++++++++++++++++++++++ 2 files changed, 140 insertions(+) create mode 100644 tests/integration/ucc.cpp diff --git a/tests/integration/CMakeLists.txt b/tests/integration/CMakeLists.txt index cbb387ae94..495eb71441 100644 --- a/tests/integration/CMakeLists.txt +++ b/tests/integration/CMakeLists.txt @@ -13,6 +13,8 @@ if (NOT SEQUANT_INTERNAL_SKIP_LONG_TESTS) "osstcc.cpp" # Equation-of-motion Coupled-Cluster "eomcc.cpp -> 2 2h2p R|2 1h2p R|2 2h1p R|2 3h1p R|2 1h3p R|2 4h2p R|3 3h3p R" + # Unitary Coupled-Cluster, both H̄ expansions (BCH and Bernoulli). + "ucc.cpp -> 2 bch 2|2 bch 3|2 bernoulli 2|2 bernoulli 3" ) if (TARGET Eigen3::Eigen) @@ -28,6 +30,8 @@ else() "srcc.cpp -> |2 t csv sf" # Equation-of-motion Coupled-Cluster (reduced test set) "eomcc.cpp -> 2 2h1p R" + # Unitary Coupled-Cluster (reduced test set: one variant per expansion) + "ucc.cpp -> 2 bch 2|2 bernoulli 2" ) if (TARGET Eigen3::Eigen) # these examples require Eigen for full functionality diff --git a/tests/integration/ucc.cpp b/tests/integration/ucc.cpp new file mode 100644 index 0000000000..95ffb486d4 --- /dev/null +++ b/tests/integration/ucc.cpp @@ -0,0 +1,136 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +// Unitary CC (UCC) equation derivation: the srcc.cpp analogue for the unitary +// ansatz, covering both H̄ expansions: the standard BCH commutator series and +// the Bernoulli expansion of 10.1063/1.5030344. +// +// CC::t() yields the whole equation set in one derivation: element 0 is the +// energy <0|H̄|0>, element R>0 the residual . Term counts are pinned +// below. +// +// Usage: ucc [N] [bch|bernoulli] [RANK] [print] +// N cluster/excitation rank of T (default 2) +// RANK commutator truncation rank of H̄ (default 2) + +using namespace sequant; +using namespace sequant::mbpt; + +namespace { + +#define runtime_assert(tf) \ + if (!(tf)) { \ + std::ostringstream oss; \ + oss << "failed assert at line " << __LINE__ << " in function " \ + << __func__; \ + throw std::runtime_error(oss.str().c_str()); \ + } + +TimerPool<32> tpool; + +using Hbar = CC::HbarExpansion; + +const std::map str2expansion = { + {"bch", Hbar::BCH}, {"bernoulli", Hbar::Bernoulli}}; + +/// pinned term count of one equation +struct TermCounts { + Hbar expansion; + std::size_t n; ///< cluster rank + std::size_t rank; ///< H̄ commutator truncation rank + std::size_t r; ///< projection manifold rank; 0 = energy + std::size_t nterms; ///< expected number of terms +}; + +// Regression pins, not independent references +const std::vector pins = { + // clang-format off + // expansion, N, rank, R, terms + {Hbar::BCH, 2, 2, 0, 20}, {Hbar::BCH, 2, 2, 1, 44}, {Hbar::BCH, 2, 2, 2, 42}, + {Hbar::BCH, 2, 3, 0, 74}, {Hbar::BCH, 2, 3, 1, 219}, {Hbar::BCH, 2, 3, 2, 267}, + {Hbar::BCH, 2, 4, 0, 307}, {Hbar::BCH, 2, 4, 1, 1100}, {Hbar::BCH, 2, 4, 2, 1433}, + {Hbar::Bernoulli, 2, 2, 0, 6}, {Hbar::Bernoulli, 2, 2, 1, 32}, {Hbar::Bernoulli, 2, 2, 2, 38}, + {Hbar::Bernoulli, 2, 3, 0, 46}, {Hbar::Bernoulli, 2, 3, 1, 141}, {Hbar::Bernoulli, 2, 3, 2, 191}, + {Hbar::Bernoulli, 2, 4, 0, 203}, {Hbar::Bernoulli, 2, 4, 1, 722}, {Hbar::Bernoulli, 2, 4, 2, 1044}, + // clang-format on +}; + +void check(Hbar expansion, std::size_t n, std::size_t rank, std::size_t r, + std::size_t nterms) { + for (const auto& p : pins) + if (expansion == p.expansion && n == p.n && rank == p.rank && r == p.r) { + if (nterms != p.nterms) + std::wcout << "MISMATCH: expected " << p.nterms << " terms, got " + << nterms << std::endl; + runtime_assert(nterms == p.nterms); + return; + } +} + +} // namespace + +int main(int argc, char* argv[]) { + std::wcout.precision(std::numeric_limits::max_digits10); + sequant::set_locale(); + + const std::size_t N = argc > 1 ? string_to(argv[1]) : 2; + const std::string expansion_str = argc > 2 ? argv[2] : "bch"; + const auto expansion = str2expansion.at(expansion_str); + const std::size_t RANK = argc > 3 ? string_to(argv[3]) : 2; + const bool print = argc > 4 && std::string(argv[4]) == "print"; + + sequant::detail::OpIdRegistrar op_id_registrar; + set_default_context({.index_space_registry_shared_ptr = make_sr_spaces(), + .vacuum = Vacuum::SingleProduct, + .metric = IndexSpaceMetric::Unit, + .spbasis = SPBasis::Spinor, + .first_dummy_index_ordinal = 100}); + TensorCanonicalizer::set_cardinal_tensor_labels(cardinal_tensor_labels()); + set_default_mbpt_context( + {.csv = mbpt::CSV::No, .op_registry_ptr = make_legacy_registry()}); + + std::cout << "SeQuant revision: " << sequant::git_revision() << "\n"; + std::cout << "Number of threads: " << sequant::num_threads() << "\n"; + + const CC cc(N, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = RANK, + .hbar_expansion = expansion}); + + tpool.clear(); + tpool.start(0); + const auto eqvec = cc.t(); + tpool.stop(0); + + std::wcout << "UCC equations [rank=" << N + << ",expansion=" << sequant::toUtf16(expansion_str) + << ",hbar_comm_rank=" << RANK << "] computed in " << tpool.read(0) + << " seconds" << std::endl; + + for (std::size_t R = 0; R < eqvec.size(); ++R) { + std::wcout << (R == 0 ? "E" : "R") << (R == 0 ? L"" : std::to_wstring(R)) + << "(expU" << N << ") has " << eqvec[R]->size() + << " terms:" << std::endl; + if (print) std::wcout << to_latex_align(eqvec[R], 20, 1) << std::endl; + check(expansion, N, RANK, R, eqvec[R]->size()); + } + + return 0; +} From fe55817ac22faba44301cba6a2365720e0cca7a4 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 23:56:19 -0400 Subject: [PATCH 12/37] =?UTF-8?q?perf(mbpt):=20collapse=20Bernoulli=20H?= =?UTF-8?q?=CC=84=20summands=20eagerly=20during=20assembly?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both accumulation sites in bernoulli.cpp built a Sum by append and left every duplicate for one final simplify. Sum::append flattens nested sums and adds up Constants but never merges like terms, so the nested commutators -- which overlap heavily by construction, the same fact that makes nest() memoize -- carried their duplicates all the way to the end. hbar() now accumulates into a HashingAccumulator, which keys summands by hash under proportional_to and merges them via Product::add_identical at insertion. The prefactor has to be folded into each summand rather than wrapped around the sum: appending Constant*Sum inserts the scaled sum as one opaque summand, since append's flatten splits a Sum but not a Product wrapping one, and nothing would collapse. Distributing it with expand() instead is shorter but materializes an intermediate Sum and gives back most of the gain (rank 4: 163.9 s vs 148.2 s). expand_to_blocks_reduced's outer loop now uses transform_sum_expr, which canonicalizes each mapped result before accumulating -- necessary here because the block assignments carry fresh temporary indices and so cannot hash-collide until canonical. It canonicalizes IN PLACE, hence expand_term now clones rather than returning one of its arguments in the two early-exit paths; that path is also parallel (std::execution::par_unseq), which is safe because Index::next_tmp_index is a static std::atomic. Derivation time, tests/integration/ucc 2 bernoulli , relwithdebinfo: rank 2 0.416 -> 0.364 s, rank 3 7.584 -> 7.216 s, rank 4 170.6 -> 148.2 s. Output is unchanged: serialize() of the rank-3 and rank-4 equations is byte-identical to the pre-change baseline (45 543 and 292 512 bytes), and repeat runs are byte-identical to each other despite the parallel index minting. Term counts alone would not have been sufficient evidence -- the use_topology bug rescaled terms while leaving counts and the VEV correct -- and to_latex() would not either, since it omits symmetry attributes. --- SeQuant/domain/mbpt/bernoulli.cpp | 80 +++++++++++++++++-------------- 1 file changed, 44 insertions(+), 36 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 5bd7690521..4dfd089126 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -184,14 +184,15 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { auto expand_term = [&](const ExprPtr& term) -> ExprPtr { // collect the residual NormalOperator's distinct general (non-base) indices const auto* nop = find_nop(term); - if (!nop) return term; // pure scalar/contraction: nothing to split + if (!nop) + return term->clone(); // pure scalar/contraction: nothing to split container::svector gens; for (const auto& op : nop->creann()) { if (!is_base_space(op.index().space()) && ranges::none_of(gens, [&](const auto& g) { return g == op.index(); })) gens.push_back(op.index()); } - if (gens.empty()) return term; + if (gens.empty()) return term->clone(); // candidate base spaces per general index: base b is a sub-block of the // general space iff its type bits are included and its quantum numbers // match (stay within the same spin sector). @@ -232,14 +233,13 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { return ExprPtr{sum}; }; + // transform_sum_expr maps in parallel, canonicalizes each result, and + // accumulates into a HashingAccumulator ExprPtr out; if (expr.is()) { - auto out_sum = std::make_shared(); - for (const auto& t : expr.as()) out_sum->append(expand_term(t)); - out = out_sum->empty() ? ex(0) : ExprPtr{out_sum}; + out = transform_sum_expr(expr.as().summands(), expand_term); } else { - // clone: expand_term may return its argument, which simplify would mutate - out = expand_term(expr)->clone(); + out = expand_term(expr); } simplify(out); return out; @@ -313,9 +313,6 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { const auto V = op::tensor::h(2); const auto T = op::tensor::T(N, skip1); const auto sigma = simplify(T - adjoint(T)); // σ = T − T† - auto c = [&](rational num, const ExprPtr& e) { - return ex(num) * e; - }; // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} // with a per-level N/R/A partition tag applied after each commutator ('A' = @@ -359,49 +356,60 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { return op; }; - // accumulate H̄ contributions via Sum::append (linear) rather than chained - // operator+ / operator+=, each of which deep-copies the accumulated Sum -- - // quadratic in the (large) term count at high rank - auto acc = std::make_shared(); - acc->append(simplify(F + V)); // H̄⁰ = F + V [Eq. (46)] + HashingAccumulator acc; + auto add = [&acc](rational num, const ExprPtr& e) { + if (e.is()) { + for (const auto& term : e.as()) { + auto scaled = ex(ExprPtrList{term}); + scaled.as().scale(num); + acc.append(std::move(scaled), /*flatten=*/false); + } + } else { + auto scaled = ex(ExprPtrList{e}); + scaled.as().scale(num); + acc.append(std::move(scaled), /*flatten=*/false); + } + }; + + add(1, simplify(F + V)); // H̄⁰ = F + V [Eq. (46)] if (rank >= 1) { // H̄¹ = [F,σ] + ½[V,σ] + ½[V_R,σ] [Eq. (47)]. F enters H̄ ONLY here // (Cancellation #1, stated just below Eq. (50)). - acc->append(wick_commutator(F, sigma)); - acc->append(c({1, 2}, nest('A', "A"))); - acc->append(c({1, 2}, nest('R', "A"))); + add(1, wick_commutator(F, sigma)); + add({1, 2}, nest('A', "A")); + add({1, 2}, nest('R', "A")); } if (rank >= 2) { // H̄² = 1/12[[V_N,σ],σ] + ¼[[V,σ]_R,σ] + ¼[[V_R,σ]_R,σ] [Eq. (48)] - acc->append(c({1, 12}, nest('N', "AA"))); - acc->append(c({1, 4}, nest('A', "RA"))); - acc->append(c({1, 4}, nest('R', "RA"))); + add({1, 12}, nest('N', "AA")); + add({1, 4}, nest('A', "RA")); + add({1, 4}, nest('R', "RA")); } if (rank >= 3) { // H̄³ = 1/24[[[V_N,σ],σ]_R,σ] + ⅛[[[V,σ]_R,σ]_R,σ] + ⅛[[[V_R,σ]_R,σ]_R,σ] // − 1/24[[[V,σ]_R,σ],σ] − 1/24[[[V_R,σ]_R,σ],σ] [Eq. (49)] - acc->append(c({1, 24}, nest('N', "ARA"))); - acc->append(c({1, 8}, nest('A', "RRA"))); - acc->append(c({1, 8}, nest('R', "RRA"))); - acc->append(c({-1, 24}, nest('A', "RAA"))); - acc->append(c({-1, 24}, nest('R', "RAA"))); + add({1, 24}, nest('N', "ARA")); + add({1, 8}, nest('A', "RRA")); + add({1, 8}, nest('R', "RRA")); + add({-1, 24}, nest('A', "RAA")); + add({-1, 24}, nest('R', "RAA")); } if (rank >= 4) { // H̄⁴ = Eq. (50), the nine order-4 terms produced by the recursion Eq. (44), // V̄^{k+1} = σ̂F + X̂⁻¹(σ̂)e^{σ̂}V − Σ_{n≠0} B_n σ̂^n V̄_R^{k}. F is absent here // (Cancellation #1). Listed in the paper's order; the outermost tag is // always A. - acc->append(c({1, 16}, nest('R', "RRRA"))); - acc->append(c({1, 16}, nest('A', "RRRA"))); - acc->append(c({1, 48}, nest('N', "ARRA"))); - acc->append(c({-1, 48}, nest('A', "RARA"))); - acc->append(c({-1, 48}, nest('R', "RARA"))); - acc->append(c({-1, 144}, nest('N', "ARAA"))); - acc->append(c({-1, 48}, nest('A', "RRAA"))); - acc->append(c({-1, 48}, nest('R', "RRAA"))); - acc->append(c({-1, 720}, nest('N', "AAAA"))); + add({1, 16}, nest('R', "RRRA")); + add({1, 16}, nest('A', "RRRA")); + add({1, 48}, nest('N', "ARRA")); + add({-1, 48}, nest('A', "RARA")); + add({-1, 48}, nest('R', "RARA")); + add({-1, 144}, nest('N', "ARAA")); + add({-1, 48}, nest('A', "RRAA")); + add({-1, 48}, nest('R', "RRAA")); + add({-1, 720}, nest('N', "AAAA")); } - ExprPtr result{std::move(acc)}; + auto result = acc.make_expr(); return simplify(result); } From 00e1a706cc3be05a74e578b9948739d1cc127cac Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:34:47 -0400 Subject: [PATCH 13/37] docs(mbpt): tighten the Bernoulli comments Trim the development narrative out of bernoulli.{cpp,hpp} and keep the reasons the code needs. The use_topology(false) comment keeps its cause (the flag defaults to ON and silently rescales terms carrying a symmetric amplitude pair on the partial-contraction path) and drops the stale wick.hpp line reference. is_N_term keeps why rank > cutoff falls to R rather than being dropped, and loses the paper quotes. Also collapse hbar's five using-declarations into one and fix the range notation in its error message. --- SeQuant/domain/mbpt/bernoulli.cpp | 84 +++++++++++++------------------ SeQuant/domain/mbpt/bernoulli.hpp | 26 +++++----- 2 files changed, 48 insertions(+), 62 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 4dfd089126..37f1c103d2 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -24,9 +24,9 @@ // Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian // H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). Because σ mixes excitation // and de-excitation the plain BCH series does not terminate; the Bernoulli -// expansion rewrites it so that Bernoulli numbers are the expansion -// coefficients, leaving the final truncation at a chosen commutator rank as the -// only approximation. H is split as F (Fock, rank-preserving) + V (fluctuation +// expansion rewrites it with Bernoulli numbers as the expansion coefficients, +// which leaves the truncation at a chosen commutator rank as the only +// approximation. H is split as F (Fock, rank-preserving) + V (fluctuation // potential), and every operator O is split into O_N (all excitation and // de-excitation operators) and O_R = O − O_N. At a converged RHF/UHF reference // two cancellations hold: F survives only in H̄¹, and the higher orders carry @@ -50,6 +50,7 @@ const sequant::NormalOperator* find_nop( using namespace sequant; if (term.is>()) return &term.as>(); + if (term.is()) { const NormalOperator* found = nullptr; for (const auto& f : term.as().factors()) @@ -62,6 +63,7 @@ const sequant::NormalOperator* find_nop( "is expected to leave at most one residual operator"); found = &f.as>(); } + return found; } return nullptr; @@ -70,23 +72,25 @@ const sequant::NormalOperator* find_nop( /// Classifies one block-resolved term as N or R (Cancellation #2). A term is N /// iff its single residual NormalOperator is a pure excitation (all creators /// pure-unoccupied AND all annihilators pure-occupied) or a pure de-excitation -/// (the mirror), with rank ≤ @p cutoff. A term with no residual NormalOperator -/// is rank-preserving, hence R. +/// (all creators pure-occupied AND all annihilators pure-unoccupied), with rank +/// ≤ @p cutoff. A term with no residual NormalOperator is rank-preserving, +/// hence R. /// -/// Rank > @p cutoff falls to R rather than being dropped. The paper's O_N -/// (above Eq. (43): "containing all the excitation operators and de-excitation -/// operators in O") carries no rank cutoff, but this mirrors pdaggerq (`nt_bra -/// > bernoulli_excitation_level -> R`), which is the convention that defines -/// qUCCSD. The cutoff is a real degree of freedom, not a formality: it moves -/// the correlation energy at the sub-mEh level. +/// Rank > @p cutoff falls to R rather than being dropped. An ]_R filter drops a +/// term only because the amplitude condition V̄_N = 0 (Eq. (43)) makes it zero, +/// and for σ truncated at rank N that condition covers rank ≤ N only. Eq. (43) +/// states O_N with no rank limit because there σ carries every rank. bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { using namespace sequant; auto isr = get_default_context().index_space_registry(); + const auto* nop = find_nop(term); if (!nop) return false; // no residual operator => rank-preserving => R + const auto ncre = ranges::distance(nop->creators()); const auto nann = ranges::distance(nop->annihilators()); if (static_cast(std::max(ncre, nann)) > cutoff) return false; + auto all_unocc = [&](auto&& ops) { return ranges::all_of(ops, [&](const auto& o) { return isr->is_pure_unoccupied(o.index().space()); @@ -97,10 +101,12 @@ bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { return isr->is_pure_occupied(o.index().space()); }); }; + const bool pure_exc = all_unocc(nop->creators()) && all_occ(nop->annihilators()); const bool pure_deexc = all_occ(nop->creators()) && all_unocc(nop->annihilators()); + return pure_exc || pure_deexc; } @@ -110,24 +116,15 @@ namespace sequant::mbpt::bernoulli { namespace detail { -/// Operator-valued Wick reduction (see header): reduces a product of -/// normal-ordered operators to a sum of normal-ordered operators, retaining -/// partial contractions so the result is an operator, not a scalar VEV. ExprPtr wick_reduce(ExprPtr expr) { simplify(expr); - // full_contractions(false) is the whole point: it yields the normal-ordered - // operator form rather than the scalar VEV. Otherwise mirrors - // mbpt::tensor::expectation_value_impl. See core/wick.hpp. FWickTheorem wick{expr}; - // use_topology must be disabled explicitly: it defaults to ON (wick.hpp: - // `bool use_topology_ = true`; the doc block at wick.hpp:154 claims otherwise - // and is stale), so not asking for it is not enough. It keeps one - // representative per symmetry-equivalent contraction class and multiplies by - // the class size, a weight bookkeeping that only holds on the - // fully-contracted path -- not on this partial-contraction one, where it - // silently rescales the terms carrying a symmetric amplitude pair. The vacuum - // expectation value stays correct either way, so the damage shows up only - // under projection. Turning it back on costs correctness, not just speed. + // use_topology defaults to ON, so it must be turned off explicitly. It keeps + // one representative per symmetry-equivalent contraction class and multiplies + // by the class size, weight bookkeeping that holds only on the + // fully-contracted path. On this partial-contraction path it silently + // rescales the terms carrying a symmetric amplitude pair, and the damage + // shows up only under projection. wick.use_topology(false).full_contractions(false); auto result = wick.compute(/*count_only=*/false, /*skip_input_canonicalization=*/true); @@ -135,16 +132,12 @@ ExprPtr wick_reduce(ExprPtr expr) { return result; } -/// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A) (see header). ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { - // Disjoin B's (bound) indices from A's before forming the product: A and B - // are independently constructed operator expressions whose summed indices are - // local to each. If they happen to share labels (e.g. a block-resolved R/N - // part, which carries definite a/i/o/g indices, commuted with sigma, which - // also uses a/i), the naive product A*B would identify two independent - // summations, corrupting the contraction. Reindexing B to globally-fresh - // temporaries makes the two index sets disjoint; canonicalization restores - // tidy labels afterward. + // A and B are built independently, so their summed indices are local to each. + // If both use the same labels (a block-resolved R/N part and sigma both carry + // a/i), the product A*B fuses two independent summations. That corrupts the + // contraction. Reindex B to fresh temporaries. Canonicalization restores tidy + // labels. container::map repl; for (const auto& idx : get_used_indices(B)) repl.emplace(idx, Index::make_tmp_index(idx.space())); @@ -154,9 +147,8 @@ ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { namespace { -/// Core of expand_to_blocks for input already in wick_reduce'd form (a -/// simplified sum of coefficient × single-NormalOperator terms). Skipping the -/// reduction is an identity: wick_reduce is idempotent (terms with a single +/// Core of expand_to_blocks for input already in wick_reduce'd form. Skipping +/// the reduction is an identity: wick_reduce is idempotent (terms with a single /// residual NormalOperator admit no further contractions). @p expr is not /// mutated. ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { @@ -284,10 +276,10 @@ ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff) { /// R part of @p expr at truncation @p cutoff (see header): the reduced operator /// minus its N part. Because expand_to_blocks is an identity -/// (N ⊎ R = expr as operators), R = expr − N holds exactly while keeping expr -/// in its compact (general-index) form -- only N is block-resolved. This is -/// equivalent to the fully block-resolved remainder, and keeping expr compact -/// makes the nested commutators that consume R operate on far fewer terms. +/// (N ⊎ R = expr as operators), R = expr − N holds exactly while expr stays in +/// its compact (general-index) form. Only N is block-resolved. The result +/// equals the fully block-resolved remainder, and the compact expr makes the +/// nested commutators that consume R operate on far fewer terms. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { auto reduced = wick_reduce(expr->clone()); return R_part_reduced(reduced, cutoff); @@ -301,13 +293,9 @@ ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { /// the next nesting". ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { if (rank > 4) - throw Exception("bernoulli::hbar: only ranks 0..4 are implemented"); + throw Exception("bernoulli::hbar: only ranks [0,4] are implemented"); - using detail::N_part; - using detail::N_part_reduced; - using detail::R_part; - using detail::R_part_reduced; - using detail::wick_commutator; + using namespace detail; const auto cutoff = N; const auto F = op::tensor::F(); const auto V = op::tensor::h(2); diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index 5670a295e0..4f924a2900 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -19,7 +19,7 @@ namespace sequant::mbpt::bernoulli { /// 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 silently -- there is no check for this. +/// out wrong. Nothing checks for this. /// /// @param N cluster/excitation rank (also the N/R rank cutoff) /// @param rank highest Bernoulli order H̄^k to include (0..4) @@ -29,12 +29,11 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1); namespace detail { -/// Operator-valued Wick reduction: 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. +/// 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. ExprPtr wick_reduce(ExprPtr expr); /// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A). NOT the bare @@ -51,17 +50,16 @@ ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B); /// N/R classifier can act on it. Idempotent on block-resolved input. ExprPtr expand_to_blocks(const ExprPtr& expr); -/// Block-resolved N part (O_N of 10.1063/1.5030344, defined above Eq. (43)): -/// the terms whose single residual -/// NormalOperator is a pure excitation or pure de-excitation of rank ≤ +/// 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 (rank-preserving remainder) part: wick_reduce(expr) minus -/// N_part(expr, cutoff). Unlike N_part the result is NOT block-resolved -- it -/// stays in compact general-index form, which is exact here because -/// expand_to_blocks is an identity, and much cheaper for the nested -/// commutators that consume R. +/// N_part(expr, cutoff). Unlike N_part the result is NOT block-resolved. It +/// stays in compact general-index form. That is exact here, because +/// expand_to_blocks is an identity, and much cheaper for the nested commutators +/// that consume R. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff); } // namespace detail From dc409713bdc3d88562e9ae1114d85e51a5ed9a00 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:00 -0400 Subject: [PATCH 14/37] =?UTF-8?q?feat(mbpt):=20per-block=20H=CC=84=20trunc?= =?UTF-8?q?ation=20for=20the=20EOM=20sigma=20equations?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CC::eom_r gains an optional block_ranks argument: a row-major K x K matrix over the projection manifolds giving each block of the secular matrix its own H̄ commutator truncation, instead of one uniform H̄ everywhere. The manifolds are indexed by ASCENDING rank, so the qUCCSD ranks {2,1,1,0} (10.1063/5.0062090 Table I, 10.1021/acs.jctc.5c01991 Table 1) serve EE, IP and EA alike. Each block is the sandwich plus an explicit -E shift on the diagonal at the block's own rank, not the commutator form : the commutator's extra - is manifold j's amplitude residual, which vanishes only when k equals the rank the amplitudes were converged against. Under the Bernoulli expansion each block's H̄ has its N part removed for the same reason. Empty block_ranks keeps the existing uniform path. That path commutes H̄ with an operator-level R, which the tensor-level Bernoulli H̄ cannot take part in, so it now throws instead of aborting inside op.ipp. The block shape and unitarity checks throw as well: SEQUANT_ASSERT compiles away under SEQUANT_ASSERT_BEHAVIOR=IGNORE, and the shape check guards an out-of-bounds read of block_ranks. --- SeQuant/domain/mbpt/models/cc.cpp | 109 +++++++++++++++++++++++++++++- SeQuant/domain/mbpt/models/cc.hpp | 34 +++++++++- 2 files changed, 139 insertions(+), 4 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 4793984a08..6b79a34d5a 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -1,7 +1,9 @@ +#include #include #include #include #include +#include #include #include #include @@ -16,6 +18,7 @@ #include #include #include +#include namespace { // alias reserved labels for readability @@ -91,7 +94,7 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { ExprPtr CC::energy(std::optional comm_rank) const { // Bernoulli: the tensor-level H̄ is already fully expanded, so there is no - // operator connectivity left to constrain -- take the plain reference + // operator connectivity left to constrain. Take the plain reference // expectation value. The energy rank defaults to the amplitude rank; // pass comm_rank explicitly for the qUCCSD [2|3] split (energy at H̄³). if (hbar_expansion_ == HbarExpansion::Bernoulli) { @@ -371,9 +374,103 @@ std::vector CC::λʼ(size_t rank, size_t order, namespace { // EOM eigenvector operators R and L use SquareRoot normalization constexpr Normalization eom_norm = Normalization::SquareRoot; + +// Per-block-truncated EOM sigma equations (qUCCSD and its IP/EA analogues: +// 10.1063/5.0062090 Table I, 10.1021/acs.jctc.5c01991 Table 1). +// +// Each block is the sandwich , plus an explicit -E shift on the +// diagonal. The commutator form would add -. That +// term is manifold j's amplitude residual. It vanishes only when the block rank +// k equals the rank the amplitudes were converged against. At the qUCCSD ranks +// it does not, and it hits DS but not SD, so M_ij != adjoint(M_ji). +// +// The shift is <0|H̄^(k_ii)|0>, at the block's own rank. Block truncation cuts +// the single operator H̄-E. Its order-m piece is H̄^m - <0|H̄^m|0>. On the qUCCSD +// DD block this leaves the doubles diagonal bare (Eq. 48 of 10.1063/5.0062090). +std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, + const std::vector& block_ranks, + size_t N) { + if (!cc.unitary()) + throw Exception("CC::eom_r: block_ranks requires a unitary ansatz"); + + std::vector> 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::reverse(manifolds.begin(), manifolds.end()); + const auto K = manifolds.size(); + if (block_ranks.size() != K * K) + throw Exception( + "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 = cc.hbar_expansion() == CC::HbarExpansion::Bernoulli; + + // One H̄ per distinct truncation order, with its N part removed. The N part + // holds the pure excitation / de-excitation intermediates H̄_ai, H̄_ab,ij of + // rank ≤ N, which are the ground-state amplitude residual. The amplitude + // equations zero them (⟨μ|H̄|Φ₀⟩ = 0, Liu & Cheng 2021 Eq. (6)), so the + // paper's off-diagonal working equations (Eqs. 41-47) carry none. That holds + // at the amplitude rank only. The off-diagonal blocks are built one rank + // lower, where ⟨μ|H̄^(k)|0⟩ ≠ 0. Keeping the N part there feeds a spurious + // off-shell residual into the S↔D coupling. + // + // Removing it from every block is exact and simpler. An N operator of rank r + // shifts the manifold rank by r, so it cannot reach a diagonal block. It also + // has no reference expectation value, so the -⟨0|H̄|0⟩ shift below stays the + // same. When every block rank equals the amplitude rank the removed terms are + // the converged residual, i.e. the step changes the equations but not the + // numbers. + container::map hbars; + for (const auto k : block_ranks) { + auto [it, fresh] = hbars.try_emplace(k); + if (!fresh) continue; // deriving H̄ twice for one rank is not cheap + it->second = cc.hbar(k); + // BCH H̄ is operator-level and has no N/R split. Block truncation under + // BCH is unpublished, so nothing validates a correction here. + 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, &cc](const ExprPtr& e) { + return tensor_level ? op::tensor::ref_av(e) + : op::ref_av(e, {.connect = {}, + .screen = cc.screen(), + .use_topology = cc.use_topology()}); + }; + + using std::min; + std::vector 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(); + for (size_t j = 0; j < K; ++j) { + const auto [kp, kh] = manifolds[j]; + const auto& hbar_ij = hbars.at(block_ranks[i * K + j]); + const auto ket = ket_of(kp, kh); + acc->append(vev(bra * hbar_ij * ket)); + // -<0|H̄^(k_ii)|0>, written as so Wick keeps E's summed + // indices disjoint from the block's external ones. + if (i == j) acc->append(ex(-1) * vev(bra * ket * hbar_ij)); + } + result.at(static_cast(min(bp, bh))) = simplify(ExprPtr{acc}); + } + return result; +} } // namespace -std::vector CC::eom_r(nₚ np, nₕ nh) const { +std::vector CC::eom_r(nₚ np, nₕ nh, + const std::vector& block_ranks) const { SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order"); if (np != nh) SEQUANT_ASSERT( @@ -384,6 +481,14 @@ std::vector CC::eom_r(nₚ np, nₕ nh) const { "hbar_comm_rank must be specified for unitary ansatz " "in CC::eom_r"); + if (!block_ranks.empty()) return eom_r_blocked(*this, np, nh, block_ranks, N); + + // the uniform path below commutes H̄ with an operator-level R, which the + // tensor-level Bernoulli H̄ cannot take part in + if (hbar_expansion_ == HbarExpansion::Bernoulli) + throw Exception( + "CC::eom_r: the Bernoulli expansion requires non-empty block_ranks"); + // construct hbar const auto hbar = this->hbar(); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index dd55202189..1e9a007fed 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -187,11 +187,41 @@ class CC { size_t rank = 1, size_t order = 1, std::optional nbatch = std::nullopt) const; + // clang-format off /// @brief derives right-side sigma equations for EOM-CC /// @param np number of particle creators in R operator /// @param nh number of hole creators in R operator - /// @return vector of right side sigma equations, element 0 is always null - [[nodiscard]] std::vector eom_r(nₚ np, nₕ nh) const; + /// @param block_ranks optional per-block H̄ commutator truncation ranks: a + /// different H̄ in each block of the secular matrix instead of one uniform + /// H̄ everywhere. For singles+doubles the matrix and its ranks are + /// | H_SS H_SD | qUCCSD: | 2 1 | + /// | H_DS H_DD | | 1 0 | + /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator + /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare f+v. + /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING + /// manifold rank, so one set of numbers serves EE, IP and EA (read S as + /// 1h/1p and D as 2h1p/1h2p: qUCCSD, IP-qUCCSD and EA-qUCCSD are all + /// `{2,1,1,0}`). Empty (the default) selects the uniform H̄ at + /// `hbar_comm_rank`, which the Bernoulli expansion does not support. + /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and + /// has nothing to truncate. + /// @throw Exception if `block_ranks` is neither empty nor `K`×`K`, if it is + /// non-empty under a non-unitary ansatz, or if it is empty under the + /// Bernoulli expansion + /// @note each block is the sandwich \f$ \langle i|\bar{H}|j \rangle \f$ + /// (Eq. 7 of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the + /// diagonal, taken at the block's own truncation rank; the returned object + /// is \f$ (\bar{H}-E)\hat{R} \f$. + /// @note under the Bernoulli expansion each block's H̄ has its N part (the + /// ground-state amplitude residual) removed. See `eom_r_blocked` in cc.cpp + /// for why. The removed terms vanish at converged amplitudes when a block + /// rank equals `hbar_comm_rank`, so this changes those blocks' equations + /// but not the numbers they evaluate to. + /// @return vector of right side sigma equations; element 0 is null iff + /// `np == nh` + // clang-format on + [[nodiscard]] std::vector eom_r( + nₚ np, nₕ nh, const std::vector& block_ranks = {}) const; /// @brief derives left-side sigma equations for EOM-CC /// @param np number of particle annihilators in L operator From 581d71eadc1973a155eeb2bc8858b6962ca98c83 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:10 -0400 Subject: [PATCH 15/37] test(mbpt): block-truncated qUCCSD EOM unit test Pin the term counts of the {2,1,1,0} EE and IP sigma equations under the Bernoulli expansion, and cover the three ways CC::eom_r rejects a block_ranks argument: a non-square matrix, a non-unitary ansatz, and an empty matrix under Bernoulli. --- tests/unit/test_mbpt_cc.cpp | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 2be3187f0a..daba3ed681 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -197,6 +197,37 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { #endif // !defined(SEQUANT_SKIP_LONG_TESTS) } + SECTION("bernoulli_quccsd_eom") { + using namespace sequant; + using namespace sequant::mbpt; + const CC cc(2, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 2, + .hbar_expansion = CC::HbarExpansion::Bernoulli}); + // qUCCSD block ranks, 10.1063/5.0062090 Table I: leading block at the + // double commutator, coupling blocks at the single, highest manifold bare. + const std::vector quccsd = {2, 1, 1, 0}; + + const auto ee = cc.eom_r(nₚ(2), nₕ(2), quccsd); + REQUIRE(ee.size() == 3); + REQUIRE(!ee[0]); + REQUIRE(size(ee[1]) == 121); + REQUIRE(size(ee[2]) == 21); + + // the same ranks drive IP: manifolds are indexed by ascending rank, so + // {1h, 2h1p} takes the place of {S, D} + const auto ip = cc.eom_r(nₚ(1), nₕ(2), quccsd); + REQUIRE(ip.size() == 2); + REQUIRE(size(ip[0]) == 32); + REQUIRE(size(ip[1]) == 11); + + // block_ranks must be a K x K matrix over the manifolds ... + REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2), {2, 1, 0}), Exception); + // ... the ansatz must be unitary ... + REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); + // ... and the Bernoulli H̄ has no uniform path to fall back on + REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2)), Exception); + } + SECTION("energy") { // CC::energy() must equal the p==0 element of CC::t() for both ansätze. const auto N = 2; From 89849b713f8bcd4a42a8de8874414945102ce0dc Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:11 -0400 Subject: [PATCH 16/37] chore: ignore the local developer setup files CMakeUserPresets.json is the documented per-developer companion to CMakePresets.json, and Notes is a symlink into a personal notes repo. --- .gitignore | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 1aac3de827..13d9e955be 100644 --- a/.gitignore +++ b/.gitignore @@ -71,4 +71,8 @@ _codeql_detected_source_root .clangd .vscode out -run \ No newline at end of file +run + +# local developer setup +CMakeUserPresets.json +Notes From a027df177f2af5b4fd72168672c7c1c60739abf6 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 4 Aug 2026 11:58:00 -0400 Subject: [PATCH 17/37] refactor(mbpt): tighten eom_r_blocked comments, restore the Bernoulli guard MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Shortens the block-truncation derivation comments in cc.cpp/cc.hpp, adds missing paper section/equation references, and switches to std::ranges::reverse. The previous version of this commit also dropped the explicit Bernoulli-empty-block_ranks throw in CC::eom_r, reasoning that CC::hbar's own ctor-time hbar_comm_rank check made it redundant. That conflated two independent preconditions: hbar_comm_rank being set says nothing about block_ranks being non-empty. Without the guard, CC::eom_r() under Bernoulli falls through to the uniform path, which commutes the tensor-level Bernoulli H̄ with an operator-level R -- exactly the mixing op.ipp's commutes_with_atom() assert exists to catch. Restored the throw; caught by running the existing bernoulli_quccsd_eom unit test, which this repo's validation history (Notes/ucc/bernoulli-ucc/PR.md) never actually re-ran after this change was first made. --- SeQuant/domain/mbpt/models/cc.cpp | 61 +++++++++++-------------------- SeQuant/domain/mbpt/models/cc.hpp | 16 +++++--- tests/unit/test_mbpt_cc.cpp | 4 +- 3 files changed, 33 insertions(+), 48 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 6b79a34d5a..d7138b76e2 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -93,12 +93,10 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { } ExprPtr CC::energy(std::optional comm_rank) const { - // Bernoulli: the tensor-level H̄ is already fully expanded, so there is no - // operator connectivity left to constrain. Take the plain reference - // expectation value. The energy rank defaults to the amplitude rank; - // pass comm_rank explicitly for the qUCCSD [2|3] split (energy at H̄³). + // 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_); + 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 @@ -113,10 +111,8 @@ std::vector CC::t(size_t pmax, size_t pmin) const { pmax = (pmax == std::numeric_limits::max() ? N : pmax); SEQUANT_ASSERT(pmax >= pmin && "pmax should be >= pmin"); - // Bernoulli: project the tensor-level H̄ (built at the amplitude - // rank = hbar_comm_rank_) onto each manifold hbar(); std::vector result(pmax + 1); @@ -375,23 +371,17 @@ namespace { // EOM eigenvector operators R and L use SquareRoot normalization constexpr Normalization eom_norm = Normalization::SquareRoot; -// Per-block-truncated EOM sigma equations (qUCCSD and its IP/EA analogues: -// 10.1063/5.0062090 Table I, 10.1021/acs.jctc.5c01991 Table 1). +// 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 Table 1. // -// Each block is the sandwich , plus an explicit -E shift on the -// diagonal. The commutator form would add -. That -// term is manifold j's amplitude residual. It vanishes only when the block rank -// k equals the rank the amplitudes were converged against. At the qUCCSD ranks -// it does not, and it hits DS but not SD, so M_ij != adjoint(M_ji). -// -// The shift is <0|H̄^(k_ii)|0>, at the block's own rank. Block truncation cuts -// the single operator H̄-E. Its order-m piece is H̄^m - <0|H̄^m|0>. On the qUCCSD -// DD block this leaves the doubles diagonal bare (Eq. 48 of 10.1063/5.0062090). +// Each block is the sandwich 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 eom_r_blocked(const CC& cc, nₚ np, nₕ nh, const std::vector& block_ranks, size_t N) { - if (!cc.unitary()) - throw Exception("CC::eom_r: block_ranks requires a unitary ansatz"); + if (!cc.unitary()) throw Exception("eom_r_blocked requires a unitary ansatz"); std::vector> manifolds; for (std::int64_t rp = np, rh = nh; rp >= 0 && rh >= 0; --rp, --rh) { @@ -399,7 +389,8 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, manifolds.emplace_back(rp, rh); if (rp == 0 || rh == 0) break; } - std::reverse(manifolds.begin(), manifolds.end()); + + std::ranges::reverse(manifolds); const auto K = manifolds.size(); if (block_ranks.size() != K * K) throw Exception( @@ -410,28 +401,17 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, // below must match it. Empty connectivity, as everywhere on the unitary path. const bool tensor_level = cc.hbar_expansion() == CC::HbarExpansion::Bernoulli; - // One H̄ per distinct truncation order, with its N part removed. The N part - // holds the pure excitation / de-excitation intermediates H̄_ai, H̄_ab,ij of - // rank ≤ N, which are the ground-state amplitude residual. The amplitude - // equations zero them (⟨μ|H̄|Φ₀⟩ = 0, Liu & Cheng 2021 Eq. (6)), so the - // paper's off-diagonal working equations (Eqs. 41-47) carry none. That holds - // at the amplitude rank only. The off-diagonal blocks are built one rank - // lower, where ⟨μ|H̄^(k)|0⟩ ≠ 0. Keeping the N part there feeds a spurious - // off-shell residual into the S↔D coupling. - // - // Removing it from every block is exact and simpler. An N operator of rank r - // shifts the manifold rank by r, so it cannot reach a diagonal block. It also - // has no reference expectation value, so the -⟨0|H̄|0⟩ shift below stays the - // same. When every block rank equals the amplitude rank the removed terms are - // the converged residual, i.e. the step changes the equations but not the - // numbers. + // One H̄ per distinct truncation order, reduced to its R part. The N part is + // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes at + // the amplitude rank only, so a block truncated below it would keep the + // residual. Dropping it everywhere is exact: an N operator of rank r shifts + // the manifold rank by r, so it never reaches a diagonal block, and it has + // no reference expectation value, so the shift below is unchanged. container::map hbars; for (const auto k : block_ranks) { auto [it, fresh] = hbars.try_emplace(k); if (!fresh) continue; // deriving H̄ twice for one rank is not cheap it->second = cc.hbar(k); - // BCH H̄ is operator-level and has no N/R split. Block truncation under - // BCH is unpublished, so nothing validates a correction here. 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) { @@ -481,6 +461,7 @@ std::vector CC::eom_r(nₚ np, nₕ nh, "hbar_comm_rank must be specified for unitary ansatz " "in CC::eom_r"); + // if block ranks are specified, dispatch and early return if (!block_ranks.empty()) return eom_r_blocked(*this, np, nh, block_ranks, N); // the uniform path below commutes H̄ with an operator-level R, which the diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 1e9a007fed..f4a487b717 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -197,21 +197,25 @@ class CC { /// | H_SS H_SD | qUCCSD: | 2 1 | /// | H_DS H_DD | | 1 0 | /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator - /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare f+v. + /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare f+v + /// (10.1063/5.0062090 Sec. II C, Eqs. (29), (41), (44), (48)). /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING /// manifold rank, so one set of numbers serves EE, IP and EA (read S as /// 1h/1p and D as 2h1p/1h2p: qUCCSD, IP-qUCCSD and EA-qUCCSD are all - /// `{2,1,1,0}`). Empty (the default) selects the uniform H̄ at - /// `hbar_comm_rank`, which the Bernoulli expansion does not support. + /// `{2,1,1,0}`, 10.1021/acs.jctc.5c01991 Table 1). Empty (the default) + /// selects the uniform H̄ at `hbar_comm_rank`, which the Bernoulli + /// expansion does not support. /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and /// has nothing to truncate. /// @throw Exception if `block_ranks` is neither empty nor `K`×`K`, if it is /// non-empty under a non-unitary ansatz, or if it is empty under the /// Bernoulli expansion /// @note each block is the sandwich \f$ \langle i|\bar{H}|j \rangle \f$ - /// (Eq. 7 of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the - /// diagonal, taken at the block's own truncation rank; the returned object - /// is \f$ (\bar{H}-E)\hat{R} \f$. + /// (Eq. (7) of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the + /// diagonal, taken at the block's own truncation rank. Eq. (10) there + /// instead splits \f$ \bar{H} = E_{gr} + {} \f$ a normal-ordered remainder + /// and forms the blocks from the remainder. The returned object is + /// \f$ (\bar{H}-E)\hat{R} \f$. /// @note under the Bernoulli expansion each block's H̄ has its N part (the /// ground-state amplitude residual) removed. See `eom_r_blocked` in cc.cpp /// for why. The removed terms vanish at converged amplitudes when a block diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index daba3ed681..d4fec07590 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -203,8 +203,8 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { const CC cc(2, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2, .hbar_expansion = CC::HbarExpansion::Bernoulli}); - // qUCCSD block ranks, 10.1063/5.0062090 Table I: leading block at the - // double commutator, coupling blocks at the single, highest manifold bare. + // qUCCSD block ranks, 10.1063/5.0062090 Sec. II C: SS at the double + // commutator (Eq. 29), SD/DS at the single (Eqs. 41, 44), DD bare (Eq. 48). const std::vector quccsd = {2, 1, 1, 0}; const auto ee = cc.eom_r(nₚ(2), nₕ(2), quccsd); From 4ff7a754f11c4c40033696b5915a17d6babbef27 Mon Sep 17 00:00:00 2001 From: Ajay Date: Wed, 5 Aug 2026 09:17:42 -0400 Subject: [PATCH 18/37] Revert .gitignore: drop unrelated local-dev-setup entries CMakeUserPresets.json and Notes are personal/local exclusions, not project-wide ignores, and don't belong in a Bernoulli-expansion feature PR. Use .git/info/exclude or a global excludesFile for these instead. --- .gitignore | 4 ---- 1 file changed, 4 deletions(-) diff --git a/.gitignore b/.gitignore index 13d9e955be..64fe8813e5 100644 --- a/.gitignore +++ b/.gitignore @@ -72,7 +72,3 @@ _codeql_detected_source_root .vscode out run - -# local developer setup -CMakeUserPresets.json -Notes From fd0c70f468db7505458c017a02511bfcac66dcc4 Mon Sep 17 00:00:00 2001 From: Ajay Date: Wed, 5 Aug 2026 10:02:54 -0400 Subject: [PATCH 19/37] docs(mbpt): clean up Bernoulli comments to defer to the paper Trim prose that paraphrased or re-derived the paper's math in favor of citing the equation it corresponds to. Fix a few inaccuracies found by re-checking against 10.1063/1.5030344 and 10.1063/5.0062090 directly: R_part was mislabeled a 'rank-preserving remainder' (it isn't, for rank >= 2 operators), the DD EOM block was missing its bare-Fock content (f_ij, f_ab), and the 'higher orders carry only R-subscripted inner commutators' claim doesn't hold for the V_N-seeded terms. Also drops the invented 'Cancellation #2' label and renames the orphaned 'Cancellation #1' to 'the F-cancellation' throughout. --- SeQuant/domain/mbpt/bernoulli.cpp | 53 ++++++++++++++++--------------- SeQuant/domain/mbpt/bernoulli.hpp | 9 +++--- SeQuant/domain/mbpt/models/cc.cpp | 11 ++++--- SeQuant/domain/mbpt/models/cc.hpp | 4 ++- tests/unit/test_mbpt_cc.cpp | 6 ++-- 5 files changed, 43 insertions(+), 40 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 37f1c103d2..de409310f7 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -22,22 +22,22 @@ #include // Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian -// H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). Because σ mixes excitation -// and de-excitation the plain BCH series does not terminate; the Bernoulli -// expansion rewrites it with Bernoulli numbers as the expansion coefficients, -// which leaves the truncation at a chosen commutator rank as the only -// approximation. H is split as F (Fock, rank-preserving) + V (fluctuation -// potential), and every operator O is split into O_N (all excitation and -// de-excitation operators) and O_R = O − O_N. At a converged RHF/UHF reference -// two cancellations hold: F survives only in H̄¹, and the higher orders carry -// only R-subscripted inner commutators. +// H̄ = e^{−σ} H e^{σ}, σ = T − T† (anti-Hermitian). For UCC the plain BCH +// series does not terminate, because σ mixes excitation and de-excitation. +// This file implements the Bernoulli-number resummation of 10.1063/1.5030344 +// that fixes that. H splits as F (Fock) + V (fluctuation potential); every +// operator O splits into O_N (its pure excitation/de-excitation part) and +// O_R = O − O_N. // -// All equation references are to 10.1063/1.5030344 (Sec. III B): superoperator -// inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, B₃=0, B₄=−1/720 -// Eq. (40); the N/R split and the UCC amplitude condition V̄_N = 0 above and at -// Eq. (43); the iterative recursion for V̄ Eq. (44); the assembly -// H̄ = Σ_k H̄^k Eq. (45); the rank-by-rank operators H̄⁰..H̄⁴ Eqs. (46)-(50). -// Cancellation #1 (F enters only H̄¹) is stated just below Eq. (50). +// Equation numbers below are all from 10.1063/1.5030344, Sec. III B: +// superoperator inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, +// B₃=0, B₄=−1/720, Eq. (40); the N/R split and the UCC amplitude condition +// V̄_N = 0, above and at Eq. (43); the H̄ recursion, Eq. (44); the sum +// H̄ = Σ_k H̄^k, Eq. (45); H̄⁰..H̄⁴, Eqs. (46)-(50). +// +// The F-cancellation: at a canonical (converged) HF reference, F has no +// occupied-virtual block (Eq. (32)), so F enters H̄ only through H̄¹ (stated +// just below Eq. (50)). Every H̄^k built below relies on this. namespace { @@ -69,17 +69,18 @@ const sequant::NormalOperator* find_nop( return nullptr; } -/// Classifies one block-resolved term as N or R (Cancellation #2). A term is N -/// iff its single residual NormalOperator is a pure excitation (all creators -/// pure-unoccupied AND all annihilators pure-occupied) or a pure de-excitation -/// (all creators pure-occupied AND all annihilators pure-unoccupied), with rank -/// ≤ @p cutoff. A term with no residual NormalOperator is rank-preserving, +/// Classifies one block-resolved term as N or R, per the O_N/O_R split above +/// Eq. (43) of 10.1063/1.5030344. A term is N iff its single residual +/// NormalOperator is a pure excitation (all creators pure-unoccupied AND all +/// annihilators pure-occupied) or a pure de-excitation (the reverse), with +/// rank ≤ @p cutoff. A term with no residual NormalOperator is rank-preserving, /// hence R. /// -/// Rank > @p cutoff falls to R rather than being dropped. An ]_R filter drops a -/// term only because the amplitude condition V̄_N = 0 (Eq. (43)) makes it zero, -/// and for σ truncated at rank N that condition covers rank ≤ N only. Eq. (43) -/// states O_N with no rank limit because there σ carries every rank. +/// Rank > @p cutoff falls to R rather than being dropped. The R filter drops a +/// term only because the amplitude condition V̄_N = 0 (Eq. (43)) makes it +/// zero, and that condition holds for rank ≤ N only, since σ here is +/// truncated at rank N. Eq. (43) itself states O_N with no rank limit, +/// because there σ carries every rank. bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { using namespace sequant; auto isr = get_default_context().index_space_registry(); @@ -362,7 +363,7 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { add(1, simplify(F + V)); // H̄⁰ = F + V [Eq. (46)] if (rank >= 1) { // H̄¹ = [F,σ] + ½[V,σ] + ½[V_R,σ] [Eq. (47)]. F enters H̄ ONLY here - // (Cancellation #1, stated just below Eq. (50)). + // (the F-cancellation, stated just below Eq. (50)). add(1, wick_commutator(F, sigma)); add({1, 2}, nest('A', "A")); add({1, 2}, nest('R', "A")); @@ -385,7 +386,7 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { if (rank >= 4) { // H̄⁴ = Eq. (50), the nine order-4 terms produced by the recursion Eq. (44), // V̄^{k+1} = σ̂F + X̂⁻¹(σ̂)e^{σ̂}V − Σ_{n≠0} B_n σ̂^n V̄_R^{k}. F is absent here - // (Cancellation #1). Listed in the paper's order; the outermost tag is + // (the F-cancellation). Listed in the paper's order; the outermost tag is // always A. add({1, 16}, nest('R', "RRRA")); add({1, 16}, nest('A', "RRRA")); diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index 4f924a2900..f7bdae34ea 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -55,11 +55,10 @@ ExprPtr expand_to_blocks(const ExprPtr& expr); /// @p cutoff. Applies expand_to_blocks first. ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff); -/// R (rank-preserving remainder) part: wick_reduce(expr) minus -/// N_part(expr, cutoff). Unlike N_part the result is NOT block-resolved. It -/// stays in compact general-index form. That is exact here, because -/// expand_to_blocks is an identity, and much cheaper for the nested commutators -/// that consume R. +/// 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. +/// This is exact because expand_to_blocks is an identity, and it is much +/// cheaper for the nested commutators that consume R. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff); } // namespace detail diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index d7138b76e2..104ba4581b 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -402,11 +402,12 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, const bool tensor_level = cc.hbar_expansion() == CC::HbarExpansion::Bernoulli; // One H̄ per distinct truncation order, reduced to its R part. The N part is - // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes at - // the amplitude rank only, so a block truncated below it would keep the - // residual. Dropping it everywhere is exact: an N operator of rank r shifts - // the manifold rank by r, so it never reaches a diagonal block, and it has - // no reference expectation value, so the shift below is unchanged. + // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes only + // at the amplitude rank. A block truncated below that rank would keep the + // residual, so drop N everywhere instead. This is exact: 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, so the −E shift below + // is unchanged either way. container::map hbars; for (const auto k : block_ranks) { auto [it, fresh] = hbars.try_emplace(k); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index f4a487b717..a6eb3aa1ab 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -197,7 +197,9 @@ class CC { /// | H_SS H_SD | qUCCSD: | 2 1 | /// | H_DS H_DD | | 1 0 | /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator - /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare f+v + /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare + /// Hamiltonian integrals (no commutators): f (Eqs. (30),(34)) plus + /// \f$ \langle ij\|kl\rangle,\ \langle ab\|cd\rangle \f$ (Eq. (48)) /// (10.1063/5.0062090 Sec. II C, Eqs. (29), (41), (44), (48)). /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING /// manifold rank, so one set of numbers serves EE, IP and EA (read S as diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index d4fec07590..177dfc1fd2 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -130,8 +130,8 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { using namespace sequant; using namespace sequant::mbpt; // Equation references are to 10.1063/1.5030344, Sec. III B. - // Cancellation #1: F appears only in H̄¹, so rank r − rank r−1 is F-free - // for r ≥ 2. + // The F-cancellation: F appears only in H̄¹, so rank r − rank r−1 is + // F-free for r ≥ 2. auto h0 = bernoulli::hbar(2, 0, false); auto h1 = bernoulli::hbar(2, 1, false); auto h2 = bernoulli::hbar(2, 2, false); @@ -214,7 +214,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE(size(ee[2]) == 21); // the same ranks drive IP: manifolds are indexed by ascending rank, so - // {1h, 2h1p} takes the place of {S, D} + // {1h, 2h1p} takes the place of {S, D} (10.1021/acs.jctc.5c01991 Table 1) const auto ip = cc.eom_r(nₚ(1), nₕ(2), quccsd); REQUIRE(ip.size() == 2); REQUIRE(size(ip[0]) == 32); From 9c9b67c686eeb0d2d2d363bb4eba6baa3053948f Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 00:56:17 -0400 Subject: [PATCH 20/37] refactor(mbpt): assert, not throw, on the block-EOM preconditions Ansatz, block_ranks shape and the Bernoulli/uniform-path mismatch are all caller errors on an internal precondition, matching how the rest of CC validates its configuration. --- SeQuant/domain/mbpt/models/cc.cpp | 16 +++++++--------- SeQuant/domain/mbpt/models/cc.hpp | 3 +-- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 104ba4581b..6b1cdc300a 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -3,7 +3,6 @@ #include #include #include -#include #include #include #include @@ -381,7 +380,7 @@ constexpr Normalization eom_norm = Normalization::SquareRoot; std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, const std::vector& block_ranks, size_t N) { - if (!cc.unitary()) throw Exception("eom_r_blocked requires a unitary ansatz"); + SEQUANT_ASSERT(cc.unitary(), "eom_r_blocked requires a unitary ansatz"); std::vector> manifolds; for (std::int64_t rp = np, rh = nh; rp >= 0 && rh >= 0; --rp, --rh) { @@ -392,10 +391,9 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, std::ranges::reverse(manifolds); const auto K = manifolds.size(); - if (block_ranks.size() != K * K) - throw Exception( - "CC::eom_r: block_ranks must be a K x K row-major matrix, " - "K = number of projection manifolds"); + SEQUANT_ASSERT(block_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. @@ -467,9 +465,9 @@ std::vector CC::eom_r(nₚ np, nₕ nh, // the uniform path below commutes H̄ with an operator-level R, which the // tensor-level Bernoulli H̄ cannot take part in - if (hbar_expansion_ == HbarExpansion::Bernoulli) - throw Exception( - "CC::eom_r: the Bernoulli expansion requires non-empty block_ranks"); + SEQUANT_ASSERT( + hbar_expansion_ != HbarExpansion::Bernoulli, + "CC::eom_r: the Bernoulli expansion requires non-empty block_ranks"); // construct hbar const auto hbar = this->hbar(); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index a6eb3aa1ab..7c6770a024 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -209,8 +209,7 @@ class CC { /// expansion does not support. /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and /// has nothing to truncate. - /// @throw Exception if `block_ranks` is neither empty nor `K`×`K`, if it is - /// non-empty under a non-unitary ansatz, or if it is empty under the + /// @pre `block_ranks` is either empty or `K`×`K`, and is non-empty under the /// Bernoulli expansion /// @note each block is the sandwich \f$ \langle i|\bar{H}|j \rangle \f$ /// (Eq. (7) of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the From 635d2162b8f0cd9f5a0d2ca99315a6bcc2130f76 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 00:56:35 -0400 Subject: [PATCH 21/37] docs(mbpt): say that the Bernoulli hbar is tensor-level and unscreened --- SeQuant/domain/mbpt/bernoulli.hpp | 4 ++++ SeQuant/domain/mbpt/models/cc.hpp | 3 +++ 2 files changed, 7 insertions(+) diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index f7bdae34ea..5a90f9d6b5 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -21,6 +21,10 @@ namespace sequant::mbpt::bernoulli { /// 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. +/// /// @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 diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 7c6770a024..127506d14a 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -125,6 +125,9 @@ class CC { /// the explicit form. For a unitary ansatz the reverse holds: H̄ is already /// self-contained, so connectivity must be left empty. See the "Using H̄ /// outside the CC class" section of the user guide. + /// @note Under `HbarExpansion::Bernoulli` the result is tensor-level, so it + /// takes `op::tensor` projectors and `op::tensor::ref_av`, not their `op` + /// counterparts, and it is unscreened: `screen` has no effect there. [[nodiscard]] ExprPtr hbar( std::optional truncation_rank = std::nullopt) const; From 737049f2bff9cc4f0aeaf2db1e4eb60ba560b77f Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:03:23 -0400 Subject: [PATCH 22/37] refactor(mbpt): take wick_reduce's input by const reference It simplified its argument in place, so the pointee the caller passed came back canonicalized; both in-library callers cloned first to avoid that. Clone inside instead. --- SeQuant/domain/mbpt/bernoulli.cpp | 7 ++++--- SeQuant/domain/mbpt/bernoulli.hpp | 3 ++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index de409310f7..823e1e9eb9 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -117,7 +117,8 @@ namespace sequant::mbpt::bernoulli { namespace detail { -ExprPtr wick_reduce(ExprPtr expr) { +ExprPtr wick_reduce(const ExprPtr& expr_in) { + auto expr = expr_in->clone(); simplify(expr); FWickTheorem wick{expr}; // use_topology defaults to ON, so it must be turned off explicitly. It keeps @@ -266,7 +267,7 @@ ExprPtr R_part_reduced(const ExprPtr& reduced, std::size_t cutoff) { /// header): after expansion every residual index is definite, so the N/R /// classifier can act on it. ExprPtr expand_to_blocks(const ExprPtr& expr_in) { - return expand_to_blocks_reduced(wick_reduce(expr_in->clone())); + return expand_to_blocks_reduced(wick_reduce(expr_in)); } /// N part of @p expr at truncation @p cutoff (see header): block-resolve, then @@ -282,7 +283,7 @@ ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff) { /// equals the fully block-resolved remainder, and the compact expr makes the /// nested commutators that consume R operate on far fewer terms. ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff) { - auto reduced = wick_reduce(expr->clone()); + auto reduced = wick_reduce(expr); return R_part_reduced(reduced, cutoff); } diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index 5a90f9d6b5..ff85c36a26 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -38,7 +38,8 @@ namespace detail { /// 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. -ExprPtr wick_reduce(ExprPtr expr); +/// @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 From 86a03c9c5eabfa19c4abb83c75f6366af6d105cd Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:04:16 -0400 Subject: [PATCH 23/37] docs(mbpt): correct the Bernoulli expansion comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The B_n listed at Eq. (40) are the paper's, i.e. B_n/n! in the textbook normalization; say so. F-cancellation needs Brillouin, not canonicality, and H̄⁰ does carry F. expand_to_blocks drops the other base spaces rather than finding them empty; the projection is what makes that harmless. --- SeQuant/domain/mbpt/bernoulli.cpp | 12 +++++++----- SeQuant/domain/mbpt/bernoulli.hpp | 18 +++++++++++------- SeQuant/domain/mbpt/models/cc.cpp | 12 +++++++----- 3 files changed, 25 insertions(+), 17 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 823e1e9eb9..a7f935460e 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -31,13 +31,15 @@ // // Equation numbers below are all from 10.1063/1.5030344, Sec. III B: // superoperator inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, -// B₃=0, B₄=−1/720, Eq. (40); the N/R split and the UCC amplitude condition -// V̄_N = 0, above and at Eq. (43); the H̄ recursion, Eq. (44); the sum +// B₃=0, B₄=−1/720, Eq. (40) -- these are Bₙ/n! in the textbook normalization, +// so they will not match a table of Bₙ; the N/R split and the UCC amplitude +// condition V̄_N = 0, above and at Eq. (43); the H̄ recursion, Eq. (44); the sum // H̄ = Σ_k H̄^k, Eq. (45); H̄⁰..H̄⁴, Eqs. (46)-(50). // -// The F-cancellation: at a canonical (converged) HF reference, F has no -// occupied-virtual block (Eq. (32)), so F enters H̄ only through H̄¹ (stated -// just below Eq. (50)). Every H̄^k built below relies on this. +// The F-cancellation: at an HF reference F has no occupied-virtual block +// (Brillouin, Eq. (32)), so H̄² and higher contain no F (stated just below +// Eq. (50)). H̄⁰ and H̄¹ do, and they are built from the full one-body operator, +// so their f_ov terms survive symbolically and vanish only on substitution. namespace { diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index ff85c36a26..28c5fff6d4 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -25,6 +25,10 @@ namespace sequant::mbpt::bernoulli { /// 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 @@ -49,10 +53,11 @@ ExprPtr wick_reduce(const ExprPtr& expr); 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), an -/// identity in the single-reference setting where the other base spaces are -/// empty. After expansion every residual index is definite so the -/// N/R classifier can act on it. Idempotent on block-resolved input. +/// 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. ExprPtr expand_to_blocks(const ExprPtr& expr); /// Block-resolved N part (O_N of 10.1063/1.5030344): the terms whose single @@ -61,9 +66,8 @@ ExprPtr expand_to_blocks(const ExprPtr& expr); 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. -/// This is exact because expand_to_blocks is an identity, and it is much -/// cheaper for the nested commutators that consume R. +/// 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 diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 6b1cdc300a..db5bdab493 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -402,12 +402,14 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, // One H̄ per distinct truncation order, reduced to its R part. The N part is // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes only // at the amplitude rank. A block truncated below that rank would keep the - // residual, so drop N everywhere instead. This is exact: 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, so the −E shift below - // is unchanged either way. + // residual, so drop N everywhere instead. The diagonal is untouched: 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, so the −E shift + // below is unchanged either way. Off the diagonal the equations do change, + // though not what they evaluate to at converged amplitudes. Bernoulli only: + // the operator-level BCH H̄ has no N/R split to take. container::map hbars; - for (const auto k : block_ranks) { + 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 = cc.hbar(k); From 1341008f37e7f975b530384ed75aa51d52fbc5f0 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 11 Aug 2026 14:42:39 -0400 Subject: [PATCH 24/37] =?UTF-8?q?fix(mbpt):=20reject=20the=20Bernoulli=20e?= =?UTF-8?q?xpansion=20in=20CC::t=CA=BC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit tʼ builds an operator-level similarity transform, which a tensor-level H̄ cannot enter; it reached the mismatch and aborted deep inside op algebra. --- SeQuant/domain/mbpt/models/cc.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index db5bdab493..cf8c5cb75b 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -253,6 +253,8 @@ std::vector CC::tʼ(size_t rank, size_t order, "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 From b87e715df3fd439f43d86c546bbb056477f98416 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:05:09 -0400 Subject: [PATCH 25/37] test(mbpt): pin coefficients, cross-check the blocked path, guard the assert tests Term counts are blind to the coefficients, so pin one projected equation in full. For a single manifold the blocked path must reproduce the uniform one, which shares no code with it. REQUIRE_THROWS_AS on a SEQUANT_ASSERT only holds in a THROW build. --- tests/unit/test_mbpt_cc.cpp | 47 ++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 177dfc1fd2..f6e3a51b04 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -190,6 +190,32 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE(size(amps[1]) == 32); REQUIRE(size(amps[2]) == 38); + // one projected equation pinned in full: term counts are blind to the + // coefficients, which is where a mis-weighted Wick reduction shows up. + // Doubles at H̄¹, the smallest such equation. + const auto R2_h1 = CC(2, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 1, + .hbar_expansion = CC::HbarExpansion::Bernoulli}) + .t() + .at(2); + REQUIRE_THAT( + R2_h1, EquivalentTo( + L"1/2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;a_3,a_4}:A-C-S " + L"* t{a_3,a_4;i_1,i_2}:A-N-S " + L"+ Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;i_1,i_2}:A-C-S " + L"+ 1/2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,i_4;i_1,i_2}:A-C-S " + L"* t{a_1,a_2;i_3,i_4}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * f{i_3;i_1}:A-C-S " + L"* t{a_1,a_2;i_2,i_3}:A-N-S " + L"- 2 Â{i_1,i_2;a_1,a_2}:A-C-S * f{a_1;a_3}:A-C-S " + L"* t{a_2,a_3;i_1,i_2}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;i_1,a_3}:A-C-S " + L"* t{a_3;i_2}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,a_1;i_1,i_2}:A-C-S " + L"* t{a_2;i_3}:A-N-S " + L"- 4 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,a_1;i_1,a_3}:A-C-S " + L"* t{a_2,a_3;i_2,i_3}:A-N-S")); + #ifndef SEQUANT_SKIP_LONG_TESTS const auto E = cc.energy(3); REQUIRE_THAT(E, !EquivalentTo(amps.at(0))); @@ -220,12 +246,21 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE(size(ip[0]) == 32); REQUIRE(size(ip[1]) == 11); - // block_ranks must be a K x K matrix over the manifolds ... - REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2), {2, 1, 0}), Exception); - // ... the ansatz must be unitary ... - REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); - // ... and the Bernoulli H̄ has no uniform path to fall back on - REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2)), Exception); + // one manifold means one block, whose sandwich minus shift is the + // commutator the uniform path builds; that path shares no code with + // eom_r_blocked, so this pins the block construction against it + const CC bch(2, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}); + REQUIRE_THAT(bch.eom_r(nₚ(1), nₕ(1), {2}).at(1), + EquivalentTo(bch.eom_r(nₚ(1), nₕ(1)).at(1))); + + if (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { + // block_ranks must be a K x K matrix over the manifolds ... + REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2), {2, 1, 0}), Exception); + // ... the ansatz must be unitary ... + REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); + // ... and the Bernoulli H̄ has no uniform path to fall back on + REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2)), Exception); + } } SECTION("energy") { From dbc5106bbb318514c771ac3dfb02efb9f241381c Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:05:27 -0400 Subject: [PATCH 26/37] feat(mbpt): allow uniform-rank EOM under the Bernoulli expansion MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit An empty block_ranks now fills every block with hbar_comm_rank instead of being rejected: the restriction was only that the operator-level commutator path cannot take a tensor-level H̄, which the blocked path handles. Index the matrix with .at() so a wrong size cannot read past the end. --- SeQuant/domain/mbpt/models/cc.cpp | 22 ++++++++++++---------- SeQuant/domain/mbpt/models/cc.hpp | 6 ++---- tests/unit/test_mbpt_cc.cpp | 8 +++++--- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index cf8c5cb75b..b0127f9e1d 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -393,7 +393,12 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, std::ranges::reverse(manifolds); const auto K = manifolds.size(); - SEQUANT_ASSERT(block_ranks.size() == K * K, + // empty means uniform truncation at hbar_comm_rank in every block + const std::vector ranks = + block_ranks.empty() + ? std::vector(K * K, cc.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"); @@ -439,7 +444,7 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, auto acc = std::make_shared(); for (size_t j = 0; j < K; ++j) { const auto [kp, kh] = manifolds[j]; - const auto& hbar_ij = hbars.at(block_ranks[i * K + 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 so Wick keeps E's summed @@ -464,14 +469,11 @@ std::vector CC::eom_r(nₚ np, nₕ nh, "hbar_comm_rank must be specified for unitary ansatz " "in CC::eom_r"); - // if block ranks are specified, dispatch and early return - if (!block_ranks.empty()) return eom_r_blocked(*this, np, nh, block_ranks, N); - - // the uniform path below commutes H̄ with an operator-level R, which the - // tensor-level Bernoulli H̄ cannot take part in - SEQUANT_ASSERT( - hbar_expansion_ != HbarExpansion::Bernoulli, - "CC::eom_r: the Bernoulli expansion requires non-empty block_ranks"); + // 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(*this, np, nh, block_ranks, N); // construct hbar const auto hbar = this->hbar(); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 127506d14a..a9722f31f0 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -208,12 +208,10 @@ class CC { /// manifold rank, so one set of numbers serves EE, IP and EA (read S as /// 1h/1p and D as 2h1p/1h2p: qUCCSD, IP-qUCCSD and EA-qUCCSD are all /// `{2,1,1,0}`, 10.1021/acs.jctc.5c01991 Table 1). Empty (the default) - /// selects the uniform H̄ at `hbar_comm_rank`, which the Bernoulli - /// expansion does not support. + /// selects the uniform H̄ at `hbar_comm_rank` everywhere. /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and /// has nothing to truncate. - /// @pre `block_ranks` is either empty or `K`×`K`, and is non-empty under the - /// Bernoulli expansion + /// @pre `block_ranks` is either empty or `K`×`K` /// @note each block is the sandwich \f$ \langle i|\bar{H}|j \rangle \f$ /// (Eq. (7) of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the /// diagonal, taken at the block's own truncation rank. Eq. (10) there diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index f6e3a51b04..a637e3b0d7 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -256,11 +256,13 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { if (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { // block_ranks must be a K x K matrix over the manifolds ... REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2), {2, 1, 0}), Exception); - // ... the ansatz must be unitary ... + // ... and the ansatz must be unitary REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); - // ... and the Bernoulli H̄ has no uniform path to fall back on - REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2)), Exception); } + + // no matrix means uniform truncation at hbar_comm_rank + REQUIRE_THAT(cc.eom_r(nₚ(2), nₕ(2)).at(1), + EquivalentTo(cc.eom_r(nₚ(2), nₕ(2), {2, 2, 2, 2}).at(1))); } SECTION("energy") { From 835a264302959aead0d3f8546c0482979a9fca70 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:16:49 -0400 Subject: [PATCH 27/37] fix(mbpt): assert the block expansion found a base space to expand into The hole/particle candidate list falls back to all base spaces when empty; if that is empty too, the assignment loop indexed an empty vector and emitted a term carrying a garbage index space. --- SeQuant/domain/mbpt/bernoulli.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index a7f935460e..8c8de22b77 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -201,6 +201,9 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { if (physical(b)) c.push_back(b); } choices.push_back(c.empty() ? c_all : c); + SEQUANT_ASSERT(!choices.back().empty(), + "bernoulli: general index spans no base space with " + "matching quantum numbers"); } // cartesian product of assignments => sum of transformed terms; // accumulate via Sum::append (linear) rather than operator+, which From d022fa6564db7c94818fd15ee950801ff91dcda5 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 6 Aug 2026 02:22:32 -0400 Subject: [PATCH 28/37] docs(mbpt): trim the Bernoulli comments Cut the debugging narrative, the restated code, and the per-rank node counts; tighten the paragraphs that stay. --- SeQuant/domain/mbpt/bernoulli.cpp | 43 ++++++++++++------------------- SeQuant/domain/mbpt/models/cc.cpp | 17 +++++------- tests/unit/test_mbpt_cc.cpp | 4 +-- 3 files changed, 25 insertions(+), 39 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 8c8de22b77..e240cfbd15 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -31,9 +31,9 @@ // // Equation numbers below are all from 10.1063/1.5030344, Sec. III B: // superoperator inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, -// B₃=0, B₄=−1/720, Eq. (40) -- these are Bₙ/n! in the textbook normalization, -// so they will not match a table of Bₙ; the N/R split and the UCC amplitude -// condition V̄_N = 0, above and at Eq. (43); the H̄ recursion, Eq. (44); the sum +// B₃=0, B₄=−1/720, Eq. (40), i.e. Bₙ/n! in the textbook normalization; the N/R +// split and the UCC amplitude condition V̄_N = 0, above and at Eq. (43); the H̄ +// recursion, Eq. (44); the sum // H̄ = Σ_k H̄^k, Eq. (45); H̄⁰..H̄⁴, Eqs. (46)-(50). // // The F-cancellation: at an HF reference F has no occupied-virtual block @@ -88,7 +88,7 @@ bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { auto isr = get_default_context().index_space_registry(); const auto* nop = find_nop(term); - if (!nop) return false; // no residual operator => rank-preserving => R + if (!nop) return false; const auto ncre = ranges::distance(nop->creators()); const auto nann = ranges::distance(nop->annihilators()); @@ -127,8 +127,7 @@ ExprPtr wick_reduce(const ExprPtr& expr_in) { // one representative per symmetry-equivalent contraction class and multiplies // by the class size, weight bookkeeping that holds only on the // fully-contracted path. On this partial-contraction path it silently - // rescales the terms carrying a symmetric amplitude pair, and the damage - // shows up only under projection. + // rescales the terms carrying a symmetric amplitude pair. wick.use_topology(false).full_contractions(false); auto result = wick.compute(/*count_only=*/false, /*skip_input_canonicalization=*/true); @@ -137,11 +136,9 @@ ExprPtr wick_reduce(const ExprPtr& expr_in) { } ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { - // A and B are built independently, so their summed indices are local to each. - // If both use the same labels (a block-resolved R/N part and sigma both carry - // a/i), the product A*B fuses two independent summations. That corrupts the - // contraction. Reindex B to fresh temporaries. Canonicalization restores tidy - // labels. + // A and B are built independently, so shared labels (both a block-resolved + // part and sigma carry a/i) would fuse two independent summations in A*B. + // Reindex B to fresh temporaries; canonicalization restores tidy labels. container::map repl; for (const auto& idx : get_used_indices(B)) repl.emplace(idx, Index::make_tmp_index(idx.space())); @@ -213,12 +210,9 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { for (;;) { container::map repl; for (std::size_t k = 0; k < gens.size(); ++k) - // fresh ordinal (not gens[k].ordinal()): reusing the general index's - // ordinal would collide with any pre-existing definite index of the - // same base space and ordinal already in the term (e.g. an a/i index - // from an amplitude in a commutator result), producing a duplicate - // index. A globally-unique temporary is disjoint by construction; - // canonicalization restores tidy labels. + // fresh ordinal (not gens[k].ordinal()): reusing it would collide with + // a definite index of the same base space already in the term, e.g. + // one an amplitude brought in. Canonicalization restores tidy labels. repl.emplace(gens[k], Index::make_tmp_index(choices[k][idx[k]])); sum->append(transform_expr(term, repl)); // increment mixed-radix counter over the assignments @@ -232,8 +226,6 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { return ExprPtr{sum}; }; - // transform_sum_expr maps in parallel, canonicalizes each result, and - // accumulates into a HashingAccumulator ExprPtr out; if (expr.is()) { out = transform_sum_expr(expr.as().summands(), expand_term); @@ -307,18 +299,15 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { const auto F = op::tensor::F(); const auto V = op::tensor::h(2); const auto T = op::tensor::T(N, skip1); - const auto sigma = simplify(T - adjoint(T)); // σ = T − T† + const auto sigma = simplify(T - adjoint(T)); // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} // with a per-level N/R/A partition tag applied after each commutator ('A' = // no filter). nest(p0, f) evaluates such a node, memoizing every prefix - // (key = p0 + tags applied so far): the terms share prefixes both within a - // rank (the 9 rank-4 terms have only 3 distinct level-1 and 6 level-2 nodes) - // and across ranks (all four ranks share the same three level-1 nodes), so - // the memo avoids recomputing them. Reusing a memoized ExprPtr is safe: - // expression composition deep-copies operands (Product/Sum append clone), so - // wick_commutator does not mutate its arguments. Commutator outputs are - // already wick_reduce'd, so the reduced-input N/R filters apply. + // (key = p0 + tags applied so far): prefixes repeat both within a rank and + // across ranks. Reuse is safe because expression composition deep-copies its + // operands. Commutator outputs are already wick_reduce'd, so the + // reduced-input N/R filters apply. container::map memo; auto nest = [&](char p0, const char* f) -> ExprPtr { SEQUANT_ASSERT((p0 == 'A' || p0 == 'N' || p0 == 'R') && diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index b0127f9e1d..d7b8e3647f 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -56,8 +56,6 @@ CC::CC(size_t n, const Options& opts) if (hbar_expansion_ == HbarExpansion::Bernoulli) { SEQUANT_ASSERT(unitary(), "CC: Bernoulli expansion requires a unitary ansatz"); - // without hbar_comm_rank CC::hbar() falls back to rank 4, silently - // selecting the most expensive (and least exercised) order SEQUANT_ASSERT(hbar_comm_rank_, "CC: Bernoulli expansion requires hbar_comm_rank"); } @@ -406,15 +404,14 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, // below must match it. Empty connectivity, as everywhere on the unitary path. const bool tensor_level = cc.hbar_expansion() == CC::HbarExpansion::Bernoulli; - // One H̄ per distinct truncation order, reduced to its R part. The N part is + // 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. A block truncated below that rank would keep the - // residual, so drop N everywhere instead. The diagonal is untouched: 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, so the −E shift - // below is unchanged either way. Off the diagonal the equations do change, - // though not what they evaluate to at converged amplitudes. Bernoulli only: - // the operator-level BCH H̄ has no N/R split to take. + // 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 the equations do change, + // though not what they evaluate to at converged amplitudes. container::map hbars; for (const auto k : ranks) { auto [it, fresh] = hbars.try_emplace(k); diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index a637e3b0d7..1a7690ef93 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -67,7 +67,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { using namespace sequant::mbpt; // [V, T2] is antisymmetric: [A,B] == -[B,A] after Wick reduction const auto V = op::tensor::h(2); - const auto T2 = op::tensor::t(2); // rank-2 excitation, tensor form + const auto T2 = op::tensor::t(2); const auto ab = bernoulli::detail::wick_commutator(V, T2); const auto ba = bernoulli::detail::wick_commutator(T2, V); REQUIRE_THAT(ab, EquivalentTo(simplify(ex(-1) * ba))); @@ -109,7 +109,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { SECTION("bernoulli_N_R_split") { using namespace sequant; using namespace sequant::mbpt; - const auto V = op::tensor::h(2); // fluctuation potential g (general) + const auto V = op::tensor::h(2); // general g const auto Vn = bernoulli::detail::N_part(V, 2); const auto Vr = bernoulli::detail::R_part(V, 2); // N ⊎ R reconstructs V. R stays in compact general-index form, so check the From d8a582aeb85576fb844dd6dd8517e816bafa92e9 Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 14:54:28 -0400 Subject: [PATCH 29/37] refactor(mbpt): make eom_r_blocked a private CC member It read six accessors off the CC it was handed and needed N threaded in by hand because N is private, so it was already a member in all but name. Taking it as one drops the cc and N parameters. Also add the include std::ranges::reverse needs; it was only ever reaching it transitively. --- SeQuant/domain/mbpt/models/cc.cpp | 27 +++++++++++++-------------- SeQuant/domain/mbpt/models/cc.hpp | 7 +++++++ 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index d7b8e3647f..3a210b7415 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -12,6 +12,7 @@ #include #include +#include #include #include #include @@ -369,6 +370,7 @@ std::vector CC::λʼ(size_t rank, size_t order, namespace { // EOM eigenvector operators R and L use SquareRoot normalization constexpr Normalization eom_norm = Normalization::SquareRoot; +} // namespace // 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, @@ -377,10 +379,9 @@ constexpr Normalization eom_norm = Normalization::SquareRoot; // Each block is the sandwich 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 eom_r_blocked(const CC& cc, nₚ np, nₕ nh, - const std::vector& block_ranks, - size_t N) { - SEQUANT_ASSERT(cc.unitary(), "eom_r_blocked requires a unitary ansatz"); +std::vector CC::eom_r_blocked( + nₚ np, nₕ nh, const std::vector& block_ranks) const { + SEQUANT_ASSERT(unitary(), "eom_r_blocked requires a unitary ansatz"); std::vector> manifolds; for (std::int64_t rp = np, rh = nh; rp >= 0 && rh >= 0; --rp, --rh) { @@ -393,16 +394,15 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, const auto K = manifolds.size(); // empty means uniform truncation at hbar_comm_rank in every block const std::vector ranks = - block_ranks.empty() - ? std::vector(K * K, cc.hbar_comm_rank().value()) - : block_ranks; + block_ranks.empty() ? std::vector(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 = cc.hbar_expansion() == CC::HbarExpansion::Bernoulli; + 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 @@ -416,7 +416,7 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, 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 = cc.hbar(k); + 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) { @@ -426,11 +426,11 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, return tensor_level ? op::tensor::r(nₚ(p), nₕ(h), eom_norm) : op::r(nₚ(p), nₕ(h), eom_norm); }; - auto vev = [tensor_level, &cc](const ExprPtr& e) { + auto vev = [tensor_level, this](const ExprPtr& e) { return tensor_level ? op::tensor::ref_av(e) : op::ref_av(e, {.connect = {}, - .screen = cc.screen(), - .use_topology = cc.use_topology()}); + .screen = screen_, + .use_topology = use_topology_}); }; using std::min; @@ -452,7 +452,6 @@ std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, } return result; } -} // namespace std::vector CC::eom_r(nₚ np, nₕ nh, const std::vector& block_ranks) const { @@ -470,7 +469,7 @@ std::vector CC::eom_r(nₚ np, nₕ nh, // 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(*this, np, nh, block_ranks, N); + return eom_r_blocked(np, nh, block_ranks); // construct hbar const auto hbar = this->hbar(); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index a9722f31f0..ae306e5825 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -245,6 +245,13 @@ class CC { std::optional pertbar_comm_rank_ = std::nullopt; HbarExpansion hbar_expansion_ = HbarExpansion::BCH; + /// @brief `eom_r`'s per-block-truncated path, taken whenever `block_ranks` is + /// non-empty or the expansion is Bernoulli + /// @param block_ranks see `eom_r`; empty means uniform `hbar_comm_rank` + /// @pre a unitary ansatz + [[nodiscard]] std::vector eom_r_blocked( + nₚ np, nₕ nh, const std::vector& block_ranks) const; + /// @return the `LSTOptions` this engine uses for every `mbpt::lst()` call /// @note The choice of commutator representation is really a question of /// whether the caller supplies operator connectivity downstream; for this From 37b7b1a02b51bb1b5cd4de59286c4dead69b6a20 Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 14:55:12 -0400 Subject: [PATCH 30/37] refactor(mbpt): require hole and particle spaces in the block expansion The candidate-space loop kept a second vector to fall back on when the registry defined no hole/particle split. Every registry in convention.cpp sets is_hole and is_particle, so that branch was unreachable. Use the throwing accessors instead and drop the fallback; the assert now names what is actually required. --- SeQuant/domain/mbpt/bernoulli.cpp | 25 ++++++++++++------------- SeQuant/domain/mbpt/bernoulli.hpp | 1 + 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index e240cfbd15..9d3715953a 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -166,10 +166,10 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { // in those are annihilated by the single-reference projection onto the // hole/particle manifolds, so dropping them changes no projected quantity. // This keeps the expansion 2-way per index instead of 4-way, which otherwise - // compounds across the nested commutators. Falls back to all base spaces if - // the registry defines no hole/particle split. - const auto& hole_t = isr->hole_space(/*nulltype_ok=*/true); - const auto& particle_t = isr->particle_space(/*nulltype_ok=*/true); + // compounds across the nested commutators. Both spaces are required; the + // accessors throw if the registry leaves either unspecified. + const auto& hole_t = isr->hole_space(); + const auto& particle_t = isr->particle_space(); auto physical = [&](const IndexSpace& b) { return hole_t.includes(b.type()) || particle_t.includes(b.type()); }; @@ -191,16 +191,15 @@ ExprPtr expand_to_blocks_reduced(const ExprPtr& expr) { // match (stay within the same spin sector). container::svector> choices; for (const auto& g : gens) { - container::svector c, c_all; + container::svector c; for (const auto& b : bases) - if (b.qns() == g.space().qns() && g.space().type().includes(b.type())) { - c_all.push_back(b); - if (physical(b)) c.push_back(b); - } - choices.push_back(c.empty() ? c_all : c); - SEQUANT_ASSERT(!choices.back().empty(), - "bernoulli: general index spans no base space with " - "matching quantum numbers"); + if (physical(b) && b.qns() == g.space().qns() && + g.space().type().includes(b.type())) + c.push_back(b); + SEQUANT_ASSERT(!c.empty(), + "bernoulli: general index spans no hole/particle base " + "space with matching quantum numbers"); + choices.push_back(std::move(c)); } // cartesian product of assignments => sum of transformed terms; // accumulate via Sum::append (linear) rather than operator+, which diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp index 28c5fff6d4..70e950c795 100644 --- a/SeQuant/domain/mbpt/bernoulli.hpp +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -58,6 +58,7 @@ ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B); /// 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 From e4155d650b926cb8322ab1d80499105ce84b79ec Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 14:55:37 -0400 Subject: [PATCH 31/37] style(mbpt): pass the assert message as SEQUANT_ASSERT's second argument SEQUANT_ASSERT(EXPR, ...) takes the message itself, so the legacy `cond && "msg"` idiom leaves the string in the stringified condition. --- SeQuant/domain/mbpt/bernoulli.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 9d3715953a..f639aa7370 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -60,7 +60,7 @@ const sequant::NormalOperator* find_nop( // The one-residual-operator invariant is load-bearing: N/R // classification reads this operator alone, so a second one would be // silently ignored and misclassify the term. - SEQUANT_ASSERT(!found && + SEQUANT_ASSERT(!found, "find_nop: term carries >1 NormalOperator; wick_reduce " "is expected to leave at most one residual operator"); found = &f.as>(); From a9f627cbaec0d24fc83db38eeec5e873eb74048d Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 14:56:21 -0400 Subject: [PATCH 32/37] refactor(mbpt): fold the Bernoulli partition tags into one lambda hbar's nest() carried the tag validation twice and the tag-to-filter cascade twice, once for the V base and once per nesting level. One `part` lambda covers both, taking a flag for whether the input is already wick_reduce'd. Rename the loop variable to `cur`; it was called `op`, which shadows the mbpt::op namespace the surrounding function uses. --- SeQuant/domain/mbpt/bernoulli.cpp | 40 +++++++++++++++---------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index f639aa7370..20c022abd7 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -300,43 +300,43 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { const auto T = op::tensor::T(N, skip1); const auto sigma = simplify(T - adjoint(T)); + // Applies one partition tag to an expression; 'A' is no filter. `reduced` + // says the input is already wick_reduce'd, which every commutator output is; + // the V base is not, so it takes the reducing form. + auto part = [&](char tag, const ExprPtr& e, bool reduced) -> ExprPtr { + SEQUANT_ASSERT(tag == 'A' || tag == 'N' || tag == 'R', + "bernoulli::hbar: partition tag must be one of A, N, R"); + if (tag == 'A') return e; + if (tag == 'N') + return reduced ? N_part_reduced(e, cutoff) : N_part(e, cutoff); + return reduced ? R_part_reduced(e, cutoff) : R_part(e, cutoff); + }; + // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} // with a per-level N/R/A partition tag applied after each commutator ('A' = // no filter). nest(p0, f) evaluates such a node, memoizing every prefix // (key = p0 + tags applied so far): prefixes repeat both within a rank and // across ranks. Reuse is safe because expression composition deep-copies its - // operands. Commutator outputs are already wick_reduce'd, so the - // reduced-input N/R filters apply. + // operands. container::map memo; auto nest = [&](char p0, const char* f) -> ExprPtr { - SEQUANT_ASSERT((p0 == 'A' || p0 == 'N' || p0 == 'R') && - "bernoulli::hbar: partition tag must be one of A, N, R"); // grow `key` in place rather than deriving it from the memo iterator: // container::map is a flat_map, whose insertions invalidate iterators std::string key{p0}; auto it = memo.find(key); - if (it == memo.end()) { - ExprPtr base = (p0 == 'N') ? N_part(V, cutoff) - : (p0 == 'R') ? R_part(V, cutoff) - : V; - it = memo.emplace(key, std::move(base)).first; - } - ExprPtr op = it->second; + if (it == memo.end()) + it = memo.emplace(key, part(p0, V, /*reduced=*/false)).first; + ExprPtr cur = it->second; for (int i = 0; f[i] != '\0'; ++i) { - SEQUANT_ASSERT((f[i] == 'A' || f[i] == 'N' || f[i] == 'R') && - "bernoulli::hbar: partition tag must be one of A, N, R"); key += f[i]; it = memo.find(key); if (it == memo.end()) { - auto cx = wick_commutator(op, sigma); - ExprPtr filtered = (f[i] == 'R') ? R_part_reduced(cx, cutoff) - : (f[i] == 'N') ? N_part_reduced(cx, cutoff) - : cx; - it = memo.emplace(key, std::move(filtered)).first; + auto cx = wick_commutator(cur, sigma); + it = memo.emplace(key, part(f[i], cx, /*reduced=*/true)).first; } - op = it->second; + cur = it->second; } - return op; + return cur; }; HashingAccumulator acc; From 72d448f01ec32b23a5128ab18f8945f188d0877a Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 14:57:19 -0400 Subject: [PATCH 33/37] docs(mbpt): correct the Bernoulli and block-rank comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Checked against the cited papers: - "F enters H̄ ONLY here" on H̄¹ contradicted the H̄⁰ = F + V line three lines above. Only the F-commutators truncate at the first power of σ, which is what 10.1063/1.5030344 states below Eq. (50). - The {2,1,1,0} ranks were credited to Table 1 of 10.1021/acs.jctc.5c01991, which lists which H̄ components enter each block and no ranks at all. The ranks are the superscripts in its Fig. 1. Cite each for what it carries, here and in cc.cpp and the unit test. - eom_r no longer presents {2,1,1,0} as "qUCCSD". It reproduces the published method only under the Bernoulli expansion, since BCH keeps the N part the paper drops. Give the numbers as the ranks that paper truncates qUCCSD at, and name its own UCCSD[2|2,1,0] notation. - Drop the H_DD integral list, which was missing , and point at the paper instead of restating it. --- SeQuant/domain/mbpt/bernoulli.cpp | 6 ++++-- SeQuant/domain/mbpt/models/cc.cpp | 3 ++- SeQuant/domain/mbpt/models/cc.hpp | 22 ++++++++++++---------- tests/unit/test_mbpt_cc.cpp | 3 ++- 4 files changed, 20 insertions(+), 14 deletions(-) diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp index 20c022abd7..1380ccf96f 100644 --- a/SeQuant/domain/mbpt/bernoulli.cpp +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -356,8 +356,10 @@ ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1) { add(1, simplify(F + V)); // H̄⁰ = F + V [Eq. (46)] if (rank >= 1) { - // H̄¹ = [F,σ] + ½[V,σ] + ½[V_R,σ] [Eq. (47)]. F enters H̄ ONLY here - // (the F-cancellation, stated just below Eq. (50)). + // H̄¹ = [F,σ] + ½[V,σ] + ½[V_R,σ] [Eq. (47)]. F-commutators enter H̄ ONLY + // here. This is the F-cancellation, stated just below Eq. (50) as "the + // terms in H̄ involving F now truncate to the first power of σ". H̄⁰ carries + // the bare F. add(1, wick_commutator(F, sigma)); add({1, 2}, nest('A', "A")); add({1, 2}, nest('R', "A")); diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 3a210b7415..bee354d38c 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -374,7 +374,8 @@ constexpr Normalization eom_norm = Normalization::SquareRoot; // 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 Table 1. +// 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 of Eq. (7). Eq. (10) writes H̄ as // E_gr + a normal-ordered remainder and builds the blocks from the remainder diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index ae306e5825..28c9f1cdc2 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -196,19 +196,20 @@ class CC { /// @param nh number of hole creators in R operator /// @param block_ranks optional per-block H̄ commutator truncation ranks: a /// different H̄ in each block of the secular matrix instead of one uniform - /// H̄ everywhere. For singles+doubles the matrix and its ranks are - /// | H_SS H_SD | qUCCSD: | 2 1 | - /// | H_DS H_DD | | 1 0 | + /// H̄ everywhere. For singles+doubles the matrix is + /// | H_SS H_SD | e.g. | 2 1 | + /// | H_DS H_DD | | 1 0 | /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare - /// Hamiltonian integrals (no commutators): f (Eqs. (30),(34)) plus - /// \f$ \langle ij\|kl\rangle,\ \langle ab\|cd\rangle \f$ (Eq. (48)) - /// (10.1063/5.0062090 Sec. II C, Eqs. (29), (41), (44), (48)). + /// Hamiltonian integrals (no commutators). Those are the ranks + /// 10.1063/5.0062090 Sec. II C truncates qUCCSD at, Eqs. (29), (41), (44) + /// and (48), which it writes UCCSD[2|2,1,0]; that section also says which + /// H̄ components each block then retains. /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING /// manifold rank, so one set of numbers serves EE, IP and EA (read S as - /// 1h/1p and D as 2h1p/1h2p: qUCCSD, IP-qUCCSD and EA-qUCCSD are all - /// `{2,1,1,0}`, 10.1021/acs.jctc.5c01991 Table 1). Empty (the default) - /// selects the uniform H̄ at `hbar_comm_rank` everywhere. + /// 1h/1p and D as 2h1p/1h2p; 10.1021/acs.jctc.5c01991 Table 1 maps the + /// IP/EA blocks onto qUCCSD's and its Fig. 1 carries their ranks). Empty + /// (the default) selects the uniform H̄ at `hbar_comm_rank` everywhere. /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and /// has nothing to truncate. /// @pre `block_ranks` is either empty or `K`×`K` @@ -222,7 +223,8 @@ class CC { /// ground-state amplitude residual) removed. See `eom_r_blocked` in cc.cpp /// for why. The removed terms vanish at converged amplitudes when a block /// rank equals `hbar_comm_rank`, so this changes those blocks' equations - /// but not the numbers they evaluate to. + /// but not the numbers they evaluate to. `BCH` has no N/R split to take, + /// so it keeps them. /// @return vector of right side sigma equations; element 0 is null iff /// `np == nh` // clang-format on diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 1a7690ef93..de60f7f036 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -240,7 +240,8 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE(size(ee[2]) == 21); // the same ranks drive IP: manifolds are indexed by ascending rank, so - // {1h, 2h1p} takes the place of {S, D} (10.1021/acs.jctc.5c01991 Table 1) + // {1h, 2h1p} takes the place of {S, D} (10.1021/acs.jctc.5c01991 Table 1 + // for the block mapping, its Fig. 1 for the ranks) const auto ip = cc.eom_r(nₚ(1), nₕ(2), quccsd); REQUIRE(ip.size() == 2); REQUIRE(size(ip[0]) == 32); From 4867eaac059314aa54a2609130ae88a37c3d725d Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 16:12:37 -0400 Subject: [PATCH 34/37] test: disable long qUCCSD unit tests --- tests/unit/test_mbpt_cc.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index de60f7f036..e4f481f941 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -223,6 +223,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { #endif // !defined(SEQUANT_SKIP_LONG_TESTS) } +#ifndef SEQUANT_SKIP_LONG_TESTS SECTION("bernoulli_quccsd_eom") { using namespace sequant; using namespace sequant::mbpt; @@ -265,6 +266,7 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE_THAT(cc.eom_r(nₚ(2), nₕ(2)).at(1), EquivalentTo(cc.eom_r(nₚ(2), nₕ(2), {2, 2, 2, 2}).at(1))); } +#endif // !defined(SEQUANT_SKIP_LONG_TESTS) SECTION("energy") { // CC::energy() must equal the p==0 element of CC::t() for both ansätze. From 4c699844faf54c6b6bab09a1caeac8babc74c7a8 Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 17 Aug 2026 17:01:14 -0400 Subject: [PATCH 35/37] style(mbpt): use SEQUANT_ASSERT's message argument in cc.cpp SEQUANT_ASSERT(EXPR, ...) takes the message itself, so the legacy `cond && "msg"` idiom folds the string into the stringified condition. Converts the twelve remaining sites in this file; conditions are unchanged. Master has been moving the same way, which is what both conflicts in the merge just before this were about. Roughly 90 `&&`-form sites remain elsewhere in SeQuant/ and tests/; converting those is a separate repo-wide sweep. --- SeQuant/domain/mbpt/models/cc.cpp | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 4d9c099d7b..4c679e27e2 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -107,7 +107,7 @@ ExprPtr CC::energy(std::optional comm_rank) const { std::vector CC::t(size_t pmax, size_t pmin) const { pmax = (pmax == std::numeric_limits::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. @@ -170,11 +170,11 @@ std::vector CC::t(size_t pmax, size_t pmin) const { } std::vector 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 @@ -270,10 +270,10 @@ ExprPtr CC::rdm(size_t rank, std::optional comm_rank) const { std::vector CC::tʼ(size_t rank, size_t order, std::optional 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()) @@ -332,15 +332,14 @@ std::vector CC::tʼ(size_t rank, size_t order, std::vector CC::λʼ(size_t rank, size_t order, std::optional 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 @@ -483,10 +482,10 @@ std::vector CC::eom_r_blocked( std::vector CC::eom_r(nₚ np, nₕ nh, const std::vector& block_ranks) const { - SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order"); + 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̄ @@ -536,9 +535,9 @@ std::vector CC::eom_r(nₚ np, nₕ nh, } std::vector 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( From 224b36a02725752ffee5e8b6300a0d545a865599 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 20 Aug 2026 23:29:37 -0400 Subject: [PATCH 36/37] fix(mbpt): tighten the Bernoulli ansatz and RDM preconditions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bernoulli required only unitary(), which admits oU; require U. CC::rdm ignored hbar_expansion_ and silently returned the BCH density, so reject it as tʼ does. --- SeQuant/domain/mbpt/models/cc.cpp | 7 +++++-- SeQuant/domain/mbpt/models/cc.hpp | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 4c679e27e2..9921bdf20e 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -55,8 +55,8 @@ CC::CC(size_t n, const Options& opts) "CC: skip_singles must be true for orbital-optimized " "ansatz"); if (hbar_expansion_ == HbarExpansion::Bernoulli) { - SEQUANT_ASSERT(unitary(), - "CC: Bernoulli expansion requires a unitary ansatz"); + SEQUANT_ASSERT(ansatz_ == Ansatz::U, + "CC: Bernoulli expansion requires the U ansatz"); SEQUANT_ASSERT(hbar_comm_rank_, "CC: Bernoulli expansion requires hbar_comm_rank"); } @@ -238,6 +238,9 @@ std::vector CC::λ() const { } ExprPtr CC::rdm(size_t rank, std::optional 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); diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 332dad45bf..8fb480f939 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -35,7 +35,7 @@ class CC { enum class HbarExpansion { /// standard Baker-Campbell-Hausdorff commutator expansion BCH, - /// Bernoulli expansion, 10.1063/1.5030344 (unitary ansatz only) + /// Bernoulli expansion, 10.1063/1.5030344 (U ansatz only) Bernoulli }; From 281e4a1a8001a1908be8b1092ae3a870eceabc8c Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 20 Aug 2026 23:29:37 -0400 Subject: [PATCH 37/37] docs(mbpt): correct the off-diagonal claim about the N part 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. --- SeQuant/domain/mbpt/models/cc.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 9921bdf20e..0f4f482951 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -440,8 +440,10 @@ std::vector CC::eom_r_blocked( // 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 the equations do change, - // though not what they evaluate to at converged amplitudes. + // 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 hbars; for (const auto k : ranks) { auto [it, fresh] = hbars.try_emplace(k);