From 0fee3e144a761a706f5de8d48e50d3615325fa70 Mon Sep 17 00:00:00 2001 From: Ajay Date: Fri, 28 Aug 2026 14:09:48 -0400 Subject: [PATCH 1/5] cc: add selectable UCC EOM assembly --- SeQuant/domain/mbpt/models/cc.cpp | 103 ++++++++++++++++-------------- SeQuant/domain/mbpt/models/cc.hpp | 48 ++++++++------ doc/user/guide/cc.rst | 9 ++- tests/unit/test_mbpt_cc.cpp | 42 +++++++++++- 4 files changed, 133 insertions(+), 69 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 144a082731..2af56c7585 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -411,34 +411,46 @@ container::svector> eom_manifolds(nₚ np, } } // namespace -// Per-block-truncated EOM sigma equations. For the qUCCSD ranks see +// UCC EOM sigma equations. For the qUCCSD ranks see // 10.1063/5.0062090 Sec. II C, Eqs. (29)-(48); for IP/EA, // 10.1021/acs.jctc.5c01991 Fig. 1. // // Eq. (10) splits the single physical E_gr from the normal-ordered components -// that build the blocks of Eq. (7). Each cumulative H̄^(k) used below still -// carries its rank-dependent scalar part, so remove that same scalar on the -// diagonal before assembling the separately truncated components. -std::vector CC::eom_r_blocked( - nₚ np, nₕ nh, const std::vector& block_ranks) const { - if (!unitary()) - throw Exception("CC::eom_r: block_ranks require a unitary ansatz"); +// that build the blocks of Eq. (7). In projected-H̄ UCC assembly, each +// cumulative H̄^(k) still carries its rank-dependent scalar part, so remove +// that same scalar on the diagonal. +std::vector CC::assemble_ucc_eom( + nₚ np, nₕ nh, const std::vector& block_ranks, + UCCEOMAssembly assembly) const { + if (assembly == UCCEOMAssembly::Commutator && + hbar_expansion_ == HbarExpansion::Bernoulli) + throw Exception("CC::eom_r: Bernoulli requires projected Hbar assembly"); + + using std::min; + if (assembly == UCCEOMAssembly::Commutator && block_ranks.empty()) { + const auto hbar_R = commutator(hbar(), R(np, nh, eom_norm)); + std::vector result(min(np, nh) + 1); + for (const auto& [rp, rh] : eom_manifolds(np, nh)) + result.at(min(rp, rh)) = ref_av(δl(nₚ(rp), nₕ(rh)) * hbar_R, {}); + return result; + } const auto manifolds = eom_manifolds(np, nh); const auto K = manifolds.size(); // `block_ranks` is read at i * K + j, so the ascending order above is what - // makes `{2,1,1,0}` mean SS, SD, DS, DD; empty means uniform hbar_comm_rank + // makes `{2,1,1,0}` mean SS, SD, DS, DD; empty means the configured H̄ rank, + // or the fourth commutator when no rank is configured. const std::vector ranks = - block_ranks.empty() ? std::vector(K * K, hbar_comm_rank().value()) - : block_ranks; + block_ranks.empty() + ? std::vector(K * K, hbar_comm_rank().value_or(4)) + : block_ranks; if (ranks.size() != K * K) throw Exception( "CC::eom_r: block_ranks must be a K x K row-major matrix, " "K = number of projection manifolds"); // Bernoulli H̄ is tensor-level, BCH H̄ operator-level; the bra/ket/vev trio - // below must match it. Connectivity is empty either way, as everywhere on the - // unitary path; only the operator-level vev forwards screen/use_topology. + // below must match it. UCC uses no connectivity. const bool tensor_level = hbar_expansion_ == HbarExpansion::Bernoulli; // One H̄ per distinct truncation order, reduced to its R part (Bernoulli only; @@ -469,7 +481,6 @@ std::vector CC::eom_r_blocked( .use_topology = use_topology_}); }; - using std::min; std::vector result(min(np, nh) + 1); for (size_t i = 0; i < K; ++i) { const auto [bp, bh] = manifolds[i]; @@ -479,13 +490,17 @@ std::vector CC::eom_r_blocked( const auto [kp, kh] = manifolds[j]; const auto& hbar_ij = hbars.at(ranks.at(i * K + j)); const auto ket = ket_of(kp, kh); - acc->append(vev(bra * hbar_ij * ket)); - // Remove the scalar part of this block's temporary H̄^(k_ii). This is - // not a block-dependent physical E_gr: it leaves the normal-ordered - // coefficients selected for this block in Eq. (10). Written as - // so Wick keeps its summed indices disjoint from the block's - // external ones. - if (i == j) acc->append(ex(-1) * vev(bra * ket * hbar_ij)); + if (assembly == UCCEOMAssembly::Commutator) { + acc->append(vev(bra * commutator(hbar_ij, ket))); + } else { + acc->append(vev(bra * hbar_ij * ket)); + // Remove the scalar part of this block's temporary H̄^(k_ii). This is + // not a block-dependent physical E_gr: it leaves the normal-ordered + // coefficients selected for this block in Eq. (10). Written as + // so Wick keeps its summed indices disjoint from the + // block's external ones. + if (i == j) acc->append(ex(-1) * vev(bra * ket * hbar_ij)); + } } result.at(static_cast(min(bp, bh))) = simplify(ExprPtr{acc}); } @@ -493,44 +508,36 @@ std::vector CC::eom_r_blocked( } std::vector CC::eom_r(nₚ np, nₕ nh, - const std::vector& block_ranks) const { + const std::vector& block_ranks, + std::optional assembly) const { SEQUANT_ASSERT(np > 0 || nh > 0, "Unsupported excitation order"); if (np != nh) SEQUANT_ASSERT( get_default_context().spbasis() != SPBasis::Spinfree, "spin-free basis does not yet support non particle-conserving cases"); - // Bernoulli always takes the blocked path: the uniform one below commutes H̄ - // with an operator-level R, which a tensor-level H̄ cannot take part in. - if (!block_ranks.empty() || hbar_expansion_ == HbarExpansion::Bernoulli) - return eom_r_blocked(np, nh, block_ranks); + const auto selected_assembly = assembly.value_or( + (!block_ranks.empty() || hbar_expansion_ == HbarExpansion::Bernoulli) + ? UCCEOMAssembly::ProjectedHbar + : UCCEOMAssembly::Commutator); + if (unitary()) + return assemble_ucc_eom(np, nh, block_ranks, selected_assembly); - // construct hbar - const auto hbar = this->hbar(); + if (!block_ranks.empty()) + throw Exception("CC::eom_r: block_ranks require a unitary ansatz"); + if (selected_assembly == UCCEOMAssembly::ProjectedHbar) + throw Exception( + "CC::eom_r: projected Hbar assembly requires a unitary ansatz"); - // construct [hbar, R] - ExprPtr hbar_R; - // for unitary ansatz, we need to compute the commutator [hbar, R], otherwise - // just hbar * R is sufficient because ref_av uses connectivity - if (this->unitary()) { - hbar_R = commutator(hbar, R(np, nh, eom_norm)); - } else { - hbar_R = hbar * R(np, nh, eom_norm); - } + const auto hbar = this->hbar(); + const auto hbar_R = hbar * R(np, nh, eom_norm); + const auto op_connect = concat(default_op_connections(), + {{L"h", L"R"}, {L"f", L"R"}, {L"g", L"R"}}); - // connectivity: empty for unitary ansatz, build otherwise - OpConnections op_connect; - if (!this->unitary()) { - // default connections + connect R with {h,f,g} - op_connect = concat(default_op_connections(), - {{L"h", L"R"}, {L"f", L"R"}, {L"g", L"R"}}); - } using std::min; - std::vector result(min(np, nh) + 1); // for EE element 0 stays null - // project with result(min(np, nh) + 1); for (const auto& [rp, rh] : eom_manifolds(np, nh)) - result.at(min(rp, rh)) = - this->ref_av(δl(nₚ(rp), nₕ(rh)) * hbar_R, op_connect); + result.at(min(rp, rh)) = ref_av(δl(nₚ(rp), nₕ(rh)) * hbar_R, op_connect); return result; } diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 83ab3145e7..b2724d5fc3 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -39,6 +39,14 @@ class CC { Bernoulli }; + /// Assembly used for the right-hand UCC EOM equations. + enum class UCCEOMAssembly { + /// Project \f$ [\bar{H}, R] \f$. + Commutator, + /// Assemble the projected \f$ \bar{H} \f$ matrix. + ProjectedHbar + }; + /// Configuration options for CC class struct Options { SEQUANT_DESIGNATED_INIT_ONLY; @@ -209,23 +217,26 @@ class CC { /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING /// manifold rank, so one set of numbers serves EE, IP and EA (read S as /// 1h/1p and D as 2h1p/1h2p; 10.1021/acs.jctc.5c01991 Fig. 1 carries the - /// IP/EA ranks). Empty (the default) fills every block with - /// `hbar_comm_rank`. Which form gets built depends on the expansion: - /// `Bernoulli` always uses the blocked form, `BCH` uses it only when a - /// matrix is given and builds a commutator otherwise. The BCH forms differ - /// by off-diagonal amplitude-residual terms unless those residuals vanish. - /// @pre if non-empty, requires a unitary ansatz; a non-unitary H̄ terminates - /// and has nothing to truncate. - /// @pre `block_ranks` is either empty or `K`×`K` + /// IP/EA ranks). For UCC, `Bernoulli` and non-empty `block_ranks` select + /// projected-H̄ assembly; otherwise `BCH` selects commutator assembly. The + /// two BCH assemblies differ by ground-state amplitude-equation residuals + /// and agree when those vanish. Traditional CC always uses its connected + /// H̄R product and does not support block ranks. + /// @param assembly optional explicit UCC assembly; when omitted, the default + /// described above is used. Traditional CC accepts only `Commutator`. + /// @throw Exception if `Commutator` is used with the Bernoulli expansion, if + /// `ProjectedHbar` or non-empty `block_ranks` is used with a traditional + /// ansatz, or if `block_ranks` is not `K`×`K` /// @return projected sigma equations in a vector of size `min(np, nh) + 1`, /// indexed by the smaller particle/hole count of each projection manifold. - /// With per-block truncation, each diagonal subtracts the scalar carried by - /// the same temporary \f$ \bar{H}^{(k)} \f$ used for that block. This leaves - /// the blockwise normal-ordered components of Eq. (10), not distinct - /// physical energy zeros. Element 0 is null iff `np == nh` + /// With projected-H̄ UCC assembly, each diagonal subtracts the scalar + /// carried by the same temporary \f$ \bar{H}^{(k)} \f$ used for that block. + /// This leaves the per-block normal-ordered components of Eq. (10), not + /// distinct physical energy zeros. Element 0 is null iff `np == nh` // clang-format on [[nodiscard]] std::vector eom_r( - nₚ np, nₕ nh, const std::vector& block_ranks = {}) const; + nₚ np, nₕ nh, const std::vector& block_ranks = {}, + std::optional assembly = std::nullopt) const; /// @brief derives left-side sigma equations for EOM-CC /// @param np number of particle annihilators in L operator @@ -271,16 +282,17 @@ class CC { std::optional pertbar_comm_rank_ = std::nullopt; HbarExpansion hbar_expansion_ = HbarExpansion::BCH; - /// @brief `eom_r`'s per-block-truncated path, taken whenever `block_ranks` is - /// non-empty or the expansion is Bernoulli - /// @param block_ranks see `eom_r`; empty means uniform `hbar_comm_rank` + /// @brief assembles the right-hand UCC EOM equations + /// @param block_ranks see `eom_r`; empty uses the configured H̄ rank, or 4 + /// @param assembly assembly to use /// @pre a unitary ansatz /// @note under the Bernoulli expansion each block's H̄ has its N part removed. /// Where a block's rank equals `hbar_comm_rank` the removed terms vanish at /// converged amplitudes, so its equations change but its values do not; /// below that rank the values change too. `BCH` keeps them. - [[nodiscard]] std::vector eom_r_blocked( - nₚ np, nₕ nh, const std::vector& block_ranks) const; + [[nodiscard]] std::vector assemble_ucc_eom( + nₚ np, nₕ nh, const std::vector& block_ranks, + UCCEOMAssembly assembly) const; /// @return the `LSTOptions` this engine uses for every `mbpt::lst()` call /// @note The choice of commutator representation is really a question of diff --git a/doc/user/guide/cc.rst b/doc/user/guide/cc.rst index 5bf1a738f7..e5eeb00155 100644 --- a/doc/user/guide/cc.rst +++ b/doc/user/guide/cc.rst @@ -78,10 +78,15 @@ Equation-of-Motion Coupled-Cluster .. code-block:: cpp - std::vector eom_r(nₚ np, nₕ nh); + std::vector eom_r( + nₚ np, nₕ nh, + const std::vector& block_ranks = {}, + std::optional assembly = std::nullopt); std::vector eom_l(nₚ np, nₕ nh); -Derives equation-of-motion coupled-cluster (EOM-CC) equations for excited states. The ``eom_r`` method generates equations for the right eigenvectors, while ``eom_l`` generates equations for the left eigenvectors. +Derives equation-of-motion coupled-cluster (EOM-CC) equations for excited states. The ``eom_r`` method generates equations for the right eigenvectors, while ``eom_l`` generates equations for the left eigenvectors. Traditional CC always uses the connected H̄R product. For UCC, passing ``CC::UCCEOMAssembly::Commutator`` or ``CC::UCCEOMAssembly::ProjectedHbar`` explicitly selects the right-hand assembly. Without that argument, a Bernoulli expansion or non-empty ``block_ranks`` selects projected-H̄ assembly; otherwise, a BCH expansion selects commutator assembly. + +For UCC, the optional ``block_ranks`` argument gives the per-block truncation orders as a row-major matrix over the EOM manifolds: nested-commutator order for BCH and :math:`\bar{H}^{k}` order for Bernoulli. Traditional CC does not support block ranks or projected-H̄ assembly; Bernoulli does not support commutator assembly. Examples -------- diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 9b18cf211a..6a9f2533ce 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -281,7 +281,8 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { // ... and the ansatz must be unitary REQUIRE_THROWS_AS(CC(2).eom_r(nₚ(2), nₕ(2), quccsd), Exception); - // Public API contract: under Bernoulli these are the same blocked request. + // Public API contract: under Bernoulli these are the same projected-H̄ + // request. // Check both returned manifolds; this does not claim BCH equivalence, since // no matrix takes the commutator path instead. const auto uniform = cc.eom_r(nₚ(2), nₕ(2)); @@ -442,6 +443,45 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { REQUIRE(size(eqs[1]) == 32); } + SECTION("EOM assembly") { + const auto tcc = CC{2}; + const auto ucc = CC(2, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}); + const auto ee_np = nₚ(2); + const auto ee_nh = nₕ(2); + const std::vector uniform_ranks = {2, 2, 2, 2}; + + REQUIRE_THAT( + tcc.eom_r(ee_np, ee_nh, {}, CC::UCCEOMAssembly::Commutator).at(2), + EquivalentTo(tcc.eom_r(ee_np, ee_nh).at(2))); + REQUIRE_THROWS_AS( + tcc.eom_r(ee_np, ee_nh, {}, CC::UCCEOMAssembly::ProjectedHbar), + Exception); + + const auto ucc_commutator = + ucc.eom_r(ee_np, ee_nh, uniform_ranks, CC::UCCEOMAssembly::Commutator); + const auto ucc_projected = ucc.eom_r(ee_np, ee_nh, uniform_ranks, + CC::UCCEOMAssembly::ProjectedHbar); + REQUIRE_THAT(ucc_commutator.at(2), + EquivalentTo(ucc.eom_r(ee_np, ee_nh).at(2))); + REQUIRE_THAT(ucc_projected.at(2), + EquivalentTo(ucc.eom_r(ee_np, ee_nh, uniform_ranks).at(2))); + REQUIRE_THAT(ucc_commutator.at(2), !EquivalentTo(ucc_projected.at(2))); + + // An asymmetric matrix pins the row-major mapping: lowering only H_SD + // changes the singles bra row, not the doubles row. + const auto ucc_row_mixed = + ucc.eom_r(ee_np, ee_nh, {2, 1, 2, 2}, CC::UCCEOMAssembly::Commutator); + REQUIRE_THAT(ucc_row_mixed.at(1), !EquivalentTo(ucc_commutator.at(1))); + REQUIRE_THAT(ucc_row_mixed.at(2), EquivalentTo(ucc_commutator.at(2))); + + const auto bern = CC(2, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 2, + .hbar_expansion = CC::HbarExpansion::Bernoulli}); + REQUIRE_THROWS_AS( + bern.eom_r(nₚ(1), nₕ(1), {}, CC::UCCEOMAssembly::Commutator), + Exception); + } + SECTION("EE-EOM-CCSD L") { const auto np = 2; const auto nh = 2; From 45f6c6f829d6f26a229335de0c659ff7b8cc3c89 Mon Sep 17 00:00:00 2001 From: Ajay Date: Fri, 28 Aug 2026 14:12:18 -0400 Subject: [PATCH 2/5] cc: format hbar and eom_r documentation as bullet points --- SeQuant/domain/mbpt/models/cc.hpp | 87 +++++++++++++++---------------- 1 file changed, 42 insertions(+), 45 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index b2724d5fc3..fa295a7aa7 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -116,26 +116,23 @@ class CC { /// the expansion; if not specified, will use the value of member /// `hbar_comm_rank`. If that is also not specified, will use 4 as the default /// value. If provided, will override all defaults. - /// @note For a non-unitary ansatz each commutator is represented as a - /// connected product, \f$ (\hat{A}\hat{B})_c \f$, equivalent to \f$ - /// [\hat{A},\hat{B}] \f$ only once operator connectivity is supplied - /// downstream; this engine always supplies it, by handing `ref_av` - /// `default_op_connections()` or a superset thereof. A unitary ansatz uses - /// explicit commutators instead, and is handed empty connectivity. - /// @warning Because of the above, the expression returned for a non-unitary - /// ansatz is NOT self-contained: evaluating it with empty connectivity, e.g. - /// `op::ref_av(P(nₚ(2)) * cc.hbar(), {.connect = {}})`, silently retains the - /// disconnected terms that the commutator would have cancelled. Pass - /// `default_op_connections()` (which `op::ref_av`/`op::vac_av` use when the - /// options argument is omitted -- note that passing `{}` instead gives empty - /// connections), or build H̄ yourself with `mbpt::lst(..., {})` if you need - /// the explicit form. For a unitary ansatz the reverse holds: H̄ is already - /// self-contained, so connectivity must be left empty. See the "Using H̄ - /// outside the CC class" section of the user guide. - /// @note Under `HbarExpansion::Bernoulli` the result is tensor-level, so it - /// takes `op::tensor` projectors and `op::tensor::ref_av`, not their `op` - /// counterparts. That path never reaches `CC::ref_av`, so `screen` and - /// `use_topology` have no effect on it. + /// @note The returned expression depends on the ansatz and expansion: + /// - A non-unitary ansatz represents each commutator as a connected + /// product, + /// \f$ (\hat{A}\hat{B})_c \f$. It is equivalent to \f$ + /// [\hat{A},\hat{B}] \f$ only when operator connectivity is supplied + /// downstream. + /// - A unitary BCH ansatz uses explicit commutators and empty connectivity. + /// - `HbarExpansion::Bernoulli` returns a tensor-level expression for use + /// with `op::tensor` projectors and `op::tensor::ref_av`. It never + /// reaches `CC::ref_av`, so `screen` and `use_topology` have no effect. + /// @warning A non-unitary H̄ is not self-contained. Evaluating it with empty + /// connectivity, e.g. `op::ref_av(P(nₚ(2)) * cc.hbar(), {.connect = {}})`, + /// retains disconnected terms. Pass `default_op_connections()` (the default + /// when the options argument is omitted), or build H̄ with + /// `mbpt::lst(..., {})` for an explicit form. A unitary H̄ is + /// self-contained, so its connectivity must be empty. See the "Using H̄ + /// outside the CC class" section of the user guide. [[nodiscard]] ExprPtr hbar( std::optional truncation_rank = std::nullopt) const; @@ -202,37 +199,37 @@ class CC { /// @brief derives right-side sigma equations for EOM-CC /// @param np number of particle creators in R operator /// @param nh number of hole creators in R operator - /// @param block_ranks optional per-block H̄ commutator truncation ranks: a - /// different H̄ in each block of the secular matrix instead of one uniform - /// H̄ everywhere. For singles+doubles the matrix is + /// @param block_ranks optional per-block H̄ truncation ranks: + /// - A `K`×`K` matrix is read row by row in ascending manifold rank. For + /// singles+doubles it is /// | H_SS H_SD | e.g. | 2 1 | /// | H_DS H_DD | | 1 0 | - /// read row by row, i.e. `{2,1,1,0}`: H_SS through the double commutator - /// [[H,σ],σ], H_SD and H_DS through the single [H,σ], H_DD the bare - /// Hamiltonian integrals (no commutators). Under `Bernoulli` a rank is the - /// Bernoulli order H̄^k instead, whose commutators are in V alone. These are - /// that paper's Bernoulli qUCCSD ranks, UCCSD[2|2,1,0] of 10.1063/5.0062090 - /// Eqs. (29), (41), (44), (48); its BCH scheme instead needs F one rank - /// above V, which one number per block cannot express. - /// `K` manifolds give a row-major `K`×`K` matrix ordered by ASCENDING - /// manifold rank, so one set of numbers serves EE, IP and EA (read S as - /// 1h/1p and D as 2h1p/1h2p; 10.1021/acs.jctc.5c01991 Fig. 1 carries the - /// IP/EA ranks). For UCC, `Bernoulli` and non-empty `block_ranks` select - /// projected-H̄ assembly; otherwise `BCH` selects commutator assembly. The - /// two BCH assemblies differ by ground-state amplitude-equation residuals - /// and agree when those vanish. Traditional CC always uses its connected - /// H̄R product and does not support block ranks. - /// @param assembly optional explicit UCC assembly; when omitted, the default - /// described above is used. Traditional CC accepts only `Commutator`. + /// Thus `{2,1,1,0}` uses [[H,σ],σ] for H_SS, [H,σ] for H_SD and + /// H_DS, and bare Hamiltonian integrals for H_DD. + /// - The same ordering serves EE, IP, and EA: read S as 1h/1p and D as + /// 2h1p/1h2p (10.1021/acs.jctc.5c01991, Fig. 1). + /// - Under `Bernoulli`, each rank is a Bernoulli order H̄^k whose + /// commutators are in V alone. These are the Bernoulli qUCCSD ranks + /// UCCSD[2|2,1,0] of 10.1063/5.0062090, Eqs. (29), (41), (44), (48). Its + /// BCH scheme needs F one rank above V, which one rank per block cannot + /// express. + /// @param assembly optional explicit UCC assembly: + /// - If omitted, `Bernoulli` and non-empty `block_ranks` select + /// `ProjectedHbar`; otherwise `BCH` selects `Commutator`. + /// - The two BCH assemblies differ by ground-state amplitude-equation + /// residuals and agree when those vanish. + /// - Traditional CC accepts only `Commutator`, uses its connected H̄R + /// product, and does not support block ranks. /// @throw Exception if `Commutator` is used with the Bernoulli expansion, if /// `ProjectedHbar` or non-empty `block_ranks` is used with a traditional /// ansatz, or if `block_ranks` is not `K`×`K` /// @return projected sigma equations in a vector of size `min(np, nh) + 1`, - /// indexed by the smaller particle/hole count of each projection manifold. - /// With projected-H̄ UCC assembly, each diagonal subtracts the scalar - /// carried by the same temporary \f$ \bar{H}^{(k)} \f$ used for that block. - /// This leaves the per-block normal-ordered components of Eq. (10), not - /// distinct physical energy zeros. Element 0 is null iff `np == nh` + /// indexed by the smaller particle/hole count of each projection manifold: + /// - Element 0 is null iff `np == nh`. + /// - With projected-H̄ UCC assembly, each diagonal subtracts the scalar + /// carried by the same temporary \f$ \bar{H}^{(k)} \f$ used for that block, + /// leaving the per-block normal-ordered components of Eq. (10), not + /// distinct physical energy zeros. // clang-format on [[nodiscard]] std::vector eom_r( nₚ np, nₕ nh, const std::vector& block_ranks = {}, From 68232d0739750f5f3ecf1e4a8d933c4a0d887b46 Mon Sep 17 00:00:00 2001 From: Ajay Date: Thu, 9 Jul 2026 01:34:35 -0400 Subject: [PATCH 3/5] cc: implement support for addtional sigma_{1} commutators in UCC --- SeQuant/domain/mbpt/models/cc.cpp | 13 ++++++++- SeQuant/domain/mbpt/models/cc.hpp | 7 +++++ tests/unit/test_mbpt_cc.cpp | 46 +++++++++++++++++++++++++++++++ 3 files changed, 65 insertions(+), 1 deletion(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 2af56c7585..4cd0a88ce4 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -45,6 +45,7 @@ CC::CC(size_t n, const Options& opts) screen_(opts.screen), use_topology_(opts.use_topology), hbar_comm_rank_(opts.hbar_comm_rank), + hbar_singles_comm_rank_(opts.hbar_singles_comm_rank), pertbar_comm_rank_(opts.pertbar_comm_rank), hbar_expansion_(opts.hbar_expansion) { if (hbar_expansion_ == HbarExpansion::Bernoulli && ansatz_ != Ansatz::U) @@ -82,7 +83,17 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { // for a non-unitary ansatz this is the cheaper connected-product form, which // is only equivalent to the commutator once the caller supplies operator // connectivity to ref_av (see lst_options() and the @warning on hbar()) - return mbpt::lst(H(), T(N, skip_singles()), truncation, lst_options()); + auto result = mbpt::lst(H(), T(N, skip_singles()), truncation, lst_options()); + + // extra singles-only commutators: wrap H̄ with the t1 similarity transform + // to order K. Same commutator form as above, since the same connectivity is + // supplied downstream. + if (hbar_singles_comm_rank_.value_or(0) > 0) { + auto opts = lst_options(); + opts.skip_clone = true; + result = mbpt::lst(result, op::t(1), *hbar_singles_comm_rank_, opts); + } + return result; } ExprPtr CC::energy(std::optional comm_rank) const { diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index fa295a7aa7..862db5e3b5 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -65,6 +65,12 @@ class CC { /// maximum order of nested commutators in H̄; must be specified if unitary /// ansatz is used std::optional hbar_comm_rank = std::nullopt; + /// order K of the additional singles-only (t1) similarity transform applied + /// on top of the standard hbar_comm_rank (R) commutator series. nullopt or + /// 0 disables it (default). Designed for the unitary ansatz; applied + /// generically. See docs/superpowers/specs/2026-07-07-ucc-extra-singles- + /// commutators-design.md + std::optional hbar_singles_comm_rank = std::nullopt; /// maximum order of nested commutators in the similarity transformed /// perturbation operator; must be specified if unitary ansatz is used in /// perturbed amplitude derivation @@ -276,6 +282,7 @@ class CC { bool screen_ = true; bool use_topology_ = true; std::optional hbar_comm_rank_ = std::nullopt; + std::optional hbar_singles_comm_rank_ = std::nullopt; std::optional pertbar_comm_rank_ = std::nullopt; HbarExpansion hbar_expansion_ = HbarExpansion::BCH; diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 6a9f2533ce..03c02d6cc6 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -411,6 +412,51 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { } } // SECTION("rdm") + SECTION("extra singles commutators") { + // hbar_singles_comm_rank wraps H̄_R with a factored t1 similarity transform + // to order K: H̄ = Σ_a (1/a!) [H̄_R, t1−t1†]_a. + const auto N = 2; + const CC::Options base{.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}; + const auto hbar_R = CC(N, base).hbar(); // K unset == the R=2 baseline + + // anti-Hermitian singles generator σ1 = t1 − t1† + const auto sigma1 = op::t(1) - adjoint(op::t(1)); + // commutator built independently of mbpt::lst (plain expression arithmetic) + const auto comm = [](const ExprPtr& A, const ExprPtr& B) { + auto r = A * B - B * A; + simplify(r); + return r; + }; + + SECTION("K=0/unset reproduces plain hbar") { + CC::Options k0 = base; + k0.hbar_singles_comm_rank = 0; + REQUIRE_THAT(CC(N, k0).hbar(), EquivalentTo(hbar_R)); + } + + SECTION("K=1: extra term is exactly [H̄_R, σ1]") { + CC::Options k1 = base; + k1.hbar_singles_comm_rank = 1; + const auto hbar_k1 = CC(N, k1).hbar(); + + // the wrap genuinely adds terms not present in the baseline + const auto extra = comm(hbar_R, sigma1); + REQUIRE(size(extra) > 0); + REQUIRE(size(hbar_k1) > size(hbar_R)); + + REQUIRE_THAT(hbar_k1, EquivalentTo(hbar_R + extra)); + } + + SECTION("K=2: increment is (1/2!)[[H̄_R, σ1], σ1]") { + CC::Options k2 = base; + k2.hbar_singles_comm_rank = 2; + const auto expected = + hbar_R + comm(hbar_R, sigma1) + + ex(rational{1, 2}) * comm(comm(hbar_R, sigma1), sigma1); + REQUIRE_THAT(CC(N, k2).hbar(), EquivalentTo(expected)); + } + } // SECTION("extra singles commutators") + SECTION("eom_cc"){SECTION("EOM-CCSD"){const auto N = 2; auto cc = CC{N}; SECTION("EE-EOM-CCSD R") { From bfa5f437697af51314f0858609526eae674bc7c5 Mon Sep 17 00:00:00 2001 From: Ajay Date: Tue, 11 Aug 2026 18:03:46 -0400 Subject: [PATCH 4/5] cc: store Options and add with() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Keeps CC::Options whole rather than decomposing it into one member per option, so with() can round-trip it losslessly — Options::skip_singles stays optional and is resolved on demand by skip_singles(), rather than being pinned to a resolved bool at construction. with() takes a single callable mutating Options and rebuilds through the normal constructor, so the result is validated while no intermediate state exists to be invalid. That matters because the ctor invariants couple the options: a chain of per-option setters would have to pass through configurations that cannot be valid. This supersedes with_hbar_comm_rank, which was removed separately. --- SeQuant/domain/mbpt/models/cc.cpp | 97 ++++++++++++++++--------------- SeQuant/domain/mbpt/models/cc.hpp | 42 ++++++++++--- tests/unit/test_mbpt_cc.cpp | 52 +++++++++++++++-- 3 files changed, 130 insertions(+), 61 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 4cd0a88ce4..2fc56c1b7b 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -37,47 +37,46 @@ namespace sequant::mbpt { CC::CC(size_t n) : CC(n, Options{}) {} -CC::CC(size_t n, const Options& opts) - : N(n), - ansatz_(opts.ansatz), - skip_singles_(opts.skip_singles.value_or(ansatz_ == Ansatz::oT || - ansatz_ == Ansatz::oU)), - screen_(opts.screen), - use_topology_(opts.use_topology), - hbar_comm_rank_(opts.hbar_comm_rank), - hbar_singles_comm_rank_(opts.hbar_singles_comm_rank), - pertbar_comm_rank_(opts.pertbar_comm_rank), - hbar_expansion_(opts.hbar_expansion) { - if (hbar_expansion_ == HbarExpansion::Bernoulli && ansatz_ != Ansatz::U) +CC::CC(size_t n, const Options& opts) : N(n), opts_(opts) { + if (opts_.hbar_expansion == HbarExpansion::Bernoulli && + opts_.ansatz != Ansatz::U) throw Exception("CC: Bernoulli expansion requires the U ansatz"); - if (unitary() && !hbar_comm_rank_) + if (unitary() && !opts_.hbar_comm_rank) throw Exception("CC: hbar_comm_rank is required for unitary ansatz"); - if (ansatz_ == Ansatz::oT || ansatz_ == Ansatz::oU) - SEQUANT_ASSERT(skip_singles_, + if (opts_.ansatz == Ansatz::oT || opts_.ansatz == Ansatz::oU) + SEQUANT_ASSERT(skip_singles(), "CC: skip_singles must be true for orbital-optimized " "ansatz"); } -CC::Ansatz CC::ansatz() const { return ansatz_; } +CC::Ansatz CC::ansatz() const { return opts_.ansatz; } bool CC::unitary() const { - return ansatz_ == Ansatz::U || ansatz_ == Ansatz::oU; + return opts_.ansatz == Ansatz::U || opts_.ansatz == Ansatz::oU; } -std::optional CC::hbar_comm_rank() const { return hbar_comm_rank_; } +std::optional CC::hbar_comm_rank() const { + return opts_.hbar_comm_rank; +} -CC::HbarExpansion CC::hbar_expansion() const { return hbar_expansion_; } +CC::HbarExpansion CC::hbar_expansion() const { return opts_.hbar_expansion; } -bool CC::skip_singles() const { return skip_singles_; } +/// resolves the `skip_singles` default: on for orbital-optimized ansätze, off +/// otherwise. Kept here rather than in the ctor so `Options` round-trips. +bool CC::skip_singles() const { + return opts_.skip_singles.value_or(opts_.ansatz == Ansatz::oT || + opts_.ansatz == Ansatz::oU); +} -bool CC::screen() const { return screen_; } +bool CC::screen() const { return opts_.screen; } -bool CC::use_topology() const { return use_topology_; } +bool CC::use_topology() const { return opts_.use_topology; } ExprPtr CC::hbar(std::optional truncation_rank) const { - const auto truncation = truncation_rank.value_or(hbar_comm_rank_.value_or(4)); + const auto truncation = + truncation_rank.value_or(opts_.hbar_comm_rank.value_or(4)); - if (hbar_expansion_ == HbarExpansion::Bernoulli) + if (opts_.hbar_expansion == HbarExpansion::Bernoulli) return bernoulli::hbar(N, truncation, skip_singles()); // for a non-unitary ansatz this is the cheaper connected-product form, which @@ -88,10 +87,10 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { // extra singles-only commutators: wrap H̄ with the t1 similarity transform // to order K. Same commutator form as above, since the same connectivity is // supplied downstream. - if (hbar_singles_comm_rank_.value_or(0) > 0) { + if (opts_.hbar_singles_comm_rank.value_or(0) > 0) { auto opts = lst_options(); opts.skip_clone = true; - result = mbpt::lst(result, op::t(1), *hbar_singles_comm_rank_, opts); + result = mbpt::lst(result, op::t(1), *opts_.hbar_singles_comm_rank, opts); } return result; } @@ -99,8 +98,8 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { ExprPtr CC::energy(std::optional comm_rank) const { // Bernoulli: the hbar expansion is at tensor level, call the tensor level // ref_av directly. No connectivity or screening. - if (hbar_expansion_ == HbarExpansion::Bernoulli) { - const auto erank = comm_rank.value_or(hbar_comm_rank_.value()); + if (opts_.hbar_expansion == HbarExpansion::Bernoulli) { + const auto erank = comm_rank.value_or(opts_.hbar_comm_rank.value()); return op::tensor::ref_av(this->hbar(erank)); } // <0|H̄|0>: reference expectation value of H̄ at the requested commutator @@ -117,7 +116,7 @@ std::vector CC::t(size_t pmax, size_t pmin) const { // Bernoulli: the hbar expansion is at tensor level, project and call the // tensor level ref_av directly. - if (hbar_expansion_ == HbarExpansion::Bernoulli) { + if (opts_.hbar_expansion == HbarExpansion::Bernoulli) { const auto hbar = this->hbar(); std::vector result(pmax + 1); for (std::int64_t p = pmax; p >= static_cast(pmin); --p) { @@ -145,7 +144,7 @@ std::vector CC::t(size_t pmax, size_t pmin) const { std::shared_ptr hbar_le_p; // keeps products that can produce excitations rank <=p - if (screen_) { // if operator level screening is on + if (opts_.screen) { // if operator level screening is on for (auto& term : *hbar) { SEQUANT_ASSERT(term->is() || term->is()); if (raises_vacuum_up_to_rank(term, p)) { @@ -179,7 +178,7 @@ std::vector CC::λ() const { SEQUANT_ASSERT(!unitary(), "there is no need for CC::λ for unitary ansatz"); // construct hbar - const auto commutator_rank = hbar_comm_rank_.value_or(4); + const auto commutator_rank = opts_.hbar_comm_rank.value_or(4); SEQUANT_ASSERT(commutator_rank >= 1, "CC::λ: hbar_comm_rank must be >= 1"); auto hbar = this->hbar(commutator_rank - 1); // -1 because of the connection with the projector @@ -212,8 +211,8 @@ std::vector CC::λ() const { std::shared_ptr lhbar_for_vev; // keeps products that can produce non-zero VEV std::shared_ptr - lhbar_le_p; // keeps products that can produce excitations rank <=p - if (screen_) { // if operator level screening is enabled + lhbar_le_p; // keeps products that can produce excitations rank <=p + if (opts_.screen) { // if operator level screening is enabled for (auto& term : *lhbar) { // pick terms from lhbar SEQUANT_ASSERT(term->is() || term->is()); @@ -244,7 +243,7 @@ std::vector CC::λ() const { } ExprPtr CC::rdm(size_t rank, std::optional comm_rank) const { - SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli, + SEQUANT_ASSERT(opts_.hbar_expansion != HbarExpansion::Bernoulli, "CC::rdm: the Bernoulli expansion is not supported yet"); // 1. replacement operator {ã^{p_1..p_r}_{p_{r+1}..p_{2r}}} (see op::ã); its @@ -262,7 +261,7 @@ ExprPtr CC::rdm(size_t rank, std::optional comm_rank) const { // connected-product form; the {ã,t} connectivity handed to ref_av below is // what makes it equivalent to the explicit commutator. const auto commutator_rank = comm_rank.value_or( - unitary() ? hbar_comm_rank_.value() : std::min(2 * rank, rank + N)); + unitary() ? opts_.hbar_comm_rank.value() : std::min(2 * rank, rank + N)); auto bar = mbpt::lst(replacer, T(N, skip_singles()), commutator_rank, lst_options()); @@ -286,22 +285,23 @@ std::vector CC::tʼ(size_t rank, size_t order, "sequant::mbpt::CC::tʼ(): only one-body perturbation " "operator is supported now"); if (unitary()) - SEQUANT_ASSERT(pertbar_comm_rank_, + SEQUANT_ASSERT(opts_.pertbar_comm_rank, "pertbar_comm_rank must be specified for unitary ansatz"); - SEQUANT_ASSERT(hbar_expansion_ != HbarExpansion::Bernoulli, + SEQUANT_ASSERT(opts_.hbar_expansion != HbarExpansion::Bernoulli, "CC::tʼ: the Bernoulli expansion is not supported yet"); // construct h1_bar // truncate h1_bar at rank 2 for one-body perturbation operator and at rank 4 // for two-body perturbation operator; unless specified otherwise const auto h1_truncate_default = rank == 1 ? 2 : 4; - const auto h1_truncate_at = pertbar_comm_rank_.value_or(h1_truncate_default); + const auto h1_truncate_at = + opts_.pertbar_comm_rank.value_or(h1_truncate_default); const auto h1_bar = mbpt::lst(Hʼ(rank, {.order = order, .nbatch = nbatch}), T(N, skip_singles()), h1_truncate_at, lst_options()); // construct [hbar, Tʼ(1)] - const auto hbar_truncate_at = hbar_comm_rank_.value_or( + const auto hbar_truncate_at = opts_.hbar_comm_rank.value_or( 3); // notice 3 instead of 4 here, this is because of the commutator with // T'(1). In case 4 is used, it will generate more terms but they // will not contribute. @@ -348,7 +348,7 @@ std::vector CC::λʼ(size_t rank, size_t order, "sequant::mbpt::CC::λʼ(): only one-body perturbation " "operator is supported now"); SEQUANT_ASSERT(!unitary(), "there is no need for CC::λʼ for unitary ansatz"); - SEQUANT_ASSERT(ansatz_ == Ansatz::T, + SEQUANT_ASSERT(opts_.ansatz == Ansatz::T, "CC::λʼ: only traditional ansatz is supported"); // construct hbar @@ -357,8 +357,8 @@ std::vector CC::λʼ(size_t rank, size_t order, // construct h1_bar // truncate h1_bar at rank 2 for one-body perturbation operator and at rank 4 // for two-body perturbation operator; unless specified otherwise - const auto h1_truncate_at = (rank == 1) ? pertbar_comm_rank_.value_or(2) - : pertbar_comm_rank_.value_or(4); + const auto h1_truncate_at = (rank == 1) ? opts_.pertbar_comm_rank.value_or(2) + : opts_.pertbar_comm_rank.value_or(4); // connected form (this path is non-unitary, see the assert above); the // op_connect built below is a superset of default_op_connections() and so // supplies the connectivity that makes it equivalent to the commutator @@ -434,7 +434,7 @@ std::vector CC::assemble_ucc_eom( nₚ np, nₕ nh, const std::vector& block_ranks, UCCEOMAssembly assembly) const { if (assembly == UCCEOMAssembly::Commutator && - hbar_expansion_ == HbarExpansion::Bernoulli) + opts_.hbar_expansion == HbarExpansion::Bernoulli) throw Exception("CC::eom_r: Bernoulli requires projected Hbar assembly"); using std::min; @@ -461,8 +461,9 @@ std::vector CC::assemble_ucc_eom( "K = number of projection manifolds"); // Bernoulli H̄ is tensor-level, BCH H̄ operator-level; the bra/ket/vev trio - // below must match it. UCC uses no connectivity. - const bool tensor_level = hbar_expansion_ == HbarExpansion::Bernoulli; + // below must match it. Connectivity is empty either way, as everywhere on the + // unitary path; only the operator-level vev forwards screen/use_topology. + const bool tensor_level = opts_.hbar_expansion == HbarExpansion::Bernoulli; // One H̄ per distinct truncation order, reduced to its R part (Bernoulli only; // BCH H̄ is operator-level and has no N/R split). The N part is the amplitude @@ -488,8 +489,8 @@ std::vector CC::assemble_ucc_eom( auto vev = [tensor_level, this](const ExprPtr& e) { return tensor_level ? op::tensor::ref_av(e) : op::ref_av(e, {.connect = {}, - .screen = screen_, - .use_topology = use_topology_}); + .screen = opts_.screen, + .use_topology = opts_.use_topology}); }; std::vector result(min(np, nh) + 1); @@ -528,7 +529,7 @@ std::vector CC::eom_r(nₚ np, nₕ nh, "spin-free basis does not yet support non particle-conserving cases"); const auto selected_assembly = assembly.value_or( - (!block_ranks.empty() || hbar_expansion_ == HbarExpansion::Bernoulli) + (!block_ranks.empty() || opts_.hbar_expansion == HbarExpansion::Bernoulli) ? UCCEOMAssembly::ProjectedHbar : UCCEOMAssembly::Commutator); if (unitary()) diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 862db5e3b5..490bf6c6e5 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include namespace sequant { @@ -105,6 +106,30 @@ class CC { /// @return the choice of H̄ expansion [[nodiscard]] HbarExpansion hbar_expansion() const; + // clang-format off + /// @brief returns a copy of this engine with @p mutate applied to its options + /// @param mutate a callable taking `Options&` + /// @return the mutated copy + /// @note The copy is built through the normal constructor, so the *result* is + /// validated; intermediate states are not, since there are none. This is why + /// the options are mutated in one callable rather than by a chain of + /// per-option setters: the ctor invariants couple the options (unitary needs + /// `hbar_comm_rank`; `oT`/`oU` need `skip_singles`), so a chain would have + /// to visit configurations that cannot be valid. + /// @code + /// auto report = cc.with([](auto& o) { + /// o.ansatz = Ansatz::U; + /// o.hbar_comm_rank = 3; + /// }); + /// @endcode + // clang-format on + template + [[nodiscard]] CC with(F&& mutate) const { + auto o = opts_; + std::forward(mutate)(o); + return CC(N, o); + } + /// @return true if singles amplitudes are excluded from \f$ \hat{T} \f$ and /// \f$ \hat{\Lambda} \f$ [[nodiscard]] bool skip_singles() const; @@ -277,14 +302,15 @@ class CC { private: size_t N; - Ansatz ansatz_ = Ansatz::T; - bool skip_singles_ = false; - bool screen_ = true; - bool use_topology_ = true; - std::optional hbar_comm_rank_ = std::nullopt; - std::optional hbar_singles_comm_rank_ = std::nullopt; - std::optional pertbar_comm_rank_ = std::nullopt; - HbarExpansion hbar_expansion_ = HbarExpansion::BCH; + /// @note Stored whole, rather than decomposed into one member per option, so + /// that `with()` can round-trip it losslessly. In particular + /// `Options::skip_singles` stays `std::optional` and is resolved on + /// demand by `skip_singles()`. Were the resolved `bool` stored instead, + /// `with()` would pin it, and `cc.with([](auto& o){ o.ansatz = Ansatz::oU; + /// })` on a non-orbital-optimized engine would carry `false` into an ansatz + /// whose ctor requires `true`, where fresh construction defaults it + /// correctly. + Options opts_; /// @brief assembles the right-hand UCC EOM equations /// @param block_ranks see `eom_r`; empty uses the configured H̄ rank, or 4 diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 03c02d6cc6..14f452643d 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -310,17 +310,22 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { } } // SECTION("energy") - SECTION("hbar_comm_rank") { + SECTION("with") { const auto N = 2; - REQUIRE(CC(N, {.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}) - .hbar_comm_rank() == 2); + const CC::Options opts{.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}; + auto bch2 = CC(N, opts); + auto bch3 = bch2.with([](auto& o) { o.hbar_comm_rank = 3; }); + REQUIRE(bch2.hbar_comm_rank() == 2); + REQUIRE(bch3.hbar_comm_rank() == 3); // rank 0 is a valid truncation (H̄ = H); only CC::λ rejects it, since it // derives at rank - 1 REQUIRE_THROWS_AS(CC(N, {.ansatz = CC::Ansatz::U}), Exception); + REQUIRE(bch2.with([](auto& o) { o.hbar_comm_rank = 0; }).hbar_comm_rank() == + 0); if (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { REQUIRE_THROWS_AS(CC(N, {.hbar_comm_rank = 0}).λ(), Exception); } - } // SECTION("hbar_comm_rank") + } // SECTION("with") SECTION("rdm") { constexpr auto N = 2; @@ -457,6 +462,42 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { } } // SECTION("extra singles commutators") + SECTION("copy-and-mutate") { + const auto N = 2; + const CC::Options base{.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}; + + SECTION("identity mutation reproduces the engine") { + const CC cc(N, base); + REQUIRE_THAT(cc.with([](auto&) {}).hbar(), EquivalentTo(cc.hbar())); + } + + SECTION("changed rank matches fresh construction") { + const CC cc(N, base); + CC::Options at3 = base; + at3.hbar_comm_rank = 3; + // N.B. cannot compare against with_hbar_comm_rank -- this commit deletes + // it + REQUIRE_THAT(cc.with([](auto& o) { o.hbar_comm_rank = 3; }).hbar(), + EquivalentTo(CC(N, at3).hbar())); + } + + SECTION("skip_singles is not pinned by the round trip") { + // base leaves skip_singles unset, so CC resolves it to false for U. + // Flipping to oU must pick up the oU default (true), not carry that + // resolved false into an ansatz whose ctor requires true. + const CC cc(N, base); + const auto ou = cc.with([](auto& o) { o.ansatz = CC::Ansatz::oU; }); + REQUIRE(ou.skip_singles()); + } + + SECTION("mutation does not disturb the source engine") { + const CC cc(N, base); + const auto hbar_before = cc.hbar(); + (void)cc.with([](auto& o) { o.hbar_comm_rank = 4; }); + REQUIRE_THAT(cc.hbar(), EquivalentTo(hbar_before)); + } + } // SECTION("copy-and-mutate") + SECTION("eom_cc"){SECTION("EOM-CCSD"){const auto N = 2; auto cc = CC{N}; SECTION("EE-EOM-CCSD R") { @@ -571,7 +612,8 @@ SECTION("ucc") { for (auto k = 0; k <= N; ++k) { REQUIRE(t_eqs[k]); } - // these are numerically verified against http://arxiv.org/abs/2503.00617 + // these are numerically verified against + // http://arxiv.org/abs/2503.00617 const auto energy_nterms = size(t_eqs[0]); if (c == 2) REQUIRE(energy_nterms == 20); if (c == 3) REQUIRE(energy_nterms == 74); From 952501e7c1ca922472d1381d4e86391505c0e9aa Mon Sep 17 00:00:00 2001 From: Ajay Date: Fri, 28 Aug 2026 15:08:18 -0400 Subject: [PATCH 5/5] cc: reject incompatible construction options --- SeQuant/domain/mbpt/models/cc.cpp | 22 ++++++------- SeQuant/domain/mbpt/models/cc.hpp | 51 +++++++++------------------- doc/user/guide/cc.rst | 8 +++-- tests/unit/test_mbpt_cc.cpp | 55 +++++++------------------------ 4 files changed, 44 insertions(+), 92 deletions(-) diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 2fc56c1b7b..4a42bf7f93 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -41,12 +41,16 @@ CC::CC(size_t n, const Options& opts) : N(n), opts_(opts) { if (opts_.hbar_expansion == HbarExpansion::Bernoulli && opts_.ansatz != Ansatz::U) throw Exception("CC: Bernoulli expansion requires the U ansatz"); + if (opts_.hbar_expansion == HbarExpansion::Bernoulli && + opts_.hbar_singles_comm_rank > 0) + throw Exception( + "CC: hbar_singles_comm_rank is not supported with Bernoulli"); if (unitary() && !opts_.hbar_comm_rank) throw Exception("CC: hbar_comm_rank is required for unitary ansatz"); - if (opts_.ansatz == Ansatz::oT || opts_.ansatz == Ansatz::oU) - SEQUANT_ASSERT(skip_singles(), - "CC: skip_singles must be true for orbital-optimized " - "ansatz"); + if ((opts_.ansatz == Ansatz::oT || opts_.ansatz == Ansatz::oU) && + !skip_singles()) + throw Exception( + "CC: skip_singles must be true for orbital-optimized ansatz"); } CC::Ansatz CC::ansatz() const { return opts_.ansatz; } @@ -61,8 +65,6 @@ std::optional CC::hbar_comm_rank() const { CC::HbarExpansion CC::hbar_expansion() const { return opts_.hbar_expansion; } -/// resolves the `skip_singles` default: on for orbital-optimized ansätze, off -/// otherwise. Kept here rather than in the ctor so `Options` round-trips. bool CC::skip_singles() const { return opts_.skip_singles.value_or(opts_.ansatz == Ansatz::oT || opts_.ansatz == Ansatz::oU); @@ -84,13 +86,11 @@ ExprPtr CC::hbar(std::optional truncation_rank) const { // connectivity to ref_av (see lst_options() and the @warning on hbar()) auto result = mbpt::lst(H(), T(N, skip_singles()), truncation, lst_options()); - // extra singles-only commutators: wrap H̄ with the t1 similarity transform - // to order K. Same commutator form as above, since the same connectivity is - // supplied downstream. - if (opts_.hbar_singles_comm_rank.value_or(0) > 0) { + // Apply the optional singles-only transform after the primary BCH series. + if (opts_.hbar_singles_comm_rank > 0) { auto opts = lst_options(); opts.skip_clone = true; - result = mbpt::lst(result, op::t(1), *opts_.hbar_singles_comm_rank, opts); + result = mbpt::lst(result, op::t(1), opts_.hbar_singles_comm_rank, opts); } return result; } diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 490bf6c6e5..73cb6dbb0c 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -66,12 +66,11 @@ class CC { /// maximum order of nested commutators in H̄; must be specified if unitary /// ansatz is used std::optional hbar_comm_rank = std::nullopt; - /// order K of the additional singles-only (t1) similarity transform applied - /// on top of the standard hbar_comm_rank (R) commutator series. nullopt or - /// 0 disables it (default). Designed for the unitary ansatz; applied - /// generically. See docs/superpowers/specs/2026-07-07-ucc-extra-singles- - /// commutators-design.md - std::optional hbar_singles_comm_rank = std::nullopt; + /// for BCH, maximum order of the additional singles-only similarity + /// transform applied after H̄; unitary ansätze use + /// \f$ \sigma_1 = T_1 - T_1^\dagger \f$, while non-unitary ansätze use + /// \f$ T_1 \f$. Zero disables the transform + size_t hbar_singles_comm_rank = 0; /// maximum order of nested commutators in the similarity transformed /// perturbation operator; must be specified if unitary ansatz is used in /// perturbed amplitude derivation @@ -89,8 +88,10 @@ class CC { /// @brief constructs CC engine with custom options /// @param n coupled cluster excitation rank /// @param opts configuration options @see CC::Options - /// @throw Exception if a unitary ansatz has no `hbar_comm_rank`, or if the - /// Bernoulli expansion is requested with an ansatz other than `Ansatz::U` + /// @throw Exception if a unitary ansatz has no `hbar_comm_rank`, if an + /// orbital-optimized ansatz includes singles, or if the Bernoulli expansion + /// is requested with an ansatz other than `Ansatz::U` or with a positive + /// `hbar_singles_comm_rank` explicit CC(size_t n, const Options& opts); /// @return the type of ansatz @@ -106,23 +107,9 @@ class CC { /// @return the choice of H̄ expansion [[nodiscard]] HbarExpansion hbar_expansion() const; - // clang-format off - /// @brief returns a copy of this engine with @p mutate applied to its options - /// @param mutate a callable taking `Options&` + /// @brief returns a validated copy with @p mutate applied to its options + /// @param mutate callable taking `Options&` /// @return the mutated copy - /// @note The copy is built through the normal constructor, so the *result* is - /// validated; intermediate states are not, since there are none. This is why - /// the options are mutated in one callable rather than by a chain of - /// per-option setters: the ctor invariants couple the options (unitary needs - /// `hbar_comm_rank`; `oT`/`oU` need `skip_singles`), so a chain would have - /// to visit configurations that cannot be valid. - /// @code - /// auto report = cc.with([](auto& o) { - /// o.ansatz = Ansatz::U; - /// o.hbar_comm_rank = 3; - /// }); - /// @endcode - // clang-format on template [[nodiscard]] CC with(F&& mutate) const { auto o = opts_; @@ -160,10 +147,10 @@ class CC { /// @warning A non-unitary H̄ is not self-contained. Evaluating it with empty /// connectivity, e.g. `op::ref_av(P(nₚ(2)) * cc.hbar(), {.connect = {}})`, /// retains disconnected terms. Pass `default_op_connections()` (the default - /// when the options argument is omitted), or build H̄ with - /// `mbpt::lst(..., {})` for an explicit form. A unitary H̄ is - /// self-contained, so its connectivity must be empty. See the "Using H̄ - /// outside the CC class" section of the user guide. + /// when the options argument is omitted), or reproduce every configured BCH + /// transform with explicit-commutator `mbpt::lst(..., {})` calls. A unitary + /// H̄ is self-contained, so its connectivity must be empty. See the "Using + /// H̄ outside the CC class" section of the user guide. [[nodiscard]] ExprPtr hbar( std::optional truncation_rank = std::nullopt) const; @@ -302,14 +289,6 @@ class CC { private: size_t N; - /// @note Stored whole, rather than decomposed into one member per option, so - /// that `with()` can round-trip it losslessly. In particular - /// `Options::skip_singles` stays `std::optional` and is resolved on - /// demand by `skip_singles()`. Were the resolved `bool` stored instead, - /// `with()` would pin it, and `cc.with([](auto& o){ o.ansatz = Ansatz::oU; - /// })` on a non-orbital-optimized engine would carry `false` into an ansatz - /// whose ctor requires `true`, where fresh construction defaults it - /// correctly. Options opts_; /// @brief assembles the right-hand UCC EOM equations diff --git a/doc/user/guide/cc.rst b/doc/user/guide/cc.rst index e5eeb00155..2f40b19d22 100644 --- a/doc/user/guide/cc.rst +++ b/doc/user/guide/cc.rst @@ -142,6 +142,8 @@ The similarity-transformed Hamiltonian is built by ``mbpt::lst``, see :ref:`mbpt :end-before: end-snippet-4 :dedent: 2 +For BCH expansions, ``Options::hbar_singles_comm_rank`` applies an additional singles-only similarity transform after the primary expansion. Unitary ansätze use :math:`\sigma_1 = T_1 - T_1^\dagger`, while non-unitary ansätze use :math:`T_1`. A zero rank disables this transform; Bernoulli expansions reject a positive rank. + .. _cc-hbar-connectivity: Using :math:`\bar{H}` outside the CC class @@ -161,10 +163,12 @@ Using :math:`\bar{H}` outside the CC class // WRONG: empty connectivity keeps terms the commutator would have cancelled auto bad = op::ref_av(op::P(nₚ(2)) * cc.hbar(), {.connect = {}}); - // also correct: ask lst for the explicit commutator form, which needs - // no connectivity + // also correct when no additional singles transform is configured: ask lst + // for the explicit commutator form, which needs no connectivity auto hbar = lst(op::H(), op::T(2), 4); +When ``hbar_singles_comm_rank`` is positive, reproducing a non-unitary :math:`\bar{H}` explicitly requires a second ``lst`` call for the configured singles transform. + Note that ``op::ref_av(expr)`` and ``op::ref_av(expr, {})`` are *not* the same call: the connectivity defaults to ``default_op_connections()`` only when the argument is omitted entirely, since ``EVOptions::connect`` is itself empty by default. The two ``ref_av`` calls above swap roles for a unitary ansatz, so it is not simply "unaffected". There :func:`CC::hbar() ` already returns explicit commutators, and imposing connectivity on top of them would drop terms that must survive; a caller must therefore pass empty connections explicitly, as the class does internally: diff --git a/tests/unit/test_mbpt_cc.cpp b/tests/unit/test_mbpt_cc.cpp index 14f452643d..533501c970 100644 --- a/tests/unit/test_mbpt_cc.cpp +++ b/tests/unit/test_mbpt_cc.cpp @@ -317,11 +317,23 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { auto bch3 = bch2.with([](auto& o) { o.hbar_comm_rank = 3; }); REQUIRE(bch2.hbar_comm_rank() == 2); REQUIRE(bch3.hbar_comm_rank() == 3); + REQUIRE( + bch2.with([](auto& o) { o.ansatz = CC::Ansatz::oU; }).skip_singles()); // rank 0 is a valid truncation (H̄ = H); only CC::λ rejects it, since it // derives at rank - 1 REQUIRE_THROWS_AS(CC(N, {.ansatz = CC::Ansatz::U}), Exception); REQUIRE(bch2.with([](auto& o) { o.hbar_comm_rank = 0; }).hbar_comm_rank() == 0); + REQUIRE_THROWS_AS(bch2.with([](auto& o) { + o.ansatz = CC::Ansatz::oU; + o.skip_singles = false; + }), + Exception); + REQUIRE_THROWS_AS(CC(N, {.ansatz = CC::Ansatz::U, + .hbar_comm_rank = 2, + .hbar_singles_comm_rank = 1, + .hbar_expansion = CC::HbarExpansion::Bernoulli}), + Exception); if (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { REQUIRE_THROWS_AS(CC(N, {.hbar_comm_rank = 0}).λ(), Exception); } @@ -433,12 +445,6 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { return r; }; - SECTION("K=0/unset reproduces plain hbar") { - CC::Options k0 = base; - k0.hbar_singles_comm_rank = 0; - REQUIRE_THAT(CC(N, k0).hbar(), EquivalentTo(hbar_R)); - } - SECTION("K=1: extra term is exactly [H̄_R, σ1]") { CC::Options k1 = base; k1.hbar_singles_comm_rank = 1; @@ -447,7 +453,6 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { // the wrap genuinely adds terms not present in the baseline const auto extra = comm(hbar_R, sigma1); REQUIRE(size(extra) > 0); - REQUIRE(size(hbar_k1) > size(hbar_R)); REQUIRE_THAT(hbar_k1, EquivalentTo(hbar_R + extra)); } @@ -462,42 +467,6 @@ TEST_CASE("mbpt_cc", "[mbpt/cc][valgrind_skip]") { } } // SECTION("extra singles commutators") - SECTION("copy-and-mutate") { - const auto N = 2; - const CC::Options base{.ansatz = CC::Ansatz::U, .hbar_comm_rank = 2}; - - SECTION("identity mutation reproduces the engine") { - const CC cc(N, base); - REQUIRE_THAT(cc.with([](auto&) {}).hbar(), EquivalentTo(cc.hbar())); - } - - SECTION("changed rank matches fresh construction") { - const CC cc(N, base); - CC::Options at3 = base; - at3.hbar_comm_rank = 3; - // N.B. cannot compare against with_hbar_comm_rank -- this commit deletes - // it - REQUIRE_THAT(cc.with([](auto& o) { o.hbar_comm_rank = 3; }).hbar(), - EquivalentTo(CC(N, at3).hbar())); - } - - SECTION("skip_singles is not pinned by the round trip") { - // base leaves skip_singles unset, so CC resolves it to false for U. - // Flipping to oU must pick up the oU default (true), not carry that - // resolved false into an ansatz whose ctor requires true. - const CC cc(N, base); - const auto ou = cc.with([](auto& o) { o.ansatz = CC::Ansatz::oU; }); - REQUIRE(ou.skip_singles()); - } - - SECTION("mutation does not disturb the source engine") { - const CC cc(N, base); - const auto hbar_before = cc.hbar(); - (void)cc.with([](auto& o) { o.hbar_comm_rank = 4; }); - REQUIRE_THAT(cc.hbar(), EquivalentTo(hbar_before)); - } - } // SECTION("copy-and-mutate") - SECTION("eom_cc"){SECTION("EOM-CCSD"){const auto N = 2; auto cc = CC{N}; SECTION("EE-EOM-CCSD R") {