diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index 144a082731..4a42bf7f93 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -37,59 +37,69 @@ 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), - 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 (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 (ansatz_ == Ansatz::oT || 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 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_; } +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 // 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()); + + // 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); + } + return result; } 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 @@ -106,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) { @@ -134,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)) { @@ -168,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 @@ -201,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()); @@ -233,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 @@ -251,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()); @@ -275,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. @@ -337,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 @@ -346,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 @@ -411,26 +422,39 @@ 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 && + opts_.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, " @@ -439,7 +463,7 @@ std::vector CC::eom_r_blocked( // 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. - const bool tensor_level = hbar_expansion_ == HbarExpansion::Bernoulli; + 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 @@ -465,11 +489,10 @@ std::vector CC::eom_r_blocked( 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}); }; - 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 +502,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 +520,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() || opts_.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..73cb6dbb0c 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 { @@ -39,6 +40,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; @@ -57,6 +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; + /// 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 @@ -74,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 @@ -91,6 +107,16 @@ class CC { /// @return the choice of H̄ expansion [[nodiscard]] HbarExpansion hbar_expansion() const; + /// @brief returns a validated copy with @p mutate applied to its options + /// @param mutate callable taking `Options&` + /// @return the mutated copy + 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; @@ -108,26 +134,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 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; @@ -194,38 +217,41 @@ 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). 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` + /// 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 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` + /// 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 = {}) 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 @@ -263,24 +289,19 @@ 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 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` + Options opts_; + + /// @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..2f40b19d22 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 -------- @@ -137,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 @@ -156,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 9b18cf211a..533501c970 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 @@ -281,7 +282,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)); @@ -308,17 +310,34 @@ 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); + 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); } - } // SECTION("hbar_comm_rank") + } // SECTION("with") SECTION("rdm") { constexpr auto N = 2; @@ -410,6 +429,44 @@ 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=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_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") { @@ -442,6 +499,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; @@ -485,7 +581,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);