From 46fc146db105dc199b963cc904331dc6437de74c Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 00:57:34 -0400 Subject: [PATCH 01/16] 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 29643264fc..24b6afabf9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -398,6 +398,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 5119afba4803dbbc192d712ccea6eb44444f0614 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:01:29 -0400 Subject: [PATCH 02/16] 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 4dc4e71c7f419c7544399978e05b970470ba9546 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:02:03 -0400 Subject: [PATCH 03/16] =?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 fb88ce80198dd6123ff2e48aa2833e462dd4cc17 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:03:02 -0400 Subject: [PATCH 04/16] =?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 e1ec613b23..e8319f9007 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_; } @@ -67,11 +79,22 @@ 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()); + return mbpt::lst(H(), T(N, skip_singles()), truncation, {.unitary = unitary()}); } 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). @@ -84,6 +107,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 1be5181818..dc2642a9b0 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -30,6 +30,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 { /// type of CC ansatz. see CC::Ansatz @@ -51,6 +58,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, @@ -73,6 +82,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; @@ -172,6 +184,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; /// @brief computes reference expectation value of an expression. Dispatches /// to `mbpt::op::ref_av()` From 93938df1c27e29f5b9ddef7c39c94e5856d0e783 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 26 Jul 2026 01:03:30 -0400 Subject: [PATCH 05/16] 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 4f9d4356c43978bff3e069ff36b727aeb42778dd Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 27 Jul 2026 01:41:58 -0400 Subject: [PATCH 06/16] 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 e78e554e1727f1fc923b0f451cb0898e8b9d4e88 Mon Sep 17 00:00:00 2001 From: Ajay Date: Mon, 27 Jul 2026 08:03:08 -0400 Subject: [PATCH 07/16] 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 0541e047f2e37c8acfb6f2be7bf0362d443d69ad Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:38 -0400 Subject: [PATCH 08/16] 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 0b8e20581412b5cc8076d14ad8a0179f2587cdab Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:46 -0400 Subject: [PATCH 09/16] =?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 dc2642a9b0..903d1823e3 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -58,7 +58,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 e51ac449616be81a8a394a98cc86294905ddfb89 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:56 -0400 Subject: [PATCH 10/16] 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 18ef06e7d39422bdc16b96101321704eca4f891a Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 13:43:56 -0400 Subject: [PATCH 11/16] 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 1f443001b8aa15917cf0b65b59914e7538a31e8c Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 28 Jul 2026 23:56:19 -0400 Subject: [PATCH 12/16] =?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 eef9044f5a05fe9370489185afe77f859afb1fed Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:34:47 -0400 Subject: [PATCH 13/16] 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 07251084f55f6f8c8969f4c598838b49da0dc171 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:00 -0400 Subject: [PATCH 14/16] =?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 | 110 +++++++++++++++++++++++++++++- SeQuant/domain/mbpt/models/cc.hpp | 34 ++++++++- 2 files changed, 140 insertions(+), 4 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index e8319f9007..b4140f0d23 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 @@ -11,11 +13,13 @@ #include #include +#include #include #include #include #include #include +#include namespace { // alias reserved labels for readability @@ -88,7 +92,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) { @@ -362,9 +366,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( @@ -375,6 +473,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 903d1823e3..ce16cd7f8b 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -168,11 +168,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 87e72ec4e5305fe431b804e76366e53d2a16b65a Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:10 -0400 Subject: [PATCH 15/16] 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 65112a011a414491a7655e76dd0ccac79b139e74 Mon Sep 17 00:00:00 2001 From: Ajay Date: Sun, 2 Aug 2026 18:35:11 -0400 Subject: [PATCH 16/16] 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