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 diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..4a6a6213a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -399,6 +399,8 @@ set(SeQuant_symb_src set(SeQuant_mbpt_src SeQuant/domain/mbpt/antisymmetrizer.cpp SeQuant/domain/mbpt/antisymmetrizer.hpp + SeQuant/domain/mbpt/bernoulli.cpp + SeQuant/domain/mbpt/bernoulli.hpp SeQuant/domain/mbpt/biorthogonalization.cpp SeQuant/domain/mbpt/biorthogonalization.hpp SeQuant/domain/mbpt/context.cpp diff --git a/SeQuant/domain/mbpt/bernoulli.cpp b/SeQuant/domain/mbpt/bernoulli.cpp new file mode 100644 index 0000000000..37f1c103d2 --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -0,0 +1,404 @@ +#include + +#include +#include +#include +#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 with Bernoulli numbers as the expansion coefficients, +// which leaves the truncation at a chosen commutator rank as the only +// approximation. H is split as F (Fock, rank-preserving) + V (fluctuation +// potential), and every operator O is split into O_N (all excitation and +// de-excitation operators) and O_R = O − O_N. At a converged RHF/UHF reference +// two cancellations hold: F survives only in H̄¹, and the higher orders carry +// only R-subscripted inner commutators. +// +// All equation references are to 10.1063/1.5030344 (Sec. III B): superoperator +// inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, B₃=0, B₄=−1/720 +// Eq. (40); the N/R split and the UCC amplitude condition V̄_N = 0 above and at +// Eq. (43); the iterative recursion for V̄ Eq. (44); the assembly +// H̄ = Σ_k H̄^k Eq. (45); the rank-by-rank operators H̄⁰..H̄⁴ Eqs. (46)-(50). +// Cancellation #1 (F enters only H̄¹) is stated just below Eq. (50). + +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 at +/// most 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 at most 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 +/// (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. 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()); + }); + }; + 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 { + +namespace detail { + +ExprPtr wick_reduce(ExprPtr expr) { + simplify(expr); + FWickTheorem wick{expr}; + // 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); + simplify(result); + return result; +} + +ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B) { + // A and B are built independently, so their summed indices are local to each. + // If both use the same labels (a block-resolved R/N part and sigma both carry + // a/i), the product A*B fuses two independent summations. That corrupts the + // contraction. Reindex B to fresh temporaries. Canonicalization restores tidy + // labels. + 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 { + +/// 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) { + 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; }); + }; + + // 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) { + 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->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->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). + 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}; + }; + + // transform_sum_expr maps in parallel, canonicalizes each result, and + // accumulates into a HashingAccumulator + ExprPtr out; + if (expr.is()) { + out = transform_sum_expr(expr.as().summands(), expand_term); + } else { + out = expand_term(expr); + } + 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 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); +} + +} // namespace detail + +/// Assembles H̄ order by order (see header), summing H̄⁰..H̄^rank of Eq. (45). +/// 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) + throw Exception("bernoulli::hbar: only ranks [0,4] are implemented"); + + using namespace detail; + 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† + + // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} + // with a per-level N/R/A partition tag applied after each commutator ('A' = + // no filter). nest(p0, f) evaluates such a node, memoizing every prefix + // (key = p0 + tags applied so far): the terms share prefixes both within a + // rank (the 9 rank-4 terms have only 3 distinct level-1 and 6 level-2 nodes) + // and across ranks (all four ranks share the same three level-1 nodes), so + // the memo avoids recomputing them. Reusing a memoized ExprPtr is safe: + // expression composition deep-copies operands (Product/Sum append clone), so + // wick_commutator does not mutate its arguments. Commutator outputs are + // already wick_reduce'd, so the reduced-input N/R filters apply. + container::map memo; + auto nest = [&](char p0, const char* f) -> ExprPtr { + SEQUANT_ASSERT((p0 == 'A' || p0 == 'N' || p0 == 'R') && + "bernoulli::hbar: partition tag must be one of A, N, R"); + // grow `key` in place rather than deriving it from the memo iterator: + // container::map is a flat_map, whose insertions invalidate iterators + std::string key{p0}; + auto it = memo.find(key); + if (it == memo.end()) { + ExprPtr base = (p0 == 'N') ? N_part(V, cutoff) + : (p0 == 'R') ? R_part(V, cutoff) + : V; + it = memo.emplace(key, std::move(base)).first; + } + ExprPtr op = it->second; + for (int i = 0; f[i] != '\0'; ++i) { + SEQUANT_ASSERT((f[i] == 'A' || f[i] == 'N' || f[i] == 'R') && + "bernoulli::hbar: partition tag must be one of A, N, R"); + key += f[i]; + it = memo.find(key); + if (it == memo.end()) { + auto cx = wick_commutator(op, sigma); + ExprPtr filtered = (f[i] == 'R') ? R_part_reduced(cx, cutoff) + : (f[i] == 'N') ? N_part_reduced(cx, cutoff) + : cx; + it = memo.emplace(key, std::move(filtered)).first; + } + op = it->second; + } + return op; + }; + + 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)). + 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)] + 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)] + 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. + 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")); + } + auto result = acc.make_expr(); + return simplify(result); +} + +} // namespace sequant::mbpt::bernoulli diff --git a/SeQuant/domain/mbpt/bernoulli.hpp b/SeQuant/domain/mbpt/bernoulli.hpp new file mode 100644 index 0000000000..4f924a2900 --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -0,0 +1,69 @@ +#ifndef SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP +#define SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP + +#include +#include + +namespace sequant::mbpt::bernoulli { + +/// Tensor-level H̄ = Σ_{k=0..rank} H̄^k in the Bernoulli expansion, for +/// σ = T−T† of rank N. +/// +/// The Bernoulli expansion rewrites the non-terminating UCC +/// similarity-transform series so that Bernoulli numbers appear as the +/// expansion coefficients; the rank-by-rank operators H̄⁰..H̄⁴ are Eqs. (46)-(50) +/// of 10.1063/1.5030344. +/// +/// @warning Single-reference only. The N/R split expands general indices over +/// the hole and particle spaces alone (see detail::expand_to_blocks), dropping +/// any other base space the registry defines. That is harmless only because the +/// single-reference projection manifolds annihilate the dropped terms. Under a +/// multireference registry they contribute, and both the N and the R part come +/// out wrong. Nothing checks for this. +/// +/// @param N cluster/excitation rank (also the N/R rank cutoff) +/// @param rank highest Bernoulli order H̄^k to include (0..4) +/// @param skip1 exclude singles from T +/// @throw Exception if @p rank > 4 +ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1); + +namespace detail { + +/// Applies Wick's theorem to @p expr retaining PARTIAL contractions, +/// reducing a product of normal-ordered operators to a sum of normal-ordered +/// operators (each = coefficient tensor × at most one residual NormalOperator; +/// fully-contracted terms carry none). Unlike the expectation-value path it +/// keeps operators rather than collapsing to a scalar VEV. +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); + +/// 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): 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. 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 + +} // namespace sequant::mbpt::bernoulli + +#endif // SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index b30ca6486e..36d5173344 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -1,8 +1,11 @@ +#include #include #include #include #include +#include #include +#include #include #include #include @@ -16,6 +19,7 @@ #include #include #include +#include namespace { // alias reserved labels for readability @@ -42,7 +46,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"); @@ -50,6 +55,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_; } @@ -66,6 +79,8 @@ CC CC::with_hbar_comm_rank(size_t rank) const { return result; } +CC::HbarExpansion CC::hbar_expansion() const { return hbar_expansion_; } + bool CC::skip_singles() const { return skip_singles_; } bool CC::screen() const { return screen_; } @@ -75,6 +90,9 @@ bool CC::use_topology() const { return use_topology_; } ExprPtr CC::hbar(std::optional truncation_rank) const { const auto truncation = truncation_rank.value_or(hbar_comm_rank_.value_or(4)); + if (hbar_expansion_ == HbarExpansion::Bernoulli) + return bernoulli::hbar(N, truncation, skip_singles()); + // for a non-unitary ansatz this is the cheaper connected-product form, which // is only equivalent to the commutator once the caller supplies operator // connectivity to ref_av (see lst_options() and the @warning on hbar()) @@ -82,6 +100,12 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { } ExprPtr CC::energy(std::optional comm_rank) const { + // Bernoulli: the hbar expansion is at tensor level, call the tensor level + // ref_av directly. No connectivity or screening. + if (hbar_expansion_ == HbarExpansion::Bernoulli) { + const auto erank = comm_rank.value_or(hbar_comm_rank_.value()); + return op::tensor::ref_av(this->hbar(erank)); + } // <0|H̄|0>: reference expectation value of H̄ at the requested commutator // truncation. No projector ⇒ this is the energy. ref_av applies the // connectivity (empty for unitary, default otherwise). @@ -94,6 +118,18 @@ 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: the hbar expansion is at tensor level, project and call the + // tensor level ref_av directly. + if (hbar_expansion_ == HbarExpansion::Bernoulli) { + const auto hbar = this->hbar(); + std::vector 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(); @@ -368,15 +404,102 @@ 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. For the qUCCSD ranks see +// 10.1063/5.0062090 Sec. II C, Eqs. (29)-(48); for the IP/EA analogues, +// 10.1021/acs.jctc.5c01991 Table 1. +// +// Each block is the sandwich of Eq. (7). Eq. (10) writes H̄ as +// E_gr + a normal-ordered remainder and builds the blocks from the remainder +// alone, so here the diagonal carries an explicit -<0|H̄|0> instead. +std::vector eom_r_blocked(const CC& cc, nₚ np, nₕ nh, + const std::vector& block_ranks, + size_t N) { + if (!cc.unitary()) throw Exception("eom_r_blocked 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::ranges::reverse(manifolds); + const auto K = manifolds.size(); + if (block_ranks.size() != K * K) + throw Exception( + "CC::eom_r: block_ranks must be a K x K row-major matrix, " + "K = number of projection manifolds"); + + // 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, reduced to its R part. The N part is + // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes at + // the amplitude rank only, so a block truncated below it would keep the + // residual. Dropping it everywhere is exact: an N operator of rank r shifts + // the manifold rank by r, so it never reaches a diagonal block, and it has + // no reference expectation value, so the shift below is unchanged. + container::map hbars; + for (const auto k : block_ranks) { + auto [it, fresh] = hbars.try_emplace(k); + if (!fresh) continue; // deriving H̄ twice for one rank is not cheap + it->second = cc.hbar(k); + 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( get_default_context().spbasis() != SPBasis::Spinfree && "spin-free basis does not yet support non particle-conserving cases"); + // if block ranks are specified, dispatch and early return + if (!block_ranks.empty()) return eom_r_blocked(*this, np, nh, block_ranks, N); + + // the uniform path below commutes H̄ with an operator-level R, which the + // tensor-level Bernoulli H̄ cannot take part in + 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 f572b338a3..09bfa18e19 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -32,6 +32,13 @@ class CC { oU }; + enum class HbarExpansion { + /// standard Baker-Campbell-Hausdorff commutator expansion + BCH, + /// Bernoulli expansion, 10.1063/1.5030344 (unitary ansatz only) + Bernoulli + }; + /// Configuration options for CC class struct Options { SEQUANT_DESIGNATED_INIT_ONLY; @@ -54,6 +61,12 @@ 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. + /// @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; }; /// @brief constructs CC engine with default options (traditional ansatz, @@ -81,6 +94,9 @@ class CC { /// over [[nodiscard]] CC with_hbar_comm_rank(size_t 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; @@ -176,11 +192,45 @@ 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 + /// (10.1063/5.0062090 Sec. II C, Eqs. (29), (41), (44), (48)). + /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING + /// manifold rank, so one set of numbers serves EE, IP and EA (read S as + /// 1h/1p and D as 2h1p/1h2p: qUCCSD, IP-qUCCSD and EA-qUCCSD are all + /// `{2,1,1,0}`, 10.1021/acs.jctc.5c01991 Table 1). Empty (the default) + /// selects the uniform H̄ at `hbar_comm_rank`, which the Bernoulli + /// expansion does not support. + /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and + /// has nothing to truncate. + /// @throw Exception if `block_ranks` is neither empty nor `K`×`K`, if it is + /// non-empty under a non-unitary ansatz, or if it is empty under the + /// Bernoulli expansion + /// @note each block is the sandwich \f$ \langle i|\bar{H}|j \rangle \f$ + /// (Eq. (7) of 10.1063/5.0062090) plus an explicit \f$ -E \f$ shift on the + /// diagonal, taken at the block's own truncation rank. Eq. (10) there + /// instead splits \f$ \bar{H} = E_{gr} + {} \f$ a normal-ordered remainder + /// and forms the blocks from the remainder. The returned object is + /// \f$ (\bar{H}-E)\hat{R} \f$. + /// @note under the Bernoulli expansion each block's H̄ has its N part (the + /// ground-state amplitude residual) removed. See `eom_r_blocked` in cc.cpp + /// for why. The removed terms vanish at converged amplitudes when a block + /// 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 @@ -224,6 +274,7 @@ class CC { bool use_topology_ = true; std::optional hbar_comm_rank_ = std::nullopt; std::optional pertbar_comm_rank_ = std::nullopt; + HbarExpansion hbar_expansion_ = HbarExpansion::BCH; /// @return the `LSTOptions` this engine uses for every `mbpt::lst()` call /// @note The choice of commutator representation is really a question of 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; +} diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 365213df81..7310e9ec00 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -6,7 +6,9 @@ #include #include #include +#include #include +#include #include #include "catch2_sequant.hpp" @@ -15,6 +17,16 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { using namespace sequant; using namespace sequant::mbpt; + auto has_tensor = [](const ExprPtr& e, const 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 @@ -51,6 +63,172 @@ 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; + // 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 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_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. + 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); + + // amplitudes through H̄² (hbar_comm_rank) + const auto amps = cc.t(); + REQUIRE(amps.size() == 3); + REQUIRE_THAT(amps[1], !EquivalentTo(ex(0))); + REQUIRE_THAT(amps[2], !EquivalentTo(ex(0))); + + 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("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 Sec. II C: SS at the double + // commutator (Eq. 29), SD/DS at the single (Eqs. 41, 44), DD bare (Eq. 48). + const std::vector quccsd = {2, 1, 1, 0}; + + const auto ee = cc.eom_r(nₚ(2), nₕ(2), quccsd); + 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;