diff --git a/CMakeLists.txt b/CMakeLists.txt index 11a020e53f..b408d9b309 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -373,6 +373,9 @@ set(SeQuant_src SeQuant/domain/mbpt/convention.hpp SeQuant/domain/mbpt/models/cc.cpp SeQuant/domain/mbpt/models/cc.hpp + SeQuant/domain/mbpt/reserved.hpp + SeQuant/domain/mbpt/op_registry.hpp + SeQuant/domain/mbpt/op_registry.cpp SeQuant/domain/mbpt/op.cpp SeQuant/domain/mbpt/op.hpp SeQuant/domain/mbpt/op.ipp diff --git a/SeQuant/domain/mbpt/context.cpp b/SeQuant/domain/mbpt/context.cpp index 6cb671ef2a..49bd20d36f 100644 --- a/SeQuant/domain/mbpt/context.cpp +++ b/SeQuant/domain/mbpt/context.cpp @@ -2,10 +2,52 @@ namespace sequant::mbpt { -Context::Context(CSV csv) noexcept : csv_(csv) {} +/// forward declaration of OpClass +enum class OpClass; + +Context::Context(Options options) + : csv_(options.csv), + op_registry_(options.op_registry_ptr + ? std::move(options.op_registry_ptr) + : (options.op_registry + ? std::make_shared( + std::move(options.op_registry.value())) + : nullptr)) {} + +Context Context::clone() const { + Context ctx(*this); + ctx.op_registry_ = std::make_shared(op_registry_->clone()); + return ctx; +} + +CSV Context::csv() const { return csv_; } + +std::shared_ptr Context::op_registry() const { + return op_registry_; +} + +std::shared_ptr Context::mutable_op_registry() { + return op_registry_; +} + +Context& Context::set(const OpRegistry& op_registry) { + op_registry_ = std::make_shared(op_registry); + return *this; +} + +Context& Context::set(std::shared_ptr op_registry) { + op_registry_ = std::move(op_registry); + return *this; +} + +Context& Context::set(CSV csv) { + csv_ = csv; + return *this; +} bool operator==(Context const& left, Context const& right) { - return left.csv() == right.csv(); + return left.csv() == right.csv() && + *(left.op_registry()) == *(right.op_registry()); } bool operator!=(Context const& left, Context const& right) { @@ -20,6 +62,10 @@ void set_default_mbpt_context(const Context& ctx) { sequant::detail::set_implicit_context(ctx); } +void set_default_mbpt_context(const Context::Options& options) { + return set_default_mbpt_context(Context(options)); +} + void reset_default_mbpt_context() { sequant::detail::reset_implicit_context(); } @@ -29,4 +75,49 @@ set_scoped_default_mbpt_context(const Context& f) { return sequant::detail::set_scoped_implicit_context(f); } +[[nodiscard]] sequant::detail::ImplicitContextResetter +set_scoped_default_mbpt_context(const Context::Options& f) { + return sequant::detail::set_scoped_implicit_context(Context(f)); +} + +std::shared_ptr make_minimal_registry() { + auto registry = std::make_shared(); + + registry + ->add(L"h", OpClass::gen) /// 1-body Hamiltonian + .add(L"g", OpClass::gen) /// 2-body Coulomb + .add(L"f", OpClass::gen) /// Fock operator + .add(L"θ", OpClass::gen) /// general fock space operator + .add(L"t", OpClass::ex) /// cluster operator + .add(L"λ", OpClass::deex) /// deexcitation cluster operator + .add(L"R", OpClass::ex) /// right-hand eigenstate + .add(L"L", OpClass::deex); /// left-hand eigenstate + + return registry; +} + +std::shared_ptr make_legacy_registry() { + auto registry = std::make_shared(); + + registry->add(L"h", OpClass::gen) + .add(L"f", OpClass::gen) + /// closed Fock operator (i.e. Fock operator due to fully-occupied + /// orbitals) + .add(L"f̃", OpClass::gen) + .add(L"g", OpClass::gen) + .add(L"θ", OpClass::gen) + .add(L"t", OpClass::ex) + .add(L"λ", OpClass::deex) + .add(L"R", OpClass::ex) + .add(L"L", OpClass::deex) + /// R12 + .add(L"F", OpClass::gen) + .add(L"GR", OpClass::gen) + .add(L"C", OpClass::gen) + /// RDM and RDM Cumulant + .add(L"γ", OpClass::gen) + .add(L"κ", OpClass::gen); + + return registry; +} } // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/context.hpp b/SeQuant/domain/mbpt/context.hpp index a1758347fb..3d2809d3c3 100644 --- a/SeQuant/domain/mbpt/context.hpp +++ b/SeQuant/domain/mbpt/context.hpp @@ -5,6 +5,7 @@ #include #include +#include namespace sequant::mbpt { @@ -20,13 +21,61 @@ class Context { constexpr static auto csv = CSV::No; }; - /// @brief Construct a Formalism object - explicit Context(CSV csv_formalism = Defaults::csv) noexcept; + struct Options { + /// whether to use cluster-specific virtuals + CSV csv = Defaults::csv; + /// shared pointer to operator registry + std::shared_ptr op_registry_ptr = nullptr; + /// optional operator registry, use if no shared pointer is provided + std::optional op_registry = std::nullopt; + }; + + /// @brief makes default options for mbpt::Context + static Options make_default_options() { return Options{}; } + + /// @brief Construct a Context object, uses default options if none are given + Context(Options options = make_default_options()); + + /// @brief destructor + ~Context() = default; + + /// @brief move constructor + Context(Context&&) = default; + + /// @brief copy constructor + Context(Context const& other) = default; + + /// @brief copy assignment + Context& operator=(Context const& other) = default; + + /// @brief clones this object and its OpRegistry + Context clone() const; + + /// @return the value of CSV in this context + CSV csv() const; - CSV csv() const { return csv_; } + /// @return a constant pointer to the OpRegistry for this context + /// @warning can be null when user did not provide one to Context (i.e., it + /// was default constructed) + std::shared_ptr op_registry() const; + + /// @return a pointer to the OpRegistry for this context + /// @warning can be null when user did not provide one to Context (i.e., it + /// was default constructed) + std::shared_ptr mutable_op_registry(); + + /// @brief sets the OpRegistry for this context + Context& set(const OpRegistry& op_registry); + + /// @brief sets the OpRegistry for this context + Context& set(std::shared_ptr op_registry); + + /// @brief sets whether to use cluster-specific virtuals + Context& set(CSV csv); private: CSV csv_ = Defaults::csv; + std::shared_ptr op_registry_; }; bool operator==(Context const& left, Context const& right); @@ -37,10 +86,39 @@ const Context& get_default_mbpt_context(); void set_default_mbpt_context(const Context& ctx); +void set_default_mbpt_context(const Context::Options& options); + void reset_default_mbpt_context(); + [[nodiscard]] sequant::detail::ImplicitContextResetter set_scoped_default_mbpt_context(const Context& ctx); +[[nodiscard]] sequant::detail::ImplicitContextResetter +set_scoped_default_mbpt_context(const Context::Options& options); + +/// predefined operator registries + +/// @brief makes a minimal operator registry with only essential operators for +/// MBPT +std::shared_ptr make_minimal_registry(); + +/// @brief make a legacy operator registry with SeQuant's old predefined +/// operators set +std::shared_ptr make_legacy_registry(); + +/// @brief converts an operator label to its OpClass using the default MBPT +/// context +/// @param op the operator label +/// @return the OpClass of the operator +inline OpClass to_op_class(const std::wstring& op) { + // check reserved labels first + if (ranges::contains(reserved::labels(), op)) { + return OpClass::gen; // all reserved labels are gen + } else { + return get_default_mbpt_context().op_registry()->to_class(op); + } +} + } // namespace sequant::mbpt #endif // SEQUANT_DOMAIN_MBPT_CONTEXT_HPP diff --git a/SeQuant/domain/mbpt/models/cc.cpp b/SeQuant/domain/mbpt/models/cc.cpp index def47f38f4..68230e8d8d 100644 --- a/SeQuant/domain/mbpt/models/cc.cpp +++ b/SeQuant/domain/mbpt/models/cc.cpp @@ -93,12 +93,12 @@ std::vector CC::λ(size_t commutator_rank) { auto lhbar = simplify((One + Λ(N)) * hbar); const auto op_connect = concat(default_op_connections(), - OpConnections{{OpType::h, OpType::A}, - {OpType::f, OpType::A}, - {OpType::g, OpType::A}, - {OpType::h, OpType::S}, - {OpType::f, OpType::S}, - {OpType::g, OpType::S}}); + OpConnections{{L"h", L"A"}, + {L"f", L"A"}, + {L"g", L"A"}, + {L"h", L"S"}, + {L"f", L"S"}, + {L"g", L"S"}}); // 2. project onto each manifold, screen, lower to tensor form and wick it std::vector result(N + 1); @@ -168,10 +168,8 @@ std::vector CC::t_pt(size_t rank, size_t order, // connect h1 with t const auto op_connect = concat(default_op_connections(), - OpConnections{{OpType::h, OpType::t_1}, - {OpType::f, OpType::t_1}, - {OpType::g, OpType::t_1}, - {OpType::h_1, OpType::t}}); + OpConnections{ + {L"h", L"t¹"}, {L"f", L"t¹"}, {L"g", L"t¹"}, {L"h¹", L"t"}}); std::vector result(N + 1); for (auto p = N; p >= 1; --p) { @@ -218,20 +216,19 @@ std::vector CC::λ_pt(size_t rank, size_t order, // projectors with {h,f,g} // h1 with t // h1 with projectors - const auto op_connect = - concat(default_op_connections(), - OpConnections{{OpType::h, OpType::t_1}, - {OpType::f, OpType::t_1}, - {OpType::g, OpType::t_1}, - {OpType::h_1, OpType::t}, - {OpType::h, OpType::A}, - {OpType::f, OpType::A}, - {OpType::g, OpType::A}, - {OpType::h, OpType::S}, - {OpType::f, OpType::S}, - {OpType::g, OpType::S}, - {OpType::h_1, OpType::A}, - {OpType::h_1, OpType::S}}); + const auto op_connect = concat(default_op_connections(), + OpConnections{{L"h", L"t¹"}, + {L"f", L"t¹"}, + {L"g", L"t¹"}, + {L"h¹", L"t"}, + {L"h", L"A"}, + {L"f", L"A"}, + {L"g", L"A"}, + {L"h", L"S"}, + {L"f", L"S"}, + {L"g", L"S"}, + {L"h¹", L"A"}, + {L"h¹", L"S"}}); std::vector result(N + 1); for (auto p = N; p >= 1; --p) { @@ -262,10 +259,9 @@ std::vector CC::eom_r(nₚ np, nₕ nh) { // connectivity: // default connections + connect R with {h,f,g} - const auto op_connect = concat(default_op_connections(), - OpConnections{{OpType::h, OpType::R}, - {OpType::f, OpType::R}, - {OpType::g, OpType::R}}); + const auto op_connect = concat( + default_op_connections(), + OpConnections{{L"h", L"R"}, {L"f", L"R"}, {L"g", L"R"}}); // initialize result vector std::vector result; @@ -305,12 +301,12 @@ std::vector CC::eom_l(nₚ np, nₕ nh) { // connectivity: // default connections + connect H with projectors const auto op_connect = concat(default_op_connections(), - OpConnections{{OpType::h, OpType::A}, - {OpType::f, OpType::A}, - {OpType::g, OpType::A}, - {OpType::h, OpType::S}, - {OpType::f, OpType::S}, - {OpType::g, OpType::S}}); + OpConnections{{L"h", L"A"}, + {L"f", L"A"}, + {L"g", L"A"}, + {L"h", L"S"}, + {L"f", L"S"}, + {L"g", L"S"}}); // initialize result vector std::vector result; diff --git a/SeQuant/domain/mbpt/models/cc.hpp b/SeQuant/domain/mbpt/models/cc.hpp index 5753428c88..157c050e93 100644 --- a/SeQuant/domain/mbpt/models/cc.hpp +++ b/SeQuant/domain/mbpt/models/cc.hpp @@ -129,7 +129,7 @@ class CC { /// @param[in] op_connect list of pairs of operators to be connected. Default /// is given by `mbpt::default_op_connections()`. auto ref_av(const ExprPtr& expr, - const OpConnections& op_connect = + const OpConnections& op_connect = default_op_connections()) const { return op::ref_av(expr, op_connect, this->use_topology(), this->screen()); } diff --git a/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp index adde256dc1..797648a632 100644 --- a/SeQuant/domain/mbpt/op.cpp +++ b/SeQuant/domain/mbpt/op.cpp @@ -46,42 +46,6 @@ std::vector cardinal_tensor_labels() { L"E"}; } -std::wstring to_wstring(OpType op) { - auto found_it = optype2label.find(op); - if (found_it != optype2label.end()) - return found_it->second; - else - throw std::invalid_argument("to_wstring(OpType op): invalid op"); -} - -OpClass to_class(OpType op) { - switch (op) { - case OpType::h: - case OpType::f: - case OpType::f̃: - case OpType::g: - case OpType::RDM: - case OpType::RDMCumulant: - case OpType::δ: - case OpType::A: - case OpType::S: - case OpType::h_1: - case OpType::θ: - return OpClass::gen; - case OpType::t: - case OpType::R: - case OpType::R12: - case OpType::t_1: - return OpClass::ex; - case OpType::λ: - case OpType::L: - case OpType::λ_1: - return OpClass::deex; - default: - throw std::invalid_argument("to_class(OpType op): invalid op"); - } -} - // Excitation type QNs will have quasiparticle annihilators in every space which // intersects with the active hole space. The active particle space will have // quasiparticle creators. The presence of additional blocks depends on whether @@ -342,7 +306,11 @@ template std::wstring to_latex(const mbpt::Operator& op) { using namespace sequant::mbpt; - auto result = L"{\\hat{" + utf_to_latex(op.label()) + L"}"; + // start with label and perturbation order (if any) + auto result = + L"{\\hat{" + + utf_to_latex(mbpt::decorate_with_pert_order(op.label(), op.order())) + + L"}"; // check if operator has adjoint label, remove if present for base label auto base_lbl = sequant::to_wstring(op.label()); @@ -354,18 +322,24 @@ std::wstring to_latex(const mbpt::Operator& op) { auto op_qns = op(); // operator action i.e. quantum number change - auto it = label2optype.find(base_lbl); // look for OpType - const bool known_optype = it != label2optype.end(); + auto registry = mbpt::get_default_mbpt_context().op_registry(); + // if it is not a reserved label, make sure it is registered + if (reserved::is_nonreserved(base_lbl)) { + SEQUANT_ASSERT(registry->contains(base_lbl) && + "to_latex(mbpt::Operator): " + "unregistered operator label"); + } + // find the `class` of Operator + OpClass opclass = mbpt::to_op_class(base_lbl); // special handling for general operators // - Ops like f and g does not need ranks, it is implied // - Ops like A, S, θ are general, but need rank information // - θ needs to be treated differently because it can have variable number of // quantum numbers - - auto skip_rank_info = [](const OpType& optype) { - return to_class(optype) == OpClass::gen && - !(optype == OpType::θ || optype == OpType::A || optype == OpType::S); + auto skip_rank_info = [opclass](const auto& label) { + return opclass == OpClass::gen && label != reserved::antisymm_label() && + label != reserved::symm_label() && label != L"θ"; }; // batch index handling @@ -385,12 +359,12 @@ std::wstring to_latex(const mbpt::Operator& op) { return str; }; - if (known_optype && skip_rank_info(it->second)) { + if (skip_rank_info(base_lbl)) { result += L"}"; // close the brace return has_batching ? add_batch_suffix(result) : result; } // specially handle θ operator - if (known_optype && it->second == OpType::θ) { + if (base_lbl == L"θ") { result += L"_{" + std::to_wstring(op_qns[0].upper()) + L"}"; result += L"}"; // close the brace return has_batching ? add_batch_suffix(result) : result; @@ -421,7 +395,8 @@ std::wstring to_latex(const mbpt::Operator& op) { // check if the Op is a projector (A or S) // projectors can have negative ranks, need special handling [[maybe_unused]] const bool is_projector = - known_optype && (it->second == OpType::A || it->second == OpType::S); + base_lbl == reserved::antisymm_label() || + base_lbl == reserved::symm_label(); // pure quasiparticle creator/annihilator? const auto qprank_cre = ncre_p.lower() + nann_h.lower(); @@ -490,13 +465,15 @@ sequant::container::svector make_batch_indices( namespace sequant::mbpt { template -OpMaker::OpMaker(OpType op) : op_(op) {} +OpMaker::OpMaker(const std::wstring& label) : label_(label) {} template -OpMaker::OpMaker(OpType op, ncre nc, nann na) { - op_ = op; +OpMaker::OpMaker(const std::wstring& label, ncre nc, nann na) { + label_ = label; SEQUANT_ASSERT(nc > 0 || na > 0); - switch (to_class(op)) { + auto registry = get_default_mbpt_context().op_registry(); + + switch (registry->to_class(label_)) { case OpClass::ex: cre_spaces_ = IndexSpaceContainer(nc, get_particle_space(Spin::any)); ann_spaces_ = IndexSpaceContainer(na, get_hole_space(Spin::any)); @@ -513,24 +490,28 @@ OpMaker::OpMaker(OpType op, ncre nc, nann na) { } template -OpMaker::OpMaker(OpType op, std::size_t rank) - : OpMaker(op, ncre(rank), nann(rank)) {} +OpMaker::OpMaker(const std::wstring& label, std::size_t rank) + : OpMaker(label, ncre(rank), nann(rank)) {} template -OpMaker::OpMaker(OpType op, ncre nc, nann na, +OpMaker::OpMaker(const std::wstring& label, ncre nc, nann na, const cre& cre_space, const ann& ann_space) { - op_ = op; + label_ = label; SEQUANT_ASSERT(nc > 0 || na > 0); cre_spaces_ = IndexSpaceContainer(nc, cre_space); ann_spaces_ = IndexSpaceContainer(na, ann_space); } template -OpMaker::OpMaker(OpType op, ncre nc, nann na, const OpParams& params) - : OpMaker(op, nc, na) { +OpMaker::OpMaker(const std::wstring& label, ncre nc, nann na, + const OpParams& params) + : OpMaker(label, nc, na) { params.validate(); + // set perturbation order + order_ = params.order; + // Handle batching indices if specified if (!params.batch_ordinals.empty()) { SEQUANT_ASSERT(ranges::is_sorted(params.batch_ordinals) && @@ -560,17 +541,21 @@ template ExprPtr OpMaker::operator()(std::optional dep, std::optional opsymm_opt) const { auto isr = get_default_context(Statistics::FermiDirac).index_space_registry(); + // if not given dep, use mbpt::Context::CSV to determine whether to use // dependent indices for pure (de)excitation ops - if (!dep && get_default_mbpt_context().csv() == mbpt::CSV::Yes) { - if (to_class(op_) == OpClass::ex) { + const auto csv = get_default_mbpt_context().csv() == mbpt::CSV::Yes; + const auto opclass = mbpt::to_op_class(label_); + + if (!dep && csv) { + if (opclass == OpClass::ex) { #ifdef SEQUANT_ASSERT_ENABLED for (auto&& s : cre_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); } #endif dep = UseDepIdx::Bra; - } else if (to_class(op_) == OpClass::deex) { + } else if (opclass == OpClass::deex) { #ifdef SEQUANT_ASSERT_ENABLED for (auto&& s : ann_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); @@ -581,14 +566,15 @@ ExprPtr OpMaker::operator()(std::optional dep, dep = UseDepIdx::None; } } + const std::wstring full_label = decorate_with_pert_order(label_, order_); // if batching indices are present, use them if (batch_indices_) { return make( cre_spaces_, ann_spaces_, batch_indices_.value(), - [this, opsymm_opt](const auto& creidxs, const auto& annidxs, - const auto& batchidxs, Symmetry opsymm) { - return ex(to_wstring(op_), bra(creidxs), ket(annidxs), + [this, opsymm_opt, full_label](const auto& creidxs, const auto& annidxs, + const auto& batchidxs, Symmetry opsymm) { + return ex(full_label, bra(creidxs), ket(annidxs), aux(batchidxs), opsymm_opt ? *opsymm_opt : opsymm); }, dep ? *dep : UseDepIdx::None); @@ -596,9 +582,9 @@ ExprPtr OpMaker::operator()(std::optional dep, // else no batching return make( cre_spaces_, ann_spaces_, - [this, opsymm_opt](const auto& creidxs, const auto& annidxs, - Symmetry opsymm) { - return ex(to_wstring(op_), bra(creidxs), ket(annidxs), + [this, opsymm_opt, full_label](const auto& creidxs, const auto& annidxs, + Symmetry opsymm) { + return ex(full_label, bra(creidxs), ket(annidxs), opsymm_opt ? *opsymm_opt : opsymm); }, dep ? *dep : UseDepIdx::None); @@ -615,20 +601,25 @@ inline namespace op { namespace tensor { ExprPtr H_(std::size_t k) { SEQUANT_ASSERT(k > 0 && k <= 2); + auto registry = get_default_mbpt_context().op_registry(); switch (k) { case 1: switch (get_default_context().vacuum()) { case Vacuum::Physical: - return OpMaker(OpType::h, 1)(); + SEQUANT_ASSERT(registry->contains(L"h")); + return OpMaker(L"h", 1)(); case Vacuum::SingleProduct: - return OpMaker(OpType::f, 1)(); + SEQUANT_ASSERT(registry->contains(L"f")); + return OpMaker(L"f", 1)(); case Vacuum::MultiProduct: - return OpMaker(OpType::f, 1)(); + SEQUANT_ASSERT(registry->contains(L"f")); + return OpMaker(L"f", 1)(); } SEQUANT_UNREACHABLE; case 2: - return OpMaker(OpType::g, 2)(); + SEQUANT_ASSERT(registry->contains(L"g")); + return OpMaker(L"g", 2)(); } SEQUANT_ABORT("Unhandled k value"); @@ -639,15 +630,17 @@ ExprPtr H(std::size_t k) { return k == 1 ? tensor::H_(1) : tensor::H_(1) + tensor::H_(2); } -ExprPtr F(bool use_tensor, IndexSpace reference_occupied) { +ExprPtr F(bool use_tensor, const IndexSpace& reference_occupied) { + auto registry = get_default_mbpt_context().op_registry(); if (use_tensor) { - return OpMaker(OpType::f, 1)(); + SEQUANT_ASSERT(registry->contains(L"f")); + return OpMaker(L"f", 1)(); } else { // explicit density matrix construction SEQUANT_ASSERT( reference_occupied); // cannot explicitly instantiate fock operator // without providing an occupied indexspace // add \bar{g}^{\kappa x}_{\lambda y} \gamma^y_x with x,y in occ_space_type - auto make_g_contribution = [](const auto occ_space) { + auto make_g_contribution = [](const auto& occ_space) { auto isr = get_default_context().index_space_registry(); return mbpt::OpMaker::make( {isr->complete_space(Spin::any)}, {isr->complete_space(Spin::any)}, @@ -659,10 +652,9 @@ ExprPtr F(bool use_tensor, IndexSpace reference_occupied) { if (opsymm == Symmetry::Antisymm) { braidxs.push_back(m1); ketidxs.push_back(m2); - return ex(to_wstring(mbpt::OpType::g), - bra(std::move(braidxs)), + return ex(L"g", bra(std::move(braidxs)), ket(std::move(ketidxs)), Symmetry::Antisymm) * - ex(to_wstring(mbpt::OpType::δ), bra{m2}, ket{m1}, + ex(kronecker_label(), bra{m2}, ket{m1}, Symmetry::Nonsymm); } else { // opsymm == Symmetry::Nonsymm auto braidx_J = braidxs; @@ -674,33 +666,36 @@ ExprPtr F(bool use_tensor, IndexSpace reference_occupied) { auto ketidxs_K = ketidxs; using std::begin; ketidxs_K.emplace(begin(ketidxs_K), m2); - return (ex(to_wstring(mbpt::OpType::g), - bra(std::move(braidx_J)), + return (ex(L"g", bra(std::move(braidx_J)), ket(std::move(ketidxs_J)), Symmetry::Nonsymm) - - ex( - to_wstring(mbpt::OpType::g), bra(std::move(braidx_K)), - ket(std::move(ketidxs_K)), Symmetry::Nonsymm)) * - ex(to_wstring(mbpt::OpType::δ), bra{m2}, ket{m1}, + ex(L"g", bra(std::move(braidx_K)), + ket(std::move(ketidxs_K)), + Symmetry::Nonsymm)) * + ex(kronecker_label(), bra{m2}, ket{m1}, Symmetry::Nonsymm); } }); }; auto isr = get_default_context().index_space_registry(); - return OpMaker(OpType::h, 1)() + + SEQUANT_ASSERT(registry->contains(L"h")); + return OpMaker(L"h", 1)() + make_g_contribution(reference_occupied); } } ExprPtr θ(std::size_t K) { - return OpMaker(OpType::θ, K)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"θ")); + return OpMaker(L"θ", K)(); } ExprPtr T_(std::size_t K) { - return OpMaker(OpType::t, K)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); + return OpMaker(L"t", K)(); } ExprPtr T(std::size_t K, bool skip1) { SEQUANT_ASSERT(K > (skip1 ? 1 : 0)); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); ExprPtr result; for (auto k = skip1 ? 2ul : 1ul; k <= K; ++k) { result += tensor::T_(k); @@ -709,12 +704,13 @@ ExprPtr T(std::size_t K, bool skip1) { } ExprPtr Λ_(std::size_t K) { - return OpMaker(OpType::λ, K)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); + return OpMaker(L"λ", K)(); } ExprPtr Λ(std::size_t K, bool skip1) { SEQUANT_ASSERT(K > (skip1 ? 1 : 0)); - + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); ExprPtr result; for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) { result = k > 1 ? result + tensor::Λ_(k) : tensor::Λ_(k); @@ -724,24 +720,26 @@ ExprPtr Λ(std::size_t K, bool skip1) { ExprPtr R_(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { - return OpMaker(OpType::R, nc, na, cre_space, - ann_space)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R")); + return OpMaker(L"R", nc, na, cre_space, ann_space)(); } ExprPtr R_(nₚ np, nₕ nh) { SEQUANT_ASSERT(np >= 0 && nh >= 0); - return OpMaker(OpType::R, ncre(np.value()), + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R")); + return OpMaker(L"R", ncre(np.value()), nann(nh.value()))(); } ExprPtr L_(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { - return OpMaker(OpType::L, nc, na, cre_space, - ann_space)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L")); + return OpMaker(L"L", nc, na, cre_space, ann_space)(); } ExprPtr L_(nₚ np, nₕ nh) { SEQUANT_ASSERT(np >= 0 && nh >= 0); - return OpMaker(OpType::L, ncre(nh.value()), + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L")); + return OpMaker(L"L", ncre(nh.value()), nann(np.value()))(); } @@ -782,8 +780,9 @@ ExprPtr A(nₚ np, nₕ nh) { if (get_default_mbpt_context().csv() == mbpt::CSV::Yes) dep = (np > 0 || nh > 0) ? OpMaker::UseDepIdx::Bra : OpMaker::UseDepIdx::Ket; - return OpMaker( - OpType::A, cre(creators), ann(annihilators))(dep, {Symmetry::Antisymm}); + return OpMaker(reserved::antisymm_label(), + cre(creators), ann(annihilators))( + dep, {Symmetry::Antisymm}); } ExprPtr S(std::int64_t K) { @@ -807,18 +806,16 @@ ExprPtr S(std::int64_t K) { if (get_default_mbpt_context().csv() == mbpt::CSV::Yes) dep = K > 0 ? OpMaker::UseDepIdx::Bra : OpMaker::UseDepIdx::Ket; - return OpMaker( - OpType::S, cre(creators), ann(annihilators))(dep, {Symmetry::Nonsymm}); + return OpMaker(reserved::symm_label(), cre(creators), + ann(annihilators))( + dep, {Symmetry::Nonsymm}); } ExprPtr H_pt(std::size_t R, const OpParams& params) { params.validate(); - SEQUANT_ASSERT(params.order == 1 && - "sequant::mbpt::tensor::H_pt: only supports first " - "order perturbation"); SEQUANT_ASSERT(R > 0 && "Operator rank must be > 0"); - return OpMaker(OpType::h_1, ncre(R), nann(R), - params)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"h")); + return OpMaker(L"h", ncre(R), nann(R), params)(); } ExprPtr T_pt_(std::size_t K, const OpParams& params) { @@ -827,13 +824,14 @@ ExprPtr T_pt_(std::size_t K, const OpParams& params) { "sequant::mbpt::tensor::T_pt_: only supports first " "order perturbation"); SEQUANT_ASSERT(K > 0 && "Operator rank must be > 0"); - return OpMaker(OpType::t_1, ncre(K), nann(K), - params)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); + return OpMaker(L"t", ncre(K), nann(K), params)(); } ExprPtr T_pt(std::size_t K, const OpParams& params) { params.validate(); if (params.skip1) SEQUANT_ASSERT(K > 1); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); ExprPtr result; for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) { result += tensor::T_pt_(k, {.order = params.order, @@ -846,17 +844,15 @@ ExprPtr T_pt(std::size_t K, const OpParams& params) { ExprPtr Λ_pt_(std::size_t K, const OpParams& params) { params.validate(); - SEQUANT_ASSERT(params.order == 1 && - "sequant::mbpt::tensor::Λ_pt_: only supports first " - "order perturbation"); SEQUANT_ASSERT(K > 0 && "Operator rank must be > 0"); - return OpMaker(OpType::λ_1, ncre(K), nann(K), - params)(); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); + return OpMaker(L"λ", ncre(K), nann(K), params)(); } ExprPtr Λ_pt(std::size_t K, const OpParams& params) { params.validate(); if (params.skip1) SEQUANT_ASSERT(K > 1); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); ExprPtr result; for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) { result += tensor::Λ_pt_(k, {.order = params.order, @@ -871,16 +867,21 @@ ExprPtr Λ_pt(std::size_t K, const OpParams& params) { ExprPtr H_(std::size_t k) { SEQUANT_ASSERT(k > 0 && k <= 2); + auto registry = get_default_mbpt_context().op_registry(); switch (k) { case 1: return ex( - [vacuum = get_default_context().vacuum()]() -> std::wstring_view { + [vacuum = get_default_context().vacuum(), + registry]() -> std::wstring_view { switch (vacuum) { case Vacuum::Physical: + SEQUANT_ASSERT(registry->contains(L"h")); return L"h"; case Vacuum::SingleProduct: + SEQUANT_ASSERT(registry->contains(L"f")); return L"f"; case Vacuum::MultiProduct: + SEQUANT_ASSERT(registry->contains(L"f")); return L"f"; } @@ -893,6 +894,7 @@ ExprPtr H_(std::size_t k) { }); case 2: + SEQUANT_ASSERT(registry->contains(L"g")); return ex([]() -> std::wstring_view { return L"g"; }, [=]() -> ExprPtr { return tensor::H_(2); }, [=](qnc_t& qns) { @@ -921,6 +923,7 @@ ExprPtr θ(std::size_t K) { ExprPtr T_(std::size_t K) { SEQUANT_ASSERT(K > 0); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); return ex([]() -> std::wstring_view { return L"t"; }, [=]() -> ExprPtr { return tensor::T_(K); }, [=](qnc_t& qns) { @@ -931,6 +934,7 @@ ExprPtr T_(std::size_t K) { ExprPtr T(std::size_t K, bool skip1) { SEQUANT_ASSERT(K > (skip1 ? 1 : 0)); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); ExprPtr result; for (auto k = skip1 ? 2ul : 1ul; k <= K; ++k) { result += T_(k); @@ -940,6 +944,7 @@ ExprPtr T(std::size_t K, bool skip1) { ExprPtr Λ_(std::size_t K) { SEQUANT_ASSERT(K > 0); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); return ex([]() -> std::wstring_view { return L"λ"; }, [=]() -> ExprPtr { return tensor::Λ_(K); }, [=](qnc_t& qns) { @@ -950,6 +955,7 @@ ExprPtr Λ_(std::size_t K) { ExprPtr Λ(std::size_t K, bool skip1) { SEQUANT_ASSERT(K > (skip1 ? 1 : 0)); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); ExprPtr result; for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) { result = k > 1 ? result + Λ_(k) : Λ_(k); @@ -957,7 +963,8 @@ ExprPtr Λ(std::size_t K, bool skip1) { return result; } -ExprPtr F(bool use_f_tensor, IndexSpace occupied_density) { +ExprPtr F(bool use_f_tensor, const IndexSpace& occupied_density) { + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"f")); if (use_f_tensor) { return ex( []() -> std::wstring_view { return L"f"; }, @@ -982,26 +989,27 @@ ExprPtr A(nₚ np, nₕ nh) { auto particle_space = get_particle_space(Spin::any); auto hole_space = get_hole_space(Spin::any); - return ex([]() -> std::wstring_view { return L"A"; }, - [=]() -> ExprPtr { return tensor::A(np, nh); }, - [=](qnc_t& qns) { - const std::size_t abs_nh = std::abs(nh); - const std::size_t abs_np = std::abs(np); - if (deexcitation) { - qnc_t op_qnc_t = generic_deexcitation_qns( - abs_np, abs_nh, particle_space, hole_space); - qns = combine(op_qnc_t, qns); - } else { - qnc_t op_qnc_t = generic_excitation_qns( - abs_np, abs_nh, particle_space, hole_space); - qns = combine(op_qnc_t, qns); - } - }); + return ex( + []() -> std::wstring_view { return reserved::antisymm_label(); }, + [=]() -> ExprPtr { return tensor::A(np, nh); }, + [=](qnc_t& qns) { + const std::size_t abs_nh = std::abs(nh); + const std::size_t abs_np = std::abs(np); + if (deexcitation) { + qnc_t op_qnc_t = generic_deexcitation_qns(abs_np, abs_nh, + particle_space, hole_space); + qns = combine(op_qnc_t, qns); + } else { + qnc_t op_qnc_t = generic_excitation_qns(abs_np, abs_nh, + particle_space, hole_space); + qns = combine(op_qnc_t, qns); + } + }); } ExprPtr S(std::int64_t K) { SEQUANT_ASSERT(K != 0); - return ex([]() -> std::wstring_view { return L"S"; }, + return ex([]() -> std::wstring_view { return reserved::symm_label(); }, [=]() -> ExprPtr { return tensor::S(K); }, [=](qnc_t& qns) { const std::size_t abs_K = std::abs(K); @@ -1029,48 +1037,20 @@ ExprPtr P(nₚ np, nₕ nh) { } } -namespace { -/// @brief Helper to create Operator from OpParams -/// @param label_gen Callable to generate operator label -/// @param tensor_gen Callable to generate tensor form -/// @param qn_action Callable to apply quantum number changes -/// @param params OpParams containing operator parameters (batching, etc.) -/// @return ExprPtr to the created Operator -ExprPtr make_op_from_params(std::function label_gen, - std::function tensor_gen, - std::function qn_action, - const OpParams& params) { - params.validate(); - if (!params.batch_ordinals.empty()) { - mbpt::check_for_batching_space(); - return ex(label_gen, tensor_gen, qn_action, params.batch_ordinals); - } else if (params.nbatch) { - mbpt::check_for_batching_space(); - return ex(label_gen, tensor_gen, qn_action, params.nbatch.value()); - } else { - return ex(label_gen, tensor_gen, qn_action); - } -} -} // anonymous namespace - ExprPtr H_pt(std::size_t R, const OpParams& params) { SEQUANT_ASSERT(R > 0); - SEQUANT_ASSERT(params.order == 1 && - "only first order perturbation is supported now"); - - return make_op_from_params( - []() -> std::wstring_view { return optype2label.at(OpType::h_1); }, - [R, params]() -> ExprPtr { return tensor::H_pt(R, params); }, - [R](qnc_t& qns) { qns = combine(general_type_qns(R), qns); }, params); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"h")); + return ex([]() -> std::wstring_view { return L"h"; }, + [R, params]() -> ExprPtr { return tensor::H_pt(R, params); }, + [R](qnc_t& qns) { qns = combine(general_type_qns(R), qns); }, + params); } ExprPtr T_pt_(std::size_t K, const OpParams& params) { SEQUANT_ASSERT(K > 0); - SEQUANT_ASSERT(params.order == 1 && - "only first order perturbation is supported now"); - - return make_op_from_params( - []() -> std::wstring_view { return optype2label.at(OpType::t_1); }, + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); + return ex( + []() -> std::wstring_view { return L"t"; }, [K, params]() -> ExprPtr { return tensor::T_pt_(K, params); }, [K](qnc_t& qns) { qns = combine(excitation_type_qns(K), qns); }, params); } @@ -1078,6 +1058,7 @@ ExprPtr T_pt_(std::size_t K, const OpParams& params) { ExprPtr T_pt(std::size_t K, const OpParams& params) { params.validate(); SEQUANT_ASSERT(K > (params.skip1 ? 1 : 0)); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"t")); ExprPtr result; for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) { result += T_pt_(k, {.order = params.order, @@ -1090,11 +1071,9 @@ ExprPtr T_pt(std::size_t K, const OpParams& params) { ExprPtr Λ_pt_(std::size_t K, const OpParams& params) { SEQUANT_ASSERT(K > 0); - SEQUANT_ASSERT(params.order == 1 && - "only first order perturbation is supported now"); - - return make_op_from_params( - []() -> std::wstring_view { return optype2label.at(OpType::λ_1); }, + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); + return ex( + []() -> std::wstring_view { return L"λ"; }, [K, params]() -> ExprPtr { return tensor::Λ_pt_(K, params); }, [K](qnc_t& qns) { qns = combine(deexcitation_type_qns(K), qns); }, params); @@ -1103,6 +1082,7 @@ ExprPtr Λ_pt_(std::size_t K, const OpParams& params) { ExprPtr Λ_pt(std::size_t K, const OpParams& params) { params.validate(); SEQUANT_ASSERT(K > (params.skip1 ? 1 : 0)); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"λ")); ExprPtr result; for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) { result += Λ_pt_(k, {.order = params.order, @@ -1115,8 +1095,9 @@ ExprPtr Λ_pt(std::size_t K, const OpParams& params) { ExprPtr R_(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R")); return ex( - []() -> std::wstring_view { return optype2label.at(OpType::R); }, + []() -> std::wstring_view { return L"R"; }, [=]() -> ExprPtr { return tensor::R_(na, nc, cre_space, ann_space); }, [=](qnc_t& qns) { // ex -> creators in particle_space, annihilators in hole_space @@ -1131,8 +1112,9 @@ ExprPtr R_(nₚ np, nₕ nh) { return R_(nann(nh), ncre(np)); } ExprPtr L_(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L")); return ex( - []() -> std::wstring_view { return optype2label.at(OpType::L); }, + []() -> std::wstring_view { return L"L"; }, [=]() -> ExprPtr { return tensor::L_(na, nc, cre_space, ann_space); }, [=](qnc_t& qns) { // deex -> creators in hole_space, annihilators in particle_space @@ -1148,6 +1130,7 @@ ExprPtr L_(nₚ np, nₕ nh) { return L_(nann(np), ncre(nh)); } ExprPtr R(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { SEQUANT_ASSERT(na > 0 || nc > 0); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"R")); ExprPtr result; std::int64_t ra = na, rc = nc; @@ -1166,6 +1149,7 @@ ExprPtr R(nₚ np, nₕ nh) { return R(nann(nh), ncre(np)); } ExprPtr L(nann na, ncre nc, const cre& cre_space, const ann& ann_space) { SEQUANT_ASSERT(na > 0 || nc > 0); + SEQUANT_ASSERT(get_default_mbpt_context().op_registry()->contains(L"L")); ExprPtr result; std::int64_t ra = na, rc = nc; @@ -1245,7 +1229,7 @@ ExprPtr expectation_value_impl(ExprPtr expr, auto isr = get_default_context().index_space_registry(); const auto spinor = get_default_context().spbasis() == SPBasis::Spinor; // convention is to use different label for spin-orbital and spin-free RDM - const auto rdm_label = spinor ? optype2label.at(OpType::RDM) : L"Γ"; + const auto rdm_label = spinor ? L"γ" : L"Γ"; // N.B. reference < vacuum is not yet supported if (isr->reference_occupied_space().intersection( diff --git a/SeQuant/domain/mbpt/op.hpp b/SeQuant/domain/mbpt/op.hpp index f980e5f321..53823a6baa 100644 --- a/SeQuant/domain/mbpt/op.hpp +++ b/SeQuant/domain/mbpt/op.hpp @@ -8,6 +8,7 @@ #include #include +#include #include #include @@ -51,6 +52,10 @@ namespace sequant { namespace mbpt { +namespace detail { +inline constexpr std::wstring_view pert_superscripts = L"⁰¹²³⁴⁵⁶⁷⁸⁹"; +} // namespace detail + DEFINE_STRONG_TYPE_FOR_INTEGER(nₚ, std::int64_t); // define nₚ DEFINE_STRONG_TYPE_FOR_INTEGER(nₕ, std::int64_t); // define nₕ @@ -69,69 +74,9 @@ inline IndexSpace make_space(const IndexSpace::Type& type) { Spin::any); } -/// enumerates the known Operator types -enum class OpType { - h, //!< 1-body Hamiltonian - f, //!< Fock operator - f̃, //!< closed Fock operator (i.e. Fock operator due to fully-occupied - //!< orbitals) - s, //!< 1-body overlap - θ, //!< general fock space operator - g, //!< 2-body Coulomb - t, //!< cluster amplitudes - λ, //!< deexcitation cluster amplitudes - A, //!< antisymmetrizer - S, //!< particle symmetrizer - L, //!< left-hand eigenstate - R, //!< right-hand eigenstate - R12, //!< geminal kernel - GR, //!< GR kernel from f12 theory - C, //!< cabs singles op - RDM, //!< RDM - RDMCumulant, //!< RDM cumulant - δ, //!< Kronecker delta (=identity) operator - h_1, //!< Hamiltonian perturbation - t_1, //!< first order perturbed excitation cluster operator - λ_1, //!< first order perturbed deexcitation cluster operator -}; - -/// maps operator types to their labels -inline const container::map optype2label{ - {OpType::h, L"h"}, - {OpType::f, L"f"}, - {OpType::f̃, L"f̃"}, - {OpType::s, overlap_label()}, - {OpType::δ, kronecker_label()}, - {OpType::g, L"g"}, - {OpType::θ, L"θ"}, - {OpType::t, L"t"}, - {OpType::λ, L"λ"}, - {OpType::A, L"A"}, - {OpType::S, L"S"}, - {OpType::L, L"L"}, - {OpType::R, L"R"}, - {OpType::R12, L"F"}, - {OpType::GR, L"GR"}, - {OpType::C, L"C"}, - {OpType::RDM, L"γ"}, - // see https://en.wikipedia.org/wiki/Cumulant - {OpType::RDMCumulant, L"κ"}, - {OpType::h_1, L"h¹"}, - {OpType::t_1, L"t¹"}, - {OpType::λ_1, L"λ¹"}}; - -/// maps operator labels to their types -inline const container::map label2optype = - ranges::views::zip(ranges::views::values(optype2label), - ranges::views::keys(optype2label)) | - ranges::to>(); - -/// Operator character relative to Fermi vacuum -enum class OpClass { ex, deex, gen }; - /// @brief Struct which holds parameters for operator construction struct OpParams { - std::size_t order = 1; ///< perturbation order (for _pt operators) + std::size_t order = 0; ///< perturbation order (for _pt operators) std::optional nbatch = std::nullopt; ///< number of batching indices container::svector batch_ordinals{}; ///< custom batching index ordinals (empty = no batching) @@ -157,8 +102,21 @@ struct OpParams { /// @return the tensor labels in the cardinal order std::vector cardinal_tensor_labels(); -std::wstring to_wstring(OpType op); -OpClass to_class(OpType op); +/// @brief decorates a base label with perturbation order as superscript +/// @param base_label the base label to decorate +/// @param pert_order the perturbation order to decorate with +/// @return the decorated label +inline std::wstring decorate_with_pert_order(std::wstring_view base_label, + int pert_order = 0) { + if (pert_order == 0) return std::wstring(base_label); + SEQUANT_ASSERT( + pert_order >= 0 && pert_order <= 9, + "decorate_with_pert_order: perturbation order out of range [0,9]"); + + std::wstring result(base_label); + result += detail::pert_superscripts[pert_order]; + return result; +} ////////////////////////////// @@ -616,9 +574,10 @@ class OpMaker { /// @param[in] cre_list list of creator indices /// @param[in] ann_list list of annihilator indices template - OpMaker(OpType op, const cre& cre_list, - const ann& ann_list) - : op_(op), + OpMaker(const std::wstring& label, const cre& cre_list, + const ann& ann_list, size_t order = 0) + : label_(label), + order_(order), cre_spaces_(cre_list.begin(), cre_list.end()), ann_spaces_(ann_list.begin(), ann_list.end()) { SEQUANT_ASSERT(ncreators() > 0 || nannihilators() > 0); @@ -627,28 +586,28 @@ class OpMaker { /// @param[in] op the operator type /// @param[in] nc number of bra indices/creators /// @param[in] na number of ket indices/annihilators - OpMaker(OpType op, ncre nc, nann na); + OpMaker(const std::wstring& label, ncre nc, nann na); /// @brief creates a particle-conserving replacement operator /// @param[in] op the operator type /// @param[in] rank particle rank of the operator (# of creators = # of /// annihilators = @p rank ) - OpMaker(OpType op, std::size_t rank); + OpMaker(const std::wstring& label, std::size_t rank); /// @param[in] op the operator type /// @param[in] nc number of bra indices/creators /// @param[in] na number of ket indices/annihilators /// @param[in] cre_space IndexSpace referred to be the creator /// @param[in] ann_space IndexSpace referred to be the annihilators - OpMaker(OpType op, ncre nc, nann na, const cre& cre_space, - const ann& ann_space); + OpMaker(const std::wstring& label, ncre nc, nann na, + const cre& cre_space, const ann& ann_space); /// @brief Creates operator from OpParams /// @param[in] op the operator type /// @param[in] nc number of bra indices/creators /// @param[in] na number of ket indices/annihilators /// @param[in] params named parameters for operator construction - OpMaker(OpType op, ncre nc, nann na, const OpParams& params); + OpMaker(const std::wstring& label, ncre nc, nann na, const OpParams& params); enum class UseDepIdx { /// bra/cre indices depend on ket @@ -838,12 +797,13 @@ class OpMaker { } protected: - OpType op_; + std::wstring label_; + size_t order_ = 0; IndexSpaceContainer cre_spaces_; IndexSpaceContainer ann_spaces_; std::optional batch_indices_ = std::nullopt; - OpMaker(OpType op); + OpMaker(const std::wstring& op); [[nodiscard]] auto ncreators() const { return cre_spaces_.size(); }; [[nodiscard]] auto nannihilators() const { return ann_spaces_.size(); }; @@ -874,31 +834,10 @@ class Operator : public Operator { std::function tensor_form_generator, std::function qn_action); - /// @brief Constructs an operator with the given label and tensor form and - /// quantum number action - /// @param label_generator a function that generates the label for the - /// operator - /// @param tensor_form_generator a function that generates the tensor form of - /// the operator - /// @param qn_action a function that modifies the quantum numbers - /// @param batch_idx_rank the rank of the batch index, must be non-zero Operator(std::function label_generator, std::function tensor_form_generator, std::function qn_action, - size_t batch_idx_rank); - - /// @brief Constructs an operator with the given label and tensor form and - /// quantum number action - /// @param label_generator a function that generates the label for the - /// operator - /// @param tensor_form_generator a function that generates the tensor form of - /// the operator - /// @param qn_action a function that modifies the quantum numbers - /// @param batch_ordinals the unique, sorted ordinals of the batch indices - Operator(std::function label_generator, - std::function tensor_form_generator, - std::function qn_action, - const container::svector& batch_ordinals); + const OpParams& params); virtual ~Operator(); @@ -939,6 +878,9 @@ class Operator : public Operator { return batch_ordinals_; } + /// @brief returns the pertubation order of this operator + [[nodiscard]] size_t order() const { return order_; } + private: std::function qn_action_; @@ -946,6 +888,8 @@ class Operator : public Operator { std::optional> batch_ordinals_ = std::nullopt; + size_t order_ = 0; + bool less_than_rank_of(const this_type& that) const; Expr::type_id_type type_id() const override; @@ -984,7 +928,8 @@ ExprPtr H(std::size_t k = 2); /// @brief Fock operator implied one-body operator, optional explicit /// construction requires user to specify the IndexSpace corresponding to all /// orbitals which may have non-zero density. -ExprPtr F(bool use_tensor = true, IndexSpace reference_occupied = {L"", 0}); +ExprPtr F(bool use_tensor = true, + const IndexSpace& reference_occupied = {L"", 0}); /// A general operator of rank \p K ExprPtr θ(std::size_t K); @@ -1118,7 +1063,8 @@ ExprPtr H(std::size_t k = 2); /// @brief Fock operator implied one-body operator, optional explicit /// construction requires user to specify the IndexSpace corresponding to all /// orbitals which may have non-zero density. -ExprPtr F(bool use_tensor = true, IndexSpace reference_occupied = {L"", 0}); +ExprPtr F(bool use_tensor = true, + const IndexSpace& reference_occupied = {L"", 0}); /// A general operator of rank \p K ExprPtr θ(std::size_t K); diff --git a/SeQuant/domain/mbpt/op.ipp b/SeQuant/domain/mbpt/op.ipp index 3a10c75564..aaf160d55d 100644 --- a/SeQuant/domain/mbpt/op.ipp +++ b/SeQuant/domain/mbpt/op.ipp @@ -26,31 +26,24 @@ template Operator::Operator( std::function label_generator, std::function tensor_form_generator, - std::function qn_action, size_t batch_idx_rank) + std::function qn_action, const OpParams& params) : Operator(std::move(label_generator), std::move(tensor_form_generator), std::move(qn_action)) { - mbpt::check_for_batching_space(); - SEQUANT_ASSERT(batch_idx_rank != 0 && - "Operator: batch_idx_rank cannot be zero"); - // make aux ordinals [1 to batch_idx_rank] - batch_ordinals_ = ranges::views::iota(1ul, 1ul + batch_idx_rank) | - ranges::to>(); -} - -template -Operator::Operator( - std::function label_generator, - std::function tensor_form_generator, - std::function qn_action, - const container::svector& batch_ordinals) - : Operator(std::move(label_generator), std::move(tensor_form_generator), - std::move(qn_action)) { - mbpt::check_for_batching_space(); - SEQUANT_ASSERT(!batch_ordinals.empty() && - "Operator: batch_ordinals cannot be empty"); - SEQUANT_ASSERT(ranges::is_sorted(batch_ordinals) && - "Operator: batch_ordinals must be sorted"); - batch_ordinals_ = batch_ordinals; + params.validate(); + order_ = params.order; + // set up batch ordinals if any + if (!params.batch_ordinals.empty()) { + mbpt::check_for_batching_space(); + SEQUANT_ASSERT(ranges::is_sorted(params.batch_ordinals) && + "Operator: batch_ordinals must be sorted"); + batch_ordinals_ = params.batch_ordinals; + } else if (params.nbatch) { + mbpt::check_for_batching_space(); + SEQUANT_ASSERT(params.nbatch != 0 && "Operator: nbatch cannot be zero"); + // make aux ordinals [1 to nbatch] + batch_ordinals_ = ranges::views::iota(1ul, 1ul + params.nbatch.value()) | + ranges::to>(); + } } template @@ -207,6 +200,7 @@ Expr::hash_type Operator::memoizing_hash() const { auto qns = (*this)(QuantumNumbers{}); auto val = sequant::hash::value(qns); sequant::hash::combine(val, std::wstring(this->label())); + sequant::hash::combine(val, this->order()); if (batch_ordinals()) { const auto ordinals = batch_ordinals().value(); sequant::hash::range(val, begin(ordinals), end(ordinals)); diff --git a/SeQuant/domain/mbpt/op_registry.cpp b/SeQuant/domain/mbpt/op_registry.cpp new file mode 100644 index 0000000000..5133ced6e7 --- /dev/null +++ b/SeQuant/domain/mbpt/op_registry.cpp @@ -0,0 +1,72 @@ +// +// Created by Ajay Melekamburath on 12/14/25. +// + +#include +#include +#include +#include + +namespace sequant::mbpt { +void OpRegistry::validate_op(const std::wstring& op) const { + if (!reserved::is_nonreserved(op)) { + throw std::runtime_error("mbpt::OpRegistry::add: operator " + + sequant::to_string(op) + " uses a reserved label"); + } + if (this->contains(op)) { + throw std::runtime_error("mbpt::OpRegistry::add: operator " + + sequant::to_string(op) + + " already exists in registry"); + } +} + +OpRegistry& OpRegistry::operator=(const OpRegistry& other) { + ops_ = other.ops_; + return *this; +} + +OpRegistry& OpRegistry::operator=(OpRegistry&& other) noexcept { + ops_ = std::move(other.ops_); + return *this; +} + +OpRegistry& OpRegistry::add(const std::wstring& op, OpClass action) { + this->validate_op(op); + ops_->emplace(op, action); + return *this; +} + +OpRegistry& OpRegistry::remove(const std::wstring& op) { + if (!this->contains(op)) { + throw std::runtime_error("mbpt::OpRegistry::remove: operator " + + sequant::to_string(op) + + " does not exist in registry"); + } + ops_->erase(op); + return *this; +} + +bool OpRegistry::contains(const std::wstring& op) const { + return ops_->contains(op); +} + +OpClass OpRegistry::to_class(const std::wstring& op) const { + if (!this->contains(op)) { + throw std::runtime_error("mbpt::OpRegistry::to_class: operator " + + sequant::to_string(op) + + " does not exist in registry"); + } + return ops_->at(op); +} + +OpRegistry OpRegistry::clone() const { + OpRegistry result(*this); + result.ops_ = std::make_shared>(*ops_); + return result; +} + +container::svector OpRegistry::ops() const { + return ranges::views::keys(*ops_) | + ranges::to>(); +} +} // namespace sequant::mbpt diff --git a/SeQuant/domain/mbpt/op_registry.hpp b/SeQuant/domain/mbpt/op_registry.hpp new file mode 100644 index 0000000000..0cb897ba3d --- /dev/null +++ b/SeQuant/domain/mbpt/op_registry.hpp @@ -0,0 +1,92 @@ +// +// Created by Ajay Melekamburath on 12/14/25. +// + +#ifndef SEQUANT_DOMAIN_MBPT_OP_REGISTRY_HPP +#define SEQUANT_DOMAIN_MBPT_OP_REGISTRY_HPP + +#include +#include +#include + +#include + +#include + +namespace sequant::mbpt { + +/// Operator character relative to Fermi vacuum +enum class OpClass { ex, deex, gen }; + +/// @brief A Registry for MBPT Operators +/// +/// A registry that keeps track of MBPT operators by their labels and +/// properties. +class OpRegistry { + public: + /// default constructor, creates an empty registry + OpRegistry() + : ops_(std::make_shared>()) {} + + /// constructs an OpRegistry from an existing map of operators and their + /// classes + OpRegistry(std::shared_ptr> ops) + : ops_(std::move(ops)) {} + + /// copy constructor + OpRegistry(const OpRegistry& other) : ops_(other.ops_) {} + + /// move constructor + OpRegistry(OpRegistry&& other) noexcept : ops_(std::move(other.ops_)) {} + + /// copy assignment operator + OpRegistry& operator=(const OpRegistry& other); + + /// move assignment operator + OpRegistry& operator=(OpRegistry&& other) noexcept; + + /// @brief clones this OpRegistry, creates a copy of ops_ + OpRegistry clone() const; + + /// @brief Adds a new operator to the registry + /// @param op the operator label + /// @param action the class of the operator + OpRegistry& add(const std::wstring& op, OpClass action); + + /// @brief Removes an operator from the registry + OpRegistry& remove(const std::wstring& op); + + /// @brief Checks if the registry contains an operator with the given label + /// @param op the operator label + /// @return true if the operator exists, false otherwise + bool contains(const std::wstring& op) const; + + /// @brief Returns the class of the operator corresponding to the given label + /// if it exists + /// @param op the operator label + /// @return the class of the operator + [[nodiscard]] OpClass to_class(const std::wstring& op) const; + + /// @brief returns a list of registered operators + [[nodiscard]] container::svector ops() const; + + /// @brief clears all registered operators + void purge() { ops_->clear(); } + + private: + std::shared_ptr> ops_; + + /// @brief Validates that the operator label is not reserved and not already + /// registered + /// @param op the operator label to validate + /// @throws std::runtime_error if the label is reserved or already exists + void validate_op(const std::wstring& op) const; + + /// @brief Equality operator for OpRegistry + friend bool operator==(const OpRegistry& reg1, const OpRegistry& reg2) { + return *reg1.ops_ == *reg2.ops_; + } +}; +} // namespace sequant::mbpt + +#endif // SEQUANT_DOMAIN_MBPT_OP_REGISTRY_HPP diff --git a/SeQuant/domain/mbpt/rdm.cpp b/SeQuant/domain/mbpt/rdm.cpp index ee90ddc9c7..f3b3423924 100644 --- a/SeQuant/domain/mbpt/rdm.cpp +++ b/SeQuant/domain/mbpt/rdm.cpp @@ -1,36 +1,45 @@ #include #include +// anonymous namespace for holding specific labels used here +namespace { +std::wstring rdm_label() { + static const std::wstring label = L"γ"; + return label; +} + +std::wstring rdm_cumulant_label() { + static const std::wstring label = L"κ"; + return label; +} +} // namespace + namespace sequant::mbpt::decompositions { ExprPtr cumu_to_density(ExprPtr ex_) { SEQUANT_ASSERT(ex_->is()); SEQUANT_ASSERT(ex_->as().rank() == 1); - SEQUANT_ASSERT(ex_->as().label() == - optype2label.at(OpType::RDMCumulant)); + SEQUANT_ASSERT(ex_->as().label() == rdm_label()); auto down_0 = ex_->as().ket()[0]; auto up_0 = ex_->as().bra()[0]; - auto density = - ex(optype2label.at(OpType::RDM), bra{up_0}, ket{down_0}); + auto density = ex(rdm_label(), bra{up_0}, ket{down_0}); return density; } sequant::ExprPtr cumu2_to_density(sequant::ExprPtr ex_) { SEQUANT_ASSERT(ex_->is()); SEQUANT_ASSERT(ex_->as().rank() == 2); - SEQUANT_ASSERT(ex_->as().label() == - optype2label.at(OpType::RDMCumulant)); + SEQUANT_ASSERT(ex_->as().label() == rdm_cumulant_label()); auto down_0 = ex_->as().ket()[0]; auto up_0 = ex_->as().bra()[0]; auto down_1 = ex_->as().ket()[1]; auto up_1 = ex_->as().bra()[1]; - const auto rdm_label = optype2label.at(OpType::RDM); - auto density2 = ex(rdm_label, bra{up_0, up_1}, ket{down_0, down_1}); - auto density_1 = ex(rdm_label, bra{up_0}, ket{down_0}); - auto density_2 = ex(rdm_label, bra{up_1}, ket{down_1}); + auto density2 = ex(rdm_label(), bra{up_0, up_1}, ket{down_0, down_1}); + auto density_1 = ex(rdm_label(), bra{up_0}, ket{down_0}); + auto density_2 = ex(rdm_label(), bra{up_1}, ket{down_1}); auto d1_d2 = antisymmetrize(density_1 * density_2); return density2 + ex(-1) * d1_d2.result; @@ -39,8 +48,7 @@ sequant::ExprPtr cumu2_to_density(sequant::ExprPtr ex_) { ExprPtr cumu3_to_density(ExprPtr ex_) { SEQUANT_ASSERT(ex_->is()); SEQUANT_ASSERT(ex_->as().rank() == 3); - SEQUANT_ASSERT(ex_->as().label() == - optype2label.at(OpType::RDMCumulant)); + SEQUANT_ASSERT(ex_->as().label() == rdm_cumulant_label()); auto down_0 = ex_->as().ket()[0]; auto up_0 = ex_->as().bra()[0]; @@ -49,14 +57,13 @@ ExprPtr cumu3_to_density(ExprPtr ex_) { auto down_2 = ex_->as().ket()[2]; auto up_2 = ex_->as().bra()[2]; - const auto rdm_label = optype2label.at(OpType::RDM); - auto cumulant2 = ex(optype2label.at(OpType::RDMCumulant), - bra{up_1, up_2}, ket{down_1, down_2}); - auto density_1 = ex(rdm_label, bra{up_0}, ket{down_0}); - auto density_2 = ex(rdm_label, bra{up_1}, ket{down_1}); - auto density_3 = ex(rdm_label, bra{up_2}, ket{down_2}); - auto density3 = - ex(rdm_label, bra{up_0, up_1, up_2}, ket{down_0, down_1, down_2}); + auto cumulant2 = + ex(rdm_cumulant_label(), bra{up_1, up_2}, ket{down_1, down_2}); + auto density_1 = ex(rdm_label(), bra{up_0}, ket{down_0}); + auto density_2 = ex(rdm_label(), bra{up_1}, ket{down_1}); + auto density_3 = ex(rdm_label(), bra{up_2}, ket{down_2}); + auto density3 = ex(rdm_label(), bra{up_0, up_1, up_2}, + ket{down_0, down_1, down_2}); auto d1_d2 = antisymmetrize(density_1 * density_2 * density_3 + density_1 * cumulant2); @@ -65,8 +72,7 @@ ExprPtr cumu3_to_density(ExprPtr ex_) { for (auto&& product : temp_result->as().summands()) { for (auto&& factor : product->as().factors()) { if (factor->is() && - (factor->as().label() == - optype2label.at(OpType::RDMCumulant)) && + (factor->as().label() == rdm_cumulant_label()) && (factor->as().rank() == 2)) { factor = cumu2_to_density(factor); } @@ -75,8 +81,7 @@ ExprPtr cumu3_to_density(ExprPtr ex_) { for (auto&& product : temp_result->as().summands()) { for (auto&& factor : product->as().factors()) { if (factor->is() && - factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + factor->as().label() == rdm_cumulant_label() && factor->as().rank() == 1) { factor = cumu_to_density(factor); } @@ -93,8 +98,7 @@ ExprPtr one_body_sub(ExprPtr ex_) { // J. Chem. Phys. 132, 234107 (2010); auto up_0 = ex_->as().creators()[0].index(); const auto a = ex(cre({up_0}), ann({down_0})); - const auto cumu1 = - ex(optype2label.at(OpType::RDMCumulant), bra{down_0}, ket{up_0}); + const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); auto result = a + (ex(-1) * cumu1); return (result); @@ -113,14 +117,12 @@ ExprPtr two_body_decomp( auto up_0 = ex_->as().creators()[0].index(); auto up_1 = ex_->as().creators()[1].index(); - const auto cumu1 = - ex(optype2label.at(OpType::RDMCumulant), bra{down_0}, ket{up_0}); - const auto cumu2 = - ex(optype2label.at(OpType::RDMCumulant), bra{down_1}, ket{up_1}); + const auto cumu1 = ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); + const auto cumu2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}); const auto a = ex(cre{up_1}, ann{down_1}); const auto a2 = ex(cre{up_0, up_1}, ann{down_0, down_1}); - const auto double_cumu = ex(optype2label.at(OpType::RDMCumulant), - bra{down_0, down_1}, ket{up_0, up_1}); + const auto double_cumu = + ex(rdm_cumulant_label(), bra{down_0, down_1}, ket{up_0, up_1}); auto term1 = cumu1 * a; auto term2 = cumu1 * cumu2; @@ -150,21 +152,19 @@ three_body_decomp(ExprPtr ex_, bool approx) { std::vector initial_upper{up_0, up_1, up_2}; const auto cumulant = - ex(optype2label.at(OpType::RDMCumulant), bra{down_0}, ket{up_0}); + ex(rdm_cumulant_label(), bra{down_0}, ket{up_0}); const auto a = ex(cre{up_1, up_2}, ann{down_1, down_2}); auto a_cumulant = cumulant * a; - auto cumulant2 = - ex(optype2label.at(OpType::RDMCumulant), bra{down_1}, ket{up_1}); - auto cumulant3 = - ex(optype2label.at(OpType::RDMCumulant), bra{down_2}, ket{up_2}); + auto cumulant2 = ex(rdm_cumulant_label(), bra{down_1}, ket{up_1}); + auto cumulant3 = ex(rdm_cumulant_label(), bra{down_2}, ket{up_2}); auto cumulant_3x = cumulant * cumulant2 * cumulant3; auto a1 = ex(cre{up_0}, ann{down_0}); auto a1_cumu1_cumu2 = a1 * cumulant2 * cumulant3; - auto two_body_cumu = ex(optype2label.at(OpType::RDMCumulant), - bra{down_1, down_2}, ket{up_1, up_2}); + auto two_body_cumu = + ex(rdm_cumulant_label(), bra{down_1, down_2}, ket{up_1, up_2}); auto a1_cumu2 = a1 * two_body_cumu; auto cumu1_cumu2 = cumulant * two_body_cumu; @@ -172,8 +172,8 @@ three_body_decomp(ExprPtr ex_, bool approx) { a1_cumu2 + cumu1_cumu2); if (!approx) { - auto cumu3 = ex(optype2label.at(OpType::RDMCumulant), - bra{down_0, down_1, down_2}, ket{up_0, up_1, up_2}); + auto cumu3 = ex(rdm_cumulant_label(), bra{down_0, down_1, down_2}, + ket{up_0, up_1, up_2}); sum_of_terms.result = cumu3 + sum_of_terms.result; } @@ -229,21 +229,18 @@ three_body_decomposition(ExprPtr ex_, int rank, bool fast) { if (product->is()) { for (auto&& factor : product->as().factors()) { if (factor->is()) { - if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() == 3) { factor = cumu3_to_density(factor); - } else if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + } else if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() == 2) { factor = cumu2_to_density(factor); - } else if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + } else if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() == 1) { factor = cumu_to_density(factor); } else { SEQUANT_ASSERT(factor->as().label() != - optype2label.at(OpType::RDMCumulant)); + rdm_cumulant_label()); } } } @@ -287,20 +284,17 @@ three_body_decomposition(ExprPtr ex_, int rank, bool fast) { if (product->is()) { for (auto&& factor : product->as().factors()) { if (factor->is()) { - if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() > 2) { factor = ex(0); - } else if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + } else if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() == 2) { factor = cumu2_to_density(factor); - } else if (factor->as().label() == - optype2label.at(OpType::RDMCumulant)) { + } else if (factor->as().label() == rdm_cumulant_label()) { factor = cumu_to_density(factor); } else { SEQUANT_ASSERT(factor->as().label() != - optype2label.at(OpType::RDMCumulant)); + rdm_cumulant_label()); } } } @@ -318,16 +312,14 @@ three_body_decomposition(ExprPtr ex_, int rank, bool fast) { if (product->is()) { for (auto&& factor : product->as().factors()) { if (factor->is()) { - if (factor->as().label() == - optype2label.at(OpType::RDMCumulant) && + if (factor->as().label() == rdm_cumulant_label() && factor->as().rank() > 1) { factor = ex(0); - } else if (factor->as().label() == - optype2label.at(OpType::RDMCumulant)) { + } else if (factor->as().label() == rdm_cumulant_label()) { factor = cumu_to_density(factor); } else { SEQUANT_ASSERT(factor->as().label() != - optype2label.at(OpType::RDMCumulant)); + rdm_cumulant_label()); } } } diff --git a/SeQuant/domain/mbpt/reserved.hpp b/SeQuant/domain/mbpt/reserved.hpp new file mode 100644 index 0000000000..e9f8c6990e --- /dev/null +++ b/SeQuant/domain/mbpt/reserved.hpp @@ -0,0 +1,40 @@ +// +// Created by Ajay Melekamburath on 12/24/25. +// + +#ifndef SEQUANT_DOMAIN_MBPT_RESERVED_HPP +#define SEQUANT_DOMAIN_MBPT_RESERVED_HPP + +namespace sequant::mbpt { + +namespace reserved { +/// @brief returns the reserved label for the antisymmetrization operator +inline const std::wstring& antisymm_label() { + static const std::wstring label = L"A"; + return label; +} + +/// @brief returns the reserved label for the symmetrization operator +inline const std::wstring& symm_label() { + static const std::wstring label = L"S"; + return label; +} + +/// @brief returns a list of all reserved operator labels +inline auto labels() { + static const std::array reserved{antisymm_label(), symm_label(), + sequant::kronecker_label(), + sequant::overlap_label()}; + return reserved; +} + +/// @brief checks if a label is not reserved +inline bool is_nonreserved(const std::wstring& label) { + return !ranges::contains(labels(), label); +} + +} // namespace reserved + +} // namespace sequant::mbpt + +#endif // SEQUANT_DOMAIN_MBPT_RESERVED_HPP diff --git a/SeQuant/domain/mbpt/vac_av.cpp b/SeQuant/domain/mbpt/vac_av.cpp index f9aba37bf1..7a533967be 100644 --- a/SeQuant/domain/mbpt/vac_av.cpp +++ b/SeQuant/domain/mbpt/vac_av.cpp @@ -154,12 +154,6 @@ ExprPtr ref_av(ExprPtr expr, const OpConnections& op_connections, skip_clone, full_contractions); } -ExprPtr ref_av(ExprPtr expr, const OpConnections& op_connections, - bool use_topology, bool screen, bool skip_clone) { - return ref_av(expr, to_label_connections(op_connections), use_topology, - screen, skip_clone); -} - ExprPtr vac_av(ExprPtr expr, const OpConnections& op_connections, bool use_topology, bool screen, bool skip_clone) { return expectation_value_impl(expr, op_connections, use_topology, screen, @@ -167,12 +161,6 @@ ExprPtr vac_av(ExprPtr expr, const OpConnections& op_connections, /* full_contractions */ true); } -ExprPtr vac_av(ExprPtr expr, const OpConnections& op_connections, - bool use_topology, bool screen, bool skip_clone) { - return vac_av(expr, to_label_connections(op_connections), use_topology, - screen, skip_clone); -} - } // namespace op } // namespace mbpt } // namespace sequant diff --git a/SeQuant/domain/mbpt/vac_av.hpp b/SeQuant/domain/mbpt/vac_av.hpp index 4830f219d6..9e837cb723 100644 --- a/SeQuant/domain/mbpt/vac_av.hpp +++ b/SeQuant/domain/mbpt/vac_av.hpp @@ -14,13 +14,12 @@ template using OpConnections = std::vector>; /// defines the default op connections -inline OpConnections default_op_connections() { - using mbpt::OpType; - static const OpConnections defaults = { - {OpType::h, OpType::t}, - {OpType::f, OpType::t}, - {OpType::f̃, OpType::t}, - {OpType::g, OpType::t}, +inline OpConnections default_op_connections() { + static const OpConnections defaults = { + {L"h", L"t"}, + {L"f", L"t"}, + {L"f̃", L"t"}, + {L"g", L"t"}, // NBs // - for unitary ansatz can also connect t^+ with Hamiltonian // - for exact (non-approximated) unitary ansatz will also need to connect @@ -30,36 +29,20 @@ inline OpConnections default_op_connections() { // as same connection ... this points out the need to have a separate // OpType for t^+ and t, and in general the contents of OpType must be // customizable - {OpType::t, OpType::h}, - {OpType::t, OpType::f}, - {OpType::t, OpType::f̃}, - {OpType::t, OpType::g}}; + {L"t", L"h"}, + {L"t", L"f"}, + {L"t", L"f̃"}, + {L"t", L"g"}}; return defaults; } /// concat 2 sets of op connections -template || - std::is_same_v>> -inline OpConnections concat(const OpConnections& connections1, - const OpConnections& connections2) { +inline OpConnections concat( + const OpConnections& connections1, + const OpConnections& connections2) { return ranges::concat_view(connections1, connections2) | ranges::to_vector; } -/// lowers representation of op connections from mbpt::OpType to labels -/// @param[in] op_connections vector of pairs of operators to be connected -/// @return vector of pairs of operator labels to be connected -inline OpConnections to_label_connections( - const OpConnections& op_connections) { - // convert mbpt::OpType to std::wstring - using mbpt::optype2label; - OpConnections op_connect_wstr; - for (const auto& [op1, op2] : op_connections) { - op_connect_wstr.emplace_back(optype2label.at(op1), optype2label.at(op2)); - } - return op_connect_wstr; -} - /// @brief lowers an expression composed of Operators to tensor form /// @param[in] expr input expression /// @return expression with all Operators lowered to tensor form @@ -101,33 +84,7 @@ ExprPtr expectation_value_impl( /// @param[in] skip_clone if true, will not clone the input expression /// @return the reference expectation value // clang_format on -ExprPtr ref_av(ExprPtr expr, const OpConnections& op_connections, - bool use_topology = true, bool screen = true, - bool skip_clone = false); - -// clang-format off -/// @brief computes the reference expectation value -/// @note equivalent to vac_av if the reference state is the Wick vacuum, -/// i.e. if `get_default_context().index_space_registry()->reference_occupied_space() == get_default_context().index_space_registry()->vacuum_occupied_space()` -/// @param[in] expr input expression -/// @param[in] op_connections list of pairs of operators to be -/// connected; connections are left-to-right, i.e., pair `{opL,opR}` declares -/// that `opL` and `opR` are to be connected when `opR` precedes `opL`, i.e. -/// `opL` is to the left of `opR`; -/// e.g., `{{OpType::h, OpType::t}}` will ensure that each -/// `OpType::h` will be connected to at least one `OpType::t` on its right -/// (preceding it). The default list of connections is given -/// by default_op_connections() -/// @param[in] use_topology if true, WickTheorem uses topological equivalence of -/// terms, default is true -/// @param[in] screen if true, expressions are screened before lowering to -/// Tensor level and calling WickTheorem, default is true -/// @param[in] skip_clone if true, will not clone the input expression -/// @return the reference expectation value -// clang-format on -ExprPtr ref_av(ExprPtr expr, - const OpConnections& op_connections = - default_op_connections(), +ExprPtr ref_av(ExprPtr expr, const OpConnections& op_connections = default_op_connections(), bool use_topology = true, bool screen = true, bool skip_clone = false); @@ -144,30 +101,7 @@ ExprPtr ref_av(ExprPtr expr, /// Tensor level and calling WickTheorem, default is true /// @param[in] skip_clone if true, will not clone the input expression /// @return the VEV -ExprPtr vac_av(ExprPtr expr, const OpConnections& op_connections, - bool use_topology = true, bool screen = true, - bool skip_clone = false); - -/// @brief computes the vacuum expectation value -/// @internal evaluates only full contractions in WickTheorem -/// @param[in] expr input expression -/// @param[in] op_connections list of pairs of operators to be -/// connected; connections are left-to-right, i.e., pair `{opL,opR}` declares -/// that `opL` and `opR` are to be connected when `opR` precedes `opL`, i.e. -/// `opL` is to the left of `opR`; -/// e.g., `{{OpType::h, OpType::t}}` will ensure that each -/// `OpType::h` will be connected to at least one `OpType::t` on its right -/// (preceding it). The default list of connections is given -/// by default_op_connections() -/// @param[in] use_topology if true, WickTheorem uses topological equivalence of -/// terms, default is true -/// @param[in] screen if true, expressions are screened before lowering to -/// Tensor level and calling WickTheorem, default is true -/// @param[in] skip_clone if true, will not clone the input expression -/// @return the VEV -ExprPtr vac_av(ExprPtr expr, - const OpConnections& op_connections = - default_op_connections(), +ExprPtr vac_av(ExprPtr expr, const OpConnections& op_connections = default_op_connections(), bool use_topology = true, bool screen = true, bool skip_clone = false); diff --git a/benchmarks/wick.cpp b/benchmarks/wick.cpp index 5141d2695a..31d2d9d0a5 100644 --- a/benchmarks/wick.cpp +++ b/benchmarks/wick.cpp @@ -150,7 +150,7 @@ VacAvPair get_mbpt_expr(std::size_t i) { static void mbpt_vac_av(benchmark::State &state, bool csv) { auto ctx = sequant::mbpt::set_scoped_default_mbpt_context( - mbpt::Context(csv ? mbpt::CSV::Yes : mbpt::CSV::No)); + mbpt::Context({.csv = csv ? CSV::Yes : CSV::No})); VacAvPair input = get_mbpt_expr(state.range(0)); diff --git a/doc/examples/synopsis/synopsis6.cpp b/doc/examples/synopsis/synopsis6.cpp index 82a1b51101..13515a5c62 100644 --- a/doc/examples/synopsis/synopsis6.cpp +++ b/doc/examples/synopsis/synopsis6.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -12,6 +13,7 @@ int main() { using namespace sequant::mbpt; set_default_context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct}); + set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()}); auto hbar = H(2) + commutator(H(2), T_(2)) + diff --git a/doc/examples/user/cc.cpp b/doc/examples/user/cc.cpp index 91fcac7ce2..409e8902a8 100644 --- a/doc/examples/user/cc.cpp +++ b/doc/examples/user/cc.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -16,6 +17,7 @@ int main() { using namespace sequant::mbpt; set_default_context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct}); + set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()}); // end-snippet-0 // start-snippet-1 diff --git a/doc/examples/user/getting_started/ccd.cpp b/doc/examples/user/getting_started/ccd.cpp index ea65012454..8e554414d1 100644 --- a/doc/examples/user/getting_started/ccd.cpp +++ b/doc/examples/user/getting_started/ccd.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -12,6 +13,7 @@ int main() { using namespace sequant::mbpt; set_default_context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct}); + set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()}); auto hbar = H(2) + commutator(H(2), T_(2)) + diff --git a/doc/examples/user/getting_started/ccd_pedestrian.cpp b/doc/examples/user/getting_started/ccd_pedestrian.cpp index cd61154f49..db52826c0c 100644 --- a/doc/examples/user/getting_started/ccd_pedestrian.cpp +++ b/doc/examples/user/getting_started/ccd_pedestrian.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -9,6 +10,7 @@ int main() { using namespace sequant::mbpt; set_default_context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct}); + set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()}); // start-snippet-1 auto t2 = ex(rational(1, 4)) * diff --git a/doc/examples/user/operator.cpp b/doc/examples/user/operator.cpp index ba7f952fdb..77a2675bdf 100644 --- a/doc/examples/user/operator.cpp +++ b/doc/examples/user/operator.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -14,6 +15,7 @@ int main() { using namespace sequant::mbpt; set_default_context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct}); + set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()}); // end-snippet-0 // start-snippet-1 diff --git a/python/src/sequant/mbpt.h b/python/src/sequant/mbpt.h index e326c78e9d..92b6f09fe6 100644 --- a/python/src/sequant/mbpt.h +++ b/python/src/sequant/mbpt.h @@ -1,12 +1,14 @@ #ifndef SEQUANT_PYTHON_MBPT_H #define SEQUANT_PYTHON_MBPT_H +#include #include #include #include #include #include +#include #include "python.h" @@ -20,9 +22,20 @@ auto make_sr_op(F f) { return op; } -template -ExprPtr VacuumAverage(const ExprPtr& e, const Args&... args) { - return sequant::mbpt::op::vac_av(e, args...); +// Overload for no op_connections +ExprPtr VacuumAverage(const ExprPtr& e) { return sequant::mbpt::op::vac_av(e); } + +// overload with string conversion +ExprPtr VacuumAverage( + const ExprPtr& e, + const std::vector>& op_connections) { + sequant::mbpt::OpConnections wop_connections; + wop_connections.reserve(op_connections.size()); + for (const auto& [op1, op2] : op_connections) { + wop_connections.emplace_back(sequant::to_wstring(op1), + sequant::to_wstring(op2)); + } + return sequant::mbpt::op::vac_av(e, wop_connections); } #define SR_OP(OP) \ @@ -33,11 +46,10 @@ inline void __init__(py::module m) { sequant::TensorCanonicalizer::register_instance( std::make_shared()); - py::enum_(m, "OpType") - .value("h", sequant::mbpt::OpType::h) - .value("f", sequant::mbpt::OpType::f) - .value("t", sequant::mbpt::OpType::t) - .export_values(); + // mbpt context setup + sequant::mbpt::Context::Options opts; + opts.op_registry_ptr = sequant::mbpt::make_legacy_registry(); + sequant::mbpt::set_default_mbpt_context(opts); m.def("F", &sequant::mbpt::F); m.def("H", &sequant::mbpt::H, @@ -48,10 +60,14 @@ inline void __init__(py::module m) { m.def(SR_OP(T)); m.def(SR_OP(T_)); - m.def("VacuumAverage", &VacuumAverage<>); + m.def("VacuumAverage", py::overload_cast(&VacuumAverage), + py::arg("expr")); m.def("VacuumAverage", - &VacuumAverage > >); + py::overload_cast< + const ExprPtr&, + const std::vector>&>( + &VacuumAverage), + py::arg("expr"), py::arg("op_connections")); } } // namespace sequant::python::mbpt diff --git a/python/test_sequant.py b/python/test_sequant.py index 5e04827e89..313735df2a 100644 --- a/python/test_sequant.py +++ b/python/test_sequant.py @@ -30,9 +30,9 @@ def _test_expressions(self): self.assertTrue((p+s).latex) def test_ccsd(self): - from _sequant.mbpt import A,H,T,T_,VacuumAverage,OpType - ccd = VacuumAverage( A(-2) * H() * T_(2) * T_(2), [(OpType.h, OpType.t)] ); # explicit list of operator connections .. - ccsd = VacuumAverage( A(-2) * H() * T(2) * T(2)); # .. is not needed since H and T are connected by default + from _sequant.mbpt import A,H,T,T_,VacuumAverage + ccd = VacuumAverage( A(-2) * H() * T_(2) * T_(2), [("h", "t")] ); # explicit list of operator connections .. + ccsd = VacuumAverage( A(-2) * H() * T(2) * T(2)); # .. is not needed since H and T are connected by default print (ccsd.latex) class String: diff --git a/tests/integration/antisymmetrizer.cpp b/tests/integration/antisymmetrizer.cpp index d3c26ed83c..5e762cc11a 100644 --- a/tests/integration/antisymmetrizer.cpp +++ b/tests/integration/antisymmetrizer.cpp @@ -9,11 +9,18 @@ using namespace sequant; +namespace { +std::wstring rdm_cumulant_label() { + static const std::wstring label = L"κ"; + return label; +} +} // namespace + void try_main() { using namespace sequant::mbpt; std::wcout << "START ANTISYMM_TEST: " << std::endl; - const auto cumulant = ex(optype2label.at(OpType::RDMCumulant), - bra{L"a_1"}, ket{L"i_1"}); + const auto cumulant = + ex(rdm_cumulant_label(), bra{L"a_1"}, ket{L"i_1"}); // const auto a =ex(L"a",bra{L"i_2", L"i_3"},ket{L"a_2", // L"a_3"}); const auto a = ex(cre{Index(L"i_2"), Index(L"i_3")}, @@ -23,10 +30,8 @@ void try_main() { antisymmetrize _a_cumulant(a_cumulant); std::wcout << to_latex_align(_a_cumulant.result) << std::endl; - auto cumulant2 = ex(optype2label.at(OpType::RDMCumulant), bra{L"a_2"}, - ket{L"i_2"}); - auto cumulant3 = ex(optype2label.at(OpType::RDMCumulant), bra{L"a_3"}, - ket{L"i_3"}); + auto cumulant2 = ex(rdm_cumulant_label(), bra{L"a_2"}, ket{L"i_2"}); + auto cumulant3 = ex(rdm_cumulant_label(), bra{L"a_3"}, ket{L"i_3"}); auto cumulant_3x = cumulant * cumulant2 * cumulant3; std::wcout << "cumulant_3x " << to_latex_align(cumulant_3x) << std::endl; antisymmetrize _cumulant_3x(cumulant_3x); @@ -38,8 +43,8 @@ void try_main() { antisymmetrize _a1_cumu1_cumu2(a1_cumu1_cumu2); std::wcout << to_latex_align(_a1_cumu1_cumu2.result) << std::endl; - auto two_body_cumu = ex(optype2label.at(OpType::RDMCumulant), - bra{L"a_2", L"a_3"}, ket{L"i_2", L"i_3"}); + auto two_body_cumu = ex(rdm_cumulant_label(), bra{L"a_2", L"a_3"}, + ket{L"i_2", L"i_3"}); auto a1_cumu2 = a1 * two_body_cumu; std::wcout << " a1 y2 " << to_latex_align(a1_cumu2) << std::endl; antisymmetrize _a1_cumu2(a1_cumu2); @@ -50,9 +55,8 @@ void try_main() { antisymmetrize _cumu1_cumu2(cumu1_cumu2); std::wcout << to_latex_align(_cumu1_cumu2.result) << std::endl; - auto cumu3 = - ex(optype2label.at(OpType::RDMCumulant), - bra{L"a_1", L"a_2", L"a_3"}, ket{L"i_1", L"i_2", L"i_3"}); + auto cumu3 = ex(rdm_cumulant_label(), bra{L"a_1", L"a_2", L"a_3"}, + ket{L"i_1", L"i_2", L"i_3"}); std::wcout << " y3 " << to_latex_align(cumu3) << std::endl; antisymmetrize _cumu3(cumu3); std::wcout << to_latex_align(_cumu3.result) << std::endl; diff --git a/tests/integration/eomcc.cpp b/tests/integration/eomcc.cpp index d70af1a58a..9ae2b7f2eb 100644 --- a/tests/integration/eomcc.cpp +++ b/tests/integration/eomcc.cpp @@ -197,6 +197,8 @@ int main(int argc, char* argv[]) { sequant::set_default_context( sequant::Context({.index_space_registry_shared_ptr = make_min_sr_spaces(), .vacuum = Vacuum::SingleProduct})); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); TensorCanonicalizer::register_instance( std::make_shared()); diff --git a/tests/integration/eval/btas/main.cpp b/tests/integration/eval/btas/main.cpp index c748e49f44..ef0b8b3be7 100644 --- a/tests/integration/eval/btas/main.cpp +++ b/tests/integration/eval/btas/main.cpp @@ -66,6 +66,8 @@ int main(int argc, char* argv[]) { CanonicalizationMethod::Complete)}); TensorCanonicalizer::register_instance( std::make_shared()); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); // for optimization tests, set occupied and unoccupied index extents { diff --git a/tests/integration/eval/ta/main.cpp b/tests/integration/eval/ta/main.cpp index 9511903609..0e59193e04 100644 --- a/tests/integration/eval/ta/main.cpp +++ b/tests/integration/eval/ta/main.cpp @@ -84,6 +84,8 @@ int main(int argc, char* argv[]) { CanonicalizationMethod::Complete)}); TensorCanonicalizer::register_instance( std::make_shared()); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); // for optimization tests, set occupied and unoccupied index extents { diff --git a/tests/integration/osstcc.cpp b/tests/integration/osstcc.cpp index 628e8be3e3..09be220ebc 100644 --- a/tests/integration/osstcc.cpp +++ b/tests/integration/osstcc.cpp @@ -35,6 +35,8 @@ int main(int argc, char* argv[]) { .canonicalization_options = CanonicalizeOptions::default_options().copy_and_set( CanonicalizationMethod::Complete)}); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); TensorCanonicalizer::register_instance( std::make_shared()); diff --git a/tests/integration/srcc.cpp b/tests/integration/srcc.cpp index 5175a2135d..0c58eacbe4 100644 --- a/tests/integration/srcc.cpp +++ b/tests/integration/srcc.cpp @@ -271,7 +271,8 @@ int main(int argc, char* argv[]) { const std::string uocc_type_str = argc > 3 ? argv[3] : "std"; const mbpt::CSV uocc_type = str2uocc.at(uocc_type_str); - auto mbpt_ctx = set_scoped_default_mbpt_context(mbpt::Context(uocc_type)); + set_default_mbpt_context( + {.csv = uocc_type, .op_registry_ptr = make_minimal_registry()}); const std::string spbasis_str = argc > 4 ? argv[4] : "so"; const SPBasis spbasis = str2spbasis.at(spbasis_str); diff --git a/tests/integration/stcc.cpp b/tests/integration/stcc.cpp index 3df305fb27..06e8be45f3 100644 --- a/tests/integration/stcc.cpp +++ b/tests/integration/stcc.cpp @@ -40,6 +40,8 @@ int main(int argc, char* argv[]) { CanonicalizationMethod::Complete)}); TensorCanonicalizer::register_instance( std::make_shared()); + mbpt::set_default_mbpt_context( + {.op_registry_ptr = mbpt::make_minimal_registry()}); #ifndef NDEBUG const size_t DEFAULT_NMAX = 3; diff --git a/tests/unit/test_main.cpp b/tests/unit/test_main.cpp index 9b682becc1..dddf6054a8 100644 --- a/tests/unit/test_main.cpp +++ b/tests/unit/test_main.cpp @@ -57,6 +57,11 @@ int main(int argc, char* argv[]) { // ... or can instead selectively set/unset particular logging flags // Logger::instance().wick_contract = true; + // MBPT Context setup + auto mbpt_ctx = mbpt::Context( + {.csv = mbpt::CSV::No, .op_registry_ptr = mbpt::make_legacy_registry()}); + sequant::mbpt::set_default_mbpt_context(mbpt_ctx); + #ifdef SEQUANT_HAS_TILEDARRAY auto& world = TA::initialize(argc, argv); TA::set_default_world(world); diff --git a/tests/unit/test_mbpt.cpp b/tests/unit/test_mbpt.cpp index abd260e615..5bbf8a3def 100644 --- a/tests/unit/test_mbpt.cpp +++ b/tests/unit/test_mbpt.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -34,6 +35,162 @@ #include "SeQuant/core/utility/debug.hpp" TEST_CASE("mbpt", "[mbpt]") { + SECTION("registry") { + using namespace sequant::mbpt; + + SECTION("empty-registry") { + OpRegistry registry; + REQUIRE(registry.ops().empty()); + REQUIRE_FALSE(registry.contains(L"T")); + } + + SECTION("add-operators") { + OpRegistry registry; + registry.add(L"T", OpClass::ex); + registry.add(L"L", OpClass::deex); + registry.add(L"F", OpClass::gen); + + REQUIRE(registry.contains(L"T")); + REQUIRE(registry.contains(L"L")); + REQUIRE(registry.contains(L"F")); + REQUIRE(registry.ops().size() == 3); + + REQUIRE(registry.to_class(L"T") == OpClass::ex); + REQUIRE(registry.to_class(L"L") == OpClass::deex); + REQUIRE(registry.to_class(L"F") == OpClass::gen); + } + + SECTION("remove-operators") { + OpRegistry registry; + registry.add(L"T", OpClass::ex); + registry.add(L"L", OpClass::deex); + + REQUIRE(registry.contains(L"T")); + registry.remove(L"T"); + REQUIRE_FALSE(registry.contains(L"T")); + REQUIRE(registry.contains(L"L")); + } + + SECTION("clone-registry") { + OpRegistry registry; + registry.add(L"T", OpClass::ex); + registry.add(L"L", OpClass::deex); + + auto cloned = registry.clone(); + REQUIRE(cloned.contains(L"T")); + REQUIRE(cloned.contains(L"L")); + REQUIRE(cloned.ops().size() == registry.ops().size()); + } + + SECTION("purge-registry") { + OpRegistry registry; + registry.add(L"T", OpClass::ex); + registry.add(L"L", OpClass::deex); + + REQUIRE_FALSE(registry.ops().empty()); + registry.purge(); + REQUIRE(registry.ops().empty()); + } + + SECTION("reserved-labels") { + OpRegistry registry; + // should not be able to add reserved labels + REQUIRE_THROWS(registry.add(reserved::antisymm_label(), OpClass::gen)); + REQUIRE_THROWS(registry.add(reserved::symm_label(), OpClass::gen)); + } + + SECTION("duplicate-operators") { + OpRegistry registry; + registry.add(L"T", OpClass::ex); + // should not be able to add duplicate + REQUIRE_THROWS(registry.add(L"T", OpClass::deex)); + } + } // SECTION("registry") + + SECTION("context") { + using namespace sequant::mbpt; + + SECTION("default-context") { + auto default_ctx = get_default_mbpt_context(); + // check default values + REQUIRE(default_ctx.csv() == CSV::No); + REQUIRE_NOTHROW(default_ctx.csv()); + REQUIRE_NOTHROW(default_ctx.op_registry()); + } + + SECTION("context-with-CSV") { + auto ctx = Context({.csv = CSV::Yes}); + REQUIRE(ctx.csv() == CSV::Yes); + + ctx.set(CSV::No); + REQUIRE(ctx.csv() == CSV::No); + } + + SECTION("context-with-registry") { + auto reg = make_minimal_registry(); + auto ctx = Context({.op_registry_ptr = reg}); + + REQUIRE(ctx.csv() == CSV::No); + REQUIRE(ctx.op_registry() != nullptr); + REQUIRE(ctx.op_registry() == reg); + } + + SECTION("set-operator-registry") { + auto ctx = Context(); + auto reg = make_legacy_registry(); + + ctx.set(reg); + REQUIRE(ctx.op_registry() != nullptr); + REQUIRE(ctx.op_registry() == reg); + } + + SECTION("clone-context") { + auto reg = std::make_shared(); + reg->add(L"T", OpClass::ex); + + auto ctx = Context({.csv = CSV::Yes, .op_registry_ptr = reg}); + auto cloned = ctx.clone(); + + REQUIRE(cloned.csv() == ctx.csv()); + REQUIRE(cloned.op_registry() != ctx.op_registry()); // Different pointer + REQUIRE(cloned.op_registry()->contains(L"T")); + } + + SECTION("context-equality") { + auto reg1 = make_legacy_registry(); + auto reg2 = make_minimal_registry(); + auto ctx1 = Context({.csv = CSV::Yes, .op_registry_ptr = reg1}); + auto ctx2 = Context({.csv = CSV::Yes, .op_registry_ptr = reg1}); + auto ctx3 = Context({.csv = CSV::No, .op_registry_ptr = reg1}); + auto ctx4 = Context({.csv = CSV::Yes, .op_registry_ptr = reg2}); + + REQUIRE(ctx1 == ctx2); + REQUIRE(ctx1 != ctx3); + REQUIRE(ctx1 != ctx4); + } + + SECTION("scoped-context") { + auto original_ctx = get_default_mbpt_context(); + auto original_csv = original_ctx.csv(); + + { + auto reg1 = make_legacy_registry(); + auto ctx = Context({.csv = CSV::Yes, .op_registry_ptr = reg1}); + auto _ = set_scoped_default_mbpt_context(ctx); + + auto current_ctx = get_default_mbpt_context(); + REQUIRE(current_ctx.csv() == CSV::Yes); + REQUIRE(current_ctx.op_registry() == reg1); + REQUIRE(current_ctx.op_registry()->contains(L"t")); + REQUIRE(current_ctx.op_registry()->contains(L"λ")); + } + + // after scope, original context should be restored + auto restored_ctx = get_default_mbpt_context(); + REQUIRE(restored_ctx.csv() == original_csv); + } + } // SECTION("context") + SECTION("nbody_operators") { using namespace sequant; @@ -531,22 +688,23 @@ TEST_CASE("mbpt", "[mbpt]") { get_default_context().index_space_registry()->retrieve(L"z")); using namespace mbpt; - REQUIRE_NOTHROW(op::H_pt(1, {.nbatch = 1})); - REQUIRE_NOTHROW(op::H_pt(2, {.nbatch = 2})); - REQUIRE_NOTHROW( - op::Λ_pt(3, {.batch_ordinals = {3, 4, 5}, .skip1 = true})); - REQUIRE_NOTHROW(op::T_pt(1, {.nbatch = 20})); + REQUIRE_NOTHROW(op::H_pt(1, {.order = 1, .nbatch = 1})); + REQUIRE_NOTHROW(op::H_pt(2, {.order = 1, .nbatch = 2})); + REQUIRE_NOTHROW(op::Λ_pt( + 3, {.order = 1, .batch_ordinals = {3, 4, 5}, .skip1 = true})); + REQUIRE_NOTHROW(op::T_pt(1, {.order = 1, .nbatch = 20})); // invalid usages #if SEQUANT_ASSERT_BEHAVIOR == SEQUANT_ASSERT_THROW // cannot set both nbatch and batch_ordinals - REQUIRE_THROWS_AS(op::H_pt(2, {.nbatch = 2, .batch_ordinals = {1, 2}}), - sequant::Exception); + REQUIRE_THROWS_AS( + op::H_pt(2, {.order = 1, .nbatch = 2, .batch_ordinals = {1, 2}}), + sequant::Exception); // all ordinals must be unique - REQUIRE_THROWS_AS(op::H_pt(2, {.batch_ordinals = {1, 2, 2}}), + REQUIRE_THROWS_AS(op::H_pt(2, {.order = 1, .batch_ordinals = {1, 2, 2}}), sequant::Exception); // ordinals must be sorted - REQUIRE_THROWS_AS(op::H_pt(1, {.batch_ordinals = {3, 2}}), + REQUIRE_THROWS_AS(op::H_pt(1, {.order = 1, .batch_ordinals = {3, 2}}), sequant::Exception); #endif @@ -554,9 +712,9 @@ TEST_CASE("mbpt", "[mbpt]") { auto h0 = op::H_pt(1); REQUIRE(to_latex(h0) == L"{\\hat{h¹}}"); - auto h1 = op::H_pt(1, {.nbatch = 1}); - auto h1_2 = op::H_pt(1, {.batch_ordinals = {1, 2}}); - auto pt1 = op::T_pt(2, {.batch_ordinals = {1}}); + auto h1 = op::H_pt(1, {.order = 1, .nbatch = 1}); + auto h1_2 = op::H_pt(1, {.order = 1, .batch_ordinals = {1, 2}}); + auto pt1 = op::T_pt(2, {.order = 1, .batch_ordinals = {1}}); auto sum0 = h0 + h1; simplify(sum0); @@ -570,8 +728,8 @@ TEST_CASE("mbpt", "[mbpt]") { simplify(sum2); // std::wcout << "sum2: " << to_latex(sum2) << std::endl; REQUIRE(to_latex(sum2) == - L"{ \\bigl({\\hat{h¹}}{[{z}_{1}]} + {\\hat{t¹}_{2}}{[{z}_{1}]} + " - L"{\\hat{t¹}_{1}}{[{z}_{1}]}\\bigr) }"); + L"{ \\bigl({\\hat{t¹}_{1}}{[{z}_{1}]} + {\\hat{h¹}}{[{z}_{1}]} + " + L"{\\hat{t¹}_{2}}{[{z}_{1}]}\\bigr) }"); auto sum3 = h1 + h1_2; simplify(sum3); @@ -694,8 +852,8 @@ SECTION("SRSO Fock") { SECTION("SRSO-PNO") { using sequant::mbpt::Context; - auto mbpt_ctx = - sequant::mbpt::set_scoped_default_mbpt_context(Context(mbpt::CSV::Yes)); + auto mbpt_ctx = sequant::mbpt::set_scoped_default_mbpt_context( + Context({.csv = CSV::Yes, .op_registry_ptr = make_minimal_registry()})); // H2**T2**T2 -> R2 SECTION("wick(H2**T2**T2 -> R2)") {