diff --git a/.gitignore b/.gitignore index 1aac3de827..64fe8813e5 100644 --- a/.gitignore +++ b/.gitignore @@ -71,4 +71,4 @@ _codeql_detected_source_root .clangd .vscode out -run \ No newline at end of file +run 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..1380ccf96f --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.cpp @@ -0,0 +1,401 @@ +#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). For UCC the plain BCH +// series does not terminate, because σ mixes excitation and de-excitation. +// This file implements the Bernoulli-number resummation of 10.1063/1.5030344 +// that fixes that. H splits as F (Fock) + V (fluctuation potential); every +// operator O splits into O_N (its pure excitation/de-excitation part) and +// O_R = O − O_N. +// +// Equation numbers below are all from 10.1063/1.5030344, Sec. III B: +// superoperator inversion Eqs. (36)-(39); Bernoulli numbers B₁=−1/2, B₂=1/12, +// B₃=0, B₄=−1/720, Eq. (40), i.e. Bₙ/n! in the textbook normalization; the N/R +// split and the UCC amplitude condition V̄_N = 0, above and at Eq. (43); the H̄ +// recursion, Eq. (44); the sum +// H̄ = Σ_k H̄^k, Eq. (45); H̄⁰..H̄⁴, Eqs. (46)-(50). +// +// The F-cancellation: at an HF reference F has no occupied-virtual block +// (Brillouin, Eq. (32)), so H̄² and higher contain no F (stated just below +// Eq. (50)). H̄⁰ and H̄¹ do, and they are built from the full one-body operator, +// so their f_ov terms survive symbolically and vanish only on substitution. + +namespace { + +/// 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, per the O_N/O_R split above +/// Eq. (43) of 10.1063/1.5030344. A term is N iff its single residual +/// NormalOperator is a pure excitation (all creators pure-unoccupied AND all +/// annihilators pure-occupied) or a pure de-excitation (the reverse), with +/// rank ≤ @p cutoff. A term with no residual NormalOperator is rank-preserving, +/// hence R. +/// +/// Rank > @p cutoff falls to R rather than being dropped. The R filter drops a +/// term only because the amplitude condition V̄_N = 0 (Eq. (43)) makes it +/// zero, and that condition holds for rank ≤ N only, since σ here is +/// truncated at rank N. Eq. (43) itself states O_N with no rank limit, +/// because there σ carries every rank. +bool is_N_term(const sequant::ExprPtr& term, std::size_t cutoff) { + using namespace sequant; + auto isr = get_default_context().index_space_registry(); + + const auto* nop = find_nop(term); + if (!nop) return false; + + 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(const ExprPtr& expr_in) { + auto expr = expr_in->clone(); + simplify(expr); + FWickTheorem wick{expr}; + // use_topology defaults to ON, so it must be turned off explicitly. It keeps + // 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. + 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 shared labels (both a block-resolved + // part and sigma carry a/i) would fuse two independent summations in A*B. + // Reindex B to fresh temporaries; canonicalization restores tidy labels. + container::map repl; + for (const auto& idx : get_used_indices(B)) + repl.emplace(idx, Index::make_tmp_index(idx.space())); + 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. Both spaces are required; the + // accessors throw if the registry leaves either unspecified. + const auto& hole_t = isr->hole_space(); + const auto& particle_t = isr->particle_space(); + auto physical = [&](const IndexSpace& b) { + return hole_t.includes(b.type()) || particle_t.includes(b.type()); + }; + + 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; + for (const auto& b : bases) + if (physical(b) && b.qns() == g.space().qns() && + g.space().type().includes(b.type())) + c.push_back(b); + SEQUANT_ASSERT(!c.empty(), + "bernoulli: general index spans no hole/particle base " + "space with matching quantum numbers"); + choices.push_back(std::move(c)); + } + // cartesian product of assignments => sum of transformed terms; + // accumulate via Sum::append (linear) rather than operator+, which + // 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 it would collide with + // a definite index of the same base space already in the term, e.g. + // one an amplitude brought in. Canonicalization restores tidy labels. + repl.emplace(gens[k], Index::make_tmp_index(choices[k][idx[k]])); + sum->append(transform_expr(term, repl)); + // increment mixed-radix counter over the assignments + 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()) { + 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)); +} + +/// 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); + 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)); + + // Applies one partition tag to an expression; 'A' is no filter. `reduced` + // says the input is already wick_reduce'd, which every commutator output is; + // the V base is not, so it takes the reducing form. + auto part = [&](char tag, const ExprPtr& e, bool reduced) -> ExprPtr { + SEQUANT_ASSERT(tag == 'A' || tag == 'N' || tag == 'R', + "bernoulli::hbar: partition tag must be one of A, N, R"); + if (tag == 'A') return e; + if (tag == 'N') + return reduced ? N_part_reduced(e, cutoff) : N_part(e, cutoff); + return reduced ? R_part_reduced(e, cutoff) : R_part(e, cutoff); + }; + + // Every term of H̄^k is a nested commutator [[..[V_{p0},σ]_{f0}..],σ]_{f_k} + // with a per-level N/R/A partition tag applied after each commutator ('A' = + // no filter). nest(p0, f) evaluates such a node, memoizing every prefix + // (key = p0 + tags applied so far): prefixes repeat both within a rank and + // across ranks. Reuse is safe because expression composition deep-copies its + // operands. + container::map memo; + auto nest = [&](char p0, const char* f) -> ExprPtr { + // 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()) + it = memo.emplace(key, part(p0, V, /*reduced=*/false)).first; + ExprPtr cur = it->second; + for (int i = 0; f[i] != '\0'; ++i) { + key += f[i]; + it = memo.find(key); + if (it == memo.end()) { + auto cx = wick_commutator(cur, sigma); + it = memo.emplace(key, part(f[i], cx, /*reduced=*/true)).first; + } + cur = it->second; + } + return cur; + }; + + 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-commutators enter H̄ ONLY + // here. This is the F-cancellation, stated just below Eq. (50) as "the + // terms in H̄ involving F now truncate to the first power of σ". H̄⁰ carries + // the bare F. + add(1, wick_commutator(F, sigma)); + add({1, 2}, nest('A', "A")); + add({1, 2}, nest('R', "A")); + } + 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 + // (the F-cancellation). Listed in the paper's order; the outermost tag is + // always A. + add({1, 16}, nest('R', "RRRA")); + add({1, 16}, nest('A', "RRRA")); + 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..70e950c795 --- /dev/null +++ b/SeQuant/domain/mbpt/bernoulli.hpp @@ -0,0 +1,78 @@ +#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. +/// +/// The result is a tensor-level expression: coefficient tensors times +/// normal-ordered operators, not `mbpt::op` operators. Nothing is screened out +/// of it, so the caller projects every term. +/// +/// @pre an HF reference: F is taken to have no occupied-virtual block, which +/// is what keeps F out of H̄² and higher. The f_ov terms of H̄⁰ and H̄¹ are +/// carried symbolically and vanish only on substitution. +/// +/// @param N cluster/excitation rank (also the N/R rank cutoff) +/// @param rank highest Bernoulli order H̄^k to include (0..4) +/// @param skip1 exclude singles from T +/// @throw Exception if @p rank > 4 +ExprPtr hbar(std::size_t N, std::size_t rank, bool skip1); + +namespace detail { + +/// Applies Wick's theorem to @p expr retaining PARTIAL contractions, +/// reducing a product of normal-ordered operators to a sum of normal-ordered +/// operators (each = coefficient tensor × at most one residual NormalOperator; +/// fully-contracted terms carry none). Unlike the expectation-value path it +/// keeps operators rather than collapsing to a scalar VEV. +/// @note @p expr is left untouched; the reduction runs on a clone. +ExprPtr wick_reduce(const ExprPtr& expr); + +/// Normal-ordered commutator [A, B] = wick_reduce(A·B − B·A). NOT the bare +/// algebraic commutator: the operator product is Wick-reduced, so contractions +/// between A and B generate the lower-rank terms the Bernoulli expansion relies +/// on. B's summed indices are reindexed to fresh temporaries first, making them +/// disjoint from A's. +ExprPtr wick_commutator(const ExprPtr& A, const ExprPtr& B); + +/// Rewrites every general (non-base) index of the residual NormalOperator as +/// the sum over the hole/particle base spaces it spans (occupied/virtual). The +/// registry's other base spaces are dropped, which changes no projected +/// quantity: the single-reference manifolds annihilate the dropped terms (see +/// the @warning on hbar). After expansion every residual index is definite so +/// the N/R classifier can act on it. Idempotent on block-resolved input. +/// @pre the registry specifies both a hole and a particle space +ExprPtr expand_to_blocks(const ExprPtr& expr); + +/// Block-resolved N part (O_N of 10.1063/1.5030344): the terms whose single +/// residual NormalOperator is a pure excitation or pure de-excitation of rank ≤ +/// @p cutoff. Applies expand_to_blocks first. +ExprPtr N_part(const ExprPtr& expr, std::size_t cutoff); + +/// R part (O_R of 10.1063/1.5030344: expr minus its N part). Unlike N_part +/// the result is NOT block-resolved; it stays in compact general-index form, +/// which is much cheaper for the nested commutators that consume R. +ExprPtr R_part(const ExprPtr& expr, std::size_t cutoff); + +} // namespace detail + +} // namespace sequant::mbpt::bernoulli + +#endif // SEQUANT_DOMAIN_MBPT_BERNOULLI_HPP diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 59ab75d49a..0f4f482951 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -1,8 +1,10 @@ +#include #include #include #include #include #include +#include #include #include #include @@ -16,6 +18,7 @@ #include #include #include +#include namespace { // alias reserved labels for readability @@ -42,7 +45,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 +54,12 @@ 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(ansatz_ == Ansatz::U, + "CC: Bernoulli expansion requires the U ansatz"); + SEQUANT_ASSERT(hbar_comm_rank_, + "CC: Bernoulli expansion requires hbar_comm_rank"); + } } CC::Ansatz CC::ansatz() const { return ansatz_; } @@ -60,6 +70,8 @@ bool CC::unitary() const { std::optional CC::hbar_comm_rank() const { return hbar_comm_rank_; } +CC::HbarExpansion CC::hbar_expansion() const { return hbar_expansion_; } + bool CC::skip_singles() const { return skip_singles_; } bool CC::screen() const { return screen_; } @@ -69,6 +81,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()) @@ -76,6 +91,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). @@ -86,7 +107,19 @@ ExprPtr CC::energy(std::optional comm_rank) const { std::vector CC::t(size_t pmax, size_t pmin) const { pmax = (pmax == std::numeric_limits::max() ? N : pmax); - SEQUANT_ASSERT(pmax >= pmin && "pmax should be >= pmin"); + SEQUANT_ASSERT(pmax >= pmin, "pmax should be >= pmin"); + + // Bernoulli: the hbar expansion is at tensor level, project and call the + // tensor level ref_av directly. + 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(); @@ -137,11 +170,11 @@ std::vector CC::t(size_t pmax, size_t pmin) const { } std::vector CC::λ() const { - SEQUANT_ASSERT(!unitary() && "there is no need for CC::λ for unitary ansatz"); + SEQUANT_ASSERT(!unitary(), "there is no need for CC::λ for unitary ansatz"); // construct hbar const auto commutator_rank = hbar_comm_rank_.value_or(4); - SEQUANT_ASSERT(commutator_rank >= 1 && "CC::λ: hbar_comm_rank must be >= 1"); + SEQUANT_ASSERT(commutator_rank >= 1, "CC::λ: hbar_comm_rank must be >= 1"); auto hbar = this->hbar(commutator_rank - 1); // -1 because of the connection with the projector @@ -205,6 +238,9 @@ std::vector CC::λ() const { } ExprPtr CC::rdm(size_t rank, std::optional comm_rank) const { + SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli, + "CC::rdm: the Bernoulli expansion is not supported yet"); + // 1. replacement operator {ã^{p_1..p_r}_{p_{r+1}..p_{2r}}} (see op::ã); its // indices are free, so they become the free indices of γ. auto replacer = op::ã(rank); @@ -237,15 +273,17 @@ ExprPtr CC::rdm(size_t rank, std::optional comm_rank) const { std::vector CC::tʼ(size_t rank, size_t order, std::optional nbatch) const { - SEQUANT_ASSERT(order == 1 && + SEQUANT_ASSERT(order == 1, "sequant::mbpt::CC::tʼ(): only first-order perturbation is " "supported now"); - SEQUANT_ASSERT(rank == 1 && + SEQUANT_ASSERT(rank == 1, "sequant::mbpt::CC::tʼ(): only one-body perturbation " "operator is supported now"); if (unitary()) SEQUANT_ASSERT(pertbar_comm_rank_, "pertbar_comm_rank must be specified for unitary ansatz"); + SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli, + "CC::tʼ: the Bernoulli expansion is not supported yet"); // construct h1_bar // truncate h1_bar at rank 2 for one-body perturbation operator and at rank 4 @@ -297,15 +335,14 @@ std::vector CC::tʼ(size_t rank, size_t order, std::vector CC::λʼ(size_t rank, size_t order, std::optional nbatch) const { - SEQUANT_ASSERT(order == 1 && + SEQUANT_ASSERT(order == 1, "sequant::mbpt::CC::λʼ(): only first-order perturbation is " "supported now"); - SEQUANT_ASSERT(rank == 1 && + SEQUANT_ASSERT(rank == 1, "sequant::mbpt::CC::λʼ(): only one-body perturbation " "operator is supported now"); - SEQUANT_ASSERT(!unitary() && - "there is no need for CC::λʼ for unitary ansatz"); - SEQUANT_ASSERT(ansatz_ == Ansatz::T && + SEQUANT_ASSERT(!unitary(), "there is no need for CC::λʼ for unitary ansatz"); + SEQUANT_ASSERT(ansatz_ == Ansatz::T, "CC::λʼ: only traditional ansatz is supported"); // construct hbar @@ -364,13 +401,104 @@ namespace { constexpr Normalization eom_norm = Normalization::SquareRoot; } // namespace -std::vector CC::eom_r(nₚ np, nₕ nh) const { - SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order"); +// Per-block-truncated EOM sigma equations. For the qUCCSD ranks see +// 10.1063/5.0062090 Sec. II C, Eqs. (29)-(48); for the IP/EA analogues, +// 10.1021/acs.jctc.5c01991 Fig. 1 (Table 1 there maps out which H̄ components +// enter each block, not the commutator ranks they are truncated at). +// +// Each block is the sandwich 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 CC::eom_r_blocked( + nₚ np, nₕ nh, const std::vector& block_ranks) const { + SEQUANT_ASSERT(unitary(), "eom_r_blocked requires a unitary ansatz"); + + std::vector> manifolds; + for (std::int64_t rp = np, rh = nh; rp >= 0 && rh >= 0; --rp, --rh) { + if (rp == 0 && rh == 0) break; + manifolds.emplace_back(rp, rh); + if (rp == 0 || rh == 0) break; + } + + std::ranges::reverse(manifolds); + const auto K = manifolds.size(); + // empty means uniform truncation at hbar_comm_rank in every block + const std::vector ranks = + block_ranks.empty() ? std::vector(K * K, hbar_comm_rank().value()) + : block_ranks; + SEQUANT_ASSERT(ranks.size() == K * K, + "CC::eom_r: block_ranks must be a K x K row-major matrix, " + "K = number of projection manifolds"); + + // Bernoulli H̄ is tensor-level, BCH H̄ operator-level; the bra/ket/vev trio + // below must match it. Empty connectivity, as everywhere on the unitary path. + const bool tensor_level = hbar_expansion_ == HbarExpansion::Bernoulli; + + // One H̄ per distinct truncation order, reduced to its R part (Bernoulli + // only: the operator-level BCH H̄ has no N/R split to take). The N part is + // the ground-state amplitude residual <Φl|H̄|Φ0>, which Eq. (6) zeroes only + // at the amplitude rank, so a block truncated below that rank would keep it. + // The diagonal is untouched either way: an N operator of rank r shifts the + // manifold rank by r, so it never lands on a diagonal block, and it has no + // reference expectation value. Off the diagonal removing it is a no-op only + // where the block rank equals hbar_comm_rank; below that rank the terms are + // off-shell, so the numbers change too. That is the point: Eqs. (41)-(47) of + // 10.1063/5.0062090 carry no such intermediate. + container::map hbars; + for (const auto k : ranks) { + auto [it, fresh] = hbars.try_emplace(k); + if (!fresh) continue; // deriving H̄ twice for one rank is not cheap + it->second = hbar(k); + if (tensor_level) it->second = bernoulli::detail::R_part(it->second, N); + } + auto bra_of = [tensor_level](std::int64_t p, std::int64_t h) { + return tensor_level ? op::tensor::δl(nₚ(p), nₕ(h)) : op::δl(nₚ(p), nₕ(h)); + }; + auto ket_of = [tensor_level](std::int64_t p, std::int64_t h) { + return tensor_level ? op::tensor::r(nₚ(p), nₕ(h), eom_norm) + : op::r(nₚ(p), nₕ(h), eom_norm); + }; + auto vev = [tensor_level, this](const ExprPtr& e) { + return tensor_level ? op::tensor::ref_av(e) + : op::ref_av(e, {.connect = {}, + .screen = screen_, + .use_topology = use_topology_}); + }; + + using std::min; + std::vector 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(ranks.at(i * K + j)); + const auto ket = ket_of(kp, kh); + acc->append(vev(bra * hbar_ij * ket)); + // -<0|H̄^(k_ii)|0>, written as so Wick keeps E's summed + // 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; +} + +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 && + get_default_context().spbasis() != SPBasis::Spinfree, "spin-free basis does not yet support non particle-conserving cases"); + // Bernoulli always takes the blocked path: the uniform one below commutes H̄ + // with an operator-level R, which a tensor-level H̄ cannot take part in. An + // empty matrix there means uniform truncation at hbar_comm_rank. + if (!block_ranks.empty() || hbar_expansion_ == HbarExpansion::Bernoulli) + return eom_r_blocked(np, nh, block_ranks); + // construct hbar const auto hbar = this->hbar(); @@ -412,9 +540,9 @@ std::vector CC::eom_r(nₚ np, nₕ nh) const { } std::vector CC::eom_l(nₚ np, nₕ nh) const { - SEQUANT_ASSERT(!unitary() && + SEQUANT_ASSERT(!unitary(), "there is no need for CC::eom_l for unitary ansatz"); - SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order"); + SEQUANT_ASSERT(np > 0 || nh > 0, "Unsupported excitation order"); if (np != nh) SEQUANT_ASSERT( diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 151d129910..8fb480f939 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 (U 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, @@ -76,6 +89,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; @@ -109,6 +125,9 @@ class CC { /// the explicit form. For a unitary ansatz the reverse holds: H̄ is already /// self-contained, so connectivity must be left empty. See the "Using H̄ /// outside the CC class" section of the user guide. + /// @note Under `HbarExpansion::Bernoulli` the result is tensor-level, so it + /// takes `op::tensor` projectors and `op::tensor::ref_av`, not their `op` + /// counterparts, and it is unscreened: `screen` has no effect there. [[nodiscard]] ExprPtr hbar( std::optional truncation_rank = std::nullopt) const; @@ -171,11 +190,46 @@ 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 is + /// | H_SS H_SD | e.g. | 2 1 | + /// | H_DS H_DD | | 1 0 | + /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator + /// [[V,σ],σ], H_SD and H_DS through the single [V,σ], H_DD the bare + /// Hamiltonian integrals (no commutators). Those are the ranks + /// 10.1063/5.0062090 Sec. II C truncates qUCCSD at, Eqs. (29), (41), (44) + /// and (48), which it writes UCCSD[2|2,1,0]; that section also says which + /// H̄ components each block then retains. + /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING + /// manifold rank, so one set of numbers serves EE, IP and EA (read S as + /// 1h/1p and D as 2h1p/1h2p; 10.1021/acs.jctc.5c01991 Table 1 maps the + /// IP/EA blocks onto qUCCSD's and its Fig. 1 carries their ranks). Empty + /// (the default) selects the uniform H̄ at `hbar_comm_rank` everywhere. + /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ is exact and + /// has nothing to truncate. + /// @pre `block_ranks` is either empty or `K`×`K` + /// @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. `BCH` has no N/R split to take, + /// so it keeps them. + /// @return vector of right side sigma equations; element 0 is null iff + /// `np == nh` + // clang-format on + [[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 @@ -219,6 +273,14 @@ 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 `eom_r`'s per-block-truncated path, taken whenever `block_ranks` is + /// non-empty or the expansion is Bernoulli + /// @param block_ranks see `eom_r`; empty means uniform `hbar_comm_rank` + /// @pre a unitary ansatz + [[nodiscard]] std::vector eom_r_blocked( + nₚ np, nₕ nh, const std::vector& block_ranks) const; /// @return the `LSTOptions` this engine uses for every `mbpt::lst()` call /// @note The choice of commutator representation is really a question of 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 3206017d29..50e52a5c12 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,212 @@ 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); + 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); // general g + const auto Vn = bernoulli::detail::N_part(V, 2); + const auto Vr = bernoulli::detail::R_part(V, 2); + // N ⊎ R reconstructs V. R stays in compact general-index form, so check the + // 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. + // The F-cancellation: F appears only in H̄¹, so rank r − rank r−1 is + // F-free for r ≥ 2. + auto h0 = bernoulli::hbar(2, 0, false); + auto h1 = bernoulli::hbar(2, 1, false); + auto h2 = bernoulli::hbar(2, 2, false); + 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); + + // one projected equation pinned in full: term counts are blind to the + // coefficients, which is where a mis-weighted Wick reduction shows up. + // Doubles at H̄¹, the smallest such equation. + const auto R2_h1 = CC(2, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 1, + .hbar_expansion = CC::HbarExpansion::Bernoulli}) + .t() + .at(2); + REQUIRE_THAT( + R2_h1, EquivalentTo( + L"1/2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;a_3,a_4}:A-C-S " + L"* t{a_3,a_4;i_1,i_2}:A-N-S " + L"+ Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;i_1,i_2}:A-C-S " + L"+ 1/2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,i_4;i_1,i_2}:A-C-S " + L"* t{a_1,a_2;i_3,i_4}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * f{i_3;i_1}:A-C-S " + L"* t{a_1,a_2;i_2,i_3}:A-N-S " + L"- 2 Â{i_1,i_2;a_1,a_2}:A-C-S * f{a_1;a_3}:A-C-S " + L"* t{a_2,a_3;i_1,i_2}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{a_1,a_2;i_1,a_3}:A-C-S " + L"* t{a_3;i_2}:A-N-S " + L"+ 2 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,a_1;i_1,i_2}:A-C-S " + L"* t{a_2;i_3}:A-N-S " + L"- 4 Â{i_1,i_2;a_1,a_2}:A-C-S * g{i_3,a_1;i_1,a_3}:A-C-S " + L"* t{a_2,a_3;i_2,i_3}:A-N-S")); + +#ifndef SEQUANT_SKIP_LONG_TESTS + const auto E = cc.energy(3); + REQUIRE_THAT(E, !EquivalentTo(amps.at(0))); + REQUIRE(size(E) == 46); +#endif // !defined(SEQUANT_SKIP_LONG_TESTS) + } + +#ifndef 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} (10.1021/acs.jctc.5c01991 Table 1 + // for the block mapping, its Fig. 1 for the ranks) + const auto ip = cc.eom_r(nₚ(1), nₕ(2), quccsd); + REQUIRE(ip.size() == 2); + REQUIRE(size(ip[0]) == 32); + REQUIRE(size(ip[1]) == 11); + + // one manifold means one block, whose sandwich minus shift is the + // commutator the uniform path builds; that path shares no code with + // eom_r_blocked, so this pins the block construction against it + const CC bch(2, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}); + REQUIRE_THAT(bch.eom_r(nₚ(1), nₕ(1), {2}).at(1), + EquivalentTo(bch.eom_r(nₚ(1), nₕ(1)).at(1))); + + if (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { + // block_ranks must be a K x K matrix over the manifolds ... + REQUIRE_THROWS_AS(cc.eom_r(nₚ(2), nₕ(2), {2, 1, 0}), Exception); + // ... and the ansatz must be unitary + REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); + } + + // no matrix means uniform truncation at hbar_comm_rank + REQUIRE_THAT(cc.eom_r(nₚ(2), nₕ(2)).at(1), + EquivalentTo(cc.eom_r(nₚ(2), nₕ(2), {2, 2, 2, 2}).at(1))); + } +#endif // !defined(SEQUANT_SKIP_LONG_TESTS) + SECTION("energy") { // CC::energy() must equal the p==0 element of CC::t() for both ansätze. const auto N = 2;