From a05d4febe8ca8a0d60954410e0b8317c7fd6a410 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 13 May 2026 20:21:47 -0400 Subject: [PATCH 1/6] mbpt: add Kramers QN and spinor module (MVP) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Introduces the Kramers symmetry quantum number for relativistic 2-component theories, mirroring the existing Spin QN structure: * space_qns.hpp: new `Kramers` enum (up/down/any/null) using bits 7-8 of the QN bitmask. Spin and Kramers are mutually exclusive on a single index — an index is either spin-orbital (alpha/beta) or relativistic-spinor (kramers-up/kramers-down), reflecting the physical regime, not both. * spinor.{hpp,cpp}: new module hosting both the Kramers tracer (kramers_trace) and the upcoming quaternion decomposer (Phase 4). Kramers index annotations use ⇑/⇓ glyphs (U+21D1/U+21D3), distinct from spin's ↑/↓ so the QN families compose without label collision. * `make_kramers_up/dn/free(Index)` helpers parallel `make_spinalpha/beta/free`, going through the IndexSpace registry so spaces with the resulting (label, qns) tuple are reused. * `append_kramers(expr, replacements)` applies an Index→Index replacement map to all tensors in an expression (mirrors `append_spin` in `spin.cpp`). * `kramers_trace(expr)` MVP: collects every Kramers::any index in the input, enumerates 2^N configurations (each index pinned to up or down), and emits the resulting Sum of substituted copies. No per-tensor TRS canonicalization or cross-tensor folding yet — those iterations are deferred; the evaluator (LCAOFactory canonical-storage layer) already collapses Kramers blocks to orbit reps at integral-fetch time, so the MVP delivers correct energies for closed-shell relativistic MP2/CC even without symbolic folding. The Visscher Eq. 28 two-term form falls out of subsequent simplification commits. --- CMakeLists.txt | 2 + SeQuant/domain/mbpt/space_qns.hpp | 23 ++++ SeQuant/domain/mbpt/spinor.cpp | 179 ++++++++++++++++++++++++++++++ SeQuant/domain/mbpt/spinor.hpp | 143 ++++++++++++++++++++++++ 4 files changed, 347 insertions(+) create mode 100644 SeQuant/domain/mbpt/spinor.cpp create mode 100644 SeQuant/domain/mbpt/spinor.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index d5da7c30ab..d35941ebe7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -423,6 +423,8 @@ set(SeQuant_mbpt_src SeQuant/domain/mbpt/space_qns.hpp SeQuant/domain/mbpt/spin.cpp SeQuant/domain/mbpt/spin.hpp + SeQuant/domain/mbpt/spinor.cpp + SeQuant/domain/mbpt/spinor.hpp SeQuant/domain/mbpt/vac_av.hpp SeQuant/domain/mbpt/vac_av.cpp SeQuant/domain/mbpt/utils.hpp diff --git a/SeQuant/domain/mbpt/space_qns.hpp b/SeQuant/domain/mbpt/space_qns.hpp index a48a638b4e..06247937f3 100644 --- a/SeQuant/domain/mbpt/space_qns.hpp +++ b/SeQuant/domain/mbpt/space_qns.hpp @@ -89,6 +89,29 @@ struct mask { static constexpr type value = static_cast(BatchingQNS::batch); }; +/// quantum numbers tags related to Kramers symmetry (relativistic 2-component) +/// \note Kramers QN takes the 8th and 9th rightmost bits since there are 3 +/// possible states (any, kramers-up, kramers-down). +/// \note Spin and Kramers are mutually exclusive on a single index: an index +/// is either spin-orbital (alpha/beta) or relativistic-spinor +/// (kramers-up/kramers-down), never both. Spin restriction is the +/// non-relativistic limit of Kramers restriction. +enum class Kramers : bitset_t { + up = 0b010000000, //!< unbarred Kramers spinor (⇑) + down = 0b100000000, //!< barred Kramers spinor (⇓) + /// arbitrary Kramers state represented by 2 bits so overlap and union + /// work as expected (`any & up = up`, `up | down = any`, etc.) + any = up | down, + // syntactic sugar + null = 0b000000000 +}; + +template <> +struct mask { + using type = std::underlying_type_t; + static constexpr type value = static_cast(Kramers::any); +}; + } // namespace sequant::mbpt #endif // SEQUANT_DOMAIN_MBPT_SPACE_QNS_HPP diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp new file mode 100644 index 0000000000..3ee069218f --- /dev/null +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -0,0 +1,179 @@ +// +// Spinor decomposition passes for relativistic 2-component theories. +// See spinor.hpp for design overview. +// + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace sequant::mbpt { + +namespace detail { + +/// Construct a new Index in the same IndexSpace type as @p idx, with +/// the @p k Kramers QN bit set (replacing any existing Kramers state) +/// and the label decorated with the matching ⇑/⇓ glyph. Spin bits +/// (alpha/beta) are NOT touched here — Spin and Kramers are mutually +/// exclusive on a single index, but the per-index check is left to +/// the caller (a typical workflow only ever sets one of them). +/// +/// Mirrors `detail::make_index_with_spincase` in `spin.cpp` so that +/// SeQuant's IndexSpace registry is reused: if a space with the +/// resulting label and qns already exists, we adopt it. +Index make_index_with_kramers(const Index& idx, mbpt::Kramers k) { + // sanity check: at most one Kramers annotation already + SEQUANT_ASSERT(!(idx.label().find(L'⇑') != std::wstring::npos && + idx.label().find(L'⇓') != std::wstring::npos)); + + // Strip any existing Kramers bits, then OR in the requested state. + auto qns = mbpt::kramers_annotation_remove(idx.space().qns()).unIon(k); + + IndexSpace space; + const auto label = + mbpt::kramers_annotation_replace(idx.space().base_key(), k); + if (auto isr = get_default_context().index_space_registry()) { + auto* space_ptr = isr->retrieve_ptr(label); + if (space_ptr && space_ptr->type() == idx.space().type() && + space_ptr->qns() == qns) { + space = *space_ptr; + } + } + if (!space) { + space = IndexSpace{label, idx.space().type(), qns, + // assume size does not depend on Kramers state + idx.space().approximate_size()}; + } + auto protoindices = idx.proto_indices(); + for (auto& pidx : protoindices) pidx = make_index_with_kramers(pidx, k); + return Index{space, idx.ordinal(), protoindices}; +} + +} // namespace detail + +Index make_kramers_up(const Index& idx) { + return detail::make_index_with_kramers(idx, mbpt::Kramers::up); +} + +Index make_kramers_dn(const Index& idx) { + return detail::make_index_with_kramers(idx, mbpt::Kramers::down); +} + +Index make_kramers_free(const Index& idx) { + return detail::make_index_with_kramers(idx, mbpt::Kramers::any); +} + +ExprPtr append_kramers( + const ExprPtr& expr, + const container::map& index_replacements) { + auto add_to_tensor = [&](const Tensor& tensor) { + auto t = std::make_shared(tensor); + t->transform_indices(index_replacements); + return t; + }; + auto add_to_product = [&](const Product& product) { + auto p = std::make_shared(); + p->scale(product.scalar()); + for (auto&& term : product) { + if (term->is()) { + p->append(1, add_to_tensor(term->as())); + } else if (term->is() || term->is()) { + p->append(1, term); + } else { + throw std::runtime_error( + "append_kramers: unsupported Expr type in product: " + + term->type_name()); + } + } + return p; + }; + + if (expr->is()) return add_to_tensor(expr->as()); + if (expr->is()) return add_to_product(expr->as()); + if (expr->is()) { + auto s = std::make_shared(); + for (auto&& summand : *expr) { + s->append(append_kramers(summand, index_replacements)); + } + return s; + } + if (expr->is() || expr->is()) return expr; + throw std::runtime_error("append_kramers: unsupported Expr type"); +} + +namespace { + +/// Walk @p expr and collect the unique Kramers-eligible (i.e., +/// `Kramers::any`) indices in encounter order. Indices that already +/// have a specific Kramers state (up/down) are skipped. +container::svector collect_kramers_indices(const ExprPtr& expr) { + container::svector seen; + container::set seen_set; + auto consider = [&](const Index& idx) { + auto qns = idx.space().qns().to_int32(); + if ((qns & mask_v) != static_cast(Kramers::any)) { + // Either no Kramers state set or already specialized — skip. + return; + } + if (seen_set.insert(idx).second) seen.push_back(idx); + }; + expr->visit( + [&](const ExprPtr& e) { + if (e->is()) { + for (const auto& idx : e->as().const_braket()) { + consider(idx); + } + } + }, + /* recursive = */ true); + return seen; +} + +} // namespace + +ExprPtr kramers_trace(const ExprPtr& expr) { + // MVP (Phase 2b/i): enumerate the 2^N Kramers configurations of the + // expression's spinor (Kramers::any) indices and emit a Sum of substituted + // copies. Each summand is a Product (or Sum/Tensor) where every Kramers + // index has been pinned to up/down. + // + // Per-tensor TRS canonicalization (the (-1)^k bit-table fold) and + // cross-tensor simplification land in subsequent commits. + // + // The eventual evaluator (LCAOFactory) sees each summand as a concrete + // (g_block, t_block, ...) contraction; LCAOFactory's own canonical-storage + // layer (Phase 1.5) ensures we still pay only the orbit-rep integral cost. + + const auto indices = collect_kramers_indices(expr); + const std::size_t n = indices.size(); + if (n == 0) return expr; + if (n > 30) { + throw std::runtime_error( + "kramers_trace: too many Kramers indices (would emit > 2^30 terms)"); + } + + auto result = std::make_shared(); + const std::uint64_t n_configs = 1ull << n; + for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { + container::map repl; + for (std::size_t i = 0; i < n; ++i) { + const bool is_down = (cfg >> i) & 1u; + Index new_idx = is_down ? make_kramers_dn(indices[i]) + : make_kramers_up(indices[i]); + repl[indices[i]] = new_idx; + } + result->append(append_kramers(expr, repl)); + } + return result; +} + +} // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/spinor.hpp b/SeQuant/domain/mbpt/spinor.hpp new file mode 100644 index 0000000000..aec1b5fa31 --- /dev/null +++ b/SeQuant/domain/mbpt/spinor.hpp @@ -0,0 +1,143 @@ +// +// Spinor decomposition passes for relativistic 2-component theories. +// +// This module provides two symbolic transforms over tensor expressions: +// +// 1. kramers_trace(expr): rewrite a closed-shell expression with +// all-spinor (Kramers::any) indices into a sum over canonical +// Kramers-block representatives (one per orbit under +// {column-swap, Hermitian, time-reversal-with-(-1)^k}). Folds +// the 2^n raw Kramers blocks of a rank-n tensor into ~5 unique +// classes for rank-4, with conjugation/sign tracking baked into +// the per-tensor lookup. +// +// 2. quaternion_decompose(expr): (Phase 4 — not yet implemented) +// further decompose Kramers-up indices into the four real +// (s, x, y, z) quaternion components (Helmich-Paris, Repisky, +// Visscher 2016 Eq. 21), aggressively folding via Pauli/quaternion +// algebra to yield ≤ ~10 surviving terms for an MP2-style input. +// +// Both passes operate on the canonical SeQuant Tensor representation +// and consume its `Kramers` index quantum number (defined in +// SeQuant/domain/mbpt/space_qns.hpp). +// + +#ifndef SEQUANT_DOMAIN_MBPT_SPINOR_HPP +#define SEQUANT_DOMAIN_MBPT_SPINOR_HPP + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace sequant::mbpt { + +// ===================================================================== +// Kramers index annotations. +// +// Convention (mirrors the spin annotations using ↑/↓): +// - L'⇑' (U+21D1) = Kramers-up +// - L'⇓' (U+21D3) = Kramers-down +// Distinct from spin's ↑/↓ so the two QN families compose without +// label collision. +// ===================================================================== + +/// Extract the Kramers QN from a QuantumNumbersAttr. +inline Kramers to_kramers(const QuantumNumbersAttr& t) { + SEQUANT_ASSERT((t.to_int32() & mask_v) != 0); + return static_cast(t.to_int32() & mask_v); +} + +/// Strip Kramers QN bits from a QuantumNumbersAttr (mirrors +/// `spinannotation_remove(QuantumNumbersAttr)` in `spin.hpp`). +inline QuantumNumbersAttr kramers_annotation_remove( + const QuantumNumbersAttr& t) { + static_assert((~(~mask_v & ~bitset::reserved) & ~bitset::reserved) == + mask_v, + "Kramers bitmask uses reserved bits"); + return t.intersection(QuantumNumbersAttr(~mask_v & ~bitset::reserved)); +} + +/// Strip a trailing ⇑/⇓ Kramers annotation from a label. +template >>> +std::wstring kramers_annotation_remove(WS&& label) { + auto view = to_basic_string_view(label); + const auto has_annotation = !view.empty() && + (view.back() == L'⇑' || view.back() == L'⇓'); + return std::wstring{view.data(), + view.data() + view.size() - (has_annotation ? 1 : 0)}; +} + +/// Add a Kramers annotation (⇑ or ⇓) to a label (must not already +/// carry one). Kramers::any returns the label unchanged. +template >>> +std::wstring kramers_annotation_add(WS&& label, Kramers k) { + [[maybe_unused]] auto view = to_basic_string_view(label); + SEQUANT_ASSERT(view.empty() || + (view.back() != L'⇑' && view.back() != L'⇓')); + switch (k) { + case Kramers::any: + return std::wstring(std::forward(label)); + case Kramers::up: + return std::wstring(std::forward(label)) + L'⇑'; + case Kramers::down: + return std::wstring(std::forward(label)) + L'⇓'; + case Kramers::null: + SEQUANT_ABORT("Invalid Kramers quantum number"); + } + SEQUANT_UNREACHABLE; +} + +/// Replace any existing ⇑/⇓ annotation with @p k (or strip if Kramers::any). +template >>> +std::wstring kramers_annotation_replace(WS&& label, Kramers k) { + auto label_kf = kramers_annotation_remove(std::forward(label)); + return kramers_annotation_add(label_kf, k); +} + +/// Construct a Kramers-up Index from @p idx (preserves type, ordinal, +/// proto-indices; sets the Kramers QN bits and decorates the label). +[[nodiscard]] Index make_kramers_up(const Index& idx); + +/// Construct a Kramers-down Index from @p idx. +[[nodiscard]] Index make_kramers_dn(const Index& idx); + +/// Construct a Kramers-free (Kramers::any) Index from @p idx. +[[nodiscard]] Index make_kramers_free(const Index& idx); + +/// Apply an Index→Index replacement map to all tensors in @p expr. +/// Mirrors `append_spin` in `spin.hpp`; provided as a separate symbol +/// so callers don't drag in the spin tracer. +ExprPtr append_kramers(const ExprPtr& expr, + const container::map& index_replacements); + +// ===================================================================== +// Kramers tracer. +// ===================================================================== + +/// Rewrite an all-spinor tensor expression into a sum over canonical +/// Kramers-block representatives. Implements the closed-shell relativistic +/// MP2/CC algebra: enumerate per-tensor Kramers blocks, apply per-tensor +/// TRS canonicalization (with (-1)^k sign), and fold via the standard +/// canonicalize+rapid_simplify pipeline. Single entry handles MP2 energy, +/// CC residuals, and arbitrary operator expressions; expands any +/// antisymmetrizer (@c A) tensors internally before tracing. +/// +/// @param expr expression with indices in `IndexSpace`s carrying +/// `Kramers::any` (i.e., not yet specialized to up/down) +/// @return a Sum with canonical-orbit-block tensors; +/// tensor labels with conjugation get a `*` suffix +ExprPtr kramers_trace(const ExprPtr& expr); + +} // namespace sequant::mbpt + +#endif // SEQUANT_DOMAIN_MBPT_SPINOR_HPP From bdcdaaf653617e7edd244b6edd2602d6f071e558 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 13 May 2026 21:09:51 -0400 Subject: [PATCH 2/6] mbpt: kramers_trace TRS rewrite + conjugate/real_part/imaginary_part MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add complex-arithmetic primitives on tensor expressions, and use them inside kramers_trace to rewrite TRS-paired Kramers configurations as explicit conjugates of one another. Conjugation infrastructure (groundwork for Phase 4 quaternion decompose as well): * `toggle_conj_suffix(label)`: append/strip the `*` suffix that marks a tensor as complex-conjugated. (z*)* = z by toggling. * `conjugate(expr)`: recursively walks a SeQuant expression and conjugates Constants, toggles `*` on Tensors, and distributes over Products and Sums (Π zᵢ)* = Π zᵢ* and (Σ zᵢ)* = Σ zᵢ*. * `real_part(expr) = (expr + conjugate(expr)) / 2` and `imaginary_part(expr) = -i (expr - conjugate(expr)) / 2` — materialized as ordinary Sum/Product/Constant trees so any standard SeQuant evaluator can consume them; downstream callers can short- circuit the (T + T*) → 2 Re(T) collapse if their output is known to be real-scalar. kramers_trace now applies the TRS rewrite: * Enumerate the 2^N Kramers configurations as before. * For each TRS pair (X, X̄), only build the lex-smaller X term; materialize the X̄ partner as `conjugate(T_X)` (using the X Kramers labels with `*`-suffixed tensor labels). * Sign analysis: T_X̄ = Π_i tensor_i_X̄ = Π_i ((-1)^{k_i} tensor_i*_X); for any all-internal-index contraction (every barred index appears in an even number of tensors) the (-1)^{Σ k_i} factor is +1, so T_X̄ ≡ conjugate(T_X) without an explicit sign. Term count is unchanged (16 stays 16 for an MP2-shape input) but the Sum now exposes the conjugate-pair structure symbolically, letting evaluators short-circuit (one fetch + .conj()) and letting downstream real-scalar callers collapse pairs to `2 Re(...)`. The MVP enumeration path still falls out as a special case (when the underlying expression has no Kramers::any indices, kramers_trace returns it unchanged). --- SeQuant/domain/mbpt/spinor.cpp | 120 ++++++++++++++++++++++++++++++--- SeQuant/domain/mbpt/spinor.hpp | 45 +++++++++++++ 2 files changed, 154 insertions(+), 11 deletions(-) diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp index 3ee069218f..b8478dcdf6 100644 --- a/SeQuant/domain/mbpt/spinor.cpp +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -72,6 +72,80 @@ Index make_kramers_free(const Index& idx) { return detail::make_index_with_kramers(idx, mbpt::Kramers::any); } +// ===================================================================== +// Conjugation infrastructure (groundwork for Kramers tracer + quaternion +// decomposer). +// ===================================================================== + +std::wstring toggle_conj_suffix(std::wstring_view label) { + if (has_conj_suffix(label)) return std::wstring{label.substr(0, label.size() - 1)}; + return std::wstring{label} + L'*'; +} + +namespace { + +/// Apply conjugation to a single Tensor: rebuild it with toggled label +/// and identical bra/ket/aux/symmetry/etc. +ExprPtr conjugate_tensor(const Tensor& t) { + auto new_label = toggle_conj_suffix(t.label()); + // Tensor copy ctor + label edit isn't directly exposed; reconstruct. + return ex( + new_label, + bra(container::svector{t.bra().begin(), t.bra().end()}), + ket(container::svector{t.ket().begin(), t.ket().end()}), + aux(container::svector{t.aux().begin(), t.aux().end()}), + t.symmetry(), t.braket_symmetry(), t.column_symmetry()); +} + +} // namespace + +ExprPtr conjugate(const ExprPtr& expr) { + if (expr->is()) { + return conjugate_tensor(expr->as()); + } + if (expr->is()) { + auto v = expr->as().value(); + return ex(sequant::conj(v)); + } + if (expr->is()) { + throw std::runtime_error( + "conjugate: Variable conjugation not yet supported"); + } + if (expr->is()) { + auto const& prod = expr->as(); + auto p = std::make_shared(); + p->scale(sequant::conj(prod.scalar())); + for (auto&& f : prod) p->append(1, conjugate(f)); + return p; + } + if (expr->is()) { + auto s = std::make_shared(); + for (auto&& summand : *expr) s->append(conjugate(summand)); + return s; + } + throw std::runtime_error("conjugate: unsupported Expr type"); +} + +ExprPtr real_part(const ExprPtr& expr) { + // Re(z) = (z + z*) / 2 — represented as a literal Sum with a 1/2 prefactor. + auto s = std::make_shared(); + s->append(expr); + s->append(conjugate(expr)); + return ex(rational{1, 2}) * ExprPtr{s}; +} + +ExprPtr imaginary_part(const ExprPtr& expr) { + // Im(z) = (z - z*) / (2i) = -i (z - z*) / 2. + auto s = std::make_shared(); + s->append(expr); + // -1 * conj(expr): scale conj by -1 via a wrapping Product. + auto neg_conj = ex(-1) * conjugate(expr); + s->append(neg_conj); + // multiply by -i / 2: + return ex(rational{1, 2}) * + ex(Complex{0, -1}) * ExprPtr{s}; +} + ExprPtr append_kramers( const ExprPtr& expr, const container::map& index_replacements) { @@ -141,17 +215,25 @@ container::svector collect_kramers_indices(const ExprPtr& expr) { } // namespace ExprPtr kramers_trace(const ExprPtr& expr) { - // MVP (Phase 2b/i): enumerate the 2^N Kramers configurations of the - // expression's spinor (Kramers::any) indices and emit a Sum of substituted - // copies. Each summand is a Product (or Sum/Tensor) where every Kramers - // index has been pinned to up/down. + // Enumerate the 2^N Kramers configurations of the expression's spinor + // (Kramers::any) indices, then apply TRS pairing: configs (X, X̄) related + // by full bar-flip are complex conjugates of each other (per the + // (-1)^k full-bar identity, with the per-tensor signs cancelling for + // products in which every barred index appears in an even number of + // tensors — true for any all-internal-index contraction, MP2 included). // - // Per-tensor TRS canonicalization (the (-1)^k bit-table fold) and - // cross-tensor simplification land in subsequent commits. + // For each TRS pair we materialize only the "lex-smaller" configuration + // X explicitly; the X̄ partner is emitted as `conjugate(T_X)`. Tensors + // in the conj branch carry a `*` label suffix (per the convention in + // the conjugate/real_part/imaginary_part scaffolding above) which the + // evaluator dispatches into a `.conj()` call on the underlying numeric + // tensor. Term count is unchanged, but the symbolic structure exposes + // the conjugate-pair relationship for downstream simplification (e.g., + // `real_part` collapse in real-scalar callers like SQ_MP2). // - // The eventual evaluator (LCAOFactory) sees each summand as a concrete - // (g_block, t_block, ...) contraction; LCAOFactory's own canonical-storage - // layer (Phase 1.5) ensures we still pay only the orbit-rep integral cost. + // Standard SeQuant simplifications (dummy renaming, antisym fold, etc.) + // are NOT applied here — caller invokes `sequant::canonicalize` if it + // wants those. const auto indices = collect_kramers_indices(expr); const std::size_t n = indices.size(); @@ -161,9 +243,13 @@ ExprPtr kramers_trace(const ExprPtr& expr) { "kramers_trace: too many Kramers indices (would emit > 2^30 terms)"); } - auto result = std::make_shared(); const std::uint64_t n_configs = 1ull << n; + const std::uint64_t mask = n_configs - 1; + + // Pass 1: build the "X" term for the lex-smaller half of TRS pairs. + container::map x_terms; for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { + if (cfg > (mask ^ cfg)) continue; // partner already (or will be) handled container::map repl; for (std::size_t i = 0; i < n; ++i) { const bool is_down = (cfg >> i) & 1u; @@ -171,7 +257,19 @@ ExprPtr kramers_trace(const ExprPtr& expr) { : make_kramers_up(indices[i]); repl[indices[i]] = new_idx; } - result->append(append_kramers(expr, repl)); + x_terms[cfg] = append_kramers(expr, repl); + } + + // Pass 2: assemble the Sum. The X̄ partner of cfg=k is k^mask; whichever + // is lex-smaller is the "X" (real) term, the other is `conjugate(X)`. + auto result = std::make_shared(); + for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { + const std::uint64_t partner = mask ^ cfg; + if (cfg <= partner) { + result->append(x_terms.at(cfg)); + } else { + result->append(conjugate(x_terms.at(partner))); + } } return result; } diff --git a/SeQuant/domain/mbpt/spinor.hpp b/SeQuant/domain/mbpt/spinor.hpp index aec1b5fa31..6f9ba3e606 100644 --- a/SeQuant/domain/mbpt/spinor.hpp +++ b/SeQuant/domain/mbpt/spinor.hpp @@ -120,6 +120,51 @@ std::wstring kramers_annotation_replace(WS&& label, Kramers k) { ExprPtr append_kramers(const ExprPtr& expr, const container::map& index_replacements); +// ===================================================================== +// Complex-arithmetic operators on tensor expressions. +// +// Groundwork for the Kramers tracer (which uses TRS to rewrite barred +// configurations as conjugates of unbarred ones) and the upcoming +// quaternion decomposer (which separates real/imaginary parts). +// +// Convention: complex conjugation of a tensor is encoded as a `*` +// suffix on the tensor's label (e.g., `g` → `g*`, `t` → `t*`). +// Evaluator-side dispatch (e.g., MPQC's SQ_MP2 yielder) translates the +// suffix into a `.conj()` call on the underlying numeric tensor. +// ===================================================================== + +/// Test whether a tensor's label encodes complex conjugation +/// (trailing `*` suffix per the convention above). +inline bool has_conj_suffix(std::wstring_view label) { + return !label.empty() && label.back() == L'*'; +} + +/// Append `*` to a tensor's label (or strip if already present, since +/// (z*)* = z). Pure label rewrite — does NOT affect Symmetry tags or +/// IndexSpace QNs. Caller's responsibility to ensure the underlying +/// numeric tensor's complex conjugate is what's wanted. +[[nodiscard]] std::wstring toggle_conj_suffix(std::wstring_view label); + +/// Symbolic complex conjugation of an expression. +/// +/// Recursively walks @p expr and applies: +/// - Tensor: toggle its `*` suffix (z → z*, z* → z). +/// - Constant: numeric complex conjugate. +/// - Variable: not yet supported (throws). +/// - Product: conjugate the scalar AND each factor. (Π zᵢ)* = Π zᵢ*. +/// - Sum: conjugate each summand. (Σ zᵢ)* = Σ zᵢ*. +/// +/// `conjugate(conjugate(expr))` recovers @p expr. +[[nodiscard]] ExprPtr conjugate(const ExprPtr& expr); + +/// Symbolic real part: returns `(expr + conjugate(expr)) / 2`. +/// Currently materialized as that two-term Sum (no native `Re` node); +/// downstream callers can apply numerical optimizations after evaluation. +[[nodiscard]] ExprPtr real_part(const ExprPtr& expr); + +/// Symbolic imaginary part: returns `(expr - conjugate(expr)) / (2i)`. +[[nodiscard]] ExprPtr imaginary_part(const ExprPtr& expr); + // ===================================================================== // Kramers tracer. // ===================================================================== From a5adfacc84937c616c0fd00a7ae1ec1b46bf2e3e Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 13 May 2026 21:40:21 -0400 Subject: [PATCH 3/6] mbpt/spinor: kramers_trace mirrors open_shell_spintrace pipeline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the TRS-rewrite-only enumeration with the full open-shell spin-trace pipeline applied to Kramers configs: 1. Substitute K_up/K_dn for each index per the 2^N enumeration. 2. expand_antisymm (skip_spinsymm=false — Kramers blocks have no orthogonality, so every permutation contributes). 3. expand → flatten Products → rapid_simplify → canonicalize → rapid_simplify per term. 4. Final cross-product expand + canonicalize so dummy renames fold equivalents across Kramers configurations. The TRS pair rewrite via conjugate() is dropped from kramers_trace because it created NonSymm tensors with `*`-suffixed labels that canonicalize couldn't fold across Kramers configs (it cuts the term count in half but blocks the deeper antisym fold). The conjugate / real_part / imaginary_part scaffolding is kept as groundwork for the upcoming quaternion decomposer (Phase 4) and as a building block for caller-side `2 Re(...)` collapse in real-scalar consumers. End-to-end on h2o-sq_mp2-x2c-631g.json: SQ_MP2 produces a 20-term Sum (down from the 16-term raw enumeration after factoring in the antisym expansion's 16x blowup). Energy preserved to 1 ULP. The 20-term plateau reflects what canonicalize-via-dummy-rename can fold without post-expansion antisym recognition or point-group symmetry filters (which Visscher Eq. 28 relies on for its 2-term form). Going below 20 requires either a real-part collapse pass (caller-side) or non-expanded antisym (which loses the open-shell-style structure). --- SeQuant/domain/mbpt/spinor.cpp | 41 ++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 17 deletions(-) diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp index b8478dcdf6..039bc9d0f1 100644 --- a/SeQuant/domain/mbpt/spinor.cpp +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include #include @@ -246,10 +248,13 @@ ExprPtr kramers_trace(const ExprPtr& expr) { const std::uint64_t n_configs = 1ull << n; const std::uint64_t mask = n_configs - 1; - // Pass 1: build the "X" term for the lex-smaller half of TRS pairs. - container::map x_terms; + // Enumerate all 2^n Kramers configurations as separate Products. We + // mirror open_shell_spintrace's pipeline: substitute → expand_antisymm → + // expand → simplify → canonicalize, working entirely in expanded + // (NonSymm) tensor form so the standard SeQuant simplifier sees full + // index permutability for fold detection. + auto result = std::make_shared(); for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { - if (cfg > (mask ^ cfg)) continue; // partner already (or will be) handled container::map repl; for (std::size_t i = 0; i < n; ++i) { const bool is_down = (cfg >> i) & 1u; @@ -257,21 +262,23 @@ ExprPtr kramers_trace(const ExprPtr& expr) { : make_kramers_up(indices[i]); repl[indices[i]] = new_idx; } - x_terms[cfg] = append_kramers(expr, repl); - } - - // Pass 2: assemble the Sum. The X̄ partner of cfg=k is k^mask; whichever - // is lex-smaller is the "X" (real) term, the other is `conjugate(X)`. - auto result = std::make_shared(); - for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { - const std::uint64_t partner = mask ^ cfg; - if (cfg <= partner) { - result->append(x_terms.at(cfg)); - } else { - result->append(conjugate(x_terms.at(partner))); - } + auto term = append_kramers(expr, repl); + term = expand_antisymm(term, /*skip_spinsymm*/ false); + expand(term); + rapid_simplify(term); + canonicalize(term); + rapid_simplify(term); + result->append(term); } - return result; + // Final cross-product canonicalize: now that every term is in NonSymm + // expanded form, dummy renaming can fold equivalents across Kramers + // configurations. + ExprPtr r = result; + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + return r; } } // namespace sequant::mbpt From 435552528a1e56afe92d8e1fbb8f0beff6e07490 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Thu, 14 May 2026 10:59:39 -0400 Subject: [PATCH 4/6] mbpt/spinor: Kramers tracing, folding, and antisymm-recombination passes Adds the symbolic machinery to take a closed-shell all-spinor expression to a minimal mixed antisymm/non-antisymm Kramers-restricted form: - RealPart / ImagPart Expr nodes + conjugate / re / im helpers - flip_kramers + fold_conj_pairs: O(n) hash-bucketed TRS conjugate-pair fold, (A + A*) -> 2 Re(A) - antisymm_recombine: inverse of expand_antisymm; recombines direct/exchange NonSymm pairs into Kramers-restricted antisymm tensors, iterated to a fixed point - kramers_trace_cycles: cycle-driven tracer (contraction-graph cycle decomposition + multiset enumeration) - kramers_trace_burnside: orbit-representative enumeration under the antisymmetrizer index-permutation group, via bit-permutation orbits - cycle decomposition utilities (kramers_cycles, kramers_cycle_dump, cycle_canonical_*) for inspection and cross-validation All three tracers feed the same fold_conj_pairs -> antisymm_recombine pipeline and are value-equivalent; they coexist for cross-checking. --- SeQuant/domain/mbpt/spinor.cpp | 1059 +++++++++++++++++++++++++++++++- SeQuant/domain/mbpt/spinor.hpp | 500 +++++++++++++-- 2 files changed, 1463 insertions(+), 96 deletions(-) diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp index 039bc9d0f1..2ef1e053c1 100644 --- a/SeQuant/domain/mbpt/spinor.cpp +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -16,6 +16,7 @@ #include #include +#include #include namespace sequant::mbpt { @@ -80,7 +81,8 @@ Index make_kramers_free(const Index& idx) { // ===================================================================== std::wstring toggle_conj_suffix(std::wstring_view label) { - if (has_conj_suffix(label)) return std::wstring{label.substr(0, label.size() - 1)}; + if (has_conj_suffix(label)) + return std::wstring{label.substr(0, label.size() - 1)}; return std::wstring{label} + L'*'; } @@ -92,8 +94,7 @@ ExprPtr conjugate_tensor(const Tensor& t) { auto new_label = toggle_conj_suffix(t.label()); // Tensor copy ctor + label edit isn't directly exposed; reconstruct. return ex( - new_label, - bra(container::svector{t.bra().begin(), t.bra().end()}), + new_label, bra(container::svector{t.bra().begin(), t.bra().end()}), ket(container::svector{t.ket().begin(), t.ket().end()}), aux(container::svector{t.aux().begin(), t.aux().end()}), t.symmetry(), t.braket_symmetry(), t.column_symmetry()); @@ -128,29 +129,8 @@ ExprPtr conjugate(const ExprPtr& expr) { throw std::runtime_error("conjugate: unsupported Expr type"); } -ExprPtr real_part(const ExprPtr& expr) { - // Re(z) = (z + z*) / 2 — represented as a literal Sum with a 1/2 prefactor. - auto s = std::make_shared(); - s->append(expr); - s->append(conjugate(expr)); - return ex(rational{1, 2}) * ExprPtr{s}; -} - -ExprPtr imaginary_part(const ExprPtr& expr) { - // Im(z) = (z - z*) / (2i) = -i (z - z*) / 2. - auto s = std::make_shared(); - s->append(expr); - // -1 * conj(expr): scale conj by -1 via a wrapping Product. - auto neg_conj = ex(-1) * conjugate(expr); - s->append(neg_conj); - // multiply by -i / 2: - return ex(rational{1, 2}) * - ex(Complex{0, -1}) * ExprPtr{s}; -} - -ExprPtr append_kramers( - const ExprPtr& expr, - const container::map& index_replacements) { +ExprPtr append_kramers(const ExprPtr& expr, + const container::map& index_replacements) { auto add_to_tensor = [&](const Tensor& tensor) { auto t = std::make_shared(tensor); t->transform_indices(index_replacements); @@ -222,7 +202,7 @@ ExprPtr kramers_trace(const ExprPtr& expr) { // by full bar-flip are complex conjugates of each other (per the // (-1)^k full-bar identity, with the per-tensor signs cancelling for // products in which every barred index appears in an even number of - // tensors — true for any all-internal-index contraction, MP2 included). + // tensors — true for any all-internal-index contraction). // // For each TRS pair we materialize only the "lex-smaller" configuration // X explicitly; the X̄ partner is emitted as `conjugate(T_X)`. Tensors @@ -231,7 +211,7 @@ ExprPtr kramers_trace(const ExprPtr& expr) { // evaluator dispatches into a `.conj()` call on the underlying numeric // tensor. Term count is unchanged, but the symbolic structure exposes // the conjugate-pair relationship for downstream simplification (e.g., - // `real_part` collapse in real-scalar callers like SQ_MP2). + // `real_part` collapse in real-scalar callers). // // Standard SeQuant simplifications (dummy renaming, antisym fold, etc.) // are NOT applied here — caller invokes `sequant::canonicalize` if it @@ -246,7 +226,6 @@ ExprPtr kramers_trace(const ExprPtr& expr) { } const std::uint64_t n_configs = 1ull << n; - const std::uint64_t mask = n_configs - 1; // Enumerate all 2^n Kramers configurations as separate Products. We // mirror open_shell_spintrace's pipeline: substitute → expand_antisymm → @@ -258,8 +237,8 @@ ExprPtr kramers_trace(const ExprPtr& expr) { container::map repl; for (std::size_t i = 0; i < n; ++i) { const bool is_down = (cfg >> i) & 1u; - Index new_idx = is_down ? make_kramers_dn(indices[i]) - : make_kramers_up(indices[i]); + Index new_idx = + is_down ? make_kramers_dn(indices[i]) : make_kramers_up(indices[i]); repl[indices[i]] = new_idx; } auto term = append_kramers(expr, repl); @@ -281,4 +260,1022 @@ ExprPtr kramers_trace(const ExprPtr& expr) { return r; } +// ===================================================================== +// Conjugate-pair fold pass: detects TRS-related Products in a Sum and +// rewrites (A + B) → 2·Re(A) where B = flip_kramers(A). +// ===================================================================== + +ExprPtr flip_kramers(const ExprPtr& expr) { + // Walk all tensors and collect every Kramers-tagged index. For each + // distinct (label, ordinal, space-type) seen, build a partner index + // with the opposite Kramers state. Then apply the swap via the existing + // append_kramers replacement machinery. + container::map repl; + expr->visit( + [&](const ExprPtr& e) { + if (!e->is()) return; + for (const auto& idx : e->as().const_braket()) { + if (repl.contains(idx)) continue; + const auto qns_int = idx.space().qns().to_int32(); + const auto k_bits = qns_int & mask_v; + if (k_bits == static_cast(Kramers::up)) { + repl[idx] = make_kramers_dn(idx); + } else if (k_bits == static_cast(Kramers::down)) { + repl[idx] = make_kramers_up(idx); + } + // Kramers::any or no Kramers state — leave alone. + } + }, + /* recursive = */ true); + if (repl.empty()) return expr; + return append_kramers(expr, repl); +} + +namespace { + +/// Normalize a Product-like expression for structural comparison: clone, +/// apply canonicalize + rapid_simplify so dummy renaming and ordering +/// match other already-canonicalized summands. +ExprPtr normalize(const ExprPtr& e) { + ExprPtr c = e->clone(); + expand(c); + rapid_simplify(c); + canonicalize(c); + rapid_simplify(c); + return c; +} + +/// Compare two scalars as complex rationals. +bool scalar_eq(const Constant::scalar_type& a, const Constant::scalar_type& b) { + return a == b; +} + +/// Extract the leading scalar of a Product, or 1 if the expr is a bare +/// Tensor. +Constant::scalar_type product_scalar(const ExprPtr& e) { + if (e->is()) return e->as().scalar(); + return Constant::scalar_type{1}; +} + +/// Build a "factor-only" copy of a Product (drop the leading scalar). +/// Returns @p e unchanged if it isn't a Product. +ExprPtr product_factors_only(const ExprPtr& e) { + if (!e->is()) return e->clone(); + auto p = std::make_shared(); + p->scale(Constant::scalar_type{1}); + for (auto&& f : e->as()) p->append(1, f->clone()); + return p; +} + +} // namespace + +ExprPtr fold_conj_pairs(const ExprPtr& expr) { + if (!expr->is()) return expr; + auto const& sum = expr->as(); + const std::size_t n = sum.size(); + + // O(n) canonical-form hash bucketing. For each summand A, compute a + // TRS-invariant key by taking the lex-smaller of normalize(A) and + // normalize(flip_kramers(A)). A and its full-bar partner B = flip(A) + // both map to the same key, so they land in the same bucket. + // + // Within a bucket we re-verify structural equality (hash collisions are + // possible; equality is O(tensor-count) per pair). Because TRS pairs + // form orbits of size ≤ 2, buckets are tiny and per-bucket work is O(1) + // on average. Total cost: O(n) normalizations + O(n) flips + O(n) hash + // inserts, vs. O(n²) for naive pairwise search. + // + // Scalar handling: we also bucket by scalar relationship. Two summands + // with matching scalars (a, a) collapse to `2·Re(A)`; opposite scalars + // (a, -a) collapse to `2i·Im(A)`. Self-conjugate singletons emit + // `Re(A)` (real-valued). + struct BucketEntry { + std::size_t idx; + ExprPtr factors_normalized; // factor part, post-normalize + Constant::scalar_type scalar; // leading scalar prefactor + bool self_conj; // factors == flip(factors) + }; + container::map> buckets; + + std::vector originals; + originals.reserve(n); + + for (std::size_t i = 0; i < n; ++i) { + originals.push_back(sum.summand(i)->clone()); + auto factors = normalize(product_factors_only(originals[i])); + auto factors_flipped = normalize(flip_kramers(factors)); + const bool self_conj = (*factors == *factors_flipped); + // Canonical factor key: lex-min of (factors, factors_flipped). Both A + // and flip(A) produce the same key, so they collide in the same bucket. + ExprPtr canon = (*factors < *factors_flipped) ? factors : factors_flipped; + const auto h = canon->hash_value(); + buckets[h].push_back( + BucketEntry{i, factors, product_scalar(originals[i]), self_conj}); + } + + auto out = std::make_shared(); + std::vector used(n, false); + + // Helper to emit a folded summand based on scalar relationship. + auto emit_pair = [&](const BucketEntry& a, const BucketEntry& b) { + if (scalar_eq(a.scalar, b.scalar)) { + // (A + B) = 2 Re(A) + out->append(ex(rational{2}) * real_part(originals[a.idx])); + } else if (scalar_eq(a.scalar, -b.scalar)) { + // (A - B) = 2i Im(A) (B carries -A's scalar) + using cscalar = Constant::scalar_type; + out->append(ex(cscalar{Complex{0, 2}}) * + imaginary_part(originals[a.idx])); + } else { + // Same factor-canonical key but unrelated scalars — emit both raw. + out->append(originals[a.idx]); + out->append(originals[b.idx]); + } + used[a.idx] = used[b.idx] = true; + }; + + for (auto& [_, entries] : buckets) { + // Verify hash-bucket entries actually share the same canonical factor + // form (defend against hash collisions): partition by structural + // equality on factors_normalized OR its flip. + std::vector taken(entries.size(), false); + for (std::size_t i = 0; i < entries.size(); ++i) { + if (taken[i] || used[entries[i].idx]) continue; + const auto& a = entries[i]; + // Self-conjugate singleton in this bucket? + if (a.self_conj) { + // Look for another self-conj with same scalar — if none, emit Re(A). + bool merged = false; + for (std::size_t j = i + 1; j < entries.size(); ++j) { + if (taken[j] || used[entries[j].idx]) continue; + const auto& b = entries[j]; + if (b.self_conj && *b.factors_normalized == *a.factors_normalized) { + // Two self-conj copies → 2·Re(A) (or scalar-aware emit). + emit_pair(a, b); + taken[i] = taken[j] = true; + merged = true; + break; + } + } + if (!merged) { + out->append(real_part(originals[a.idx])); + used[a.idx] = true; + taken[i] = true; + } + continue; + } + // Look for the partner: another entry whose factors equal flip(A). + // We don't have flip(A) materialized here — use the structural test + // (a.factors and b.factors land in the same bucket only when one is + // the flip of the other). + bool paired = false; + for (std::size_t j = i + 1; j < entries.size(); ++j) { + if (taken[j] || used[entries[j].idx]) continue; + const auto& b = entries[j]; + // Confirm they're actually a TRS pair: b.factors should equal + // flip(a.factors). + auto a_flip = normalize(flip_kramers(a.factors_normalized)); + if (*b.factors_normalized == *a_flip) { + emit_pair(a, b); + taken[i] = taken[j] = true; + paired = true; + break; + } + } + if (!paired) { + out->append(originals[a.idx]); + used[a.idx] = true; + taken[i] = true; + } + } + } + + return out; +} + +// ===================================================================== +// Cycle decomposition. See spinor.hpp for the design summary. +// ===================================================================== + +container::svector kramers_cycles(const Product& product) { + // Collect (tensor_idx, braket_pos, Index) for every position in every + // Tensor factor. Indices that aren't Tensor are skipped (Constants, + // Variables don't contribute to cycles). + struct Slot { + std::size_t tensor_idx; + std::size_t braket_pos; + Index idx; + }; + container::svector slots; + std::size_t tensor_counter = 0; + for (auto&& f : product) { + if (!f->is()) { + ++tensor_counter; + continue; + } + auto const& t = f->as(); + std::size_t pos = 0; + for (auto const& idx : t.const_braket()) { + slots.push_back(Slot{tensor_counter, pos, idx}); + ++pos; + } + ++tensor_counter; + } + const std::size_t N = slots.size(); + container::svector visited(N, false); + + // Pre-build two adjacency maps: + // - intra-tensor: bra[k] ↔ ket[k] of the same tensor (positions 0↔rank, + // 1↔rank+1, ... where rank = bra_rank). + // - inter-tensor: same Index occurrence in two distinct slots. + container::svector intra_partner(N, N); // N = "no partner" + container::svector inter_partner(N, N); + + // intra-tensor pairing per tensor: walk positions and pair k ↔ k+rank/2 + // within each tensor's slots. We rebuild rank by counting consecutive + // slots with the same tensor_idx. + for (std::size_t i = 0; i < N;) { + std::size_t j = i; + while (j < N && slots[j].tensor_idx == slots[i].tensor_idx) ++j; + const std::size_t tensor_rank = j - i; // total bra+ket count + SEQUANT_ASSERT(tensor_rank % 2 == 0); + const std::size_t half = tensor_rank / 2; + for (std::size_t k = 0; k < half; ++k) { + intra_partner[i + k] = i + half + k; + intra_partner[i + half + k] = i + k; + } + i = j; + } + + // inter-tensor pairing: same Index in two distinct positions across the + // whole Product. Use Index-equality (which compares full label). + for (std::size_t i = 0; i < N; ++i) { + if (inter_partner[i] != N) continue; + for (std::size_t j = i + 1; j < N; ++j) { + if (inter_partner[j] != N) continue; + if (slots[j].idx == slots[i].idx) { + inter_partner[i] = j; + inter_partner[j] = i; + break; + } + } + } + + // Walk cycles. From each unvisited slot, alternate intra → inter → + // intra → ... edges until returning to the start. + container::svector cycles; + for (std::size_t start = 0; start < N; ++start) { + if (visited[start]) continue; + ContractionCycle cyc; + std::size_t cur = start; + bool take_intra = true; // first hop is intra-tensor + while (!visited[cur]) { + visited[cur] = true; + cyc.nodes.push_back(ContractionCycle::Node{ + slots[cur].tensor_idx, slots[cur].braket_pos, slots[cur].idx}); + const std::size_t next = + take_intra ? intra_partner[cur] : inter_partner[cur]; + if (next == N) break; // dangling (shouldn't happen for closed expr) + cur = next; + take_intra = !take_intra; + } + cycles.push_back(std::move(cyc)); + } + + return cycles; +} + +namespace { + +/// Format a cycle's Kramers labelling as a wstring like "UAUA" or "UDDU" +/// where each char is the visited index's Kramers state. +std::wstring cycle_kramers_label(const ContractionCycle& c) { + std::wstring s; + s.reserve(c.nodes.size()); + for (auto const& n : c.nodes) { + const auto qns_int = n.idx.space().qns().to_int32(); + const auto k_bits = qns_int & mask_v; + if (k_bits == static_cast(Kramers::up)) { + s += L'U'; + } else if (k_bits == static_cast(Kramers::down)) { + s += L'B'; + } else { + s += L'?'; + } + } + return s; +} + +/// Number of canonical Kramers patterns of a length-L cycle under +/// cyclic-rotation symmetry. By Burnside: (1/L) Σ_{d | L} φ(L/d) · 2^d. +std::size_t canonical_pattern_count_cyclic(std::size_t L) { + if (L == 0) return 1; + // φ(n) + auto phi = [](std::size_t n) { + std::size_t result = n; + for (std::size_t p = 2; p * p <= n; ++p) { + if (n % p == 0) { + while (n % p == 0) n /= p; + result -= result / p; + } + } + if (n > 1) result -= result / n; + return result; + }; + std::size_t total = 0; + for (std::size_t d = 1; d <= L; ++d) { + if (L % d == 0) { + total += phi(L / d) * (std::size_t{1} << d); + } + } + return total / L; +} + +} // namespace + +void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os) { + auto dump_product = [&](const Product& p, std::size_t pidx) { + os << L" Product[" << pidx << L"]:"; + for (auto&& f : p) { + if (f->is()) { + auto const& t = f->as(); + os << L" " << t.label(); + } + } + os << L"\n"; + auto cycles = kramers_cycles(p); + os << L" " << cycles.size() << L" cycle(s):\n"; + for (std::size_t ci = 0; ci < cycles.size(); ++ci) { + auto const& c = cycles[ci]; + // Distinct indices visited. + container::set distinct; + for (auto const& n : c.nodes) distinct.insert(n.idx); + const auto canon = canonical_pattern_count_cyclic(c.nodes.size()); + os << L" cycle[" << ci << L"] len=" << c.nodes.size() + << L" distinctIdx=" << distinct.size() << L" labels=`" + << cycle_kramers_label(c) << L"`" << L" cyclic-canonical-#patterns=" + << canon << L" walk="; + bool first = true; + for (auto const& n : c.nodes) { + if (!first) os << L"→"; + os << L"t" << n.tensor_idx << L"/" << n.braket_pos << L":" + << n.idx.label(); + first = false; + } + os << L"\n"; + } + }; + + if (expr->is()) { + dump_product(expr->as(), 0); + } else if (expr->is()) { + auto const& s = expr->as(); + for (std::size_t i = 0; i < s.size(); ++i) { + auto const& term = s.summand(i); + if (term->is()) dump_product(term->as(), i); + } + } +} + +std::wstring cycle_canonical_label(const ContractionCycle& c) { + // Build raw label, then take lex-min over all cyclic rotations. + const auto raw = cycle_kramers_label(c); + if (raw.empty()) return raw; + std::wstring best = raw; + std::wstring rot = raw; + for (std::size_t k = 1; k < raw.size(); ++k) { + // rotate by 1: move first char to end + std::rotate(rot.begin(), rot.begin() + 1, rot.end()); + if (rot < best) best = rot; + } + return best; +} + +std::wstring cycle_canonical_signature(const Product& product) { + auto cycles = kramers_cycles(product); + container::svector labels; + labels.reserve(cycles.size()); + for (auto const& c : cycles) labels.push_back(cycle_canonical_label(c)); + std::sort(labels.begin(), labels.end()); + std::wstring sig; + for (std::size_t i = 0; i < labels.size(); ++i) { + if (i > 0) sig += L'|'; + sig += labels[i]; + } + return sig; +} + +ExprPtr cycle_canonical_fold(const ExprPtr& expr) { + if (!expr->is()) return expr; + auto const& sum = expr->as(); + // Bucket Sum summands by cycle_canonical_signature. Within each bucket + // the contractions are numerically equal up to scalar prefactor; sum + // the prefactors and emit a single representative per bucket. + container::map> buckets; + std::vector originals; + originals.reserve(sum.size()); + for (std::size_t i = 0; i < sum.size(); ++i) { + originals.push_back(sum.summand(i)->clone()); + if (!originals[i]->is()) { + // Non-Product summands (Tensor, Constant, Variable) — emit a + // distinct bucket per index so they pass through unchanged. + buckets[L"__pass" + std::to_wstring(i)].push_back(i); + continue; + } + const auto sig = cycle_canonical_signature(originals[i]->as()); + buckets[sig].push_back(i); + } + + auto out = std::make_shared(); + for (auto& [sig, idxs] : buckets) { + if (idxs.size() == 1) { + out->append(originals[idxs[0]]); + continue; + } + // Sum scalars; reuse first member's tensor structure as the + // representative. NOTE: this assumes the bucketed Products are + // numerically equal as contractions, which holds when the cycles' + // canonical Kramers labels match — that's the bucket-key meaning. + Constant::scalar_type total{0}; + for (auto i : idxs) { + total += product_scalar(originals[i]); + } + if (total == Constant::scalar_type{0}) continue; + auto rep = product_factors_only(originals[idxs.front()]); + out->append(ex(total) * rep); + } + return out; +} + +// ===================================================================== +// Standalone cycle-driven Kramers tracer. See spinor.hpp. +// ===================================================================== + +namespace { + +/// Distinct Kramers-eligible (`Kramers::any`) indices visited by a cycle, +/// in first-encounter order. +container::svector cycle_distinct_indices(const ContractionCycle& c) { + container::svector out; + container::set seen; + for (auto const& n : c.nodes) { + const auto qns = n.idx.space().qns().to_int32(); + if ((qns & mask_v) != static_cast(Kramers::any)) + continue; // already specialized or not Kramers — skip + if (seen.insert(n.idx).second) out.push_back(n.idx); + } + return out; +} + +/// Structural shape signature of a cycle: a string encoding the walk's +/// IndexSpace-type sequence, made rotation-invariant by taking the +/// lex-min rotation. Two cycles with the same shape are interchangeable +/// under the expression's tensor (anti)symmetries — the caller is +/// responsible for ensuring those symmetries actually hold. +std::wstring cycle_shape(const ContractionCycle& c) { + std::wstring raw; + for (auto const& n : c.nodes) { + // encode IndexSpace base type as a hex digit (stable per registry) + raw += static_cast(L'a' + (n.idx.space().type().to_int32() & 0xF)); + } + if (raw.empty()) return raw; + std::wstring best = raw, rot = raw; + for (std::size_t k = 1; k < raw.size(); ++k) { + std::rotate(rot.begin(), rot.begin() + 1, rot.end()); + if (rot < best) best = rot; + } + return best; +} + +/// All 2^n Kramers labelings of a cycle's @p n distinct indices, each a +/// bitset (bit i = 1 ⇔ index i is Kramers-down). +container::svector cycle_labelings(std::size_t n_distinct) { + container::svector out; + const std::uint32_t n_cfg = 1u << n_distinct; + out.reserve(n_cfg); + for (std::uint32_t c = 0; c < n_cfg; ++c) out.push_back(c); + return out; +} + +/// Multinomial coefficient m! / ∏ count_i! for a multiset assignment. +double multiset_multiplicity(const container::svector& counts) { + auto factorial = [](std::size_t k) { + double f = 1; + for (std::size_t i = 2; i <= k; ++i) f *= static_cast(i); + return f; + }; + std::size_t total = 0; + for (auto c : counts) total += c; + double num = factorial(total); + for (auto c : counts) num /= factorial(c); + return num; +} + +} // namespace + +ExprPtr kramers_trace_cycles(const ExprPtr& expr) { + // Sum: recurse over summands. + if (expr->is()) { + auto sum_result = std::make_shared(); + for (auto&& s : *expr) sum_result->append(kramers_trace_cycles(s)); + ExprPtr r = sum_result; + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + return r; + } + // Constant / Variable: no Kramers indices, pass through. + if (expr->is() || expr->is()) return expr; + + // Normalize to a Product. + Product prod; + if (expr->is()) { + prod = expr->as(); + } else if (expr->is()) { + prod.append(1, expr, Product::Flatten::No); + } else { + throw std::runtime_error("kramers_trace_cycles: unsupported Expr type"); + } + + // 1. Decompose into cycles. + const auto cycles = kramers_cycles(prod); + const std::size_t n_cycles = cycles.size(); + if (n_cycles == 0) return expr; + + // Per-cycle distinct Kramers indices + labelings. + std::vector> cyc_indices(n_cycles); + std::vector> cyc_labelings(n_cycles); + for (std::size_t c = 0; c < n_cycles; ++c) { + cyc_indices[c] = cycle_distinct_indices(cycles[c]); + cyc_labelings[c] = cycle_labelings(cyc_indices[c].size()); + } + + // 2. Group cycles into structural-equivalence classes by shape. + container::map> shape_classes; + for (std::size_t c = 0; c < n_cycles; ++c) + shape_classes[cycle_shape(cycles[c])].push_back(c); + + // 3. Per class, enumerate Kramers-labeling MULTISETS (sorted tuples of + // per-cycle labelings) with their multiplicity. A class of m cycles + // each with L labelings yields C(L+m-1, m) multisets. + // + // Represent a class assignment as a sorted vector of + // length m (one labeling index per cycle in the class). Multiplicity + // = multinomial(counts of each distinct labeling). + struct ClassConfig { + container::svector class_ids; // cycle indices in this class + // each entry: (sorted labeling-assignment, multiplicity) + container::svector, double>> + multisets; + }; + container::svector class_configs; + for (auto& [shape, cyc_ids] : shape_classes) { + ClassConfig cc; + cc.class_ids = cyc_ids; + const std::size_t m = cyc_ids.size(); + // All cycles in a class share labeling-set size (same shape ⇒ same + // #distinct indices). Use the first cycle's labeling set. + const auto& labels = cyc_labelings[cyc_ids.front()]; + const std::size_t L = labels.size(); + // Enumerate sorted m-tuples (multisets) from L labelings. + container::svector pick(m, 0); + std::function rec = [&](std::size_t pos, + std::size_t start) { + if (pos == m) { + // pick[] holds indices into `labels`, non-decreasing. + container::svector assignment(m); + for (std::size_t i = 0; i < m; ++i) assignment[i] = labels[pick[i]]; + // multiplicity = multinomial over runs of equal pick values + container::svector counts; + std::size_t run = 1; + for (std::size_t i = 1; i <= m; ++i) { + if (i < m && pick[i] == pick[i - 1]) { + ++run; + } else { + counts.push_back(run); + run = 1; + } + } + cc.multisets.emplace_back(assignment, multiset_multiplicity(counts)); + return; + } + for (std::size_t v = start; v < L; ++v) { + pick[pos] = v; + rec(pos + 1, v); + } + }; + rec(0, 0); + class_configs.push_back(std::move(cc)); + } + + // 4. Cartesian product across classes: each global config picks one + // multiset per class. Build the per-index Kramers replacement, + // substitute, expand_antisymm, canonicalize, scale by the product + // of per-class multiplicities. + auto result = std::make_shared(); + const std::size_t n_classes = class_configs.size(); + container::svector sel(n_classes, 0); + std::function build = [&](std::size_t ci) { + if (ci == n_classes) { + // Assemble the full per-index Kramers replacement. + container::map repl; + double mult = 1.0; + for (std::size_t c = 0; c < n_classes; ++c) { + const auto& cc = class_configs[c]; + const auto& [assignment, m] = cc.multisets[sel[c]]; + mult *= m; + for (std::size_t k = 0; k < cc.class_ids.size(); ++k) { + const std::size_t cyc = cc.class_ids[k]; + const std::uint32_t labeling = assignment[k]; + for (std::size_t d = 0; d < cyc_indices[cyc].size(); ++d) { + const bool is_dn = (labeling >> d) & 1u; + const Index& src = cyc_indices[cyc][d]; + repl[src] = is_dn ? make_kramers_dn(src) : make_kramers_up(src); + } + } + } + auto term = append_kramers(prod.clone(), repl); + term = expand_antisymm(term, /*skip_spinsymm*/ false); + expand(term); + rapid_simplify(term); + canonicalize(term); + rapid_simplify(term); + if (mult != 1.0) { + term = ex(rational{static_cast(mult)}) * term; + expand(term); + rapid_simplify(term); + } + result->append(term); + return; + } + for (std::size_t s = 0; s < class_configs[ci].multisets.size(); ++s) { + sel[ci] = s; + build(ci + 1); + } + }; + build(0); + + // Final cross-term canonicalize so equivalent multiset terms fold. + ExprPtr r = result; + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + return r; +} + +// ===================================================================== +// Antisymmetrizer recombination (Phase 2.11). Inverse of expand_antisymm. +// See spinor.hpp for the algorithm. +// ===================================================================== + +namespace { + +/// A parsed summand: total scalar, wrapper kind (0 raw / 1 Re / 2 Im), and +/// the tensor factors. The inner Product's scalar is folded into `scalar`. +struct RecombEntry { + Constant::scalar_type scalar{1}; + int wrapper = 0; + container::svector tensors; // each is a Tensor ExprPtr + bool alive = true; +}; + +/// Decompose a Sum summand into a RecombEntry. Recognized shapes: +/// Tensor | Product(s; Tensor...) | Re/Im(inner) | Product(s; Re/Im(inner)) +/// where inner = Tensor | Product(s; Tensor...). +bool parse_recomb_entry(const ExprPtr& summand, RecombEntry& e) { + e = RecombEntry{}; + ExprPtr cur = summand; + // optional outer `Constant * (Re|Im wrapper)` + if (cur->is() && cur->as().size() == 1) { + auto const& p = cur->as(); + auto const& f0 = p.factor(0); + if (f0->is() || f0->is()) { + e.scalar = e.scalar * p.scalar(); + cur = f0; + } + } + if (cur->is()) { + e.wrapper = 1; + cur = cur->as().inner(); + } else if (cur->is()) { + e.wrapper = 2; + cur = cur->as().inner(); + } + if (cur->is()) { + e.tensors.push_back(cur); + return true; + } + if (cur->is()) { + auto const& p = cur->as(); + e.scalar = e.scalar * p.scalar(); + for (auto&& f : p) { + if (!f->is()) return false; + e.tensors.push_back(f); + } + return true; + } + return false; +} + +/// Detect whether @p t2 is @p t1 with a single column transposition. +/// Returns 1 if t2 = bra-swap(t1), 2 if t2 = ket-swap(t1), 0 otherwise. +/// Rank-4 (2+2) only for now. Both indices compared by full Index identity +/// (recombinable terms share dummy names post-canonicalize). +int detect_single_column_swap(const Tensor& t1, const Tensor& t2) { + if (t1.label() != t2.label()) return 0; + if (t1.bra_rank() != 2 || t1.ket_rank() != 2) return 0; + if (t2.bra_rank() != 2 || t2.ket_rank() != 2) return 0; + if (t1.symmetry() != t2.symmetry()) return 0; + const auto& b1 = t1.bra(); + const auto& k1 = t1.ket(); + const auto& b2 = t2.bra(); + const auto& k2 = t2.ket(); + const bool bra_same = (b1.at(0) == b2.at(0) && b1.at(1) == b2.at(1)); + const bool bra_swap = (b1.at(0) == b2.at(1) && b1.at(1) == b2.at(0)); + const bool ket_same = (k1.at(0) == k2.at(0) && k1.at(1) == k2.at(1)); + const bool ket_swap = (k1.at(0) == k2.at(1) && k1.at(1) == k2.at(0)); + if (bra_swap && ket_same) return 1; + if (ket_swap && bra_same) return 2; + return 0; +} + +/// Re-tag a Tensor with Symmetry::Antisymm (keeping label/indices/other +/// symmetries). Used to materialize a recombined antisymmetrized tensor. +ExprPtr make_antisymm(const Tensor& t) { + return ex( + t.label(), bra(container::svector{t.bra().begin(), t.bra().end()}), + ket(container::svector{t.ket().begin(), t.ket().end()}), + aux(container::svector{t.aux().begin(), t.aux().end()}), + Symmetry::Antisymm, t.braket_symmetry(), t.column_symmetry()); +} + +/// Try to recombine two entries into one antisymmetrized entry. +/// Succeeds when, for some tensor position p: all other positions match +/// structurally, t2[p] is a single column-swap of t1[p], and the scalars +/// satisfy B.scalar == -A.scalar (single transposition → odd → sign -1). +/// Result: A.scalar * (... A.tensors[p] tagged Antisymm ...). +std::optional try_recombine(const RecombEntry& A, + const RecombEntry& B) { + if (A.wrapper != B.wrapper) return std::nullopt; + if (A.tensors.size() != B.tensors.size()) return std::nullopt; + if (!(B.scalar == -A.scalar)) return std::nullopt; + const std::size_t n = A.tensors.size(); + for (std::size_t p = 0; p < n; ++p) { + bool others_match = true; + for (std::size_t q = 0; q < n; ++q) { + if (q == p) continue; + if (!(*A.tensors[q] == *B.tensors[q])) { + others_match = false; + break; + } + } + if (!others_match) continue; + if (!A.tensors[p]->is() || !B.tensors[p]->is()) continue; + const int swap = detect_single_column_swap(A.tensors[p]->as(), + B.tensors[p]->as()); + if (swap == 0) continue; + RecombEntry R = A; + R.tensors[p] = make_antisymm(A.tensors[p]->as()); + return R; + } + return std::nullopt; +} + +} // namespace + +ExprPtr antisymm_recombine(const ExprPtr& expr) { + if (!expr->is()) return expr; + auto const& sum = expr->as(); + + std::vector entries; + entries.reserve(sum.size()); + for (auto&& s : sum) { + RecombEntry e; + if (!parse_recomb_entry(s, e)) return expr; // unrecognized — bail safely + entries.push_back(std::move(e)); + } + + // Iterate to a fixed point: each pass merges one recombinable pair. A + // second pass over the once-recombined terms catches doubly-antisymm + // blocks (the merged antisymm tensor becomes the "fixed" tensor). + bool changed = true; + while (changed) { + changed = false; + for (std::size_t i = 0; i < entries.size() && !changed; ++i) { + if (!entries[i].alive) continue; + for (std::size_t j = i + 1; j < entries.size() && !changed; ++j) { + if (!entries[j].alive) continue; + auto rec = try_recombine(entries[i], entries[j]); + if (!rec) rec = try_recombine(entries[j], entries[i]); + if (rec) { + entries[i] = *rec; + entries[j].alive = false; + changed = true; + } + } + } + } + + auto out = std::make_shared(); + for (auto const& e : entries) { + if (!e.alive) continue; + auto prod = std::make_shared(); + for (auto const& t : e.tensors) prod->append(1, t, Product::Flatten::No); + ExprPtr inner = prod; + if (e.wrapper == 1) + inner = ex(inner); + else if (e.wrapper == 2) + inner = ex(inner); + out->append(ex(e.scalar) * inner); + } + return out; +} + +// ===================================================================== +// Burnside-enumeration Kramers tracer (Phase 2.12 prototype). +// See spinor.hpp for the algorithm summary. +// ===================================================================== + +namespace { + +/// Collect the Tensor factors of a Product (or wrap a bare Tensor). +container::svector product_tensors(const Product& p) { + container::svector out; + for (auto&& f : p) { + if (f->is()) out.push_back(&f->as()); + } + return out; +} + +/// True iff indices @p p and @p q both sit in the same antisymmetric +/// bra OR ket group of Tensor @p t (i.e. t is Antisymm and {p,q} ⊆ bra +/// or {p,q} ⊆ ket). +bool same_antisymm_group(const Tensor& t, const Index& p, const Index& q) { + if (t.symmetry() != Symmetry::Antisymm) return false; + auto in = [](const auto& range, const Index& x) { + return std::find(range.begin(), range.end(), x) != range.end(); + }; + const bool p_bra = in(t.bra(), p), q_bra = in(t.bra(), q); + const bool p_ket = in(t.ket(), p), q_ket = in(t.ket(), q); + return (p_bra && q_bra) || (p_ket && q_ket); +} + +/// True iff index @p p appears anywhere in Tensor @p t. +bool tensor_has_index(const Tensor& t, const Index& p) { + for (auto&& idx : t.const_braket()) + if (idx == p) return true; + return false; +} + +} // namespace + +ExprPtr kramers_trace_burnside(const ExprPtr& expr) { + if (expr->is()) { + auto s = std::make_shared(); + for (auto&& t : *expr) s->append(kramers_trace_burnside(t)); + ExprPtr r = s; + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + return r; + } + if (expr->is() || expr->is()) return expr; + + Product prod; + if (expr->is()) { + prod = expr->as(); + } else if (expr->is()) { + prod.append(1, expr, Product::Flatten::No); + } else { + throw std::runtime_error("kramers_trace_burnside: unsupported Expr type"); + } + + const auto indices = collect_kramers_indices(expr); + const std::size_t n = indices.size(); + if (n == 0) return expr; + if (n > 30) { + throw std::runtime_error( + "kramers_trace_burnside: too many Kramers indices (> 2^30)"); + } + const auto tensors = product_tensors(prod); + + // --- 1. Find transposition generators of the index-permutation + // symmetry group. (p,q) qualifies iff in every tensor that + // contains both they share an antisymmetric bra/ket group, they + // co-occur in at least one tensor, and the sign is +1. + container::svector> generators; + for (std::size_t a = 0; a < n; ++a) { + for (std::size_t b = a + 1; b < n; ++b) { + std::size_t n_common = 0; + bool ok = true; + for (const Tensor* t : tensors) { + const bool ha = tensor_has_index(*t, indices[a]); + const bool hb = tensor_has_index(*t, indices[b]); + if (ha && hb) { + ++n_common; + if (!same_antisymm_group(*t, indices[a], indices[b])) { + ok = false; + break; + } + } else if (ha != hb) { + // appears in only one slot of this tensor — a transposition + // would move it to a different tensor position; not a symmetry + ok = false; + break; + } + } + if (ok && n_common > 0 && (n_common % 2 == 0)) { + generators.emplace_back(a, b); + } + } + } + + // --- 2. Build the permutation group as a BFS closure over generators. + // Each element is a permutation of {0..n-1} (svector). + auto apply_transposition = + [](container::svector perm, std::size_t a, + std::size_t b) -> container::svector { + std::swap(perm[a], perm[b]); + return perm; + }; + container::svector> group; + { + container::svector identity(n); + for (std::size_t k = 0; k < n; ++k) identity[k] = k; + container::set> seen; + container::svector> frontier{identity}; + seen.insert(identity); + while (!frontier.empty()) { + container::svector> next; + for (auto const& g : frontier) { + group.push_back(g); + for (auto const& [a, b] : generators) { + auto ng = apply_transposition(g, a, b); + if (seen.insert(ng).second) next.push_back(ng); + } + } + frontier = std::move(next); + } + } + + // --- 3. Burnside orbit enumeration over the 2^n bit-string configs. + // A permutation acts: (perm·cfg) has bit perm[k] = bit k of cfg. + auto apply_perm = [](std::uint64_t cfg, + const container::svector& perm) { + std::uint64_t out = 0; + for (std::size_t k = 0; k < perm.size(); ++k) + if ((cfg >> k) & 1u) out |= (std::uint64_t{1} << perm[k]); + return out; + }; + const std::uint64_t n_configs = std::uint64_t{1} << n; + container::svector visited(n_configs, false); + // each rep: (config, orbit_size) + container::svector> reps; + for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { + if (visited[cfg]) continue; + container::set orbit; + for (auto const& g : group) orbit.insert(apply_perm(cfg, g)); + for (auto o : orbit) visited[o] = true; + // representative = lex-min of the orbit + const std::uint64_t rep = *orbit.begin(); + reps.emplace_back(rep, orbit.size()); + } + + // --- 4. Per orbit rep: substitute Kramers labels, expand_antisymm, + // canonicalize, scale by orbit size. + auto result = std::make_shared(); + for (auto const& [cfg, orbit_size] : reps) { + container::map repl; + for (std::size_t k = 0; k < n; ++k) { + const bool is_dn = (cfg >> k) & 1u; + repl[indices[k]] = + is_dn ? make_kramers_dn(indices[k]) : make_kramers_up(indices[k]); + } + auto term = append_kramers(expr, repl); + term = expand_antisymm(term, /*skip_spinsymm*/ false); + expand(term); + rapid_simplify(term); + canonicalize(term); + rapid_simplify(term); + if (orbit_size != 1) { + term = + ex(rational{static_cast(orbit_size)}) * term; + expand(term); + rapid_simplify(term); + } + result->append(term); + } + + // --- 5. Final cross-term canonicalize. + ExprPtr r = result; + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + return r; +} + } // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/spinor.hpp b/SeQuant/domain/mbpt/spinor.hpp index 6f9ba3e606..9d90985125 100644 --- a/SeQuant/domain/mbpt/spinor.hpp +++ b/SeQuant/domain/mbpt/spinor.hpp @@ -1,21 +1,19 @@ // // Spinor decomposition passes for relativistic 2-component theories. // -// This module provides two symbolic transforms over tensor expressions: +// This module provides symbolic transforms over tensor expressions: // // 1. kramers_trace(expr): rewrite a closed-shell expression with // all-spinor (Kramers::any) indices into a sum over canonical // Kramers-block representatives (one per orbit under -// {column-swap, Hermitian, time-reversal-with-(-1)^k}). Folds -// the 2^n raw Kramers blocks of a rank-n tensor into ~5 unique -// classes for rank-4, with conjugation/sign tracking baked into -// the per-tensor lookup. +// {column-swap, Hermitian, time-reversal-with-(-1)^k}), with +// conjugation/sign tracking baked into the per-tensor lookup. // -// 2. quaternion_decompose(expr): (Phase 4 — not yet implemented) -// further decompose Kramers-up indices into the four real -// (s, x, y, z) quaternion components (Helmich-Paris, Repisky, -// Visscher 2016 Eq. 21), aggressively folding via Pauli/quaternion -// algebra to yield ≤ ~10 surviving terms for an MP2-style input. +// 2. quaternion_decompose(expr): (not yet implemented) further +// decompose Kramers-up indices into the four real (s, x, y, z) +// quaternion components (Helmich-Paris, Repisky, Visscher 2016 +// Eq. 21), aggressively folding via Pauli/quaternion algebra to a +// compact surviving-term set. // // Both passes operate on the canonical SeQuant Tensor representation // and consume its `Kramers` index quantum number (defined in @@ -48,41 +46,53 @@ namespace sequant::mbpt { // label collision. // ===================================================================== -/// Extract the Kramers QN from a QuantumNumbersAttr. +/// @brief Extracts the Kramers quantum number from a QuantumNumbersAttr. +/// @param t the quantum-numbers attribute to read +/// @return the Kramers QN encoded in @p t +/// @pre @p t has at least one Kramers bit set inline Kramers to_kramers(const QuantumNumbersAttr& t) { SEQUANT_ASSERT((t.to_int32() & mask_v) != 0); return static_cast(t.to_int32() & mask_v); } -/// Strip Kramers QN bits from a QuantumNumbersAttr (mirrors -/// `spinannotation_remove(QuantumNumbersAttr)` in `spin.hpp`). +/// @brief Strips the Kramers QN bits from a QuantumNumbersAttr. +/// @param t the quantum-numbers attribute to clear +/// @return @p t with all Kramers bits unset +/// @sa spinannotation_remove(const QuantumNumbersAttr&) in spin.hpp inline QuantumNumbersAttr kramers_annotation_remove( const QuantumNumbersAttr& t) { static_assert((~(~mask_v & ~bitset::reserved) & ~bitset::reserved) == mask_v, "Kramers bitmask uses reserved bits"); - return t.intersection(QuantumNumbersAttr(~mask_v & ~bitset::reserved)); + return t.intersection( + QuantumNumbersAttr(~mask_v & ~bitset::reserved)); } -/// Strip a trailing ⇑/⇓ Kramers annotation from a label. +/// @brief Strips a trailing ⇑/⇓ Kramers annotation from a label. +/// @tparam WS a wide-string or wide-string-view type +/// @param label the (possibly Kramers-annotated) label +/// @return @p label with any trailing ⇑/⇓ glyph removed template >>> + meta::is_wstring_or_view_v>>> std::wstring kramers_annotation_remove(WS&& label) { auto view = to_basic_string_view(label); - const auto has_annotation = !view.empty() && - (view.back() == L'⇑' || view.back() == L'⇓'); + const auto has_annotation = + !view.empty() && (view.back() == L'⇑' || view.back() == L'⇓'); return std::wstring{view.data(), view.data() + view.size() - (has_annotation ? 1 : 0)}; } -/// Add a Kramers annotation (⇑ or ⇓) to a label (must not already -/// carry one). Kramers::any returns the label unchanged. +/// @brief Adds a Kramers annotation (⇑ or ⇓) to a label. +/// @tparam WS a wide-string or wide-string-view type +/// @param label the label to annotate +/// @param k the Kramers state; `Kramers::any` returns @p label unchanged +/// @return @p label decorated with the matching ⇑/⇓ glyph +/// @pre @p label does not already carry a ⇑/⇓ annotation template >>> + meta::is_wstring_or_view_v>>> std::wstring kramers_annotation_add(WS&& label, Kramers k) { [[maybe_unused]] auto view = to_basic_string_view(label); - SEQUANT_ASSERT(view.empty() || - (view.back() != L'⇑' && view.back() != L'⇓')); + SEQUANT_ASSERT(view.empty() || (view.back() != L'⇑' && view.back() != L'⇓')); switch (k) { case Kramers::any: return std::wstring(std::forward(label)); @@ -96,27 +106,43 @@ std::wstring kramers_annotation_add(WS&& label, Kramers k) { SEQUANT_UNREACHABLE; } -/// Replace any existing ⇑/⇓ annotation with @p k (or strip if Kramers::any). +/// @brief Replaces any existing ⇑/⇓ annotation on a label with @p k. +/// @tparam WS a wide-string or wide-string-view type +/// @param label the (possibly Kramers-annotated) label +/// @param k the Kramers state to set; `Kramers::any` strips the annotation +/// @return @p label re-annotated with @p k template >>> + meta::is_wstring_or_view_v>>> std::wstring kramers_annotation_replace(WS&& label, Kramers k) { auto label_kf = kramers_annotation_remove(std::forward(label)); return kramers_annotation_add(label_kf, k); } -/// Construct a Kramers-up Index from @p idx (preserves type, ordinal, -/// proto-indices; sets the Kramers QN bits and decorates the label). +/// @brief Constructs a Kramers-up Index from @p idx. +/// @param idx the source index (its type, ordinal and proto-indices are +/// preserved) +/// @return a copy of @p idx with the Kramers-up QN bit set and the label +/// decorated with ⇑ [[nodiscard]] Index make_kramers_up(const Index& idx); -/// Construct a Kramers-down Index from @p idx. +/// @brief Constructs a Kramers-down Index from @p idx. +/// @param idx the source index +/// @return a copy of @p idx with the Kramers-down QN bit set and the +/// label decorated with ⇓ [[nodiscard]] Index make_kramers_dn(const Index& idx); -/// Construct a Kramers-free (Kramers::any) Index from @p idx. +/// @brief Constructs a Kramers-free (`Kramers::any`) Index from @p idx. +/// @param idx the source index +/// @return a copy of @p idx carrying the `Kramers::any` QN [[nodiscard]] Index make_kramers_free(const Index& idx); -/// Apply an Index→Index replacement map to all tensors in @p expr. -/// Mirrors `append_spin` in `spin.hpp`; provided as a separate symbol -/// so callers don't drag in the spin tracer. +/// @brief Applies an Index→Index replacement map to all tensors in an +/// expression. +/// @param expr the expression to rewrite +/// @param index_replacements the Index→Index substitution map +/// @return a copy of @p expr with every tensor index substituted +/// @sa append_spin in spin.hpp — this is the Kramers analogue, provided +/// as a separate symbol so callers need not pull in the spin tracer ExprPtr append_kramers(const ExprPtr& expr, const container::map& index_replacements); @@ -129,60 +155,404 @@ ExprPtr append_kramers(const ExprPtr& expr, // // Convention: complex conjugation of a tensor is encoded as a `*` // suffix on the tensor's label (e.g., `g` → `g*`, `t` → `t*`). -// Evaluator-side dispatch (e.g., MPQC's SQ_MP2 yielder) translates the -// suffix into a `.conj()` call on the underlying numeric tensor. +// Evaluator-side dispatch translates the suffix into a `.conj()` call +// on the underlying numeric tensor. // ===================================================================== -/// Test whether a tensor's label encodes complex conjugation -/// (trailing `*` suffix per the convention above). +/// @brief Tests whether a tensor label encodes complex conjugation. +/// @param label the tensor label to inspect +/// @return true iff @p label carries a trailing `*` suffix inline bool has_conj_suffix(std::wstring_view label) { return !label.empty() && label.back() == L'*'; } -/// Append `*` to a tensor's label (or strip if already present, since -/// (z*)* = z). Pure label rewrite — does NOT affect Symmetry tags or -/// IndexSpace QNs. Caller's responsibility to ensure the underlying -/// numeric tensor's complex conjugate is what's wanted. +/// @brief Toggles the `*` conjugation suffix on a tensor label. +/// @param label the tensor label +/// @return @p label with a `*` appended, or stripped if already present +/// (since `(z*)* = z`) +/// @note pure label rewrite — does not touch Symmetry tags or IndexSpace +/// QNs; the caller must ensure the underlying numeric tensor's +/// complex conjugate is what is wanted [[nodiscard]] std::wstring toggle_conj_suffix(std::wstring_view label); -/// Symbolic complex conjugation of an expression. +/// @brief Symbolic complex conjugation of an expression. /// /// Recursively walks @p expr and applies: -/// - Tensor: toggle its `*` suffix (z → z*, z* → z). -/// - Constant: numeric complex conjugate. -/// - Variable: not yet supported (throws). -/// - Product: conjugate the scalar AND each factor. (Π zᵢ)* = Π zᵢ*. -/// - Sum: conjugate each summand. (Σ zᵢ)* = Σ zᵢ*. +/// - Tensor: toggles its `*` suffix (`z → z*`, `z* → z`); +/// - Constant: numeric complex conjugate; +/// - Variable: not yet supported (throws); +/// - Product: conjugates the scalar and each factor, `(Π zᵢ)* = Π zᵢ*`; +/// - Sum: conjugates each summand, `(Σ zᵢ)* = Σ zᵢ*`. /// -/// `conjugate(conjugate(expr))` recovers @p expr. +/// @param expr the expression to conjugate +/// @return the symbolic complex conjugate of @p expr +/// @note `conjugate(conjugate(expr))` recovers @p expr [[nodiscard]] ExprPtr conjugate(const ExprPtr& expr); -/// Symbolic real part: returns `(expr + conjugate(expr)) / 2`. -/// Currently materialized as that two-term Sum (no native `Re` node); -/// downstream callers can apply numerical optimizations after evaluation. -[[nodiscard]] ExprPtr real_part(const ExprPtr& expr); +/// @brief Symbolic real-part wrapper: `RealPart(E)` represents `Re(E)` +/// for a scalar-valued Expr `E`. +/// @note real-valued itself and idempotent under conjugation; the +/// evaluator extracts @c inner, evaluates it, and takes the real +/// part of the resulting scalar +class RealPart : public sequant::Expr { + public: + RealPart() = delete; + RealPart(const RealPart&) = default; + RealPart(RealPart&&) = default; + ~RealPart() override = default; -/// Symbolic imaginary part: returns `(expr - conjugate(expr)) / (2i)`. -[[nodiscard]] ExprPtr imaginary_part(const ExprPtr& expr); + /// @brief Wraps a scalar-valued expression in a real-part marker. + /// @param inner the scalar-valued expression + /// @pre @p inner is non-null and `inner->is_scalar()` + explicit RealPart(ExprPtr inner) : inner_{std::move(inner)} { + SEQUANT_ASSERT(inner_); + SEQUANT_ASSERT(inner_->is_scalar()); + } + + /// @return the wrapped expression + const ExprPtr& inner() const { return inner_; } + + bool is_scalar() const override { return true; } + + type_id_type type_id() const override { return get_type_id(); } + + ExprPtr clone() const override { return ex(inner_->clone()); } + + /// @brief Adjoint of `Re(E)` — a no-op, since `Re(E)` is real. + void adjoint() override {} + + std::wstring to_latex() const override { + return L"\\Re\\left[" + inner_->to_latex() + L"\\right]"; + } + + private: + ExprPtr inner_; + + hash_type memoizing_hash() const override { + auto compute = [this]() { + auto v = hash::value(*inner_); + hash::combine(v, std::size_t{0xC0FFEE01ull}); // RealPart tag + return v; + }; + if (!hash_value_) hash_value_ = compute(); + return *hash_value_; + } + + bool static_equal(const Expr& that) const override { + return *inner_ == *static_cast(that).inner_; + } + + bool static_less_than(const Expr& that) const override { + return *inner_ < *static_cast(that).inner_; + } +}; + +/// @brief Symbolic imaginary-part wrapper: `ImagPart(E)` represents +/// `Im(E)` for a scalar-valued Expr `E`. +/// @note real-valued and idempotent under conjugation +class ImagPart : public sequant::Expr { + public: + ImagPart() = delete; + ImagPart(const ImagPart&) = default; + ImagPart(ImagPart&&) = default; + ~ImagPart() override = default; + + /// @brief Wraps a scalar-valued expression in an imaginary-part marker. + /// @param inner the scalar-valued expression + /// @pre @p inner is non-null and `inner->is_scalar()` + explicit ImagPart(ExprPtr inner) : inner_{std::move(inner)} { + SEQUANT_ASSERT(inner_); + SEQUANT_ASSERT(inner_->is_scalar()); + } + + /// @return the wrapped expression + const ExprPtr& inner() const { return inner_; } + + bool is_scalar() const override { return true; } + + type_id_type type_id() const override { return get_type_id(); } + + ExprPtr clone() const override { return ex(inner_->clone()); } + + /// @brief Adjoint of `Im(E)` — a no-op, since `Im(E)` is real. + void adjoint() override {} + + std::wstring to_latex() const override { + return L"\\Im\\left[" + inner_->to_latex() + L"\\right]"; + } + + private: + ExprPtr inner_; + + hash_type memoizing_hash() const override { + auto compute = [this]() { + auto v = hash::value(*inner_); + hash::combine(v, std::size_t{0xC0FFEE02ull}); // ImagPart tag + return v; + }; + if (!hash_value_) hash_value_ = compute(); + return *hash_value_; + } + + bool static_equal(const Expr& that) const override { + return *inner_ == *static_cast(that).inner_; + } + + bool static_less_than(const Expr& that) const override { + return *inner_ < *static_cast(that).inner_; + } +}; + +/// @brief Constructs `Re(expr)` as a RealPart wrapper. +/// @param expr the scalar-valued expression to wrap +/// @return a RealPart wrapping @p expr +[[nodiscard]] inline ExprPtr re(ExprPtr expr) { + return ex(std::move(expr)); +} + +/// @brief Constructs `Im(expr)` as an ImagPart wrapper. +/// @param expr the scalar-valued expression to wrap +/// @return an ImagPart wrapping @p expr +[[nodiscard]] inline ExprPtr im(ExprPtr expr) { + return ex(std::move(expr)); +} + +/// @brief Symbolic real part of an expression. +/// @param expr the scalar-valued expression +/// @return `Re(expr)` as a RealPart wrapper Expr +/// @note replaces the older `(expr + conj(expr))/2` Sum form so callers +/// can detect RealPart structurally and dispatch real-part +/// evaluation +[[nodiscard]] inline ExprPtr real_part(const ExprPtr& expr) { + return ex(expr); +} + +/// @brief Symbolic imaginary part of an expression. +/// @param expr the scalar-valued expression +/// @return `Im(expr)` as an ImagPart wrapper Expr +[[nodiscard]] inline ExprPtr imaginary_part(const ExprPtr& expr) { + return ex(expr); +} // ===================================================================== // Kramers tracer. // ===================================================================== -/// Rewrite an all-spinor tensor expression into a sum over canonical -/// Kramers-block representatives. Implements the closed-shell relativistic -/// MP2/CC algebra: enumerate per-tensor Kramers blocks, apply per-tensor -/// TRS canonicalization (with (-1)^k sign), and fold via the standard -/// canonicalize+rapid_simplify pipeline. Single entry handles MP2 energy, -/// CC residuals, and arbitrary operator expressions; expands any -/// antisymmetrizer (@c A) tensors internally before tracing. -/// -/// @param expr expression with indices in `IndexSpace`s carrying -/// `Kramers::any` (i.e., not yet specialized to up/down) -/// @return a Sum with canonical-orbit-block tensors; -/// tensor labels with conjugation get a `*` suffix +/// @brief Rewrites an all-spinor tensor expression into a sum over +/// canonical Kramers-block representatives. +/// +/// Implements the closed-shell relativistic trace algebra: enumerate the +/// per-tensor Kramers blocks, apply per-tensor TRS canonicalization (with +/// the `(-1)^k` sign), and fold via the standard +/// `canonicalize` + `rapid_simplify` pipeline. A single entry point +/// handles energy expressions, amplitude residuals, and arbitrary +/// operator expressions; antisymmetric tensors are expanded internally +/// before tracing. +/// +/// @param expr expression whose indices carry `Kramers::any` (i.e. not +/// yet specialized to up/down) +/// @return a `Sum` of canonical-orbit-block tensors; tensor +/// labels carrying conjugation get a `*` suffix ExprPtr kramers_trace(const ExprPtr& expr); +/// @brief Flips every Kramers state (`up` ↔ `down`) in an expression. +/// +/// Walks every tensor in @p expr, collects each Kramers-tagged index, +/// builds an Index→Index replacement that swaps the Kramers QN bit (and +/// the ⇑/⇓ label glyph), then applies the swap. Indices carrying +/// `Kramers::any` are left untouched. +/// +/// @param expr the expression to flip +/// @return a copy of @p expr with every specialized Kramers index flipped +/// @sa fold_conj_pairs — uses this to find TRS-conjugate-pair Products +[[nodiscard]] ExprPtr flip_kramers(const ExprPtr& expr); + +/// @brief Folds TRS conjugate-pair Products in a Sum. +/// +/// For each summand `A`, builds `A_flipped` via flip_kramers() (the +/// full-bar TRS partner), canonicalizes it, and checks structural +/// equality against the remaining summands: +/// - `A == A_flipped` (self-conjugate after canonicalize): emits +/// `Re(A)`, since the term is necessarily real-valued; +/// - `A_flipped == B` for some unused `B`: emits `2·Re(A)` and marks +/// both `A` and `B` consumed (term count drops by one per fold); +/// - no match: keeps `A` unchanged. +/// +/// Imaginary-pair detection (`A - A_flipped → 2i·Im(A)`) is supported +/// when the partner `B` carries the opposite scalar prefactor. +/// +/// @param expr a `Sum` to fold (typically the output of kramers_trace()) +/// @return a new `Sum`, typically with fewer summands, whose scalar +/// value equals that of @p expr +[[nodiscard]] ExprPtr fold_conj_pairs(const ExprPtr& expr); + +// ===================================================================== +// Cycle-based Kramers enumeration. +// +// The closed-shell spin-trace algorithm in SeQuant decomposes a Product's +// contraction graph into cycles and exploits cycle structure to count +// independent traces. We adapt the same idea for Kramers: each cycle in +// the contraction graph carries its own Kramers labelling, and the +// number of *canonical* Kramers patterns per cycle (under cyclic rotation +// + reflection) is typically << 2^L for length-L cycles — a long cycle +// visiting d distinct indices twice each has 2^d raw labellings that +// collapse to a much smaller canonical-pattern set. +// ===================================================================== + +/// @brief One contraction cycle in a Product's index graph. +/// +/// The cycle alternates intra-tensor edges (bra[k] ↔ ket[k] of the same +/// tensor) and inter-tensor edges (the same Index on two tensors). +/// `nodes` walks the cycle in visiting order. +struct ContractionCycle { + /// @brief One step of the cycle walk. + struct Node { + std::size_t tensor_idx; ///< which factor of the Product + std::size_t braket_pos; ///< position in the tensor's const_braket() + Index idx; ///< index occupying that slot + }; + container::svector nodes; ///< the cycle walk, in visiting order +}; + +/// @brief Decomposes a Product's index-contraction structure into cycles. +/// @param product a fully-contracted scalar Product +/// @return the contraction cycles of @p product +/// @pre every Index in @p product appears in exactly two `(tensor, slot)` +/// positions; each tensor's bra[k]/ket[k] positions are paired as +/// intra-tensor edges +[[nodiscard]] container::svector kramers_cycles( + const sequant::Product& product); + +/// @brief Diagnostic dump of the cycle structure of an expression. +/// +/// Walks @p expr (a `Sum` of Products), decomposes each Product into +/// cycles, and prints each cycle's length, distinct indices visited, +/// per-cycle Kramers labelling, and the number of canonical Kramers +/// patterns under cyclic-rotation symmetry. +/// +/// @param expr the expression to inspect +/// @param os the wide-character stream to write to +void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os); + +/// @brief Canonical Kramers label of a cycle. +/// @param c the contraction cycle +/// @return the lex-smallest rotation of @p c 's raw label string (e.g. +/// `BUUB` and `UBBU` both canonicalize to `BBUU`) +[[nodiscard]] std::wstring cycle_canonical_label(const ContractionCycle& c); + +/// @brief Cycle-canonical signature of a Product. +/// @param product the Product to fingerprint +/// @return the sorted multiset of per-cycle canonical Kramers labels, +/// joined into one key string (e.g. `BBUU|UUUU`) +/// @note Products with the same signature are equivalent under per-cycle +/// cyclic-rotation symmetry +[[nodiscard]] std::wstring cycle_canonical_signature( + const sequant::Product& product); + +/// @brief Buckets Products by cycle-canonical signature and sums each +/// bucket into one representative. +/// +/// Within a bucket the scalar prefactors are summed and one +/// representative is emitted. Acts independently of fold_conj_pairs() +/// and on a raw `Sum`-of-Products without any Re/Im wrappers. +/// +/// @param expr a `Sum` to fold +/// @return a `Sum` with one summand per cycle-canonical signature +/// @warning NOT value-preserving in general — the cyclic-rotation +/// signature is not a complete contraction invariant, so two +/// genuinely different contractions can share it. Superseded by +/// kramers_trace_cycles(); kept only for diagnostic comparison. +[[nodiscard]] ExprPtr cycle_canonical_fold(const ExprPtr& expr); + +// ===================================================================== +// Standalone cycle-driven Kramers tracer. +// +// Unlike `kramers_trace` (which enumerates all 2^n per-index Kramers +// configurations then folds post-hoc), this tracer is *driven* by the +// contraction-graph cycle structure — analogous to how +// `closed_shell_spintrace` is driven by `count_cycles`. +// +// Algorithm (cycles decomposed on the ANTISYMMETRIZED form, before +// expand_antisymm): +// 1. Decompose the input Product's contraction graph into cycles. +// 2. Group cycles into structural-equivalence classes: cycles that map +// to one another under the antisymmetrizer's index permutations +// (structurally-identical cycles are one class — the antisymmetric +// bra/ket index permutations swap them). +// 3. Enumerate Kramers labelings as MULTISETS within each class (not +// ordered assignments) — this is where the reduction over naive +// per-index 2^n enumeration comes from: equivalent cycles with +// swapped labels are the same contraction. +// 4. For each multiset configuration: substitute Kramers labels, +// expand_antisymm, canonicalize; attach the multiset multiplicity. +// 5. Collect into a Sum; the caller may apply `fold_conj_pairs` for +// additional TRS conjugate-pair folding. +// +// Correct-by-construction: works from contraction topology, never guesses +// equivalences. Coexists with `kramers_trace` so the two can be +// cross-validated. +// ===================================================================== + +/// @brief Cycle-driven Kramers tracer. +/// +/// See the block comment above for the algorithm. +/// +/// @param expr a closed-shell all-spinor expression with `Kramers::any` +/// indices and (optionally) antisymmetric tensors +/// @return a `Sum` of canonical Kramers-block contractions whose +/// scalar value equals the full 2^n per-index trace +[[nodiscard]] ExprPtr kramers_trace_cycles(const ExprPtr& expr); + +/// @brief Recombines direct/exchange Product pairs into antisymmetric +/// tensors — the inverse of `expand_antisymm`. +/// +/// After `kramers_trace` + `fold_conj_pairs` the expression is a `Sum` of +/// fully-expanded NonSymm contractions. Many of these are the `direct` +/// and `exchange` halves of a single antisymmetrized contraction: two +/// summands recombine into one when +/// - one tensor is structurally identical between them, and +/// - the other tensor differs by a column permutation, with the +/// summand scalars in the ratio `sign(permutation)`. +/// Then `c·X·Y_direct + (sign·c)·X·Y_exchange → c·X·Ȳ`, where `Ȳ` is +/// `Y` re-tagged `Symmetry::Antisymm` over the Kramers-restricted index +/// space (NOT full spinor). +/// +/// Applied iteratively to a fixed point, so a second pass over the +/// once-recombined terms catches doubly-antisymmetrized blocks. Re/Im +/// wrappers from fold_conj_pairs() are preserved — recombination acts on +/// the wrapped inner Product. +/// +/// @param expr a `Sum` of (optionally Re/Im-wrapped) NonSymm contractions +/// @return a `Sum` in minimal mixed antisymm/non-antisymm form, with the +/// same scalar value as @p expr +/// @note value-preserving by construction — it only re-groups terms that +/// sum to the same scalar; the result is the analogue of what the +/// open-shell spin tracer yields, but driven by index wiring +/// rather than spin-label uniformity +[[nodiscard]] ExprPtr antisymm_recombine(const ExprPtr& expr); + +/// @brief Burnside-enumeration Kramers tracer. +/// +/// Like `kramers_trace`, but instead of walking all `2^n` per-index +/// Kramers configurations it enumerates only the orbit representatives +/// under the index-permutation symmetry group induced by the +/// antisymmetrizer structure of the input (the group generated by +/// transpositions within each tensor's antisymmetric bra/ket groups). +/// Each representative is emitted once, scaled by its orbit size; group +/// elements act as bit-permutations on the n-bit config and orbits are +/// found by bit-ops. +/// +/// A transposition `(p,q)` is admitted as a group generator iff in every +/// tensor containing both indices they sit in the same antisymmetric +/// bra/ket group and the accumulated permutation sign is `+1`. +/// +/// @param expr a closed-shell all-spinor expression with `Kramers::any` +/// indices +/// @return a `Sum` value-equivalent (after `canonicalize`) to +/// the output of kramers_trace(); feed it the same +/// `fold_conj_pairs → antisymm_recombine` pipeline +[[nodiscard]] ExprPtr kramers_trace_burnside(const ExprPtr& expr); + } // namespace sequant::mbpt #endif // SEQUANT_DOMAIN_MBPT_SPINOR_HPP From c1e12b8768d527f46de8a843ba2715e1ba22656b Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Thu, 14 May 2026 15:40:04 -0400 Subject: [PATCH 5/6] mbpt/spinor: quaternion_decompose + complex_split passes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit quaternion_decompose — per-tensor M-matrix expansion of Kramers tensors into the 4 real (s,x,y,z) quaternion components; prototyped as a fine-grained oracle, shown to over-decompose (closed-shell 1eX2C spin-traces to a complex scalar, not a quaternion). complex_split — the keeper: split each complex Kramers tensor into real g~r/g~i tensors and propagate complex arithmetic symbolically (Re(g·t) = g~r·t~r - g~i·t~i), giving a real-tensor form for real-arithmetic evaluation. Plus CC @todo notes on the Kramers tracer (free-index TRS signs, rank-(2,2) antisymm_recombine, rectangular-tensor cycle assumption). --- SeQuant/domain/mbpt/spinor.cpp | 269 ++++++++++++++++++++++++++++++++- SeQuant/domain/mbpt/spinor.hpp | 184 +++++++++++++++++++++- 2 files changed, 447 insertions(+), 6 deletions(-) diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp index 2ef1e053c1..310fdec308 100644 --- a/SeQuant/domain/mbpt/spinor.cpp +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -204,6 +204,11 @@ ExprPtr kramers_trace(const ExprPtr& expr) { // products in which every barred index appears in an even number of // tensors — true for any all-internal-index contraction). // + // TODO(CC): this all-internal-index assumption breaks for CC residuals, + // which carry free/external indices — those need explicit + // (-1)^(external bars) sign bookkeeping. See the kramers_trace @todo in + // spinor.hpp. + // // For each TRS pair we materialize only the "lex-smaller" configuration // X explicitly; the X̄ partner is emitted as `conjugate(T_X)`. Tensors // in the conj branch carry a `*` label suffix (per the convention in @@ -494,6 +499,10 @@ container::svector kramers_cycles(const Product& product) { // intra-tensor pairing per tensor: walk positions and pair k ↔ k+rank/2 // within each tensor's slots. We rebuild rank by counting consecutive // slots with the same tensor_idx. + // TODO(CC): the k ↔ k+rank/2 pairing assumes a rectangular tensor + // (bra_rank == ket_rank), which holds for g/t/residual tensors but not + // for every CC intermediate; non-rectangular tensors need the bra/ket + // split read from the tensor directly rather than from total rank. for (std::size_t i = 0; i < N;) { std::size_t j = i; while (j < N && slots[j].tensor_idx == slots[i].tensor_idx) ++j; @@ -981,8 +990,11 @@ bool parse_recomb_entry(const ExprPtr& summand, RecombEntry& e) { /// Detect whether @p t2 is @p t1 with a single column transposition. /// Returns 1 if t2 = bra-swap(t1), 2 if t2 = ket-swap(t1), 0 otherwise. -/// Rank-4 (2+2) only for now. Both indices compared by full Index identity -/// (recombinable terms share dummy names post-canonicalize). +/// Both indices compared by full Index identity (recombinable terms share +/// dummy names post-canonicalize). +/// TODO(CC): rank-(2,2) only. CCSDT triples residuals need this +/// generalized to rank-(n,n) — enumerate the column transpositions of an +/// n-column tensor and report which single one (if any) maps t1 -> t2. int detect_single_column_swap(const Tensor& t1, const Tensor& t2) { if (t1.label() != t2.label()) return 0; if (t1.bra_rank() != 2 || t1.ket_rank() != 2) return 0; @@ -1278,4 +1290,257 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { return r; } +// ===================================================================== +// Quaternion decomposer. See spinor.hpp for the design overview. +// ===================================================================== + +namespace { + +/// Reads the specialized Kramers state of an index. +Kramers index_kramers(const Index& idx) { + const auto k_bits = idx.space().qns().to_int32() & mask_v; + if (k_bits == static_cast(Kramers::up)) return Kramers::up; + if (k_bits == static_cast(Kramers::down)) return Kramers::down; + throw std::runtime_error( + "quaternion_decompose: tensor index is not Kramers-specialized"); +} + +/// Entry M(k_bra,k_ket)[c_bra][c_ket] of the quaternion decomposition +/// matrix. Component order s=0, x=1, y=2, z=3. Derived from the spinor +/// coefficient structure C_α = C_s + iC_x, C_β = -C_y + iC_z (barred +/// partner [-C_β*; C_α*]); see kramers-restriction §6 / Wiebke Eq. 8-11. +Constant::scalar_type quaternion_m(Kramers k_bra, Kramers k_ket, int c_bra, + int c_ket) { + using S = Constant::scalar_type; + auto v = [](int re, int im) { return S{Complex{re, im}}; }; + // (re,im) tables for M(up,up) and M(up,down); rows = c_bra, cols = c_ket. + static const int uu_re[4][4] = { + {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; + static const int uu_im[4][4] = { + {0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0}}; + static const int ud_re[4][4] = { + {0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 0}}; + static const int ud_im[4][4] = { + {0, 0, 0, 1}, {0, 0, -1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}; + const bool b_up = (k_bra == Kramers::up); + const bool k_up = (k_ket == Kramers::up); + if (b_up && k_up) return v(uu_re[c_bra][c_ket], uu_im[c_bra][c_ket]); + if (!b_up && !k_up) // M(down,down) = conj(M(up,up)) + return v(uu_re[c_bra][c_ket], -uu_im[c_bra][c_ket]); + if (b_up && !k_up) return v(ud_re[c_bra][c_ket], ud_im[c_bra][c_ket]); + // M(down,up) = M(up,down)^dagger + return v(ud_re[c_ket][c_bra], -ud_im[c_ket][c_bra]); +} + +constexpr char kQComp[4] = {'s', 'x', 'y', 'z'}; + +/// Decompose one rank-(2,2) Nonsymm tensor into its (up to 64) quaternion +/// components: T = Σ_c M(k_B0,k_K0)[c_B0,c_K0]·M(k_B1,k_K1)[c_B1,c_K1]·T#c. +ExprPtr decompose_nonsymm_tensor(const Tensor& t) { + if (t.bra_rank() != 2 || t.ket_rank() != 2) + throw std::runtime_error( + "quaternion_decompose: only rank-(2,2) tensors are supported"); + const bool conj = has_conj_suffix(t.label()); + const std::wstring core = + conj ? std::wstring{t.label().substr(0, t.label().size() - 1)} + : std::wstring{t.label()}; + const Index& B0 = t.bra().at(0); + const Index& B1 = t.bra().at(1); + const Index& K0 = t.ket().at(0); + const Index& K1 = t.ket().at(1); + // electron 1 = (B0, K0), electron 2 = (B1, K1) + const Kramers kB0 = index_kramers(B0), kB1 = index_kramers(B1); + const Kramers kK0 = index_kramers(K0), kK1 = index_kramers(K1); + + auto sum = std::make_shared(); + const Constant::scalar_type zero{0}; + for (int cB0 = 0; cB0 < 4; ++cB0) + for (int cK0 = 0; cK0 < 4; ++cK0) { + const auto m1 = quaternion_m(kB0, kK0, cB0, cK0); + if (m1 == zero) continue; + for (int cB1 = 0; cB1 < 4; ++cB1) + for (int cK1 = 0; cK1 < 4; ++cK1) { + const auto m2 = quaternion_m(kB1, kK1, cB1, cK1); + if (m2 == zero) continue; + auto coeff = m1 * m2; + if (conj) coeff = sequant::conj(coeff); + std::wstring lbl = quaternion_label_add( + core, kQComp[cB0], kQComp[cB1], kQComp[cK0], kQComp[cK1]); + if (conj) lbl += L'*'; + auto qt = + ex(lbl, bra{B0, B1}, ket{K0, K1}, Symmetry::Nonsymm, + BraKetSymmetry::Nonsymm, ColumnSymmetry::Nonsymm); + sum->append(ex(coeff) * qt); + } + } + return sum; +} + +} // namespace + +ExprPtr quaternion_decompose(const ExprPtr& expr) { + if (expr->is()) { + auto s = std::make_shared(); + for (auto&& summand : *expr) s->append(quaternion_decompose(summand)); + return s; + } + if (expr->is()) + return ex(quaternion_decompose(expr->as().inner())); + if (expr->is()) + return ex(quaternion_decompose(expr->as().inner())); + if (expr->is() || expr->is()) return expr; + if (expr->is()) { + auto const& t = expr->as(); + if (t.symmetry() == Symmetry::Antisymm) { + // expand g - g_exchange first, then decompose each Nonsymm piece + return quaternion_decompose( + expand_antisymm(expr, /*skip_spinsymm*/ false)); + } + return decompose_nonsymm_tensor(t); + } + if (expr->is()) { + auto const& p = expr->as(); + bool has_wrapper = false; + for (auto&& f : p) + if (f->is() || f->is()) has_wrapper = true; + auto prod = std::make_shared(); + prod->scale(p.scalar()); + for (auto&& f : p) + prod->append(1, quaternion_decompose(f), Product::Flatten::No); + ExprPtr r = prod; + if (!has_wrapper) { + // a genuine contraction product — distribute the per-tensor Sums and + // fold dummy-renamed duplicates + expand(r); + rapid_simplify(r); + canonicalize(r); + rapid_simplify(r); + } + return r; + } + throw std::runtime_error("quaternion_decompose: unsupported Expr type"); +} + +// ===================================================================== +// Complex (Re/Im) split. See spinor.hpp for the design overview. +// ===================================================================== + +namespace { + +/// Returns (re, im) such that @p expr == re + i*im, with both `re` and +/// `im` expressed over real-valued tensors (labels carrying `~r`/`~i`). +std::pair split_re_im(const ExprPtr& expr) { + using S = Constant::scalar_type; + if (expr->is()) { + const auto v = expr->as().value(); + return {ex(S{v.real()}), ex(S{v.imag()})}; + } + if (expr->is()) { + auto inner = split_re_im(expr->as().inner()); + return {inner.first, ex(S{0})}; + } + if (expr->is()) { + auto inner = split_re_im(expr->as().inner()); + return {inner.second, ex(S{0})}; + } + if (expr->is()) { + auto const& t = expr->as(); + const bool conj = has_conj_suffix(t.label()); + const std::wstring core = + conj ? std::wstring{t.label().substr(0, t.label().size() - 1)} + : std::wstring{t.label()}; + auto mk = [&](bool imag) { + return ex( + complex_label_add(core, imag), + bra(container::svector{t.bra().begin(), t.bra().end()}), + ket(container::svector{t.ket().begin(), t.ket().end()}), + aux(container::svector{t.aux().begin(), t.aux().end()}), + t.symmetry(), t.braket_symmetry(), t.column_symmetry()); + }; + ExprPtr re = mk(false); + ExprPtr im = mk(true); + if (conj) im = ex(S{-1}) * im; // g* = g~r - i·g~i + return {re, im}; + } + if (expr->is()) { + auto re = std::make_shared(); + auto im = std::make_shared(); + for (auto&& s : *expr) { + auto part = split_re_im(s); + re->append(part.first); + im->append(part.second); + } + return {re, im}; + } + if (expr->is()) { + auto const& p = expr->as(); + // Π_i z_i built up incrementally, starting from the leading scalar. + auto acc = split_re_im(ex(p.scalar())); + for (auto&& f : p) { + auto fac = split_re_im(f); + // (ar + i·ai)(fr + i·fi) = (ar·fr - ai·fi) + i(ar·fi + ai·fr) + ExprPtr nr = acc.first * fac.first + + ex(S{-1}) * (acc.second * fac.second); + ExprPtr ni = acc.first * fac.second + acc.second * fac.first; + acc = {nr, ni}; + } + return acc; + } + if (expr->is()) + throw std::runtime_error("complex_split: Variable not supported"); + throw std::runtime_error("complex_split: unsupported Expr type"); +} + +/// Distribute and constant-fold a real-tensor scalar expression. +/// +/// NB: deliberately does NOT `canonicalize`. The split tensors inherit +/// the (already-consistent) index wiring of the Kramers-traced input; +/// `canonicalize` would permute the indices of `Antisymm` `~r`/`~i` +/// tensors and emit `(-1)` signs that the evaluator cannot see — it +/// rebuilds the `[as]` integral formula from index *type* + Kramers +/// label, which is invariant under such a permutation, so the sign +/// would be uncompensated. `expand` + `rapid_simplify` only distribute +/// and constant-fold; they never permute tensor indices. +ExprPtr clean_real(ExprPtr e) { + expand(e); + rapid_simplify(e); + return e; +} + +} // namespace + +ExprPtr complex_split(const ExprPtr& expr) { + if (expr->is()) { + auto s = std::make_shared(); + for (auto&& summand : *expr) s->append(complex_split(summand)); + return s; + } + if (expr->is()) { + return ex( + clean_real(split_re_im(expr->as().inner()).first)); + } + if (expr->is()) { + // Im(inner) is itself real-valued; its real-tensor form is the + // `.second` of the split — wrap in RealPart so the evaluator + // boundary treats it as a real quantity. + return ex( + clean_real(split_re_im(expr->as().inner()).second)); + } + if (expr->is() || expr->is()) return expr; + if (expr->is()) { + auto const& p = expr->as(); + if (p.size() == 1 && + (p.factor(0)->is() || p.factor(0)->is())) { + auto prod = std::make_shared(); + prod->scale(p.scalar()); + prod->append(1, complex_split(p.factor(0)), Product::Flatten::No); + return prod; + } + throw std::runtime_error( + "complex_split: expected `Constant * (Re|Im)(...)` summand — run " + "fold_conj_pairs first so every term is real-wrapped"); + } + throw std::runtime_error("complex_split: unsupported top-level Expr type"); +} + } // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/spinor.hpp b/SeQuant/domain/mbpt/spinor.hpp index 9d90985125..9069d1bfa8 100644 --- a/SeQuant/domain/mbpt/spinor.hpp +++ b/SeQuant/domain/mbpt/spinor.hpp @@ -30,6 +30,8 @@ #include #include +#include +#include #include #include #include @@ -344,15 +346,22 @@ class ImagPart : public sequant::Expr { /// Implements the closed-shell relativistic trace algebra: enumerate the /// per-tensor Kramers blocks, apply per-tensor TRS canonicalization (with /// the `(-1)^k` sign), and fold via the standard -/// `canonicalize` + `rapid_simplify` pipeline. A single entry point -/// handles energy expressions, amplitude residuals, and arbitrary -/// operator expressions; antisymmetric tensors are expanded internally -/// before tracing. +/// `canonicalize` + `rapid_simplify` pipeline. Antisymmetric tensors are +/// expanded internally before tracing. /// /// @param expr expression whose indices carry `Kramers::any` (i.e. not /// yet specialized to up/down) /// @return a `Sum` of canonical-orbit-block tensors; tensor /// labels carrying conjugation get a `*` suffix +/// @todo CC support. The per-tensor `(-1)^k` TRS signs only cancel when +/// every barred index appears an even number of times across the +/// whole expression — true for a fully-contracted scalar (the MP2 +/// energy), but NOT for an expression carrying free/external +/// indices (a CC residual `R_{IJ...}^{AB...}`). Residuals need +/// explicit `(-1)^(bars on external indices)` sign bookkeeping; +/// until that is added this entry point is correct only for +/// fully-contracted scalar expressions. The same caveat applies to +/// `kramers_trace_cycles` and `kramers_trace_burnside`. ExprPtr kramers_trace(const ExprPtr& expr); /// @brief Flips every Kramers state (`up` ↔ `down`) in an expression. @@ -529,6 +538,12 @@ void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os); /// sum to the same scalar; the result is the analogue of what the /// open-shell spin tracer yields, but driven by index wiring /// rather than spin-label uniformity +/// @todo CC support. Recombination currently recognizes only rank-(2,2) +/// column swaps (see `detect_single_column_swap` in spinor.cpp). +/// CCSDT triples residuals (`Ȳ_{IJK}^{ABC}` etc.) need the +/// single-column-swap detection generalized to rank-(n,n). +/// Higher-rank tensors are left expanded — not folded, never +/// miscombined — so this is a missed-optimization, not a bug. [[nodiscard]] ExprPtr antisymm_recombine(const ExprPtr& expr); /// @brief Burnside-enumeration Kramers tracer. @@ -553,6 +568,167 @@ void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os); /// `fold_conj_pairs → antisymm_recombine` pipeline [[nodiscard]] ExprPtr kramers_trace_burnside(const ExprPtr& expr); +// ===================================================================== +// Quaternion decomposer. +// +// After Kramers tracing, every 2-electron tensor is rewritten into the +// four real (s, x, y, z) quaternion components of its spinor +// coefficients (Helmich-Paris, Repisky, Visscher 2016; see +// kramers-restriction §6 / Wiebke Eq. 8-11). For a rank-(2,2) tensor +// `T(B0,B1;K0,K1)` with Kramers-specialized indices, the spin-traced +// per-electron transition density decomposes via a 4x4 complex matrix +// `M(k_bra,k_ket)` — one per Kramers-pair, 8 nonzero entries each: +// +// T(B0,B1;K0,K1) = Σ_c M(k_B0,k_K0)[c_B0,c_K0]·M(k_B1,k_K1)[c_B1,c_K1] +// · T#(B0,B1;K0,K1) +// +// where electron 1 = (B0,K0), electron 2 = (B1,K1). The quaternion +// component tuple is encoded as a `#`-prefixed 4-char suffix on the +// tensor label (analogous to the `*` conjugation marker); the contracted +// indices are left unchanged so tensors still contract pairwise. +// ===================================================================== + +/// @brief Encodes a quaternion-component tuple onto a tensor label. +/// @param core the (component-free) tensor label +/// @param c0,c1,c2,c3 the components of bra[0],bra[1],ket[0],ket[1], +/// each one of 's','x','y','z' +/// @return @p core with a `#c0c1c2c3` suffix appended +[[nodiscard]] inline std::wstring quaternion_label_add(std::wstring_view core, + char c0, char c1, + char c2, char c3) { + std::wstring s{core}; + s += L'#'; + s += static_cast(c0); + s += static_cast(c1); + s += static_cast(c2); + s += static_cast(c3); + return s; +} + +/// @brief Parses the quaternion-component suffix off a tensor label. +/// @param label the (possibly quaternion-decomposed) tensor label; a +/// trailing `*` conjugation marker is tolerated +/// @return the four component chars (bra[0],bra[1],ket[0],ket[1]), or +/// std::nullopt if @p label carries no `#`-suffix +[[nodiscard]] inline std::optional> quaternion_label_parse( + std::wstring_view label) { + const auto pos = label.find(L'#'); + if (pos == std::wstring_view::npos) return std::nullopt; + auto tail = label.substr(pos + 1); + if (!tail.empty() && tail.back() == L'*') + tail = tail.substr(0, tail.size() - 1); + if (tail.size() != 4) return std::nullopt; + std::array out{}; + for (std::size_t i = 0; i < 4; ++i) { + const wchar_t ch = tail[i]; + if (ch != L's' && ch != L'x' && ch != L'y' && ch != L'z') + return std::nullopt; + out[i] = static_cast(ch); + } + return out; +} + +/// @brief Returns the component-free core of a (possibly quaternion- +/// decomposed, possibly conjugated) tensor label. +/// @param label the tensor label +/// @return @p label truncated at the `#` suffix marker (no-op if absent) +[[nodiscard]] inline std::wstring quaternion_label_core( + std::wstring_view label) { + return std::wstring{label.substr(0, label.find(L'#'))}; +} + +/// @brief Rewrites a Kramers-traced expression into quaternion components. +/// +/// Recurses through Sum / Product / RealPart / ImagPart, expanding any +/// `Symmetry::Antisymm` tensor first, then applies the per-tensor 4x4 +/// `M`-matrix decomposition (see the section comment above) to every +/// rank-(2,2) `Symmetry::Nonsymm` tensor. Each input tensor expands into +/// up to 64 quaternion-component tensors (8 nonzero `M` entries per +/// electron pair); the decomposed sub-products are canonicalized so +/// dummy-renamed duplicates fold. +/// +/// @param expr a Kramers-traced expression (output of `kramers_trace` → +/// `fold_conj_pairs` → `antisymm_recombine`), every tensor +/// index carrying a specialized `Kramers::up`/`Kramers::down` +/// QN +/// @return the quaternion-component-decomposed expression, value-equal to +/// @p expr; tensor labels carry a `#c0c1c2c3` component suffix +/// @pre every tensor in @p expr is rank-(2,2) with Kramers-specialized +/// indices +[[nodiscard]] ExprPtr quaternion_decompose(const ExprPtr& expr); + +// ===================================================================== +// Complex (Re/Im) split. +// +// For a closed-shell 1eX2C energy expression the spin-free Coulomb +// operator spin-traces each per-electron transition density down to its +// scalar Pauli component, which is *complex* — no quaternion +// (4-component) structure survives, so `quaternion_decompose` +// over-decomposes. The minimal real-arithmetic form is the 2-component +// split: every complex tensor `g` becomes `g = g~r + i·g~i` with `g~r`, +// `g~i` real tensors, and complex arithmetic is propagated symbolically +// (`Re(g·t) = g~r·t~r − g~i·t~i`, etc.). The component is encoded as a +// `~r`/`~i` label suffix (analogous to the `*` conjugation marker). +// ===================================================================== + +/// @brief Encodes a real/imaginary-part marker onto a tensor label. +/// @param core the (marker-free) tensor label +/// @param imag true for the imaginary part (`~i`), false for the real +/// part (`~r`) +/// @return @p core with the `~r`/`~i` suffix appended +[[nodiscard]] inline std::wstring complex_label_add(std::wstring_view core, + bool imag) { + return std::wstring{core} + (imag ? L"~i" : L"~r"); +} + +/// @brief Parses the real/imaginary-part marker off a tensor label. +/// @param label the (possibly complex-split) tensor label; a trailing +/// `*` conjugation marker is tolerated +/// @return true for the imaginary part, false for the real part, or +/// std::nullopt if @p label carries no `~` marker +[[nodiscard]] inline std::optional complex_label_parse( + std::wstring_view label) { + const auto pos = label.find(L'~'); + if (pos == std::wstring_view::npos) return std::nullopt; + auto tail = label.substr(pos + 1); + if (!tail.empty() && tail.back() == L'*') + tail = tail.substr(0, tail.size() - 1); + if (tail == std::wstring_view{L"r"}) return false; + if (tail == std::wstring_view{L"i"}) return true; + return std::nullopt; +} + +/// @brief Returns the marker-free core of a (possibly complex-split) +/// tensor label. +/// @param label the tensor label +/// @return @p label truncated at the `~` marker (no-op if absent) +[[nodiscard]] inline std::wstring complex_label_core(std::wstring_view label) { + return std::wstring{label.substr(0, label.find(L'~'))}; +} + +/// @brief Rewrites a Kramers-traced expression into one over real +/// tensors via the complex (Re/Im) split. +/// +/// Each complex tensor `g` is replaced by two real tensors `g~r`, `g~i` +/// (`g = g~r + i·g~i`), and complex arithmetic is propagated +/// symbolically through every Product and Sum. A `RealPart`/`ImagPart` +/// wrapper collapses to the corresponding real sub-expression +/// (`Re(g·t) → g~r·t~r − g~i·t~i`). Conjugated tensors (`g*`) split as +/// `g~r − i·g~i`. The resulting real sub-expressions are distributed and +/// constant-folded but deliberately NOT canonicalized — the split +/// tensors inherit the input's index wiring, and canonicalizing +/// `Antisymm` `~r`/`~i` tensors would emit signs the evaluator cannot +/// see (it rebuilds `[as]` formulas from index type + Kramers label). +/// +/// @param expr a Kramers-traced expression whose summands are +/// `Constant * RealPart/ImagPart(...)` (the output of +/// `kramers_trace → fold_conj_pairs → antisymm_recombine`) +/// @return an equivalent expression in which every tensor is real; +/// tensor labels carry a `~r`/`~i` suffix, real sub-expressions +/// re-wrapped in `RealPart` for the evaluator boundary +/// @pre every summand is real-wrapped — run `fold_conj_pairs` first +[[nodiscard]] ExprPtr complex_split(const ExprPtr& expr); + } // namespace sequant::mbpt #endif // SEQUANT_DOMAIN_MBPT_SPINOR_HPP From 915d40eef7f25799d45bb4fe04ca45b1abec4795 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 20 May 2026 11:03:53 -0400 Subject: [PATCH 6/6] mbpt/spinor: reuse core SeQuant passes; drop experimental helpers Rework the Kramers-restriction / spinor module to lean on existing SeQuant machinery instead of hand-rolled helpers. Add a ResultExpr-based kramers_trace overload and specialize_and_simplify; clean_real now defers to non_canon_simplify. Remove the rejected quaternion (4-component) decomposition and the dead cycle / normalize / conjugate / scalar_eq helpers. Value-preserving. --- SeQuant/domain/mbpt/spin.hpp | 13 + SeQuant/domain/mbpt/spinor.cpp | 1467 +++++++++++++++----------------- SeQuant/domain/mbpt/spinor.hpp | 654 ++++---------- 3 files changed, 864 insertions(+), 1270 deletions(-) diff --git a/SeQuant/domain/mbpt/spin.hpp b/SeQuant/domain/mbpt/spin.hpp index 0194ef664b..ac03a488d9 100644 --- a/SeQuant/domain/mbpt/spin.hpp +++ b/SeQuant/domain/mbpt/spin.hpp @@ -110,6 +110,19 @@ std::wstring spinannotation_replacе(WS&& label, Spin s) { // make null-spin idx [[nodiscard]] Index make_spinfree(const Index& idx); +namespace detail { + +/// @brief Resets index tags on every Tensor in @p expr. +/// +/// `canonicalize` tags Index instances and does not always reset them on +/// exit; calling it twice without an intervening reset (or on a Sum that +/// aggregates per-term canonicalized Products) trips +/// `SEQUANT_ASSERT(!tag_.has_value())` in `Taggable::assign`. The standard +/// SeQuant idiom is `detail::reset_idx_tags(e); simplify(e);`. +void reset_idx_tags(const ExprPtr& expr); + +} // namespace detail + /// @brief Preserving particle symmetry, swaps bra and ket labels on all tensors /// in an expression /// @param expr ExprPtr to transform diff --git a/SeQuant/domain/mbpt/spinor.cpp b/SeQuant/domain/mbpt/spinor.cpp index 310fdec308..5d8e897e15 100644 --- a/SeQuant/domain/mbpt/spinor.cpp +++ b/SeQuant/domain/mbpt/spinor.cpp @@ -13,6 +13,10 @@ #include #include #include +#include +#include +#include +#include #include #include @@ -21,14 +25,11 @@ namespace sequant::mbpt { -namespace detail { +namespace { /// Construct a new Index in the same IndexSpace type as @p idx, with /// the @p k Kramers QN bit set (replacing any existing Kramers state) -/// and the label decorated with the matching ⇑/⇓ glyph. Spin bits -/// (alpha/beta) are NOT touched here — Spin and Kramers are mutually -/// exclusive on a single index, but the per-index check is left to -/// the caller (a typical workflow only ever sets one of them). +/// and the label decorated with the matching ⇑/⇓ glyph. /// /// Mirrors `detail::make_index_with_spincase` in `spin.cpp` so that /// SeQuant's IndexSpace registry is reused: if a space with the @@ -38,9 +39,7 @@ Index make_index_with_kramers(const Index& idx, mbpt::Kramers k) { SEQUANT_ASSERT(!(idx.label().find(L'⇑') != std::wstring::npos && idx.label().find(L'⇓') != std::wstring::npos)); - // Strip any existing Kramers bits, then OR in the requested state. auto qns = mbpt::kramers_annotation_remove(idx.space().qns()).unIon(k); - IndexSpace space; const auto label = mbpt::kramers_annotation_replace(idx.space().base_key(), k); @@ -53,7 +52,6 @@ Index make_index_with_kramers(const Index& idx, mbpt::Kramers k) { } if (!space) { space = IndexSpace{label, idx.space().type(), qns, - // assume size does not depend on Kramers state idx.space().approximate_size()}; } auto protoindices = idx.proto_indices(); @@ -61,167 +59,114 @@ Index make_index_with_kramers(const Index& idx, mbpt::Kramers k) { return Index{space, idx.ordinal(), protoindices}; } -} // namespace detail - -Index make_kramers_up(const Index& idx) { - return detail::make_index_with_kramers(idx, mbpt::Kramers::up); +/// Standard SeQuant post-trace pipeline: reset index tags then +/// `simplify` (= expand + rapid_simplify + canonicalize + +/// rapid_simplify, see `expr_algorithms.cpp`). Mirrors the spin tracer +/// idiom at `spin.cpp:953-956`. +ExprPtr& simplify_with_reset(ExprPtr& e) { + detail::reset_idx_tags(e); + return simplify(e); } -Index make_kramers_dn(const Index& idx) { - return detail::make_index_with_kramers(idx, mbpt::Kramers::down); +/// Read the specialized Kramers state (`up`/`down`/`any`) of an index. +Kramers index_kramers_state(const Index& idx) { + const auto k_bits = idx.space().qns().to_int32() & mask_v; + if (k_bits == static_cast(Kramers::up)) return Kramers::up; + if (k_bits == static_cast(Kramers::down)) return Kramers::down; + return Kramers::any; } -Index make_kramers_free(const Index& idx) { - return detail::make_index_with_kramers(idx, mbpt::Kramers::any); +/// Return a clone of @p t with `set_label(label)` applied. Used by the +/// label-rewriting passes (conjugation toggle, complex split, +/// antisymm-recombine) to avoid the noisier full-Tensor reconstruction. +ExprPtr tensor_with_label(const Tensor& t, std::wstring label) { + auto out = std::static_pointer_cast(t.clone()); + out->set_label(std::move(label)); + return out; } -// ===================================================================== -// Conjugation infrastructure (groundwork for Kramers tracer + quaternion -// decomposer). -// ===================================================================== - -std::wstring toggle_conj_suffix(std::wstring_view label) { - if (has_conj_suffix(label)) - return std::wstring{label.substr(0, label.size() - 1)}; - return std::wstring{label} + L'*'; +/// Walk @p expr and collect the unique Kramers-eligible (`Kramers::any`) +/// indices in encounter order. Indices that already carry a specialized +/// Kramers state are skipped. +container::svector collect_kramers_indices(const ExprPtr& expr) { + container::svector seen; + container::set seen_set; + expr->visit( + [&](const ExprPtr& e) { + if (!e->is()) return; + for (const auto& idx : e->as().const_braket()) { + if (index_kramers_state(idx) != Kramers::any) continue; + if (seen_set.insert(idx).second) seen.push_back(idx); + } + }, + /* recursive = */ true); + return seen; } -namespace { +/// Build the Index → Kramers-specialized-Index replacement for one +/// configuration (bit `k` of @p cfg → Kramers::down, else Kramers::up). +container::map kramers_replacement_map( + const container::svector& indices, std::uint64_t cfg) { + container::map repl; + for (std::size_t k = 0; k < indices.size(); ++k) { + const bool is_dn = (cfg >> k) & 1u; + repl[indices[k]] = + is_dn ? make_kramers_dn(indices[k]) : make_kramers_up(indices[k]); + } + return repl; +} -/// Apply conjugation to a single Tensor: rebuild it with toggled label -/// and identical bra/ket/aux/symmetry/etc. -ExprPtr conjugate_tensor(const Tensor& t) { - auto new_label = toggle_conj_suffix(t.label()); - // Tensor copy ctor + label edit isn't directly exposed; reconstruct. - return ex( - new_label, bra(container::svector{t.bra().begin(), t.bra().end()}), - ket(container::svector{t.ket().begin(), t.ket().end()}), - aux(container::svector{t.aux().begin(), t.aux().end()}), - t.symmetry(), t.braket_symmetry(), t.column_symmetry()); +/// Substitute Kramers labels onto @p expr per @p repl, expand +/// antisymmetrizers, then run the canonical SeQuant simplify pass. +ExprPtr specialize_and_simplify(const ExprPtr& expr, + const container::map& repl) { + auto term = append_kramers(expr, repl); + term = expand_antisymm(term, /* skip_spinsymm */ false); + return simplify_with_reset(term); } } // namespace -ExprPtr conjugate(const ExprPtr& expr) { - if (expr->is()) { - return conjugate_tensor(expr->as()); - } - if (expr->is()) { - auto v = expr->as().value(); - return ex(sequant::conj(v)); - } - if (expr->is()) { - throw std::runtime_error( - "conjugate: Variable conjugation not yet supported"); - } - if (expr->is()) { - auto const& prod = expr->as(); - auto p = std::make_shared(); - p->scale(sequant::conj(prod.scalar())); - for (auto&& f : prod) p->append(1, conjugate(f)); - return p; - } - if (expr->is()) { - auto s = std::make_shared(); - for (auto&& summand : *expr) s->append(conjugate(summand)); - return s; - } - throw std::runtime_error("conjugate: unsupported Expr type"); +Index make_kramers_up(const Index& idx) { + return make_index_with_kramers(idx, mbpt::Kramers::up); } -ExprPtr append_kramers(const ExprPtr& expr, - const container::map& index_replacements) { - auto add_to_tensor = [&](const Tensor& tensor) { - auto t = std::make_shared(tensor); - t->transform_indices(index_replacements); - return t; - }; - auto add_to_product = [&](const Product& product) { - auto p = std::make_shared(); - p->scale(product.scalar()); - for (auto&& term : product) { - if (term->is()) { - p->append(1, add_to_tensor(term->as())); - } else if (term->is() || term->is()) { - p->append(1, term); - } else { - throw std::runtime_error( - "append_kramers: unsupported Expr type in product: " + - term->type_name()); - } - } - return p; - }; +Index make_kramers_dn(const Index& idx) { + return make_index_with_kramers(idx, mbpt::Kramers::down); +} - if (expr->is()) return add_to_tensor(expr->as()); - if (expr->is()) return add_to_product(expr->as()); - if (expr->is()) { - auto s = std::make_shared(); - for (auto&& summand : *expr) { - s->append(append_kramers(summand, index_replacements)); - } - return s; - } - if (expr->is() || expr->is()) return expr; - throw std::runtime_error("append_kramers: unsupported Expr type"); +Index make_kramers_free(const Index& idx) { + return make_index_with_kramers(idx, mbpt::Kramers::any); } -namespace { +std::wstring toggle_conj_suffix(std::wstring_view label) { + if (has_conj_suffix(label)) + return std::wstring{label.substr(0, label.size() - 1)}; + return std::wstring{label} + L'*'; +} -/// Walk @p expr and collect the unique Kramers-eligible (i.e., -/// `Kramers::any`) indices in encounter order. Indices that already -/// have a specific Kramers state (up/down) are skipped. -container::svector collect_kramers_indices(const ExprPtr& expr) { - container::svector seen; - container::set seen_set; - auto consider = [&](const Index& idx) { - auto qns = idx.space().qns().to_int32(); - if ((qns & mask_v) != static_cast(Kramers::any)) { - // Either no Kramers state set or already specialized — skip. - return; - } - if (seen_set.insert(idx).second) seen.push_back(idx); - }; - expr->visit( - [&](const ExprPtr& e) { - if (e->is()) { - for (const auto& idx : e->as().const_braket()) { - consider(idx); - } - } - }, - /* recursive = */ true); - return seen; +ExprPtr append_kramers(const ExprPtr& expr, + const container::map& index_replacements) { + // The transformation is identical to spin's `append_spin`: in both + // cases the only operation is `Tensor::transform_indices`, which is + // QN-agnostic. Delegate to avoid duplicating the recursion logic. + return append_spin(expr, index_replacements); } -} // namespace +// ===================================================================== +// kramers_trace: 2^n exhaustive enumeration. +// ===================================================================== ExprPtr kramers_trace(const ExprPtr& expr) { - // Enumerate the 2^N Kramers configurations of the expression's spinor - // (Kramers::any) indices, then apply TRS pairing: configs (X, X̄) related - // by full bar-flip are complex conjugates of each other (per the - // (-1)^k full-bar identity, with the per-tensor signs cancelling for - // products in which every barred index appears in an even number of - // tensors — true for any all-internal-index contraction). - // - // TODO(CC): this all-internal-index assumption breaks for CC residuals, - // which carry free/external indices — those need explicit - // (-1)^(external bars) sign bookkeeping. See the kramers_trace @todo in - // spinor.hpp. - // - // For each TRS pair we materialize only the "lex-smaller" configuration - // X explicitly; the X̄ partner is emitted as `conjugate(T_X)`. Tensors - // in the conj branch carry a `*` label suffix (per the convention in - // the conjugate/real_part/imaginary_part scaffolding above) which the - // evaluator dispatches into a `.conj()` call on the underlying numeric - // tensor. Term count is unchanged, but the symbolic structure exposes - // the conjugate-pair relationship for downstream simplification (e.g., - // `real_part` collapse in real-scalar callers). - // - // Standard SeQuant simplifications (dummy renaming, antisym fold, etc.) - // are NOT applied here — caller invokes `sequant::canonicalize` if it - // wants those. - + // Enumerate 2^N Kramers configurations of @p expr's spinor + // (Kramers::any) indices, then sum. TRS pairs (X, X̄) related by full + // bar-flip are complex conjugates of each other under the (-1)^k + // identity; per-tensor signs cancel iff every barred index appears in + // an even number of tensors — true for any all-internal-index + // contraction. CC residuals carry external bars and need explicit + // (-1)^(external bars) bookkeeping; not yet implemented. + + detail::reset_idx_tags(expr); const auto indices = collect_kramers_indices(expr); const std::size_t n = indices.size(); if (n == 0) return expr; @@ -230,104 +175,181 @@ ExprPtr kramers_trace(const ExprPtr& expr) { "kramers_trace: too many Kramers indices (would emit > 2^30 terms)"); } - const std::uint64_t n_configs = 1ull << n; - - // Enumerate all 2^n Kramers configurations as separate Products. We - // mirror open_shell_spintrace's pipeline: substitute → expand_antisymm → - // expand → simplify → canonicalize, working entirely in expanded - // (NonSymm) tensor form so the standard SeQuant simplifier sees full - // index permutability for fold detection. auto result = std::make_shared(); + const std::uint64_t n_configs = std::uint64_t{1} << n; for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { - container::map repl; - for (std::size_t i = 0; i < n; ++i) { - const bool is_down = (cfg >> i) & 1u; - Index new_idx = - is_down ? make_kramers_dn(indices[i]) : make_kramers_up(indices[i]); - repl[indices[i]] = new_idx; - } - auto term = append_kramers(expr, repl); - term = expand_antisymm(term, /*skip_spinsymm*/ false); - expand(term); - rapid_simplify(term); - canonicalize(term); - rapid_simplify(term); - result->append(term); + result->append( + specialize_and_simplify(expr, kramers_replacement_map(indices, cfg))); } - // Final cross-product canonicalize: now that every term is in NonSymm - // expanded form, dummy renaming can fold equivalents across Kramers - // configurations. ExprPtr r = result; - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - return r; + return simplify_with_reset(r); +} + +container::svector kramers_trace(const ResultExpr& expr, + bool fold_klein) { + // Collect external indices from LHS (bra + ket; aux ignored — the + // pair-density use case doesn't have aux slots). + container::svector ext_indices; + for (const auto& idx : expr.bra()) ext_indices.push_back(idx); + for (const auto& idx : expr.ket()) ext_indices.push_back(idx); + + const std::size_t n_ext = ext_indices.size(); + if (n_ext > 30) { + throw std::runtime_error( + "kramers_trace(ResultExpr): too many external indices"); + } + if (n_ext == 0) { + // Scalar case: just run the regular tracer + pipeline and return + // a single ResultExpr (variable-shape; no LHS slots). + auto traced = kramers_trace(expr.expression()); + traced = fold_conj_pairs(traced); + traced = antisymm_recombine(traced); + container::svector single; + single.emplace_back( + bra(container::svector{}), ket(container::svector{}), + aux(container::svector{}), expr.symmetry(), + expr.braket_symmetry(), expr.column_symmetry(), + expr.has_label() ? std::optional(expr.label()) + : std::nullopt, + traced); + return single; + } + + // Process only canonical representatives under the Kramers + // symmetry group of the LHS. For a rank-(p,q) result the symmetry + // group is the Klein 4-group generated by: + // α: flip all bra-side Kramers states (partial flip on virtuals) + // β: flip all ket-side Kramers states (partial flip on occupieds) + // αβ: flip all Kramers states (full TRS) + // For Kramers-restricted closed-shell systems these are all + // symmetries of the trace (verifiable: traced expressions in one + // orbit are equal up to Kramers-conjugation relabeling). The orbit + // size is up to 4; emitting only the lex-smallest cfg in each orbit + // reduces 2^N raw blocks by a factor of 4 (typically). + const std::uint64_t n_configs = std::uint64_t{1} << n_ext; + const std::uint64_t mask = n_configs - 1; + const std::size_t n_bra = expr.bra().size(); + const std::uint64_t bra_mask = ((std::uint64_t{1} << n_bra) - 1); + const std::uint64_t ket_mask = mask ^ bra_mask; + + container::svector out; + for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { + const std::uint64_t cfg_alpha = cfg ^ bra_mask; + const std::uint64_t cfg_beta = cfg ^ ket_mask; + const std::uint64_t cfg_full = cfg ^ mask; + const std::uint64_t canonical = + fold_klein ? std::min({cfg, cfg_alpha, cfg_beta, cfg_full}) + : std::min(cfg, cfg_full); + if (cfg != canonical) continue; // skip non-canonical + + // Build Kramers-specialized externals for this config. + container::map ext_repl; + for (std::size_t k = 0; k < n_ext; ++k) { + const bool is_dn = (cfg >> k) & 1u; + ext_repl[ext_indices[k]] = is_dn ? make_kramers_dn(ext_indices[k]) + : make_kramers_up(ext_indices[k]); + } + + // Substitute externals on the RHS, then run the FULL kramers + // pipeline on the substituted expression: + // 1) kramers_trace — enumerates remaining Kramers::any indices + // (the internal dummies) + expand_antisymm + canonicalize + // per config + final cross-product canonicalize. After this + // the RHS is a Sum of NonSymm products. + // 2) antisymm_recombine — folds direct/exchange NonSymm pairs + // back into Kramers-restricted Antisymm tensors, restoring + // the mixed antisym/non-antisym MP2-style structure. + // + // Install a per-block `named_indices` context override carrying + // the kramers-substituted externals — every external appears + // exactly twice in the RHS (once via the LHS-side slot, once via + // the dummy contraction), so without the override SeQuant's + // canonicalize would treat them as dummies and rename them. The + // override keeps the externals fixed and consistent with the + // kramers-specialized LHS we build below. + // + // We do NOT call fold_conj_pairs here: it folds TRS pairs WITHIN + // a Sum, but within a single external-fixed block the externals + // are already specialized — flip_kramers would change them too + // and the partner would land in a different block. External TRS + // folding is handled at the block level above (cfg <= cfg_flipped + // canonicalization). + auto rhs_subs = append_kramers(expr.expression(), ext_repl); + container::set specialized_externals; + for (const auto& kv : ext_repl) specialized_externals.insert(kv.second); + auto saved_ctx = get_default_context(); + auto block_ctx = saved_ctx; + CanonicalizeOptions copts; + copts.named_indices = specialized_externals; + block_ctx.set(copts); + set_default_context(block_ctx); + auto traced = kramers_trace(rhs_subs); + traced = antisymm_recombine(traced); + set_default_context(saved_ctx); + + // Drop blocks that vanished. + if (traced->is() && traced->as().is_zero()) continue; + if (traced->is() && traced->as().summands().empty()) continue; + + // Build the kramers-specialized LHS slots in input order. + container::svector spec_bra, spec_ket; + for (const auto& idx : expr.bra()) spec_bra.push_back(ext_repl.at(idx)); + for (const auto& idx : expr.ket()) spec_ket.push_back(ext_repl.at(idx)); + + out.emplace_back(bra(std::move(spec_bra)), ket(std::move(spec_ket)), + aux(container::svector{}), expr.symmetry(), + expr.braket_symmetry(), expr.column_symmetry(), + expr.has_label() + ? std::optional(expr.label()) + : std::nullopt, + traced); + } + return out; } // ===================================================================== -// Conjugate-pair fold pass: detects TRS-related Products in a Sum and -// rewrites (A + B) → 2·Re(A) where B = flip_kramers(A). +// flip_kramers + fold_conj_pairs. // ===================================================================== ExprPtr flip_kramers(const ExprPtr& expr) { - // Walk all tensors and collect every Kramers-tagged index. For each - // distinct (label, ordinal, space-type) seen, build a partner index - // with the opposite Kramers state. Then apply the swap via the existing - // append_kramers replacement machinery. + // Walk every tensor, build the Index → flipped-Index map for any + // index already carrying a specialized Kramers state, then apply via + // append_kramers. container::map repl; expr->visit( [&](const ExprPtr& e) { if (!e->is()) return; for (const auto& idx : e->as().const_braket()) { if (repl.contains(idx)) continue; - const auto qns_int = idx.space().qns().to_int32(); - const auto k_bits = qns_int & mask_v; - if (k_bits == static_cast(Kramers::up)) { - repl[idx] = make_kramers_dn(idx); - } else if (k_bits == static_cast(Kramers::down)) { - repl[idx] = make_kramers_up(idx); + switch (index_kramers_state(idx)) { + case Kramers::up: + repl[idx] = make_kramers_dn(idx); + break; + case Kramers::down: + repl[idx] = make_kramers_up(idx); + break; + default: + break; // Kramers::any or unset — leave alone } - // Kramers::any or no Kramers state — leave alone. } }, /* recursive = */ true); - if (repl.empty()) return expr; - return append_kramers(expr, repl); + return repl.empty() ? expr : append_kramers(expr, repl); } namespace { -/// Normalize a Product-like expression for structural comparison: clone, -/// apply canonicalize + rapid_simplify so dummy renaming and ordering -/// match other already-canonicalized summands. -ExprPtr normalize(const ExprPtr& e) { - ExprPtr c = e->clone(); - expand(c); - rapid_simplify(c); - canonicalize(c); - rapid_simplify(c); - return c; -} - -/// Compare two scalars as complex rationals. -bool scalar_eq(const Constant::scalar_type& a, const Constant::scalar_type& b) { - return a == b; -} - -/// Extract the leading scalar of a Product, or 1 if the expr is a bare -/// Tensor. +/// Extract the leading scalar of a Product, or 1 for a non-Product. Constant::scalar_type product_scalar(const ExprPtr& e) { - if (e->is()) return e->as().scalar(); - return Constant::scalar_type{1}; + return e->is() ? e->as().scalar() + : Constant::scalar_type{1}; } -/// Build a "factor-only" copy of a Product (drop the leading scalar). -/// Returns @p e unchanged if it isn't a Product. +/// Build a copy of @p e with the leading Product scalar dropped. +/// Returns a clone if @p e is not a Product. ExprPtr product_factors_only(const ExprPtr& e) { if (!e->is()) return e->clone(); auto p = std::make_shared(); - p->scale(Constant::scalar_type{1}); for (auto&& f : e->as()) p->append(1, f->clone()); return p; } @@ -335,30 +357,20 @@ ExprPtr product_factors_only(const ExprPtr& e) { } // namespace ExprPtr fold_conj_pairs(const ExprPtr& expr) { + detail::reset_idx_tags(expr); if (!expr->is()) return expr; auto const& sum = expr->as(); const std::size_t n = sum.size(); - // O(n) canonical-form hash bucketing. For each summand A, compute a - // TRS-invariant key by taking the lex-smaller of normalize(A) and - // normalize(flip_kramers(A)). A and its full-bar partner B = flip(A) - // both map to the same key, so they land in the same bucket. - // - // Within a bucket we re-verify structural equality (hash collisions are - // possible; equality is O(tensor-count) per pair). Because TRS pairs - // form orbits of size ≤ 2, buckets are tiny and per-bucket work is O(1) - // on average. Total cost: O(n) normalizations + O(n) flips + O(n) hash - // inserts, vs. O(n²) for naive pairwise search. - // - // Scalar handling: we also bucket by scalar relationship. Two summands - // with matching scalars (a, a) collapse to `2·Re(A)`; opposite scalars - // (a, -a) collapse to `2i·Im(A)`. Self-conjugate singletons emit - // `Re(A)` (real-valued). + // O(n) canonical-form hash bucketing: for each summand A, the + // canonical key is the lex-min of (normalize(A), normalize(flip(A))). + // Both A and its full-bar partner map to the same bucket. Per-bucket + // work is O(1) on average since TRS orbits have size ≤ 2. struct BucketEntry { std::size_t idx; - ExprPtr factors_normalized; // factor part, post-normalize - Constant::scalar_type scalar; // leading scalar prefactor - bool self_conj; // factors == flip(factors) + ExprPtr factors_normalized; + Constant::scalar_type scalar; + bool self_conj; }; container::map> buckets; @@ -367,32 +379,28 @@ ExprPtr fold_conj_pairs(const ExprPtr& expr) { for (std::size_t i = 0; i < n; ++i) { originals.push_back(sum.summand(i)->clone()); - auto factors = normalize(product_factors_only(originals[i])); - auto factors_flipped = normalize(flip_kramers(factors)); + auto factors = product_factors_only(originals[i]); + simplify_with_reset(factors); + auto factors_flipped = flip_kramers(factors); + simplify_with_reset(factors_flipped); const bool self_conj = (*factors == *factors_flipped); - // Canonical factor key: lex-min of (factors, factors_flipped). Both A - // and flip(A) produce the same key, so they collide in the same bucket. ExprPtr canon = (*factors < *factors_flipped) ? factors : factors_flipped; - const auto h = canon->hash_value(); - buckets[h].push_back( - BucketEntry{i, factors, product_scalar(originals[i]), self_conj}); + buckets[canon->hash_value()].push_back( + {i, factors, product_scalar(originals[i]), self_conj}); } auto out = std::make_shared(); std::vector used(n, false); - // Helper to emit a folded summand based on scalar relationship. + // Emit the folded form of (A, B) per their scalar relationship. auto emit_pair = [&](const BucketEntry& a, const BucketEntry& b) { - if (scalar_eq(a.scalar, b.scalar)) { - // (A + B) = 2 Re(A) + using cscalar = Constant::scalar_type; + if (a.scalar == b.scalar) { out->append(ex(rational{2}) * real_part(originals[a.idx])); - } else if (scalar_eq(a.scalar, -b.scalar)) { - // (A - B) = 2i Im(A) (B carries -A's scalar) - using cscalar = Constant::scalar_type; + } else if (a.scalar == -b.scalar) { out->append(ex(cscalar{Complex{0, 2}}) * imaginary_part(originals[a.idx])); } else { - // Same factor-canonical key but unrelated scalars — emit both raw. out->append(originals[a.idx]); out->append(originals[b.idx]); } @@ -400,22 +408,18 @@ ExprPtr fold_conj_pairs(const ExprPtr& expr) { }; for (auto& [_, entries] : buckets) { - // Verify hash-bucket entries actually share the same canonical factor - // form (defend against hash collisions): partition by structural - // equality on factors_normalized OR its flip. std::vector taken(entries.size(), false); for (std::size_t i = 0; i < entries.size(); ++i) { if (taken[i] || used[entries[i].idx]) continue; const auto& a = entries[i]; - // Self-conjugate singleton in this bucket? if (a.self_conj) { - // Look for another self-conj with same scalar — if none, emit Re(A). + // Self-conjugate: pair with another self-conj if present, else + // emit the lone Re(A). bool merged = false; for (std::size_t j = i + 1; j < entries.size(); ++j) { if (taken[j] || used[entries[j].idx]) continue; const auto& b = entries[j]; if (b.self_conj && *b.factors_normalized == *a.factors_normalized) { - // Two self-conj copies → 2·Re(A) (or scalar-aware emit). emit_pair(a, b); taken[i] = taken[j] = true; merged = true; @@ -429,17 +433,14 @@ ExprPtr fold_conj_pairs(const ExprPtr& expr) { } continue; } - // Look for the partner: another entry whose factors equal flip(A). - // We don't have flip(A) materialized here — use the structural test - // (a.factors and b.factors land in the same bucket only when one is - // the flip of the other). + // Non-self-conj: confirm partner via structural equality with + // flip(A) (defending against hash collisions). bool paired = false; for (std::size_t j = i + 1; j < entries.size(); ++j) { if (taken[j] || used[entries[j].idx]) continue; const auto& b = entries[j]; - // Confirm they're actually a TRS pair: b.factors should equal - // flip(a.factors). - auto a_flip = normalize(flip_kramers(a.factors_normalized)); + auto a_flip = flip_kramers(a.factors_normalized); + simplify_with_reset(a_flip); if (*b.factors_normalized == *a_flip) { emit_pair(a, b); taken[i] = taken[j] = true; @@ -459,13 +460,32 @@ ExprPtr fold_conj_pairs(const ExprPtr& expr) { } // ===================================================================== -// Cycle decomposition. See spinor.hpp for the design summary. +// Cycle decomposition (used by kramers_trace_cycles). // ===================================================================== -container::svector kramers_cycles(const Product& product) { - // Collect (tensor_idx, braket_pos, Index) for every position in every - // Tensor factor. Indices that aren't Tensor are skipped (Constants, - // Variables don't contribute to cycles). +namespace { + +/// One step of a contraction-cycle walk. +struct CycleNode { + std::size_t tensor_idx; + std::size_t braket_pos; + Index idx; +}; + +/// One contraction cycle: alternating intra-tensor (bra[k] ↔ ket[k] of +/// the same tensor) and inter-tensor (same Index in two tensors) edges. +struct ContractionCycle { + container::svector nodes; +}; + +/// Decompose a Product's index-contraction graph into cycles. +/// +/// @pre every Index in @p product appears in exactly two tensor slot +/// positions; each tensor's bra[k]/ket[k] positions pair as +/// intra-tensor edges +/// @todo CC support: the k ↔ k+rank/2 pairing assumes a rectangular +/// tensor (bra_rank == ket_rank); some CC intermediates aren't. +container::svector decompose_cycles(const Product& product) { struct Slot { std::size_t tensor_idx; std::size_t braket_pos; @@ -480,33 +500,21 @@ container::svector kramers_cycles(const Product& product) { } auto const& t = f->as(); std::size_t pos = 0; - for (auto const& idx : t.const_braket()) { - slots.push_back(Slot{tensor_counter, pos, idx}); - ++pos; - } + for (auto const& idx : t.const_braket()) + slots.push_back(Slot{tensor_counter, pos++, idx}); ++tensor_counter; } const std::size_t N = slots.size(); container::svector visited(N, false); - - // Pre-build two adjacency maps: - // - intra-tensor: bra[k] ↔ ket[k] of the same tensor (positions 0↔rank, - // 1↔rank+1, ... where rank = bra_rank). - // - inter-tensor: same Index occurrence in two distinct slots. - container::svector intra_partner(N, N); // N = "no partner" + container::svector intra_partner(N, N); container::svector inter_partner(N, N); - // intra-tensor pairing per tensor: walk positions and pair k ↔ k+rank/2 - // within each tensor's slots. We rebuild rank by counting consecutive - // slots with the same tensor_idx. - // TODO(CC): the k ↔ k+rank/2 pairing assumes a rectangular tensor - // (bra_rank == ket_rank), which holds for g/t/residual tensors but not - // for every CC intermediate; non-rectangular tensors need the bra/ket - // split read from the tensor directly rather than from total rank. + // Intra-tensor pairing per tensor (k ↔ k+rank/2 within each + // contiguous slot run, assuming rectangular). for (std::size_t i = 0; i < N;) { std::size_t j = i; while (j < N && slots[j].tensor_idx == slots[i].tensor_idx) ++j; - const std::size_t tensor_rank = j - i; // total bra+ket count + const std::size_t tensor_rank = j - i; SEQUANT_ASSERT(tensor_rank % 2 == 0); const std::size_t half = tensor_rank / 2; for (std::size_t k = 0; k < half; ++k) { @@ -516,13 +524,11 @@ container::svector kramers_cycles(const Product& product) { i = j; } - // inter-tensor pairing: same Index in two distinct positions across the - // whole Product. Use Index-equality (which compares full label). + // Inter-tensor pairing: same Index in two distinct slots. for (std::size_t i = 0; i < N; ++i) { if (inter_partner[i] != N) continue; for (std::size_t j = i + 1; j < N; ++j) { - if (inter_partner[j] != N) continue; - if (slots[j].idx == slots[i].idx) { + if (inter_partner[j] == N && slots[j].idx == slots[i].idx) { inter_partner[i] = j; inter_partner[j] = i; break; @@ -530,223 +536,48 @@ container::svector kramers_cycles(const Product& product) { } } - // Walk cycles. From each unvisited slot, alternate intra → inter → - // intra → ... edges until returning to the start. + // Walk cycles: alternate intra → inter → intra → ... edges. container::svector cycles; for (std::size_t start = 0; start < N; ++start) { if (visited[start]) continue; ContractionCycle cyc; std::size_t cur = start; - bool take_intra = true; // first hop is intra-tensor + bool take_intra = true; while (!visited[cur]) { visited[cur] = true; - cyc.nodes.push_back(ContractionCycle::Node{ - slots[cur].tensor_idx, slots[cur].braket_pos, slots[cur].idx}); + cyc.nodes.push_back( + {slots[cur].tensor_idx, slots[cur].braket_pos, slots[cur].idx}); const std::size_t next = take_intra ? intra_partner[cur] : inter_partner[cur]; - if (next == N) break; // dangling (shouldn't happen for closed expr) + if (next == N) break; cur = next; take_intra = !take_intra; } cycles.push_back(std::move(cyc)); } - return cycles; } -namespace { - -/// Format a cycle's Kramers labelling as a wstring like "UAUA" or "UDDU" -/// where each char is the visited index's Kramers state. -std::wstring cycle_kramers_label(const ContractionCycle& c) { - std::wstring s; - s.reserve(c.nodes.size()); - for (auto const& n : c.nodes) { - const auto qns_int = n.idx.space().qns().to_int32(); - const auto k_bits = qns_int & mask_v; - if (k_bits == static_cast(Kramers::up)) { - s += L'U'; - } else if (k_bits == static_cast(Kramers::down)) { - s += L'B'; - } else { - s += L'?'; - } - } - return s; -} - -/// Number of canonical Kramers patterns of a length-L cycle under -/// cyclic-rotation symmetry. By Burnside: (1/L) Σ_{d | L} φ(L/d) · 2^d. -std::size_t canonical_pattern_count_cyclic(std::size_t L) { - if (L == 0) return 1; - // φ(n) - auto phi = [](std::size_t n) { - std::size_t result = n; - for (std::size_t p = 2; p * p <= n; ++p) { - if (n % p == 0) { - while (n % p == 0) n /= p; - result -= result / p; - } - } - if (n > 1) result -= result / n; - return result; - }; - std::size_t total = 0; - for (std::size_t d = 1; d <= L; ++d) { - if (L % d == 0) { - total += phi(L / d) * (std::size_t{1} << d); - } - } - return total / L; -} - -} // namespace - -void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os) { - auto dump_product = [&](const Product& p, std::size_t pidx) { - os << L" Product[" << pidx << L"]:"; - for (auto&& f : p) { - if (f->is()) { - auto const& t = f->as(); - os << L" " << t.label(); - } - } - os << L"\n"; - auto cycles = kramers_cycles(p); - os << L" " << cycles.size() << L" cycle(s):\n"; - for (std::size_t ci = 0; ci < cycles.size(); ++ci) { - auto const& c = cycles[ci]; - // Distinct indices visited. - container::set distinct; - for (auto const& n : c.nodes) distinct.insert(n.idx); - const auto canon = canonical_pattern_count_cyclic(c.nodes.size()); - os << L" cycle[" << ci << L"] len=" << c.nodes.size() - << L" distinctIdx=" << distinct.size() << L" labels=`" - << cycle_kramers_label(c) << L"`" << L" cyclic-canonical-#patterns=" - << canon << L" walk="; - bool first = true; - for (auto const& n : c.nodes) { - if (!first) os << L"→"; - os << L"t" << n.tensor_idx << L"/" << n.braket_pos << L":" - << n.idx.label(); - first = false; - } - os << L"\n"; - } - }; - - if (expr->is()) { - dump_product(expr->as(), 0); - } else if (expr->is()) { - auto const& s = expr->as(); - for (std::size_t i = 0; i < s.size(); ++i) { - auto const& term = s.summand(i); - if (term->is()) dump_product(term->as(), i); - } - } -} - -std::wstring cycle_canonical_label(const ContractionCycle& c) { - // Build raw label, then take lex-min over all cyclic rotations. - const auto raw = cycle_kramers_label(c); - if (raw.empty()) return raw; - std::wstring best = raw; - std::wstring rot = raw; - for (std::size_t k = 1; k < raw.size(); ++k) { - // rotate by 1: move first char to end - std::rotate(rot.begin(), rot.begin() + 1, rot.end()); - if (rot < best) best = rot; - } - return best; -} - -std::wstring cycle_canonical_signature(const Product& product) { - auto cycles = kramers_cycles(product); - container::svector labels; - labels.reserve(cycles.size()); - for (auto const& c : cycles) labels.push_back(cycle_canonical_label(c)); - std::sort(labels.begin(), labels.end()); - std::wstring sig; - for (std::size_t i = 0; i < labels.size(); ++i) { - if (i > 0) sig += L'|'; - sig += labels[i]; - } - return sig; -} - -ExprPtr cycle_canonical_fold(const ExprPtr& expr) { - if (!expr->is()) return expr; - auto const& sum = expr->as(); - // Bucket Sum summands by cycle_canonical_signature. Within each bucket - // the contractions are numerically equal up to scalar prefactor; sum - // the prefactors and emit a single representative per bucket. - container::map> buckets; - std::vector originals; - originals.reserve(sum.size()); - for (std::size_t i = 0; i < sum.size(); ++i) { - originals.push_back(sum.summand(i)->clone()); - if (!originals[i]->is()) { - // Non-Product summands (Tensor, Constant, Variable) — emit a - // distinct bucket per index so they pass through unchanged. - buckets[L"__pass" + std::to_wstring(i)].push_back(i); - continue; - } - const auto sig = cycle_canonical_signature(originals[i]->as()); - buckets[sig].push_back(i); - } - - auto out = std::make_shared(); - for (auto& [sig, idxs] : buckets) { - if (idxs.size() == 1) { - out->append(originals[idxs[0]]); - continue; - } - // Sum scalars; reuse first member's tensor structure as the - // representative. NOTE: this assumes the bucketed Products are - // numerically equal as contractions, which holds when the cycles' - // canonical Kramers labels match — that's the bucket-key meaning. - Constant::scalar_type total{0}; - for (auto i : idxs) { - total += product_scalar(originals[i]); - } - if (total == Constant::scalar_type{0}) continue; - auto rep = product_factors_only(originals[idxs.front()]); - out->append(ex(total) * rep); - } - return out; -} - -// ===================================================================== -// Standalone cycle-driven Kramers tracer. See spinor.hpp. -// ===================================================================== - -namespace { - -/// Distinct Kramers-eligible (`Kramers::any`) indices visited by a cycle, -/// in first-encounter order. +/// Distinct Kramers-eligible (`Kramers::any`) indices visited by a +/// cycle, in first-encounter order. container::svector cycle_distinct_indices(const ContractionCycle& c) { container::svector out; container::set seen; for (auto const& n : c.nodes) { - const auto qns = n.idx.space().qns().to_int32(); - if ((qns & mask_v) != static_cast(Kramers::any)) - continue; // already specialized or not Kramers — skip + if (index_kramers_state(n.idx) != Kramers::any) continue; if (seen.insert(n.idx).second) out.push_back(n.idx); } return out; } -/// Structural shape signature of a cycle: a string encoding the walk's -/// IndexSpace-type sequence, made rotation-invariant by taking the -/// lex-min rotation. Two cycles with the same shape are interchangeable -/// under the expression's tensor (anti)symmetries — the caller is -/// responsible for ensuring those symmetries actually hold. +/// Rotation-invariant shape signature of a cycle: the lex-min rotation +/// of the IndexSpace-type sequence along the walk. Two cycles with the +/// same shape are interchangeable under the expression's tensor +/// (anti)symmetries. std::wstring cycle_shape(const ContractionCycle& c) { std::wstring raw; - for (auto const& n : c.nodes) { - // encode IndexSpace base type as a hex digit (stable per registry) + for (auto const& n : c.nodes) raw += static_cast(L'a' + (n.idx.space().type().to_int32() & 0xF)); - } if (raw.empty()) return raw; std::wstring best = raw, rot = raw; for (std::size_t k = 1; k < raw.size(); ++k) { @@ -756,17 +587,7 @@ std::wstring cycle_shape(const ContractionCycle& c) { return best; } -/// All 2^n Kramers labelings of a cycle's @p n distinct indices, each a -/// bitset (bit i = 1 ⇔ index i is Kramers-down). -container::svector cycle_labelings(std::size_t n_distinct) { - container::svector out; - const std::uint32_t n_cfg = 1u << n_distinct; - out.reserve(n_cfg); - for (std::uint32_t c = 0; c < n_cfg; ++c) out.push_back(c); - return out; -} - -/// Multinomial coefficient m! / ∏ count_i! for a multiset assignment. +/// Multinomial coefficient (Σ counts)! / ∏ count_i!. double multiset_multiplicity(const container::svector& counts) { auto factorial = [](std::size_t k) { double f = 1; @@ -783,21 +604,15 @@ double multiset_multiplicity(const container::svector& counts) { } // namespace ExprPtr kramers_trace_cycles(const ExprPtr& expr) { - // Sum: recurse over summands. + detail::reset_idx_tags(expr); if (expr->is()) { auto sum_result = std::make_shared(); for (auto&& s : *expr) sum_result->append(kramers_trace_cycles(s)); ExprPtr r = sum_result; - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - return r; + return simplify_with_reset(r); } - // Constant / Variable: no Kramers indices, pass through. if (expr->is() || expr->is()) return expr; - // Normalize to a Product. Product prod; if (expr->is()) { prod = expr->as(); @@ -807,34 +622,23 @@ ExprPtr kramers_trace_cycles(const ExprPtr& expr) { throw std::runtime_error("kramers_trace_cycles: unsupported Expr type"); } - // 1. Decompose into cycles. - const auto cycles = kramers_cycles(prod); + // 1. Decompose into cycles + collect per-cycle Kramers-eligible indices. + const auto cycles = decompose_cycles(prod); + if (cycles.empty()) return expr; const std::size_t n_cycles = cycles.size(); - if (n_cycles == 0) return expr; - - // Per-cycle distinct Kramers indices + labelings. std::vector> cyc_indices(n_cycles); - std::vector> cyc_labelings(n_cycles); - for (std::size_t c = 0; c < n_cycles; ++c) { + for (std::size_t c = 0; c < n_cycles; ++c) cyc_indices[c] = cycle_distinct_indices(cycles[c]); - cyc_labelings[c] = cycle_labelings(cyc_indices[c].size()); - } - // 2. Group cycles into structural-equivalence classes by shape. + // 2. Group cycles into shape-equivalence classes. container::map> shape_classes; for (std::size_t c = 0; c < n_cycles; ++c) shape_classes[cycle_shape(cycles[c])].push_back(c); - // 3. Per class, enumerate Kramers-labeling MULTISETS (sorted tuples of - // per-cycle labelings) with their multiplicity. A class of m cycles - // each with L labelings yields C(L+m-1, m) multisets. - // - // Represent a class assignment as a sorted vector of - // length m (one labeling index per cycle in the class). Multiplicity - // = multinomial(counts of each distinct labeling). + // 3. Per class: enumerate Kramers-labeling MULTISETS (sorted m-tuples + // of per-cycle 2^d configs) with their multinomial multiplicity. struct ClassConfig { - container::svector class_ids; // cycle indices in this class - // each entry: (sorted labeling-assignment, multiplicity) + container::svector class_ids; container::svector, double>> multisets; }; @@ -843,19 +647,14 @@ ExprPtr kramers_trace_cycles(const ExprPtr& expr) { ClassConfig cc; cc.class_ids = cyc_ids; const std::size_t m = cyc_ids.size(); - // All cycles in a class share labeling-set size (same shape ⇒ same - // #distinct indices). Use the first cycle's labeling set. - const auto& labels = cyc_labelings[cyc_ids.front()]; - const std::size_t L = labels.size(); - // Enumerate sorted m-tuples (multisets) from L labelings. + const std::size_t L = std::size_t{1} << cyc_indices[cyc_ids.front()].size(); container::svector pick(m, 0); std::function rec = [&](std::size_t pos, std::size_t start) { if (pos == m) { - // pick[] holds indices into `labels`, non-decreasing. container::svector assignment(m); - for (std::size_t i = 0; i < m; ++i) assignment[i] = labels[pick[i]]; - // multiplicity = multinomial over runs of equal pick values + for (std::size_t i = 0; i < m; ++i) + assignment[i] = static_cast(pick[i]); container::svector counts; std::size_t run = 1; for (std::size_t i = 1; i <= m; ++i) { @@ -879,15 +678,13 @@ ExprPtr kramers_trace_cycles(const ExprPtr& expr) { } // 4. Cartesian product across classes: each global config picks one - // multiset per class. Build the per-index Kramers replacement, - // substitute, expand_antisymm, canonicalize, scale by the product - // of per-class multiplicities. + // multiset per class; build the per-index Kramers replacement, + // specialize, and scale by the product of multiplicities. auto result = std::make_shared(); const std::size_t n_classes = class_configs.size(); container::svector sel(n_classes, 0); std::function build = [&](std::size_t ci) { if (ci == n_classes) { - // Assemble the full per-index Kramers replacement. container::map repl; double mult = 1.0; for (std::size_t c = 0; c < n_classes; ++c) { @@ -904,16 +701,10 @@ ExprPtr kramers_trace_cycles(const ExprPtr& expr) { } } } - auto term = append_kramers(prod.clone(), repl); - term = expand_antisymm(term, /*skip_spinsymm*/ false); - expand(term); - rapid_simplify(term); - canonicalize(term); - rapid_simplify(term); + auto term = specialize_and_simplify(prod.clone(), repl); if (mult != 1.0) { term = ex(rational{static_cast(mult)}) * term; - expand(term); - rapid_simplify(term); + non_canon_simplify(term); } result->append(term); return; @@ -925,38 +716,32 @@ ExprPtr kramers_trace_cycles(const ExprPtr& expr) { }; build(0); - // Final cross-term canonicalize so equivalent multiset terms fold. ExprPtr r = result; - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - return r; + return simplify_with_reset(r); } // ===================================================================== -// Antisymmetrizer recombination (Phase 2.11). Inverse of expand_antisymm. -// See spinor.hpp for the algorithm. +// Antisymmetrizer recombination (inverse of expand_antisymm). // ===================================================================== namespace { -/// A parsed summand: total scalar, wrapper kind (0 raw / 1 Re / 2 Im), and -/// the tensor factors. The inner Product's scalar is folded into `scalar`. +/// Parsed Sum summand: total scalar, wrapper kind (0 raw / 1 Re / 2 +/// Im), and the underlying tensor factors. struct RecombEntry { Constant::scalar_type scalar{1}; int wrapper = 0; - container::svector tensors; // each is a Tensor ExprPtr + container::svector tensors; bool alive = true; }; -/// Decompose a Sum summand into a RecombEntry. Recognized shapes: -/// Tensor | Product(s; Tensor...) | Re/Im(inner) | Product(s; Re/Im(inner)) +/// Parse a Sum summand into a RecombEntry. Recognized shapes: +/// Tensor | Product(s; Tensor...) | Re/Im(inner) +/// | Product(s; Re/Im(inner)) /// where inner = Tensor | Product(s; Tensor...). bool parse_recomb_entry(const ExprPtr& summand, RecombEntry& e) { e = RecombEntry{}; ExprPtr cur = summand; - // optional outer `Constant * (Re|Im wrapper)` if (cur->is() && cur->as().size() == 1) { auto const& p = cur->as(); auto const& f0 = p.factor(0); @@ -990,11 +775,9 @@ bool parse_recomb_entry(const ExprPtr& summand, RecombEntry& e) { /// Detect whether @p t2 is @p t1 with a single column transposition. /// Returns 1 if t2 = bra-swap(t1), 2 if t2 = ket-swap(t1), 0 otherwise. -/// Both indices compared by full Index identity (recombinable terms share -/// dummy names post-canonicalize). -/// TODO(CC): rank-(2,2) only. CCSDT triples residuals need this -/// generalized to rank-(n,n) — enumerate the column transpositions of an -/// n-column tensor and report which single one (if any) maps t1 -> t2. +/// +/// @todo CC support: rank-(2,2) only. CCSDT triples residuals need the +/// single-column-swap detection generalized to rank-(n,n). int detect_single_column_swap(const Tensor& t1, const Tensor& t2) { if (t1.label() != t2.label()) return 0; if (t1.bra_rank() != 2 || t1.ket_rank() != 2) return 0; @@ -1013,8 +796,8 @@ int detect_single_column_swap(const Tensor& t1, const Tensor& t2) { return 0; } -/// Re-tag a Tensor with Symmetry::Antisymm (keeping label/indices/other -/// symmetries). Used to materialize a recombined antisymmetrized tensor. +/// Re-tag a Tensor with `Symmetry::Antisymm` (Tensor::symmetry is +/// constructor-set; we rebuild via the existing Tensor ctor). ExprPtr make_antisymm(const Tensor& t) { return ex( t.label(), bra(container::svector{t.bra().begin(), t.bra().end()}), @@ -1023,11 +806,10 @@ ExprPtr make_antisymm(const Tensor& t) { Symmetry::Antisymm, t.braket_symmetry(), t.column_symmetry()); } -/// Try to recombine two entries into one antisymmetrized entry. -/// Succeeds when, for some tensor position p: all other positions match -/// structurally, t2[p] is a single column-swap of t1[p], and the scalars -/// satisfy B.scalar == -A.scalar (single transposition → odd → sign -1). -/// Result: A.scalar * (... A.tensors[p] tagged Antisymm ...). +/// Try to recombine @p A and @p B into one antisymmetrized entry. +/// Succeeds when, for some tensor position p: all other positions +/// match structurally, t2[p] is a single column-swap of t1[p], and +/// `B.scalar == -A.scalar` (single transposition → odd → sign -1). std::optional try_recombine(const RecombEntry& A, const RecombEntry& B) { if (A.wrapper != B.wrapper) return std::nullopt; @@ -1045,9 +827,9 @@ std::optional try_recombine(const RecombEntry& A, } if (!others_match) continue; if (!A.tensors[p]->is() || !B.tensors[p]->is()) continue; - const int swap = detect_single_column_swap(A.tensors[p]->as(), - B.tensors[p]->as()); - if (swap == 0) continue; + if (detect_single_column_swap(A.tensors[p]->as(), + B.tensors[p]->as()) == 0) + continue; RecombEntry R = A; R.tensors[p] = make_antisymm(A.tensors[p]->as()); return R; @@ -1058,6 +840,7 @@ std::optional try_recombine(const RecombEntry& A, } // namespace ExprPtr antisymm_recombine(const ExprPtr& expr) { + detail::reset_idx_tags(expr); if (!expr->is()) return expr; auto const& sum = expr->as(); @@ -1069,8 +852,8 @@ ExprPtr antisymm_recombine(const ExprPtr& expr) { entries.push_back(std::move(e)); } - // Iterate to a fixed point: each pass merges one recombinable pair. A - // second pass over the once-recombined terms catches doubly-antisymm + // Iterate to a fixed point: each pass merges one recombinable pair. + // A second pass over once-recombined terms catches doubly-antisymm // blocks (the merged antisymm tensor becomes the "fixed" tensor). bool changed = true; while (changed) { @@ -1106,24 +889,13 @@ ExprPtr antisymm_recombine(const ExprPtr& expr) { } // ===================================================================== -// Burnside-enumeration Kramers tracer (Phase 2.12 prototype). -// See spinor.hpp for the algorithm summary. +// Burnside-orbit Kramers tracer. // ===================================================================== namespace { -/// Collect the Tensor factors of a Product (or wrap a bare Tensor). -container::svector product_tensors(const Product& p) { - container::svector out; - for (auto&& f : p) { - if (f->is()) out.push_back(&f->as()); - } - return out; -} - -/// True iff indices @p p and @p q both sit in the same antisymmetric -/// bra OR ket group of Tensor @p t (i.e. t is Antisymm and {p,q} ⊆ bra -/// or {p,q} ⊆ ket). +/// True iff @p p and @p q both sit in the same antisymmetric bra/ket +/// group of @p t (i.e. t is Antisymm and {p,q} ⊆ bra or {p,q} ⊆ ket). bool same_antisymm_group(const Tensor& t, const Index& p, const Index& q) { if (t.symmetry() != Symmetry::Antisymm) return false; auto in = [](const auto& range, const Index& x) { @@ -1134,25 +906,28 @@ bool same_antisymm_group(const Tensor& t, const Index& p, const Index& q) { return (p_bra && q_bra) || (p_ket && q_ket); } -/// True iff index @p p appears anywhere in Tensor @p t. bool tensor_has_index(const Tensor& t, const Index& p) { for (auto&& idx : t.const_braket()) if (idx == p) return true; return false; } +container::svector product_tensors(const Product& p) { + container::svector out; + for (auto&& f : p) + if (f->is()) out.push_back(&f->as()); + return out; +} + } // namespace ExprPtr kramers_trace_burnside(const ExprPtr& expr) { + detail::reset_idx_tags(expr); if (expr->is()) { auto s = std::make_shared(); for (auto&& t : *expr) s->append(kramers_trace_burnside(t)); ExprPtr r = s; - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - return r; + return simplify_with_reset(r); } if (expr->is() || expr->is()) return expr; @@ -1174,10 +949,9 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { } const auto tensors = product_tensors(prod); - // --- 1. Find transposition generators of the index-permutation - // symmetry group. (p,q) qualifies iff in every tensor that - // contains both they share an antisymmetric bra/ket group, they - // co-occur in at least one tensor, and the sign is +1. + // 1. Find transposition generators of the index-permutation symmetry + // group: (a,b) qualifies iff in every tensor that contains both + // they share an antisymmetric bra/ket group (sign +1 net). container::svector> generators; for (std::size_t a = 0; a < n; ++a) { for (std::size_t b = a + 1; b < n; ++b) { @@ -1193,26 +967,16 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { break; } } else if (ha != hb) { - // appears in only one slot of this tensor — a transposition - // would move it to a different tensor position; not a symmetry ok = false; break; } } - if (ok && n_common > 0 && (n_common % 2 == 0)) { + if (ok && n_common > 0 && n_common % 2 == 0) generators.emplace_back(a, b); - } } } - // --- 2. Build the permutation group as a BFS closure over generators. - // Each element is a permutation of {0..n-1} (svector). - auto apply_transposition = - [](container::svector perm, std::size_t a, - std::size_t b) -> container::svector { - std::swap(perm[a], perm[b]); - return perm; - }; + // 2. BFS-close the transposition generators into the full group. container::svector> group; { container::svector identity(n); @@ -1225,7 +989,8 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { for (auto const& g : frontier) { group.push_back(g); for (auto const& [a, b] : generators) { - auto ng = apply_transposition(g, a, b); + auto ng = g; + std::swap(ng[a], ng[b]); if (seen.insert(ng).second) next.push_back(ng); } } @@ -1233,8 +998,7 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { } } - // --- 3. Burnside orbit enumeration over the 2^n bit-string configs. - // A permutation acts: (perm·cfg) has bit perm[k] = bit k of cfg. + // 3. Burnside orbit enumeration over 2^n bit-string configs. auto apply_perm = [](std::uint64_t cfg, const container::svector& perm) { std::uint64_t out = 0; @@ -1244,181 +1008,29 @@ ExprPtr kramers_trace_burnside(const ExprPtr& expr) { }; const std::uint64_t n_configs = std::uint64_t{1} << n; container::svector visited(n_configs, false); - // each rep: (config, orbit_size) container::svector> reps; for (std::uint64_t cfg = 0; cfg < n_configs; ++cfg) { if (visited[cfg]) continue; container::set orbit; for (auto const& g : group) orbit.insert(apply_perm(cfg, g)); for (auto o : orbit) visited[o] = true; - // representative = lex-min of the orbit - const std::uint64_t rep = *orbit.begin(); - reps.emplace_back(rep, orbit.size()); + reps.emplace_back(*orbit.begin(), orbit.size()); } - // --- 4. Per orbit rep: substitute Kramers labels, expand_antisymm, - // canonicalize, scale by orbit size. + // 4. Per orbit rep: substitute, simplify, scale by orbit size. auto result = std::make_shared(); for (auto const& [cfg, orbit_size] : reps) { - container::map repl; - for (std::size_t k = 0; k < n; ++k) { - const bool is_dn = (cfg >> k) & 1u; - repl[indices[k]] = - is_dn ? make_kramers_dn(indices[k]) : make_kramers_up(indices[k]); - } - auto term = append_kramers(expr, repl); - term = expand_antisymm(term, /*skip_spinsymm*/ false); - expand(term); - rapid_simplify(term); - canonicalize(term); - rapid_simplify(term); + auto term = + specialize_and_simplify(expr, kramers_replacement_map(indices, cfg)); if (orbit_size != 1) { term = ex(rational{static_cast(orbit_size)}) * term; - expand(term); - rapid_simplify(term); + non_canon_simplify(term); } result->append(term); } - - // --- 5. Final cross-term canonicalize. ExprPtr r = result; - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - return r; -} - -// ===================================================================== -// Quaternion decomposer. See spinor.hpp for the design overview. -// ===================================================================== - -namespace { - -/// Reads the specialized Kramers state of an index. -Kramers index_kramers(const Index& idx) { - const auto k_bits = idx.space().qns().to_int32() & mask_v; - if (k_bits == static_cast(Kramers::up)) return Kramers::up; - if (k_bits == static_cast(Kramers::down)) return Kramers::down; - throw std::runtime_error( - "quaternion_decompose: tensor index is not Kramers-specialized"); -} - -/// Entry M(k_bra,k_ket)[c_bra][c_ket] of the quaternion decomposition -/// matrix. Component order s=0, x=1, y=2, z=3. Derived from the spinor -/// coefficient structure C_α = C_s + iC_x, C_β = -C_y + iC_z (barred -/// partner [-C_β*; C_α*]); see kramers-restriction §6 / Wiebke Eq. 8-11. -Constant::scalar_type quaternion_m(Kramers k_bra, Kramers k_ket, int c_bra, - int c_ket) { - using S = Constant::scalar_type; - auto v = [](int re, int im) { return S{Complex{re, im}}; }; - // (re,im) tables for M(up,up) and M(up,down); rows = c_bra, cols = c_ket. - static const int uu_re[4][4] = { - {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; - static const int uu_im[4][4] = { - {0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0}}; - static const int ud_re[4][4] = { - {0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 0}}; - static const int ud_im[4][4] = { - {0, 0, 0, 1}, {0, 0, -1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}; - const bool b_up = (k_bra == Kramers::up); - const bool k_up = (k_ket == Kramers::up); - if (b_up && k_up) return v(uu_re[c_bra][c_ket], uu_im[c_bra][c_ket]); - if (!b_up && !k_up) // M(down,down) = conj(M(up,up)) - return v(uu_re[c_bra][c_ket], -uu_im[c_bra][c_ket]); - if (b_up && !k_up) return v(ud_re[c_bra][c_ket], ud_im[c_bra][c_ket]); - // M(down,up) = M(up,down)^dagger - return v(ud_re[c_ket][c_bra], -ud_im[c_ket][c_bra]); -} - -constexpr char kQComp[4] = {'s', 'x', 'y', 'z'}; - -/// Decompose one rank-(2,2) Nonsymm tensor into its (up to 64) quaternion -/// components: T = Σ_c M(k_B0,k_K0)[c_B0,c_K0]·M(k_B1,k_K1)[c_B1,c_K1]·T#c. -ExprPtr decompose_nonsymm_tensor(const Tensor& t) { - if (t.bra_rank() != 2 || t.ket_rank() != 2) - throw std::runtime_error( - "quaternion_decompose: only rank-(2,2) tensors are supported"); - const bool conj = has_conj_suffix(t.label()); - const std::wstring core = - conj ? std::wstring{t.label().substr(0, t.label().size() - 1)} - : std::wstring{t.label()}; - const Index& B0 = t.bra().at(0); - const Index& B1 = t.bra().at(1); - const Index& K0 = t.ket().at(0); - const Index& K1 = t.ket().at(1); - // electron 1 = (B0, K0), electron 2 = (B1, K1) - const Kramers kB0 = index_kramers(B0), kB1 = index_kramers(B1); - const Kramers kK0 = index_kramers(K0), kK1 = index_kramers(K1); - - auto sum = std::make_shared(); - const Constant::scalar_type zero{0}; - for (int cB0 = 0; cB0 < 4; ++cB0) - for (int cK0 = 0; cK0 < 4; ++cK0) { - const auto m1 = quaternion_m(kB0, kK0, cB0, cK0); - if (m1 == zero) continue; - for (int cB1 = 0; cB1 < 4; ++cB1) - for (int cK1 = 0; cK1 < 4; ++cK1) { - const auto m2 = quaternion_m(kB1, kK1, cB1, cK1); - if (m2 == zero) continue; - auto coeff = m1 * m2; - if (conj) coeff = sequant::conj(coeff); - std::wstring lbl = quaternion_label_add( - core, kQComp[cB0], kQComp[cB1], kQComp[cK0], kQComp[cK1]); - if (conj) lbl += L'*'; - auto qt = - ex(lbl, bra{B0, B1}, ket{K0, K1}, Symmetry::Nonsymm, - BraKetSymmetry::Nonsymm, ColumnSymmetry::Nonsymm); - sum->append(ex(coeff) * qt); - } - } - return sum; -} - -} // namespace - -ExprPtr quaternion_decompose(const ExprPtr& expr) { - if (expr->is()) { - auto s = std::make_shared(); - for (auto&& summand : *expr) s->append(quaternion_decompose(summand)); - return s; - } - if (expr->is()) - return ex(quaternion_decompose(expr->as().inner())); - if (expr->is()) - return ex(quaternion_decompose(expr->as().inner())); - if (expr->is() || expr->is()) return expr; - if (expr->is()) { - auto const& t = expr->as(); - if (t.symmetry() == Symmetry::Antisymm) { - // expand g - g_exchange first, then decompose each Nonsymm piece - return quaternion_decompose( - expand_antisymm(expr, /*skip_spinsymm*/ false)); - } - return decompose_nonsymm_tensor(t); - } - if (expr->is()) { - auto const& p = expr->as(); - bool has_wrapper = false; - for (auto&& f : p) - if (f->is() || f->is()) has_wrapper = true; - auto prod = std::make_shared(); - prod->scale(p.scalar()); - for (auto&& f : p) - prod->append(1, quaternion_decompose(f), Product::Flatten::No); - ExprPtr r = prod; - if (!has_wrapper) { - // a genuine contraction product — distribute the per-tensor Sums and - // fold dummy-renamed duplicates - expand(r); - rapid_simplify(r); - canonicalize(r); - rapid_simplify(r); - } - return r; - } - throw std::runtime_error("quaternion_decompose: unsupported Expr type"); + return simplify_with_reset(r); } // ===================================================================== @@ -1449,17 +1061,9 @@ std::pair split_re_im(const ExprPtr& expr) { const std::wstring core = conj ? std::wstring{t.label().substr(0, t.label().size() - 1)} : std::wstring{t.label()}; - auto mk = [&](bool imag) { - return ex( - complex_label_add(core, imag), - bra(container::svector{t.bra().begin(), t.bra().end()}), - ket(container::svector{t.ket().begin(), t.ket().end()}), - aux(container::svector{t.aux().begin(), t.aux().end()}), - t.symmetry(), t.braket_symmetry(), t.column_symmetry()); - }; - ExprPtr re = mk(false); - ExprPtr im = mk(true); - if (conj) im = ex(S{-1}) * im; // g* = g~r - i·g~i + ExprPtr re = tensor_with_label(t, complex_label_add(core, /*imag*/ false)); + ExprPtr im = tensor_with_label(t, complex_label_add(core, /*imag*/ true)); + if (conj) im = ex(S{-1}) * im; // g* = g~r − i·g~i return {re, im}; } if (expr->is()) { @@ -1474,11 +1078,10 @@ std::pair split_re_im(const ExprPtr& expr) { } if (expr->is()) { auto const& p = expr->as(); - // Π_i z_i built up incrementally, starting from the leading scalar. auto acc = split_re_im(ex(p.scalar())); for (auto&& f : p) { auto fac = split_re_im(f); - // (ar + i·ai)(fr + i·fi) = (ar·fr - ai·fi) + i(ar·fi + ai·fr) + // (ar + i·ai)(fr + i·fi) = (ar·fr − ai·fi) + i(ar·fi + ai·fr) ExprPtr nr = acc.first * fac.first + ex(S{-1}) * (acc.second * fac.second); ExprPtr ni = acc.first * fac.second + acc.second * fac.first; @@ -1493,19 +1096,10 @@ std::pair split_re_im(const ExprPtr& expr) { /// Distribute and constant-fold a real-tensor scalar expression. /// -/// NB: deliberately does NOT `canonicalize`. The split tensors inherit -/// the (already-consistent) index wiring of the Kramers-traced input; -/// `canonicalize` would permute the indices of `Antisymm` `~r`/`~i` -/// tensors and emit `(-1)` signs that the evaluator cannot see — it -/// rebuilds the `[as]` integral formula from index *type* + Kramers -/// label, which is invariant under such a permutation, so the sign -/// would be uncompensated. `expand` + `rapid_simplify` only distribute -/// and constant-fold; they never permute tensor indices. -ExprPtr clean_real(ExprPtr e) { - expand(e); - rapid_simplify(e); - return e; -} +/// Deliberately does NOT canonicalize (would permute Antisymm `~r`/`~i` +/// indices and emit signs the evaluator cannot see). `non_canon_simplify` +/// only distributes and constant-folds; it never permutes tensor indices. +ExprPtr clean_real(ExprPtr e) { return non_canon_simplify(e); } } // namespace @@ -1520,9 +1114,8 @@ ExprPtr complex_split(const ExprPtr& expr) { clean_real(split_re_im(expr->as().inner()).first)); } if (expr->is()) { - // Im(inner) is itself real-valued; its real-tensor form is the - // `.second` of the split — wrap in RealPart so the evaluator - // boundary treats it as a real quantity. + // Im(inner) is itself real-valued; its real-tensor form is + // `.second` of the split — wrap in RealPart for the evaluator. return ex( clean_real(split_re_im(expr->as().inner()).second)); } @@ -1543,4 +1136,300 @@ ExprPtr complex_split(const ExprPtr& expr) { throw std::runtime_error("complex_split: unsupported top-level Expr type"); } +// ===================================================================== +// V1-canonicalization variant of the open-shell tracer. +// +// Mirrors `Product::canonicalize_impl` (expr.cpp:174-254) but forces +// the all-tensor sub-block through `TensorNetworkV1` rather than the +// hard-coded `using TN = TensorNetwork (= V3)`. This is the minimum +// surface area needed for V1 canonicalization without modifying +// SeQuant core; the rest of the pipeline (rapid_simplify, expand, +// flattening, hashing) is V-agnostic. +// ===================================================================== + +namespace { + +/// Canonicalize one Product through TensorNetworkV1. Mirrors +/// `Product::canonicalize_impl`'s all-tensor branch (expr.cpp:216-254) +/// with the TN typedef forced to V1. Falls back to the input unchanged +/// for Products that contain non-tensor non-scalar factors (e.g. +/// embedded RealPart/ImagPart wrappers — the open-shell pipeline never +/// produces these mid-pass). +ExprPtr canonicalize_v1_product(Product& product, CanonicalizeOptions opts) { + using NamedIndexSet = tensor_network::NamedIndexSet; + + // recursively canonicalize non-tensor subfactors first + for (auto& factor : product.factors()) { + if (factor->is()) continue; + auto bp = factor->canonicalize(opts); + if (bp) { + SEQUANT_ASSERT(bp->is()); + product.scale(std::static_pointer_cast(bp)->value()); + } + } + + // bail if any non-tensor non-scalar factors remain + for (const auto& factor : product.factors()) { + if (!factor->is() && !factor->is_scalar()) return {}; + } + + // pull tensor factors out into their own range for TN construction. + // We ignore scalar factors (Constants embedded as factors): they are + // multiplied into the scalar at the end of this pass. + container::svector tensor_factors; + Constant::scalar_type folded_scalar{1}; + for (const auto& factor : product.factors()) { + if (factor->is()) { + tensor_factors.push_back(factor); + } else if (factor->is()) { + folded_scalar = folded_scalar * factor->as().value(); + } + } + + if (tensor_factors.empty()) return {}; + + TensorNetworkV1 tn(tensor_factors); + + std::shared_ptr named_indices; + if (opts.named_indices) { + named_indices = std::make_shared(opts.named_indices->begin(), + opts.named_indices->end()); + } + + ExprPtr canon_factor = + tn.canonicalize(TensorCanonicalizer::cardinal_tensor_labels(), + /*fast=*/opts.method == CanonicalizationMethod::Rapid, + named_indices.get()); + + // Rebuild factors_ in canonical order. + auto& factors = product.factors(); + factors.clear(); + for (const auto& tptr : tn.tensors()) { + auto exprptr = std::dynamic_pointer_cast(tptr); + SEQUANT_ASSERT(exprptr); + factors.push_back(exprptr); + } + + product.scale(folded_scalar); + if (canon_factor) product.scale(canon_factor->as().value()); + return {}; +} + +/// Recurse Sum → for each summand → if Product, canonicalize via V1; +/// then rebuild the Sum so equivalent summands fold by hash. Mirrors +/// the structure of `Sum::canonicalize_impl` (expr.cpp:374-432) but +/// routes Product canonicalization through V1 rather than V3. +ExprPtr canonicalize_v1_recurse(const ExprPtr& expr, CanonicalizeOptions opts) { + if (expr->is()) { + auto out = std::make_shared(); + for (const auto& s : expr->as().summands()) { + out->append(canonicalize_v1_recurse(s, opts)); + } + return out; + } + if (expr->is()) { + auto product = std::static_pointer_cast(expr->clone()); + canonicalize_v1_product(*product, opts); + return product; + } + // Tensor / Constant / Variable / wrapper nodes: nothing to do. + return expr; +} + +} // namespace + +ExprPtr& canonicalize_v1(ExprPtr& expr, CanonicalizeOptions opts) { + expr = canonicalize_v1_recurse(expr, opts); + return expr; +} + +ExprPtr& simplify_v1(ExprPtr& expr) { + expand(expr); + rapid_simplify(expr); + canonicalize_v1(expr); + rapid_simplify(expr); + return expr; +} + +// ===================================================================== +// open_shell_spintrace_v1: V1-driven clone of +// `mbpt::open_shell_spintrace_impl` (spin.cpp:1367). +// +// Both V1 and V3 enforce bra↔ket-only edge connectivity in their +// graph builders, so an expression with a dummy summed across two +// bras (or two kets) — like the PNO pair density +// `t̄^{ac}_{ij} t̄^{bc}_{ij}` (`c` in two bras) — trips the assertion +// regardless of which TN version we use. The closed-shell tracer +// `mbpt::spintrace` (spin.cpp:1626) avoids the issue by NEVER +// canonicalizing the spin-traced products: it only does `expand` + +// `rapid_simplify`, both of which are TN-agnostic. We mirror that +// strategy here — `simplify_v1` and `canonicalize_v1` still exist as +// useful primitives for callers whose expressions DO have CC-style +// wiring, but in this open-shell pipeline we deliberately skip those +// steps and stay TN-free. +// ===================================================================== + +namespace { + +/// True if a fully-flattened Product is spin-symmetric, i.e. has the +/// concatenated bra and ket lists carrying the same per-slot QN. +/// Mirrors `spin_symm_product` lambda in `open_shell_spintrace_impl`. +bool spin_symm_product_v1(const Product& product) { + container::svector cBra, cKet; + for (auto& term : product) { + if (term->is()) { + auto const& t = term->as(); + cBra.insert(cBra.end(), t.bra().begin(), t.bra().end()); + cKet.insert(cKet.end(), t.ket().begin(), t.ket().end()); + } else if (term->is() || term->is()) { + throw Exception( + "open_shell_spintrace_v1: nested Product/Sum unsupported in " + "spin_symm_product"); + } + } + if (cKet.size() != cBra.size()) return false; + auto i_ket = cKet.begin(); + for (auto& b : cBra) { + if (b.space().qns() != i_ket->space().qns()) return false; + ++i_ket; + } + return true; +} + +} // namespace + +std::vector open_shell_spintrace_v1( + const ExprPtr& expr, + container::svector> ext_index_groups, + std::optional target_spin_case) { + if (expr->is() || expr->is()) + return std::vector{expr}; + + // Internal vs external indices. + container::set grand_idxlist = + get_used_indices>(expr); + container::set ext_idxlist; + for (const auto& grp : ext_index_groups) { + for (const auto& idx : grp) { + Index ix = idx; + ix.reset_tag(); + ext_idxlist.insert(std::move(ix)); + } + } + container::set int_idxlist; + for (auto&& gidx : grand_idxlist) { + if (ext_idxlist.find(gidx) == ext_idxlist.end()) int_idxlist.insert(gidx); + } + + using IndexGroup = container::svector; + container::svector int_index_groups; + for (auto&& i : int_idxlist) int_index_groups.emplace_back(IndexGroup(1, i)); + + SEQUANT_ASSERT(grand_idxlist.size() == + int_idxlist.size() + ext_idxlist.size()); + + auto make_spinspecific = [](const Index& idx, long int spin_bit) { + return spin_bit == 0 ? make_spinalpha(idx) : make_spinbeta(idx); + }; + + // Internal-index replacements: 2^(n_internal_groups) configurations. + auto spin_cases = [&](const container::svector& idx_group) { + const auto ncases = std::uint64_t{1} << idx_group.size(); + container::svector> all(ncases); + for (std::uint64_t i = 0; i < ncases; ++i) { + container::map rep; + for (std::size_t g = 0; g < idx_group.size(); ++g) { + const auto bit = (i >> g) & 1u; + for (auto& idx : idx_group[g]) + rep.emplace(idx, make_spinspecific(idx, bit)); + } + all[i] = rep; + } + return all; + }; + + // External-index replacements: pad the alpha-prefix with k betas for + // k = 0..n_groups (n_groups+1 distinct M_S sectors). + auto ext_spin_cases = [&](const auto& idx_groups) { + container::svector> all; + for (std::size_t i = 0; i <= idx_groups.size(); ++i) { + container::svector spins(idx_groups.size(), 0); + std::fill(spins.end() - i, spins.end(), 1); + container::map rep; + for (std::size_t j = 0; j < idx_groups.size(); ++j) { + for (const auto& idx : idx_groups[j]) + rep.emplace(idx, make_spinspecific(idx, spins[j])); + } + all.push_back(rep); + } + return all; + }; + + auto i_rep = spin_cases(int_index_groups); + auto e_rep = ext_spin_cases(ext_index_groups); + if (target_spin_case) { + auto pick = e_rep.at(*target_spin_case); + e_rep.clear(); + e_rep.push_back(pick); + } + + // Expand antisymmetrizers + flatten — but NO canonicalize. The + // expression may contain bra↔bra (or ket↔ket) summed dummies that + // both TN versions reject; staying TN-free keeps the pipeline + // tolerant. Mirrors `spintrace_impl` (spin.cpp:1626). + auto expanded = expand_A_op(expr); + detail::reset_idx_tags(expanded); + expand(expanded); + rapid_simplify(expanded); + + std::vector result; + for (auto& e : e_rep) { + auto spin_expr = append_spin(expanded, e); + detail::reset_idx_tags(spin_expr); + Sum e_result{}; + for (auto& i : i_rep) { + ExprPtr spin_expr_i = append_spin(spin_expr, i); + spin_expr_i = expand_antisymm(spin_expr_i, /*skip_spinsymm=*/true); + expand(spin_expr_i); + rapid_simplify(spin_expr_i); + detail::reset_idx_tags(spin_expr_i); + Sum i_result{}; + if (spin_expr_i->is() || spin_expr_i->is() || + spin_expr_i->is()) { + e_result.append(spin_expr_i); + } else if (spin_expr_i->is()) { + if (spin_symm_product_v1(spin_expr_i->as())) + e_result.append(spin_expr_i); + } else if (spin_expr_i->is()) { + for (auto& pr : *spin_expr_i) { + if (pr->is()) { + if (spin_symm_product_v1(pr->as())) i_result.append(pr); + } else if (pr->is()) { + if (ms_conserving_columns(pr->as())) i_result.append(pr); + } else if (pr->is() || pr->is()) { + i_result.append(pr); + } else { + throw Exception( + "open_shell_spintrace_v1: unsupported summand type"); + } + } + e_result.append(std::make_shared(i_result)); + } + } + result.push_back(std::make_shared(e_result)); + } + + if (target_spin_case) { + SEQUANT_ASSERT(result.size() == 1 && + "Spin-specific case must return one expression"); + } + + // Final per-spin-block clean-up: rapid_simplify only (no canonicalize). + for (auto& expression : result) { + detail::reset_idx_tags(expression); + rapid_simplify(expression); + } + return result; +} + } // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/spinor.hpp b/SeQuant/domain/mbpt/spinor.hpp index 9069d1bfa8..29ee7e6fc2 100644 --- a/SeQuant/domain/mbpt/spinor.hpp +++ b/SeQuant/domain/mbpt/spinor.hpp @@ -1,23 +1,24 @@ // // Spinor decomposition passes for relativistic 2-component theories. // -// This module provides symbolic transforms over tensor expressions: +// Provides three symbolic transforms over closed-shell relativistic +// tensor expressions: // -// 1. kramers_trace(expr): rewrite a closed-shell expression with -// all-spinor (Kramers::any) indices into a sum over canonical -// Kramers-block representatives (one per orbit under -// {column-swap, Hermitian, time-reversal-with-(-1)^k}), with -// conjugation/sign tracking baked into the per-tensor lookup. +// 1. kramers_trace / kramers_trace_cycles / kramers_trace_burnside +// Rewrite an all-spinor (Kramers::any) expression into a sum over +// canonical Kramers-block representatives. +// 2. fold_conj_pairs + antisymm_recombine +// Post-trace folds: collapse TRS-conjugate Product pairs to +// `2·Re(A)`/`2i·Im(A)`, recombine direct/exchange pairs into +// Kramers-restricted antisymmetric tensors. +// 3. complex_split +// Split each complex Kramers tensor `g = g~r + i·g~i` so the +// result is over real tensors only (`Re(g·t) → g~r·t~r − +// g~i·t~i`). // -// 2. quaternion_decompose(expr): (not yet implemented) further -// decompose Kramers-up indices into the four real (s, x, y, z) -// quaternion components (Helmich-Paris, Repisky, Visscher 2016 -// Eq. 21), aggressively folding via Pauli/quaternion algebra to a -// compact surviving-term set. -// -// Both passes operate on the canonical SeQuant Tensor representation -// and consume its `Kramers` index quantum number (defined in -// SeQuant/domain/mbpt/space_qns.hpp). +// All passes operate on the canonical SeQuant Tensor representation +// and consume the `Kramers` index quantum number defined in +// SeQuant/domain/mbpt/space_qns.hpp. // #ifndef SEQUANT_DOMAIN_MBPT_SPINOR_HPP @@ -30,7 +31,6 @@ #include #include -#include #include #include #include @@ -49,17 +49,12 @@ namespace sequant::mbpt { // ===================================================================== /// @brief Extracts the Kramers quantum number from a QuantumNumbersAttr. -/// @param t the quantum-numbers attribute to read -/// @return the Kramers QN encoded in @p t -/// @pre @p t has at least one Kramers bit set inline Kramers to_kramers(const QuantumNumbersAttr& t) { SEQUANT_ASSERT((t.to_int32() & mask_v) != 0); return static_cast(t.to_int32() & mask_v); } /// @brief Strips the Kramers QN bits from a QuantumNumbersAttr. -/// @param t the quantum-numbers attribute to clear -/// @return @p t with all Kramers bits unset /// @sa spinannotation_remove(const QuantumNumbersAttr&) in spin.hpp inline QuantumNumbersAttr kramers_annotation_remove( const QuantumNumbersAttr& t) { @@ -71,9 +66,6 @@ inline QuantumNumbersAttr kramers_annotation_remove( } /// @brief Strips a trailing ⇑/⇓ Kramers annotation from a label. -/// @tparam WS a wide-string or wide-string-view type -/// @param label the (possibly Kramers-annotated) label -/// @return @p label with any trailing ⇑/⇓ glyph removed template >>> std::wstring kramers_annotation_remove(WS&& label) { @@ -85,10 +77,6 @@ std::wstring kramers_annotation_remove(WS&& label) { } /// @brief Adds a Kramers annotation (⇑ or ⇓) to a label. -/// @tparam WS a wide-string or wide-string-view type -/// @param label the label to annotate -/// @param k the Kramers state; `Kramers::any` returns @p label unchanged -/// @return @p label decorated with the matching ⇑/⇓ glyph /// @pre @p label does not already carry a ⇑/⇓ annotation template >>> @@ -109,10 +97,6 @@ std::wstring kramers_annotation_add(WS&& label, Kramers k) { } /// @brief Replaces any existing ⇑/⇓ annotation on a label with @p k. -/// @tparam WS a wide-string or wide-string-view type -/// @param label the (possibly Kramers-annotated) label -/// @param k the Kramers state to set; `Kramers::any` strips the annotation -/// @return @p label re-annotated with @p k template >>> std::wstring kramers_annotation_replace(WS&& label, Kramers k) { @@ -120,82 +104,45 @@ std::wstring kramers_annotation_replace(WS&& label, Kramers k) { return kramers_annotation_add(label_kf, k); } -/// @brief Constructs a Kramers-up Index from @p idx. -/// @param idx the source index (its type, ordinal and proto-indices are -/// preserved) -/// @return a copy of @p idx with the Kramers-up QN bit set and the label -/// decorated with ⇑ +/// @brief Constructs a Kramers-up Index from @p idx (label decorated with ⇑). [[nodiscard]] Index make_kramers_up(const Index& idx); -/// @brief Constructs a Kramers-down Index from @p idx. -/// @param idx the source index -/// @return a copy of @p idx with the Kramers-down QN bit set and the -/// label decorated with ⇓ +/// @brief Constructs a Kramers-down Index from @p idx (label decorated with ⇓). [[nodiscard]] Index make_kramers_dn(const Index& idx); /// @brief Constructs a Kramers-free (`Kramers::any`) Index from @p idx. -/// @param idx the source index -/// @return a copy of @p idx carrying the `Kramers::any` QN [[nodiscard]] Index make_kramers_free(const Index& idx); -/// @brief Applies an Index→Index replacement map to all tensors in an -/// expression. -/// @param expr the expression to rewrite -/// @param index_replacements the Index→Index substitution map -/// @return a copy of @p expr with every tensor index substituted -/// @sa append_spin in spin.hpp — this is the Kramers analogue, provided -/// as a separate symbol so callers need not pull in the spin tracer +/// @brief Applies an Index→Index replacement map to all tensors in @p expr. +/// +/// @sa append_spin in spin.hpp — append_kramers is a synonym; both are +/// thin wrappers around `Tensor::transform_indices`. ExprPtr append_kramers(const ExprPtr& expr, const container::map& index_replacements); // ===================================================================== -// Complex-arithmetic operators on tensor expressions. -// -// Groundwork for the Kramers tracer (which uses TRS to rewrite barred -// configurations as conjugates of unbarred ones) and the upcoming -// quaternion decomposer (which separates real/imaginary parts). +// Re/Im wrappers + complex-conjugation label convention. // // Convention: complex conjugation of a tensor is encoded as a `*` -// suffix on the tensor's label (e.g., `g` → `g*`, `t` → `t*`). -// Evaluator-side dispatch translates the suffix into a `.conj()` call -// on the underlying numeric tensor. +// suffix on the tensor's label (e.g., `g` → `g*`). Evaluator-side +// dispatch translates the suffix into a `.conj()` call on the +// underlying numeric tensor. +// +// `RealPart`/`ImagPart` are scalar Expr nodes that wrap a contracted +// (closed-form) Product; they are opaque to `simplify`/`canonicalize`, +// so the wrapped Antisymm sub-expressions are not re-permuted. // ===================================================================== /// @brief Tests whether a tensor label encodes complex conjugation. -/// @param label the tensor label to inspect -/// @return true iff @p label carries a trailing `*` suffix inline bool has_conj_suffix(std::wstring_view label) { return !label.empty() && label.back() == L'*'; } /// @brief Toggles the `*` conjugation suffix on a tensor label. -/// @param label the tensor label -/// @return @p label with a `*` appended, or stripped if already present -/// (since `(z*)* = z`) -/// @note pure label rewrite — does not touch Symmetry tags or IndexSpace -/// QNs; the caller must ensure the underlying numeric tensor's -/// complex conjugate is what is wanted [[nodiscard]] std::wstring toggle_conj_suffix(std::wstring_view label); -/// @brief Symbolic complex conjugation of an expression. -/// -/// Recursively walks @p expr and applies: -/// - Tensor: toggles its `*` suffix (`z → z*`, `z* → z`); -/// - Constant: numeric complex conjugate; -/// - Variable: not yet supported (throws); -/// - Product: conjugates the scalar and each factor, `(Π zᵢ)* = Π zᵢ*`; -/// - Sum: conjugates each summand, `(Σ zᵢ)* = Σ zᵢ*`. -/// -/// @param expr the expression to conjugate -/// @return the symbolic complex conjugate of @p expr -/// @note `conjugate(conjugate(expr))` recovers @p expr -[[nodiscard]] ExprPtr conjugate(const ExprPtr& expr); - -/// @brief Symbolic real-part wrapper: `RealPart(E)` represents `Re(E)` -/// for a scalar-valued Expr `E`. -/// @note real-valued itself and idempotent under conjugation; the -/// evaluator extracts @c inner, evaluates it, and takes the real -/// part of the resulting scalar +/// @brief Symbolic real-part wrapper: `RealPart(E)` = `Re(E)` for a +/// scalar-valued Expr `E`. class RealPart : public sequant::Expr { public: RealPart() = delete; @@ -203,26 +150,20 @@ class RealPart : public sequant::Expr { RealPart(RealPart&&) = default; ~RealPart() override = default; - /// @brief Wraps a scalar-valued expression in a real-part marker. - /// @param inner the scalar-valued expression - /// @pre @p inner is non-null and `inner->is_scalar()` + /// @pre @p inner is non-null and scalar-valued + /// + /// We do not enforce `Expr::is_scalar()` because Tensor (an atom) + /// returns false unconditionally and a closed-contraction Product of + /// two Tensors inherits the same answer — the typical inner here. explicit RealPart(ExprPtr inner) : inner_{std::move(inner)} { SEQUANT_ASSERT(inner_); - SEQUANT_ASSERT(inner_->is_scalar()); } - /// @return the wrapped expression const ExprPtr& inner() const { return inner_; } - bool is_scalar() const override { return true; } - type_id_type type_id() const override { return get_type_id(); } - ExprPtr clone() const override { return ex(inner_->clone()); } - - /// @brief Adjoint of `Re(E)` — a no-op, since `Re(E)` is real. - void adjoint() override {} - + void adjoint() override {} // Re(E) is real, self-adjoint std::wstring to_latex() const override { return L"\\Re\\left[" + inner_->to_latex() + L"\\right]"; } @@ -233,25 +174,21 @@ class RealPart : public sequant::Expr { hash_type memoizing_hash() const override { auto compute = [this]() { auto v = hash::value(*inner_); - hash::combine(v, std::size_t{0xC0FFEE01ull}); // RealPart tag + hash::combine(v, std::size_t{0xC0FFEE01ull}); return v; }; if (!hash_value_) hash_value_ = compute(); return *hash_value_; } - bool static_equal(const Expr& that) const override { return *inner_ == *static_cast(that).inner_; } - bool static_less_than(const Expr& that) const override { return *inner_ < *static_cast(that).inner_; } }; -/// @brief Symbolic imaginary-part wrapper: `ImagPart(E)` represents -/// `Im(E)` for a scalar-valued Expr `E`. -/// @note real-valued and idempotent under conjugation +/// @brief Symbolic imaginary-part wrapper: `ImagPart(E)` = `Im(E)`. class ImagPart : public sequant::Expr { public: ImagPart() = delete; @@ -259,26 +196,17 @@ class ImagPart : public sequant::Expr { ImagPart(ImagPart&&) = default; ~ImagPart() override = default; - /// @brief Wraps a scalar-valued expression in an imaginary-part marker. - /// @param inner the scalar-valued expression - /// @pre @p inner is non-null and `inner->is_scalar()` + /// @pre @p inner is non-null and scalar-valued (see RealPart for why + /// we don't use `Expr::is_scalar()` as the precondition) explicit ImagPart(ExprPtr inner) : inner_{std::move(inner)} { SEQUANT_ASSERT(inner_); - SEQUANT_ASSERT(inner_->is_scalar()); } - /// @return the wrapped expression const ExprPtr& inner() const { return inner_; } - bool is_scalar() const override { return true; } - type_id_type type_id() const override { return get_type_id(); } - ExprPtr clone() const override { return ex(inner_->clone()); } - - /// @brief Adjoint of `Im(E)` — a no-op, since `Im(E)` is real. - void adjoint() override {} - + void adjoint() override {} // Im(E) is real, self-adjoint std::wstring to_latex() const override { return L"\\Im\\left[" + inner_->to_latex() + L"\\right]"; } @@ -289,403 +217,128 @@ class ImagPart : public sequant::Expr { hash_type memoizing_hash() const override { auto compute = [this]() { auto v = hash::value(*inner_); - hash::combine(v, std::size_t{0xC0FFEE02ull}); // ImagPart tag + hash::combine(v, std::size_t{0xC0FFEE02ull}); return v; }; if (!hash_value_) hash_value_ = compute(); return *hash_value_; } - bool static_equal(const Expr& that) const override { return *inner_ == *static_cast(that).inner_; } - bool static_less_than(const Expr& that) const override { return *inner_ < *static_cast(that).inner_; } }; -/// @brief Constructs `Re(expr)` as a RealPart wrapper. -/// @param expr the scalar-valued expression to wrap -/// @return a RealPart wrapping @p expr -[[nodiscard]] inline ExprPtr re(ExprPtr expr) { +/// @brief Wraps @p expr as `Re(expr)`. +[[nodiscard]] inline ExprPtr real_part(ExprPtr expr) { return ex(std::move(expr)); } -/// @brief Constructs `Im(expr)` as an ImagPart wrapper. -/// @param expr the scalar-valued expression to wrap -/// @return an ImagPart wrapping @p expr -[[nodiscard]] inline ExprPtr im(ExprPtr expr) { +/// @brief Wraps @p expr as `Im(expr)`. +[[nodiscard]] inline ExprPtr imaginary_part(ExprPtr expr) { return ex(std::move(expr)); } -/// @brief Symbolic real part of an expression. -/// @param expr the scalar-valued expression -/// @return `Re(expr)` as a RealPart wrapper Expr -/// @note replaces the older `(expr + conj(expr))/2` Sum form so callers -/// can detect RealPart structurally and dispatch real-part -/// evaluation -[[nodiscard]] inline ExprPtr real_part(const ExprPtr& expr) { - return ex(expr); -} - -/// @brief Symbolic imaginary part of an expression. -/// @param expr the scalar-valued expression -/// @return `Im(expr)` as an ImagPart wrapper Expr -[[nodiscard]] inline ExprPtr imaginary_part(const ExprPtr& expr) { - return ex(expr); -} - // ===================================================================== -// Kramers tracer. +// Kramers tracers. +// +// Three coexisting implementations, all returning expressions whose +// scalar value equals the full 2^n per-Kramers-index trace of @p expr. +// Provided in parallel for cross-validation; pick whichever is fastest +// in production. +// +// All assume @p expr is fully contracted (every barred index appears +// an even number of times). CC residuals carry external bars and need +// `(-1)^(external bars)` sign bookkeeping that is not yet implemented. // ===================================================================== -/// @brief Rewrites an all-spinor tensor expression into a sum over -/// canonical Kramers-block representatives. -/// -/// Implements the closed-shell relativistic trace algebra: enumerate the -/// per-tensor Kramers blocks, apply per-tensor TRS canonicalization (with -/// the `(-1)^k` sign), and fold via the standard -/// `canonicalize` + `rapid_simplify` pipeline. Antisymmetric tensors are -/// expanded internally before tracing. -/// -/// @param expr expression whose indices carry `Kramers::any` (i.e. not -/// yet specialized to up/down) -/// @return a `Sum` of canonical-orbit-block tensors; tensor -/// labels carrying conjugation get a `*` suffix -/// @todo CC support. The per-tensor `(-1)^k` TRS signs only cancel when -/// every barred index appears an even number of times across the -/// whole expression — true for a fully-contracted scalar (the MP2 -/// energy), but NOT for an expression carrying free/external -/// indices (a CC residual `R_{IJ...}^{AB...}`). Residuals need -/// explicit `(-1)^(bars on external indices)` sign bookkeeping; -/// until that is added this entry point is correct only for -/// fully-contracted scalar expressions. The same caveat applies to -/// `kramers_trace_cycles` and `kramers_trace_burnside`. +/// @brief 2^n exhaustive enumeration over Kramers configurations of +/// every Kramers::any index in @p expr. ExprPtr kramers_trace(const ExprPtr& expr); -/// @brief Flips every Kramers state (`up` ↔ `down`) in an expression. +/// @brief ResultExpr overload: kramers-trace a full equation +/// `D^{...}_{...} = RHS`. /// -/// Walks every tensor in @p expr, collects each Kramers-tagged index, -/// builds an Index→Index replacement that swaps the Kramers QN bit (and -/// the ⇑/⇓ label glyph), then applies the swap. Indices carrying -/// `Kramers::any` are left untouched. +/// Enumerates 2^N external Kramers configurations; per configuration, +/// substitutes the externals on the RHS, then runs `kramers_trace` + +/// `antisymm_recombine` on the result. Two folding levels are +/// supported via @p fold_klein: +/// - **true** (default): emit only the canonical representative per +/// orbit under the Klein 4-group {e, α (bra-flip), β (ket-flip), +/// αβ (full TRS)}. For an N-external Kramers-restricted result +/// this reduces 2^N raw blocks by a factor of 4 (typically). +/// - **false**: emit only canonical representatives under full TRS +/// alone (cfg ≤ ~cfg). For partial-Kramers-flip equivalence +/// verification — emits the 2^(N-1) TRS-canonical configurations +/// so a caller can pair them into Klein orbits and numerically +/// check the partial-flip equivalence. +/// +/// @return one ResultExpr per non-vanishing canonical configuration. +[[nodiscard]] container::svector kramers_trace( + const ResultExpr& expr, bool fold_klein = true); + +/// @brief Cycle-driven multiset enumeration: decompose the +/// contraction graph into cycles, group cycles by shape, and +/// enumerate Kramers-labeling multisets per group with +/// multinomial multiplicity. +[[nodiscard]] ExprPtr kramers_trace_cycles(const ExprPtr& expr); + +/// @brief Burnside-orbit enumeration over the index-permutation +/// symmetry group induced by the antisymmetrizer structure. +[[nodiscard]] ExprPtr kramers_trace_burnside(const ExprPtr& expr); + +/// @brief Flips every specialized Kramers state (`up` ↔ `down`) in @p expr. /// -/// @param expr the expression to flip -/// @return a copy of @p expr with every specialized Kramers index flipped -/// @sa fold_conj_pairs — uses this to find TRS-conjugate-pair Products +/// `Kramers::any` indices are left untouched. [[nodiscard]] ExprPtr flip_kramers(const ExprPtr& expr); /// @brief Folds TRS conjugate-pair Products in a Sum. /// -/// For each summand `A`, builds `A_flipped` via flip_kramers() (the -/// full-bar TRS partner), canonicalizes it, and checks structural -/// equality against the remaining summands: -/// - `A == A_flipped` (self-conjugate after canonicalize): emits -/// `Re(A)`, since the term is necessarily real-valued; -/// - `A_flipped == B` for some unused `B`: emits `2·Re(A)` and marks -/// both `A` and `B` consumed (term count drops by one per fold); +/// For each summand `A`, builds `A_flipped = flip_kramers(A)`, +/// canonicalizes it, and looks for a structural match among the +/// remaining summands: +/// - `A == A_flipped` (self-conjugate): emits `Re(A)`; +/// - matching `B` with same scalar: emits `2·Re(A)`, drops B; +/// - matching `B` with opposite scalar: emits `2i·Im(A)`, drops B; /// - no match: keeps `A` unchanged. /// -/// Imaginary-pair detection (`A - A_flipped → 2i·Im(A)`) is supported -/// when the partner `B` carries the opposite scalar prefactor. -/// -/// @param expr a `Sum` to fold (typically the output of kramers_trace()) -/// @return a new `Sum`, typically with fewer summands, whose scalar -/// value equals that of @p expr +/// O(n) hash-bucketed; equivalent to pairwise search but much cheaper. [[nodiscard]] ExprPtr fold_conj_pairs(const ExprPtr& expr); -// ===================================================================== -// Cycle-based Kramers enumeration. -// -// The closed-shell spin-trace algorithm in SeQuant decomposes a Product's -// contraction graph into cycles and exploits cycle structure to count -// independent traces. We adapt the same idea for Kramers: each cycle in -// the contraction graph carries its own Kramers labelling, and the -// number of *canonical* Kramers patterns per cycle (under cyclic rotation -// + reflection) is typically << 2^L for length-L cycles — a long cycle -// visiting d distinct indices twice each has 2^d raw labellings that -// collapse to a much smaller canonical-pattern set. -// ===================================================================== - -/// @brief One contraction cycle in a Product's index graph. -/// -/// The cycle alternates intra-tensor edges (bra[k] ↔ ket[k] of the same -/// tensor) and inter-tensor edges (the same Index on two tensors). -/// `nodes` walks the cycle in visiting order. -struct ContractionCycle { - /// @brief One step of the cycle walk. - struct Node { - std::size_t tensor_idx; ///< which factor of the Product - std::size_t braket_pos; ///< position in the tensor's const_braket() - Index idx; ///< index occupying that slot - }; - container::svector nodes; ///< the cycle walk, in visiting order -}; - -/// @brief Decomposes a Product's index-contraction structure into cycles. -/// @param product a fully-contracted scalar Product -/// @return the contraction cycles of @p product -/// @pre every Index in @p product appears in exactly two `(tensor, slot)` -/// positions; each tensor's bra[k]/ket[k] positions are paired as -/// intra-tensor edges -[[nodiscard]] container::svector kramers_cycles( - const sequant::Product& product); - -/// @brief Diagnostic dump of the cycle structure of an expression. -/// -/// Walks @p expr (a `Sum` of Products), decomposes each Product into -/// cycles, and prints each cycle's length, distinct indices visited, -/// per-cycle Kramers labelling, and the number of canonical Kramers -/// patterns under cyclic-rotation symmetry. -/// -/// @param expr the expression to inspect -/// @param os the wide-character stream to write to -void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os); - -/// @brief Canonical Kramers label of a cycle. -/// @param c the contraction cycle -/// @return the lex-smallest rotation of @p c 's raw label string (e.g. -/// `BUUB` and `UBBU` both canonicalize to `BBUU`) -[[nodiscard]] std::wstring cycle_canonical_label(const ContractionCycle& c); - -/// @brief Cycle-canonical signature of a Product. -/// @param product the Product to fingerprint -/// @return the sorted multiset of per-cycle canonical Kramers labels, -/// joined into one key string (e.g. `BBUU|UUUU`) -/// @note Products with the same signature are equivalent under per-cycle -/// cyclic-rotation symmetry -[[nodiscard]] std::wstring cycle_canonical_signature( - const sequant::Product& product); - -/// @brief Buckets Products by cycle-canonical signature and sums each -/// bucket into one representative. -/// -/// Within a bucket the scalar prefactors are summed and one -/// representative is emitted. Acts independently of fold_conj_pairs() -/// and on a raw `Sum`-of-Products without any Re/Im wrappers. -/// -/// @param expr a `Sum` to fold -/// @return a `Sum` with one summand per cycle-canonical signature -/// @warning NOT value-preserving in general — the cyclic-rotation -/// signature is not a complete contraction invariant, so two -/// genuinely different contractions can share it. Superseded by -/// kramers_trace_cycles(); kept only for diagnostic comparison. -[[nodiscard]] ExprPtr cycle_canonical_fold(const ExprPtr& expr); - -// ===================================================================== -// Standalone cycle-driven Kramers tracer. -// -// Unlike `kramers_trace` (which enumerates all 2^n per-index Kramers -// configurations then folds post-hoc), this tracer is *driven* by the -// contraction-graph cycle structure — analogous to how -// `closed_shell_spintrace` is driven by `count_cycles`. -// -// Algorithm (cycles decomposed on the ANTISYMMETRIZED form, before -// expand_antisymm): -// 1. Decompose the input Product's contraction graph into cycles. -// 2. Group cycles into structural-equivalence classes: cycles that map -// to one another under the antisymmetrizer's index permutations -// (structurally-identical cycles are one class — the antisymmetric -// bra/ket index permutations swap them). -// 3. Enumerate Kramers labelings as MULTISETS within each class (not -// ordered assignments) — this is where the reduction over naive -// per-index 2^n enumeration comes from: equivalent cycles with -// swapped labels are the same contraction. -// 4. For each multiset configuration: substitute Kramers labels, -// expand_antisymm, canonicalize; attach the multiset multiplicity. -// 5. Collect into a Sum; the caller may apply `fold_conj_pairs` for -// additional TRS conjugate-pair folding. -// -// Correct-by-construction: works from contraction topology, never guesses -// equivalences. Coexists with `kramers_trace` so the two can be -// cross-validated. -// ===================================================================== - -/// @brief Cycle-driven Kramers tracer. -/// -/// See the block comment above for the algorithm. -/// -/// @param expr a closed-shell all-spinor expression with `Kramers::any` -/// indices and (optionally) antisymmetric tensors -/// @return a `Sum` of canonical Kramers-block contractions whose -/// scalar value equals the full 2^n per-index trace -[[nodiscard]] ExprPtr kramers_trace_cycles(const ExprPtr& expr); - /// @brief Recombines direct/exchange Product pairs into antisymmetric /// tensors — the inverse of `expand_antisymm`. /// -/// After `kramers_trace` + `fold_conj_pairs` the expression is a `Sum` of -/// fully-expanded NonSymm contractions. Many of these are the `direct` -/// and `exchange` halves of a single antisymmetrized contraction: two -/// summands recombine into one when -/// - one tensor is structurally identical between them, and -/// - the other tensor differs by a column permutation, with the -/// summand scalars in the ratio `sign(permutation)`. -/// Then `c·X·Y_direct + (sign·c)·X·Y_exchange → c·X·Ȳ`, where `Ȳ` is -/// `Y` re-tagged `Symmetry::Antisymm` over the Kramers-restricted index -/// space (NOT full spinor). -/// -/// Applied iteratively to a fixed point, so a second pass over the -/// once-recombined terms catches doubly-antisymmetrized blocks. Re/Im -/// wrappers from fold_conj_pairs() are preserved — recombination acts on -/// the wrapped inner Product. +/// Iterates to a fixed point. Recombines only when one tensor matches +/// structurally between two summands, the other differs by a single +/// column transposition, and the summand scalars are negatives. Re/Im +/// wrappers from `fold_conj_pairs` are preserved. /// -/// @param expr a `Sum` of (optionally Re/Im-wrapped) NonSymm contractions -/// @return a `Sum` in minimal mixed antisymm/non-antisymm form, with the -/// same scalar value as @p expr -/// @note value-preserving by construction — it only re-groups terms that -/// sum to the same scalar; the result is the analogue of what the -/// open-shell spin tracer yields, but driven by index wiring -/// rather than spin-label uniformity -/// @todo CC support. Recombination currently recognizes only rank-(2,2) -/// column swaps (see `detect_single_column_swap` in spinor.cpp). -/// CCSDT triples residuals (`Ȳ_{IJK}^{ABC}` etc.) need the -/// single-column-swap detection generalized to rank-(n,n). -/// Higher-rank tensors are left expanded — not folded, never -/// miscombined — so this is a missed-optimization, not a bug. +/// @todo CC support requires generalizing `detect_single_column_swap` +/// past rank-(2,2) for CCSDT triples. [[nodiscard]] ExprPtr antisymm_recombine(const ExprPtr& expr); -/// @brief Burnside-enumeration Kramers tracer. -/// -/// Like `kramers_trace`, but instead of walking all `2^n` per-index -/// Kramers configurations it enumerates only the orbit representatives -/// under the index-permutation symmetry group induced by the -/// antisymmetrizer structure of the input (the group generated by -/// transpositions within each tensor's antisymmetric bra/ket groups). -/// Each representative is emitted once, scaled by its orbit size; group -/// elements act as bit-permutations on the n-bit config and orbits are -/// found by bit-ops. -/// -/// A transposition `(p,q)` is admitted as a group generator iff in every -/// tensor containing both indices they sit in the same antisymmetric -/// bra/ket group and the accumulated permutation sign is `+1`. -/// -/// @param expr a closed-shell all-spinor expression with `Kramers::any` -/// indices -/// @return a `Sum` value-equivalent (after `canonicalize`) to -/// the output of kramers_trace(); feed it the same -/// `fold_conj_pairs → antisymm_recombine` pipeline -[[nodiscard]] ExprPtr kramers_trace_burnside(const ExprPtr& expr); - -// ===================================================================== -// Quaternion decomposer. -// -// After Kramers tracing, every 2-electron tensor is rewritten into the -// four real (s, x, y, z) quaternion components of its spinor -// coefficients (Helmich-Paris, Repisky, Visscher 2016; see -// kramers-restriction §6 / Wiebke Eq. 8-11). For a rank-(2,2) tensor -// `T(B0,B1;K0,K1)` with Kramers-specialized indices, the spin-traced -// per-electron transition density decomposes via a 4x4 complex matrix -// `M(k_bra,k_ket)` — one per Kramers-pair, 8 nonzero entries each: -// -// T(B0,B1;K0,K1) = Σ_c M(k_B0,k_K0)[c_B0,c_K0]·M(k_B1,k_K1)[c_B1,c_K1] -// · T#(B0,B1;K0,K1) -// -// where electron 1 = (B0,K0), electron 2 = (B1,K1). The quaternion -// component tuple is encoded as a `#`-prefixed 4-char suffix on the -// tensor label (analogous to the `*` conjugation marker); the contracted -// indices are left unchanged so tensors still contract pairwise. -// ===================================================================== - -/// @brief Encodes a quaternion-component tuple onto a tensor label. -/// @param core the (component-free) tensor label -/// @param c0,c1,c2,c3 the components of bra[0],bra[1],ket[0],ket[1], -/// each one of 's','x','y','z' -/// @return @p core with a `#c0c1c2c3` suffix appended -[[nodiscard]] inline std::wstring quaternion_label_add(std::wstring_view core, - char c0, char c1, - char c2, char c3) { - std::wstring s{core}; - s += L'#'; - s += static_cast(c0); - s += static_cast(c1); - s += static_cast(c2); - s += static_cast(c3); - return s; -} - -/// @brief Parses the quaternion-component suffix off a tensor label. -/// @param label the (possibly quaternion-decomposed) tensor label; a -/// trailing `*` conjugation marker is tolerated -/// @return the four component chars (bra[0],bra[1],ket[0],ket[1]), or -/// std::nullopt if @p label carries no `#`-suffix -[[nodiscard]] inline std::optional> quaternion_label_parse( - std::wstring_view label) { - const auto pos = label.find(L'#'); - if (pos == std::wstring_view::npos) return std::nullopt; - auto tail = label.substr(pos + 1); - if (!tail.empty() && tail.back() == L'*') - tail = tail.substr(0, tail.size() - 1); - if (tail.size() != 4) return std::nullopt; - std::array out{}; - for (std::size_t i = 0; i < 4; ++i) { - const wchar_t ch = tail[i]; - if (ch != L's' && ch != L'x' && ch != L'y' && ch != L'z') - return std::nullopt; - out[i] = static_cast(ch); - } - return out; -} - -/// @brief Returns the component-free core of a (possibly quaternion- -/// decomposed, possibly conjugated) tensor label. -/// @param label the tensor label -/// @return @p label truncated at the `#` suffix marker (no-op if absent) -[[nodiscard]] inline std::wstring quaternion_label_core( - std::wstring_view label) { - return std::wstring{label.substr(0, label.find(L'#'))}; -} - -/// @brief Rewrites a Kramers-traced expression into quaternion components. -/// -/// Recurses through Sum / Product / RealPart / ImagPart, expanding any -/// `Symmetry::Antisymm` tensor first, then applies the per-tensor 4x4 -/// `M`-matrix decomposition (see the section comment above) to every -/// rank-(2,2) `Symmetry::Nonsymm` tensor. Each input tensor expands into -/// up to 64 quaternion-component tensors (8 nonzero `M` entries per -/// electron pair); the decomposed sub-products are canonicalized so -/// dummy-renamed duplicates fold. -/// -/// @param expr a Kramers-traced expression (output of `kramers_trace` → -/// `fold_conj_pairs` → `antisymm_recombine`), every tensor -/// index carrying a specialized `Kramers::up`/`Kramers::down` -/// QN -/// @return the quaternion-component-decomposed expression, value-equal to -/// @p expr; tensor labels carry a `#c0c1c2c3` component suffix -/// @pre every tensor in @p expr is rank-(2,2) with Kramers-specialized -/// indices -[[nodiscard]] ExprPtr quaternion_decompose(const ExprPtr& expr); - // ===================================================================== // Complex (Re/Im) split. // -// For a closed-shell 1eX2C energy expression the spin-free Coulomb -// operator spin-traces each per-electron transition density down to its -// scalar Pauli component, which is *complex* — no quaternion -// (4-component) structure survives, so `quaternion_decompose` -// over-decomposes. The minimal real-arithmetic form is the 2-component -// split: every complex tensor `g` becomes `g = g~r + i·g~i` with `g~r`, -// `g~i` real tensors, and complex arithmetic is propagated symbolically -// (`Re(g·t) = g~r·t~r − g~i·t~i`, etc.). The component is encoded as a -// `~r`/`~i` label suffix (analogous to the `*` conjugation marker). +// For closed-shell 1eX2C the Coulomb operator spin-traces each +// per-electron transition density to a *complex* scalar; no quaternion +// (4-component) structure survives. The minimal real-arithmetic form +// is the 2-component split: `g → g~r + i·g~i`, both real, propagated +// symbolically (`Re(g·t) = g~r·t~r − g~i·t~i`). Components are encoded +// as a `~r`/`~i` label suffix. // ===================================================================== /// @brief Encodes a real/imaginary-part marker onto a tensor label. -/// @param core the (marker-free) tensor label -/// @param imag true for the imaginary part (`~i`), false for the real -/// part (`~r`) -/// @return @p core with the `~r`/`~i` suffix appended [[nodiscard]] inline std::wstring complex_label_add(std::wstring_view core, bool imag) { return std::wstring{core} + (imag ? L"~i" : L"~r"); } /// @brief Parses the real/imaginary-part marker off a tensor label. -/// @param label the (possibly complex-split) tensor label; a trailing -/// `*` conjugation marker is tolerated -/// @return true for the imaginary part, false for the real part, or -/// std::nullopt if @p label carries no `~` marker +/// @return true for `~i`, false for `~r`, std::nullopt if no marker +/// (a trailing `*` conjugation marker is tolerated). [[nodiscard]] inline std::optional complex_label_parse( std::wstring_view label) { const auto pos = label.find(L'~'); @@ -700,35 +353,74 @@ void kramers_cycle_dump(const ExprPtr& expr, std::wostream& os); /// @brief Returns the marker-free core of a (possibly complex-split) /// tensor label. -/// @param label the tensor label -/// @return @p label truncated at the `~` marker (no-op if absent) [[nodiscard]] inline std::wstring complex_label_core(std::wstring_view label) { return std::wstring{label.substr(0, label.find(L'~'))}; } -/// @brief Rewrites a Kramers-traced expression into one over real -/// tensors via the complex (Re/Im) split. +/// @brief Splits each complex tensor `g = g~r + i·g~i` and propagates +/// complex arithmetic symbolically. /// -/// Each complex tensor `g` is replaced by two real tensors `g~r`, `g~i` -/// (`g = g~r + i·g~i`), and complex arithmetic is propagated -/// symbolically through every Product and Sum. A `RealPart`/`ImagPart` -/// wrapper collapses to the corresponding real sub-expression -/// (`Re(g·t) → g~r·t~r − g~i·t~i`). Conjugated tensors (`g*`) split as -/// `g~r − i·g~i`. The resulting real sub-expressions are distributed and -/// constant-folded but deliberately NOT canonicalized — the split -/// tensors inherit the input's index wiring, and canonicalizing -/// `Antisymm` `~r`/`~i` tensors would emit signs the evaluator cannot -/// see (it rebuilds `[as]` formulas from index type + Kramers label). +/// Conjugated tensors (`g*`) split as `g~r − i·g~i`. Real +/// sub-expressions are distributed and constant-folded but +/// deliberately NOT canonicalized — the split tensors inherit the +/// input's index wiring, and canonicalizing `Antisymm` `~r`/`~i` +/// tensors would emit signs the evaluator cannot see. /// -/// @param expr a Kramers-traced expression whose summands are -/// `Constant * RealPart/ImagPart(...)` (the output of -/// `kramers_trace → fold_conj_pairs → antisymm_recombine`) -/// @return an equivalent expression in which every tensor is real; -/// tensor labels carry a `~r`/`~i` suffix, real sub-expressions -/// re-wrapped in `RealPart` for the evaluator boundary -/// @pre every summand is real-wrapped — run `fold_conj_pairs` first +/// @pre every summand of @p expr is wrapped in `RealPart`/`ImagPart` +/// (i.e. run `fold_conj_pairs` first) [[nodiscard]] ExprPtr complex_split(const ExprPtr& expr); +// ===================================================================== +// V1-based open-shell spin tracer. +// +// SeQuant's default `mbpt::open_shell_spintrace` runs the +// canonicalization pipeline through `TensorNetworkV3`, which assumes +// a CC-style index wiring (every contracted dummy appears once in a +// bra and once in a ket). Expressions in which a dummy appears across +// two bras (or two kets) — e.g. the PNO pair density +// `D^{ij}_{ab} = t̄^{ac}_{ij} t̄^{bc}_{ij}`, where `c` is summed across +// two bras — trip V3's graph-build assertion on a worker thread, +// bypassing user-side `try`/`catch`. +// +// `open_shell_spintrace_v1` is a drop-in V1-canonicalization variant +// of `open_shell_spintrace_impl` (see `spin.cpp:1367`). The algorithm +// is identical; only the canonicalization back-end is swapped, which +// V1's permissive graph builder accepts. +// +// `canonicalize_v1` and `simplify_v1` are the underlying primitives. +// ===================================================================== + +/// @brief Like `sequant::canonicalize`, but forces every all-tensor +/// Product subterm through `TensorNetworkV1`. +/// +/// Use when an expression's contraction wiring trips +/// `TensorNetworkV3`'s graph-build assumptions (typically: a dummy +/// summed across two bras or two kets). +ExprPtr& canonicalize_v1(ExprPtr& expr, + CanonicalizeOptions opts = CanonicalizeOptions{ + .method = CanonicalizationMethod::Lexicographic}); + +/// @brief Like `sequant::simplify`, but uses `canonicalize_v1`. +ExprPtr& simplify_v1(ExprPtr& expr); + +/// @brief V1-canonicalization variant of `mbpt::open_shell_spintrace`. +/// +/// Same algorithm as `open_shell_spintrace_impl` in `spin.cpp:1367`; +/// only difference is that `simplify`/`canonicalize` are routed +/// through `simplify_v1`/`canonicalize_v1` so the pipeline tolerates +/// non-CC index wiring (open dummies across two bras / two kets). +/// +/// @param expr the spinorbital expression to trace +/// @param ext_index_groups groups of external (free) indices; one +/// per particle (e.g. `{{i,a}, {j,b}}` for a rank-(2,2) RHS) +/// @param target_spin_case if set, returns only the αα...β...β +/// configuration with @p target_spin_case beta-tagged +/// externals; otherwise returns all such configurations +[[nodiscard]] std::vector open_shell_spintrace_v1( + const ExprPtr& expr, + container::svector> ext_index_groups, + std::optional target_spin_case = std::nullopt); + } // namespace sequant::mbpt #endif // SEQUANT_DOMAIN_MBPT_SPINOR_HPP