.. _program_listing_file_SeQuant_domain_mbpt_op.cpp:

Program Listing for File op.cpp
===============================

|exhale_lsh| :ref:`Return to documentation for file <file_SeQuant_domain_mbpt_op.cpp>` (``SeQuant/domain/mbpt/op.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include <SeQuant/domain/mbpt/context.hpp>
   #include <SeQuant/domain/mbpt/op.hpp>
   
   #include <SeQuant/core/math.hpp>
   #include <SeQuant/core/op.hpp>
   #include <SeQuant/core/tensor.hpp>
   #include <SeQuant/core/wick.hpp>
   
   #include <stdexcept>
   
   namespace sequant::mbpt {
   
   std::vector<std::wstring> cardinal_tensor_labels() {
     return {L"κ", L"γ", L"Γ", L"A", L"S", L"P", L"L",  L"λ", L"λ¹",
             L"h", L"f", L"f̃", L"g", L"θ", L"t", L"t¹", L"R", L"F",
             L"X", L"μ", L"V", L"Ṽ", L"B", L"U", L"GR", L"C", overlap_label(),
             L"a", L"ã", L"b", L"b̃", 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
   // the corresponding active hole or active particle space is a base space.
   qns_t excitation_type_qns(std::size_t K, const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, K};
       result[1] = {0ul, K};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       // are the active qp spaces base spaces?
       bool aps_base = isr->is_base(isr->particle_space(SQN));
       bool ahs_base = isr->is_base(isr->hole_space(SQN));
   
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // ex -> creators in particle space
         if (includes(isr->particle_space(SQN).type(), base_space.type())) {
           result[i * 2] = {aps_base ? K : 0ul, K};
         }
         // ex -> annihilators in hole space
         if (includes(isr->hole_space(SQN).type(), base_space.type())) {
           result[i * 2 + 1] = {ahs_base ? K : 0ul, K};
         }
       }
     }
     return result;
   }
   
   qns_t interval_excitation_type_qns(std::size_t K,
                                      const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, K};
       result[1] = {0ul, K};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
   
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // ex -> creators in particle space
         if (includes(isr->particle_space(SQN).type(), base_space.type())) {
           result[i * 2] = {0ul, K};
         }
         // ex -> annihilators in hole space
         if (includes(isr->hole_space(SQN).type(), base_space.type())) {
           result[i * 2 + 1] = {0ul, K};
         }
       }
     }
     return result;
   }
   
   qns_t deexcitation_type_qns(std::size_t K,
                               const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, K};
       result[1] = {0ul, K};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       bool aps_base = isr->is_base(isr->particle_space(SQN));
       bool ahs_base = isr->is_base(isr->hole_space(SQN));
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // deex -> annihilators in particle space
         if (includes(isr->particle_space(SQN).type(), base_space.type())) {
           result[i * 2 + 1] = {aps_base ? K : 0ul, K};
         }
         // deex -> creators in hole space
         if (includes(isr->hole_space(SQN).type(), base_space.type())) {
           result[i * 2] = {ahs_base ? K : 0ul, K};
         }
       }
     }
     return result;
   }
   
   qns_t interval_deexcitation_type_qns(std::size_t K,
                                        const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, K};
       result[1] = {0ul, K};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // deex -> annihilators in particle space
         if (includes(isr->particle_space(SQN).type(), base_space.type())) {
           result[i * 2 + 1] = {0ul, K};
         }
         // deex -> creators in hole space
         if (includes(isr->hole_space(SQN).type(), base_space.type())) {
           result[i * 2] = {0ul, K};
         }
       }
     }
     return result;
   }
   
   qns_t general_type_qns(std::size_t K) {
     qnc_t result;
     for (int i = 0; i < result.size(); i++) {
       result[i] = {0, K};
     }
     return result;
   }
   
   qns_t generic_excitation_qns(std::size_t particle_rank, std::size_t hole_rank,
                                IndexSpace particle_space, IndexSpace hole_space,
                                const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, hole_rank};
       result[1] = {0ul, particle_rank};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       bool aps_base = isr->is_base(isr->particle_space(SQN));
       bool ahs_base = isr->is_base(isr->hole_space(SQN));
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // ex -> creators in particle space
         if (includes(particle_space.type(), base_space.type())) {
           result[i * 2] = {aps_base ? particle_rank : 0ul,
                            particle_rank};  // creators
         }
         // ex -> annihilators in hole space
         if (includes(hole_space.type(), base_space.type())) {
           result[i * 2 + 1] = {ahs_base ? hole_rank : 0ul,
                                hole_rank};  // annihilators
         }
       }
     }
     return result;
   }
   
   qns_t generic_deexcitation_qns(std::size_t particle_rank, std::size_t hole_rank,
                                  IndexSpace particle_space, IndexSpace hole_space,
                                  const IndexSpace::QuantumNumbers SQN) {
     qnc_t result;
     if (get_default_context().vacuum() == Vacuum::Physical) {
       result[0] = {0ul, particle_rank};
       result[1] = {0ul, hole_rank};
     } else {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       bool aps_base = isr->is_base(isr->particle_space(SQN));
       bool ahs_base = isr->is_base(isr->hole_space(SQN));
       for (int i = 0; i < base_spaces.size(); i++) {
         const auto& base_space = base_spaces[i];
         result[i * 2] = {0ul, 0ul};
         result[i * 2 + 1] = {0ul, 0ul};
         if (base_space.qns() != SQN) continue;
   
         // deex -> creators in hole space
         if (includes(hole_space.type(), base_space.type())) {
           result[i * 2] = {ahs_base ? hole_rank : 0ul, hole_rank};  // creators
         }
         // deex -> annihilators in particle space
         if (includes(particle_space.type(), base_space.type())) {
           result[i * 2 + 1] = {aps_base ? particle_rank : 0ul,
                                particle_rank};  // annihilators
         }
       }
     }
     return result;
   }
   
   // By counting the number of contractions between indices of proper type, we can
   // know the quantum numbers of a combined result.
   qns_t combine(qns_t a, qns_t b) {
     assert(a.size() == b.size());
     qns_t result;
   
     if (get_default_context().vacuum() == Vacuum::Physical) {
       qns_t result;
       const auto ncontr = qninterval_t{0, std::min(b[0].upper(), a[1].upper())};
       const auto nc = nonnegative(a[0] + b[0] - ncontr);
       const auto na = nonnegative(a[1] + b[1] - ncontr);
       result[0] = nc;
       result[1] = na;
       return result;
     } else if (get_default_context().vacuum() == Vacuum::SingleProduct) {
       auto isr = get_default_context().index_space_registry();
       const auto& base_spaces = isr->base_spaces();
       for (auto i = 0; i < base_spaces.size(); i++) {
         auto cre = i * 2;
         auto ann = (i * 2) + 1;
         auto base_is_fermi_occupied = isr->is_pure_occupied(
             base_spaces[i]);  // need to distinguish particle and hole
                               // contractions.
         auto ncontr_space =
             base_is_fermi_occupied
                 ? qninterval_t{0, std::min(b[ann].upper(), a[cre].upper())}
                 : qninterval_t{0, std::min(b[cre].upper(), a[ann].upper())};
         auto nc_space = nonnegative(b[cre] + a[cre] - ncontr_space);
         auto na_space = nonnegative(b[ann] + a[ann] - ncontr_space);
         result[cre] = nc_space;
         result[ann] = na_space;
       }
       return result;
     } else {
       throw "unsupported vacuum context.";
     }
   }
   
   // must be defined including op.ipp since it's used there
   template <>
   bool is_vacuum<qns_t>(qns_t qns) {
     return qns == qns_t{};
   }
   
   }  // namespace sequant::mbpt
   
   namespace sequant {
   
   mbpt::qns_t adjoint(mbpt::qns_t qns) {
     mbpt::qns_t new_qnst;
     new_qnst.resize(qns.size());
     auto is_even = [](int i) { return i % 2 == 0; };
     for (int i = 0; i < qns.size(); i++) {
       if (is_even(i)) {
         new_qnst[i + 1] = qns[i];
       } else {
         new_qnst[i - 1] = qns[i];
       }
     }
     return new_qnst;
   }
   
   template <Statistics S>
   std::wstring to_latex(const mbpt::Operator<mbpt::qns_t, S>& op) {
     using namespace sequant::mbpt;
   
     auto result = L"{\\hat{" + utf_to_latex(op.label()) + L"}";
   
     // check if operator has adjoint label, remove if present for base label
     auto base_lbl = sequant::to_wstring(op.label());
     bool is_adjoint = false;
     if (base_lbl.back() == adjoint_label) {
       is_adjoint = true;
       base_lbl.pop_back();
     }
   
     auto it = label2optype.find(base_lbl);
     OpType optype = OpType::invalid;
     if (it != label2optype.end()) {  // handle special cases
       optype = it->second;
       if (to_class(optype) == OpClass::gen) {
         if (optype == OpType::θ) {  // special case for θ
           result += L"_{" + std::to_wstring(op()[0].upper()) + L"}";
         }
         result += L"}";
         return result;
       }
     }
     std::wstring baseline_char = is_adjoint ? L"^" : L"_";
     if (get_default_context().vacuum() == Vacuum::Physical) {
       if (op()[0] == op()[1]) {  // particle conserving
         result += L"_{" + std::to_wstring(op()[0].lower()) + L"}";
       } else {  // non-particle conserving
         result += L"_{" + std::to_wstring(op()[1].lower()) + L"}^{" +
                   std::to_wstring(op()[0].lower()) + L"}";
       }
     } else {  // single product vacuum
       auto nann_p = is_adjoint ? op().ncre_particles() : op().nann_particles();
       auto ncre_h = is_adjoint ? op().nann_holes() : op().ncre_holes();
       auto ncre_p = is_adjoint ? op().nann_particles() : op().ncre_particles();
       auto nann_h = is_adjoint ? op().ncre_holes() : op().nann_holes();
   
       if (!is_definite(nann_p) || !is_definite(ncre_h) || !is_definite(ncre_p) ||
           !is_definite(nann_h)) {
         throw std::invalid_argument(
             "to_latex(const Operator<qns_t, S>& op): "
             "can only handle generic operators with definite cre/ann numbers");
       }
   
       // pure quasiparticle creator/annihilator?
       const auto qprank_cre = ncre_p.lower() + nann_h.lower();
       const auto qprank_ann = nann_p.lower() + ncre_h.lower();
       const auto qppure = qprank_cre == 0 || qprank_ann == 0;
       if (qppure) {
         if (qprank_cre) {
           if (ncre_p.lower() == nann_h.lower()) {  // q-particle conserving
             result +=
                 baseline_char + L"{" + std::to_wstring(nann_h.lower()) + L"}";
           } else {  // particle non-conserving
             result += baseline_char + L"{" + std::to_wstring(nann_h.lower()) +
                       L"," + std::to_wstring(ncre_p.lower()) + L"}";
           }
         } else {
           if (ncre_h.lower() == nann_p.lower()) {  // q-particle conserving
             result +=
                 baseline_char + L"{" + std::to_wstring(ncre_h.lower()) + L"}";
           } else {  // q-particle non-conserving
             result += baseline_char + L"{" + std::to_wstring(ncre_h.lower()) +
                       L"," + std::to_wstring(nann_p.lower()) + L"}";
           }
         }
       } else {  // not pure qp creator/annihilator
         result += L"_{" + std::to_wstring(nann_h.lower()) + L"," +
                   std::to_wstring(ncre_p.lower()) + L"}^{" +
                   std::to_wstring(ncre_h.lower()) + L"," +
                   std::to_wstring(nann_p.lower()) + L"}";
       }
     }
     result += L"}";
     return result;
   }
   
   }  // namespace sequant
   
   #include <SeQuant/domain/mbpt/op.ipp>
   
   namespace sequant::mbpt {
   
   template <Statistics S>
   OpMaker<S>::OpMaker(OpType op) : op_(op) {}
   
   template <Statistics S>
   OpMaker<S>::OpMaker(OpType op, ncre nc, nann na) {
     op_ = op;
     assert(nc > 0 || na > 0);
     switch (to_class(op)) {
       case OpClass::ex:
         cre_spaces_ = decltype(cre_spaces_)(nc, get_particle_space(Spin::any));
         ann_spaces_ = decltype(ann_spaces_)(na, get_hole_space(Spin::any));
         break;
       case OpClass::deex:
         cre_spaces_ = decltype(cre_spaces_)(nc, get_hole_space(Spin::any));
         ann_spaces_ = decltype(ann_spaces_)(na, get_particle_space(Spin::any));
         break;
       case OpClass::gen:
         cre_spaces_ = decltype(cre_spaces_)(nc, get_complete_space(Spin::any));
         ann_spaces_ = decltype(ann_spaces_)(na, get_complete_space(Spin::any));
         break;
     }
   }
   
   template <Statistics S>
   OpMaker<S>::OpMaker(OpType op, std::size_t rank)
       : OpMaker(op, ncre(rank), nann(rank)) {}
   
   template <Statistics S>
   OpMaker<S>::OpMaker(OpType op, ncre nc, nann na,
                       const cre<IndexSpace>& cre_space,
                       const ann<IndexSpace>& ann_space) {
     op_ = op;
     assert(nc > 0 || na > 0);
     cre_spaces_ = decltype(cre_spaces_)(nc, cre_space);
     ann_spaces_ = decltype(ann_spaces_)(na, ann_space);
   }
   
   template <Statistics S>
   ExprPtr OpMaker<S>::operator()(std::optional<UseDepIdx> dep,
                                  std::optional<Symmetry> 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) {
         for (auto&& s : cre_spaces_) {
           assert(isr->contains_unoccupied(s));
         }
         dep = UseDepIdx::Bra;
       } else if (to_class(op_) == OpClass::deex) {
         for (auto&& s : ann_spaces_) {
           assert(isr->contains_unoccupied(s));
         }
         dep = UseDepIdx::Ket;
       } else {
         dep = UseDepIdx::None;
       }
     }
   
     return make(
         cre_spaces_, ann_spaces_,
         [this, opsymm_opt](const auto& creidxs, const auto& annidxs,
                            Symmetry opsymm) {
           return ex<Tensor>(to_wstring(op_), bra(creidxs), ket(annidxs),
                             opsymm_opt ? *opsymm_opt : opsymm);
         },
         dep ? *dep : UseDepIdx::None);
   }
   
   template class OpMaker<Statistics::FermiDirac>;
   template class OpMaker<Statistics::BoseEinstein>;
   
   template class Operator<qns_t, Statistics::FermiDirac>;
   template class Operator<qns_t, Statistics::BoseEinstein>;
   
   inline namespace op {
   
   namespace tensor {
   ExprPtr H_(std::size_t k) {
     assert(k > 0 && k <= 2);
     switch (k) {
       case 1:
         switch (get_default_context().vacuum()) {
           case Vacuum::Physical:
             return OpMaker<Statistics::FermiDirac>(OpType::h, 1)();
           case Vacuum::SingleProduct:
             return OpMaker<Statistics::FermiDirac>(OpType::f, 1)();
           case Vacuum::MultiProduct:
             return OpMaker<Statistics::FermiDirac>(OpType::f, 1)();
           default:
             abort();
         }
   
       case 2:
         return OpMaker<Statistics::FermiDirac>(OpType::g, 2)();
   
       default:
         abort();
     }
   }
   
   ExprPtr H(std::size_t k) {
     assert(k > 0 && k <= 2);
     return k == 1 ? tensor::H_(1) : tensor::H_(1) + tensor::H_(2);
   }
   
   ExprPtr F(bool use_tensor, IndexSpace reference_occupied) {
     if (use_tensor) {
       return OpMaker<Statistics::FermiDirac>(OpType::f, 1)();
     } else {                       // explicit density matrix construction
       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 isr = get_default_context().index_space_registry();
         return mbpt::OpMaker<Statistics::FermiDirac>::make(
             {isr->complete_space(Spin::any)}, {isr->complete_space(Spin::any)},
             [=](auto braidxs, auto ketidxs, Symmetry opsymm) {
               auto m1 = Index::make_tmp_index(occ_space);
               auto m2 = Index::make_tmp_index(occ_space);
               assert(opsymm == Symmetry::antisymm || opsymm == Symmetry::nonsymm);
               if (opsymm == Symmetry::antisymm) {
                 braidxs.push_back(m1);
                 ketidxs.push_back(m2);
                 return ex<Tensor>(to_wstring(mbpt::OpType::g),
                                   bra(std::move(braidxs)),
                                   ket(std::move(ketidxs)), Symmetry::antisymm) *
                        ex<Tensor>(to_wstring(mbpt::OpType::δ), bra{m2}, ket{m1},
                                   Symmetry::nonsymm);
               } else {  // opsymm == Symmetry::nonsymm
                 auto braidx_J = braidxs;
                 braidx_J.push_back(m1);
                 auto ketidxs_J = ketidxs;
                 ketidxs_J.push_back(m2);
                 auto braidx_K = braidxs;
                 braidx_K.push_back(m1);
                 auto ketidxs_K = ketidxs;
                 using std::begin;
                 ketidxs_K.emplace(begin(ketidxs_K), m2);
                 return (ex<Tensor>(to_wstring(mbpt::OpType::g),
                                    bra(std::move(braidx_J)),
                                    ket(std::move(ketidxs_J)), Symmetry::nonsymm) -
                         ex<Tensor>(
                             to_wstring(mbpt::OpType::g), bra(std::move(braidx_K)),
                             ket(std::move(ketidxs_K)), Symmetry::nonsymm)) *
                        ex<Tensor>(to_wstring(mbpt::OpType::δ), bra{m2}, ket{m1},
                                   Symmetry::nonsymm);
               }
             });
       };
       auto isr = get_default_context().index_space_registry();
       return OpMaker<Statistics::FermiDirac>(OpType::h, 1)() +
              make_g_contribution(reference_occupied);
     }
   }
   
   ExprPtr θ(std::size_t K) {
     return OpMaker<Statistics::FermiDirac>(OpType::θ, K)();
   }
   
   ExprPtr T_(std::size_t K) {
     return OpMaker<Statistics::FermiDirac>(OpType::t, K)();
   }
   
   ExprPtr T(std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = skip1 ? 2ul : 1ul; k <= K; ++k) {
       result += tensor::T_(k);
     }
     return result;
   }
   
   ExprPtr Λ_(std::size_t K) {
     return OpMaker<Statistics::FermiDirac>(OpType::λ, K)();
   }
   
   ExprPtr Λ(std::size_t K) {
     assert(K > 0);
   
     ExprPtr result;
     for (auto k = 1ul; k <= K; ++k) {
       result = k > 1 ? result + tensor::Λ_(k) : tensor::Λ_(k);
     }
     return result;
   }
   
   ExprPtr R_(nann na, ncre nc, const cre<IndexSpace>& cre_space,
              const ann<IndexSpace>& ann_space) {
     return OpMaker<Statistics::FermiDirac>(OpType::R, nc, na, cre_space,
                                            ann_space)();
   }
   ExprPtr R_(nₚ np, nₕ nh) {
     assert(np >= 0 && nh >= 0);
     return OpMaker<Statistics::FermiDirac>(OpType::R, ncre(np.value()),
                                            nann(nh.value()))();
   }
   
   ExprPtr L_(nann na, ncre nc, const cre<IndexSpace>& cre_space,
              const ann<IndexSpace>& ann_space) {
     return OpMaker<Statistics::FermiDirac>(OpType::L, nc, na, cre_space,
                                            ann_space)();
   }
   
   ExprPtr L_(nₚ np, nₕ nh) {
     assert(np >= 0 && nh >= 0);
     return OpMaker<Statistics::FermiDirac>(OpType::L, ncre(nh.value()),
                                            nann(np.value()))();
   }
   
   ExprPtr P(nₚ np, nₕ nh) {
     if (np != nh)
       assert(
           get_default_context().spbasis() != SPBasis::spinfree &&
           "Spinfree basis does not support non-particle conserving projectors");
     return get_default_context().spbasis() == SPBasis::spinfree
                ? tensor::S(-nh /* nh == np */)
                : tensor::A(-np, -nh);
   }
   
   ExprPtr A(nₚ np, nₕ nh) {
     assert(!(np == 0 && nh == 0));
     // if one of them is not zero, nh and np should have the same sign
     if (np != 0 && nh != 0) {
       assert((np > 0 && nh > 0) || (np < 0 && nh < 0));
     }
   
     container::svector<IndexSpace> creators;
     container::svector<IndexSpace> annihilators;
     if (nh > 0)  // ex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, nh))
         annihilators.emplace_back(get_hole_space(Spin::any));
     else if (nh < 0)  // deex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, -nh))
         creators.emplace_back(get_hole_space(Spin::any));
     if (np > 0)  // ex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, np))
         creators.emplace_back(get_particle_space(Spin::any));
     else if (np < 0)  // deex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, -np))
         annihilators.emplace_back(get_particle_space(Spin::any));
     // don't populate if rank is zero
   
     std::optional<OpMaker<Statistics::FermiDirac>::UseDepIdx> dep;
     if (get_default_mbpt_context().csv() == mbpt::CSV::Yes)
       dep = (np > 0 || nh > 0) ? OpMaker<Statistics::FermiDirac>::UseDepIdx::Bra
                                : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
     return OpMaker<Statistics::FermiDirac>(
         OpType::A, cre(creators), ann(annihilators))(dep, {Symmetry::antisymm});
   }
   
   ExprPtr S(std::int64_t K) {
     assert(K != 0);
     container::svector<IndexSpace> creators;
     container::svector<IndexSpace> annihilators;
     if (K > 0)  // ex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, K))
         annihilators.emplace_back(get_hole_space(Spin::any));
     else  // deex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, -K))
         creators.emplace_back(get_hole_space(Spin::any));
     if (K > 0)  // ex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, K))
         creators.emplace_back(get_particle_space(Spin::any));
     else  // deex
       for ([[maybe_unused]] auto i : ranges::views::iota(0, -K))
         annihilators.emplace_back(get_particle_space(Spin::any));
   
     std::optional<OpMaker<Statistics::FermiDirac>::UseDepIdx> dep;
     if (get_default_mbpt_context().csv() == mbpt::CSV::Yes)
       dep = K > 0 ? OpMaker<Statistics::FermiDirac>::UseDepIdx::Bra
                   : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
     return OpMaker<Statistics::FermiDirac>(
         OpType::S, cre(creators), ann(annihilators))(dep, {Symmetry::nonsymm});
   }
   
   ExprPtr H_pt(std::size_t order, std::size_t R) {
     assert(order == 1 &&
            "sequant::sr::H_pt(): only supports first order perturbation");
     assert(R > 0);
     return OpMaker<Statistics::FermiDirac>(OpType::h_1, R)();
   }
   
   ExprPtr T_pt_(std::size_t order, std::size_t K) {
     assert(order == 1 &&
            "sequant::sr::T_pt_(): only supports first order perturbation");
     return OpMaker<Statistics::FermiDirac>(OpType::t_1, K)();
   }
   
   ExprPtr T_pt(std::size_t order, std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) {
       result = k > 1 ? result + tensor::T_pt_(order, k) : tensor::T_pt_(order, k);
     }
     return result;
   }
   
   ExprPtr Λ_pt_(std::size_t order, std::size_t K) {
     assert(order == 1 &&
            "sequant::sr::Λ_pt_(): only supports first order perturbation");
     return OpMaker<Statistics::FermiDirac>(OpType::λ_1, K)();
   }
   
   ExprPtr Λ_pt(std::size_t order, std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) {
       result = k > 1 ? result + tensor::Λ_pt_(order, k) : tensor::Λ_pt_(order, k);
     }
     return result;
   }
   
   }  // namespace tensor
   
   ExprPtr H_(std::size_t k) {
     assert(k > 0 && k <= 2);
     switch (k) {
       case 1:
         return ex<op_t>(
             [vacuum = get_default_context().vacuum()]() -> std::wstring_view {
               switch (vacuum) {
                 case Vacuum::Physical:
                   return L"h";
                 case Vacuum::SingleProduct:
                   return L"f";
                 case Vacuum::MultiProduct:
                   return L"f";
                 default:
                   abort();
               }
             },
             [=]() -> ExprPtr { return tensor::H_(1); },
             [=](qnc_t& qns) {
               qnc_t op_qnc_t = general_type_qns(1);
               qns = combine(op_qnc_t, qns);
             });
   
       case 2:
         return ex<op_t>([]() -> std::wstring_view { return L"g"; },
                         [=]() -> ExprPtr { return tensor::H_(2); },
                         [=](qnc_t& qns) {
                           qnc_t op_qnc_t = general_type_qns(2);
                           qns = combine(op_qnc_t, qns);
                         });
   
       default:
         abort();
     }
   }
   
   ExprPtr H(std::size_t k) {
     assert(k > 0 && k <= 2);
     return k == 1 ? H_(1) : H_(1) + H_(2);
   }
   
   ExprPtr θ(std::size_t K) {
     assert(K > 0);
     return ex<op_t>([]() -> std::wstring_view { return L"θ"; },
                     [=]() -> ExprPtr { return tensor::θ(K); },
                     [=](qnc_t& qns) {
                       qnc_t op_qnc_t = general_type_qns(K);
                       qns = combine(op_qnc_t, qns);
                     });
   }
   
   ExprPtr T_(std::size_t K) {
     assert(K > 0);
     return ex<op_t>([]() -> std::wstring_view { return L"t"; },
                     [=]() -> ExprPtr { return tensor::T_(K); },
                     [=](qnc_t& qns) {
                       qnc_t op_qnc_t = excitation_type_qns(K);
                       qns = combine(op_qnc_t, qns);
                     });
   }
   
   ExprPtr T(std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = skip1 ? 2ul : 1ul; k <= K; ++k) {
       result += T_(k);
     }
     return result;
   }
   
   ExprPtr Λ_(std::size_t K) {
     assert(K > 0);
     return ex<op_t>([]() -> std::wstring_view { return L"λ"; },
                     [=]() -> ExprPtr { return tensor::Λ_(K); },
                     [=](qnc_t& qns) {
                       qnc_t op_qnc_t = deexcitation_type_qns(K);
                       qns = combine(op_qnc_t, qns);
                     });
   }
   
   ExprPtr Λ(std::size_t K) {
     assert(K > 0);
     ExprPtr result;
     for (auto k = 1ul; k <= K; ++k) {
       result = k > 1 ? result + Λ_(k) : Λ_(k);
     }
     return result;
   }
   
   ExprPtr F(bool use_f_tensor, IndexSpace occupied_density) {
     if (use_f_tensor) {
       return ex<op_t>(
           []() -> std::wstring_view { return L"f"; },
           [=]() -> ExprPtr { return tensor::F(true, occupied_density); },
           [=](qnc_t& qns) {
             qnc_t op_qnc_t = general_type_qns(1);
             qns = combine(op_qnc_t, qns);
           });
     } else {
       throw "non-tensor use at operator level not yet supported";
     }
   }
   
   ExprPtr A(nₚ np, nₕ nh) {
     assert(!(nh == 0 && np == 0));
     // if one of them is not zero, nh and np should have the same sign
     if (nh != 0 && np != 0) {
       assert((nh > 0 && np > 0) || (nh < 0 && np < 0));
     }
     // if np or nh is negative, it's a deexcitation operator
     const auto deexcitation = (np < 0 || nh < 0);
   
     auto particle_space = get_particle_space(Spin::any);
     auto hole_space = get_hole_space(Spin::any);
     return ex<op_t>([]() -> 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);
                       }
                     });
   }
   
   ExprPtr S(std::int64_t K) {
     assert(K != 0);
     return ex<op_t>([]() -> std::wstring_view { return L"S"; },
                     [=]() -> ExprPtr { return tensor::S(K); },
                     [=](qnc_t& qns) {
                       const std::size_t abs_K = std::abs(K);
                       if (K < 0) {
                         qnc_t op_qnc_t = deexcitation_type_qns(abs_K);
                         qns = combine(op_qnc_t, qns);
                       } else {
                         qnc_t op_qnc_t = excitation_type_qns(abs_K);
                         qns = combine(op_qnc_t, qns);
                       }
                     });
   }
   
   ExprPtr P(nₚ np, nₕ nh) {
     if (get_default_context().spbasis() == SPBasis::spinfree) {
       assert(nh == np &&
              "Only particle number conserving cases are supported with spinfree "
              "basis for now");
       const auto K = np;  // K = np = nh
       return S(-K);
     } else {
       assert(get_default_context().spbasis() == SPBasis::spinorbital);
       return A(-np, -nh);
     }
   }
   
   ExprPtr H_pt(std::size_t order, std::size_t R) {
     assert(R > 0);
     assert(order == 1 && "only first order perturbation is supported now");
     return ex<op_t>(
         []() -> std::wstring_view { return optype2label.at(OpType::h_1); },
         [=]() -> ExprPtr { return tensor::H_pt(order, R); },
         [=](qnc_t& qns) { qns = combine(general_type_qns(R), qns); });
   }
   
   ExprPtr T_pt_(std::size_t order, std::size_t K) {
     assert(K > 0);
     assert(order == 1 && "only first order perturbation is supported now");
     return ex<op_t>(
         []() -> std::wstring_view { return optype2label.at(OpType::t_1); },
         [=]() -> ExprPtr { return tensor::T_pt_(order, K); },
         [=](qnc_t& qns) { qns = combine(excitation_type_qns(K), qns); });
   }
   
   ExprPtr T_pt(std::size_t order, std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) {
       result = k > 1 ? result + T_pt_(order, k) : T_pt_(order, k);
     }
     return result;
   }
   
   ExprPtr Λ_pt_(std::size_t order, std::size_t K) {
     assert(K > 0);
     assert(order == 1 && "only first order perturbation is supported now");
     return ex<op_t>(
         []() -> std::wstring_view { return optype2label.at(OpType::λ_1); },
         [=]() -> ExprPtr { return tensor::Λ_pt_(order, K); },
         [=](qnc_t& qns) { qns = combine(deexcitation_type_qns(K), qns); });
   }
   
   ExprPtr Λ_pt(std::size_t order, std::size_t K, bool skip1) {
     assert(K > (skip1 ? 1 : 0));
     ExprPtr result;
     for (auto k = (skip1 ? 2ul : 1ul); k <= K; ++k) {
       result = k > 1 ? result + Λ_pt_(order, k) : Λ_pt_(order, k);
     }
     return result;
   }
   
   ExprPtr R_(nann na, ncre nc, const cre<IndexSpace>& cre_space,
              const ann<IndexSpace>& ann_space) {
     return ex<op_t>(
         []() -> std::wstring_view { return optype2label.at(OpType::R); },
         [=]() -> ExprPtr { return tensor::R_(na, nc, cre_space, ann_space); },
         [=](qnc_t& qns) {
           // ex -> creators in particle_space, annihilators in hole_space
           qns = combine(
               generic_excitation_qns(/*particle_rank*/ nc, /*hole_rank*/ na,
                                      cre_space, ann_space),
               qns);
         });
   }
   
   ExprPtr R_(nₚ np, nₕ nh) { return R_(nann(nh), ncre(np)); }
   
   ExprPtr L_(nann na, ncre nc, const cre<IndexSpace>& cre_space,
              const ann<IndexSpace>& ann_space) {
     return ex<op_t>(
         []() -> std::wstring_view { return optype2label.at(OpType::L); },
         [=]() -> ExprPtr { return tensor::L_(na, nc, cre_space, ann_space); },
         [=](qnc_t& qns) {
           // deex -> creators in hole_space, annihilators in particle_space
           qns = combine(
               generic_deexcitation_qns(
                   /*particle_rank*/ na, /*hole_rank*/ nc, ann_space, cre_space),
               qns);
         });
   }
   
   ExprPtr L_(nₚ np, nₕ nh) { return L_(nann(np), ncre(nh)); }
   
   ExprPtr R(nann na, ncre nc, const cre<IndexSpace>& cre_space,
             const ann<IndexSpace>& ann_space) {
     assert(na > 0 || nc > 0);
     ExprPtr result;
   
     for (std::int64_t ra = na, rc = nc; ra > 0 || rc > 0; --ra, --rc) {
       result += R_(nann(ra), ncre(rc), cre_space, ann_space);
     }
     return result;
   }
   
   ExprPtr R(nₚ np, nₕ nh) {
     assert(np >= 0 && nh >= 0);
     return R(nann(nh), ncre(np));
   }
   
   ExprPtr L(nann na, ncre nc, const cre<IndexSpace>& cre_space,
             const ann<IndexSpace>& ann_space) {
     assert(na > 0 || nc > 0);
     ExprPtr result;
   
     for (std::int64_t ra = na, rc = nc; ra > 0 || rc > 0; --ra, --rc) {
       result += L_(nann(ra), ncre(rc), cre_space, ann_space);
     }
     return result;
   }
   
   ExprPtr L(nₚ np, nₕ nh) { return L(nann(np), ncre(nh)); }
   
   bool can_change_qns(const ExprPtr& op_or_op_product, const qns_t target_qns,
                       const qns_t source_qns = {}) {
     qns_t qns = source_qns;
     if (op_or_op_product.is<Product>()) {
       const auto& op_product = op_or_op_product.as<Product>();
       for (auto& op_ptr : ranges::views::reverse(op_product.factors())) {
         assert(op_ptr->template is<op_t>());
         const auto& op = op_ptr->template as<op_t>();
         qns = op(qns);
       }
       return qns.overlaps_with(target_qns);
     } else if (op_or_op_product.is<op_t>()) {
       const auto& op = op_or_op_product.as<op_t>();
       qns = op();
       return qns.overlaps_with(target_qns);
     } else
       throw std::invalid_argument(
           "sequant::mbpt::sr::contains_rank(op_or_op_product): op_or_op_product "
           "must be mbpt::sr::op_t or Product thereof");
   }
   
   bool raises_vacuum_up_to_rank(const ExprPtr& op_or_op_product,
                                 const unsigned long k) {
     assert(op_or_op_product.is<op_t>() || op_or_op_product.is<Product>());
   
     return can_change_qns(op_or_op_product, interval_excitation_type_qns(k));
   }
   
   bool lowers_rank_or_lower_to_vacuum(const ExprPtr& op_or_op_product,
                                       const unsigned long k) {
     assert(op_or_op_product.is<op_t>() || op_or_op_product.is<Product>());
     return can_change_qns(op_or_op_product, qns_t{},
                           interval_excitation_type_qns(k));
   }
   
   bool raises_vacuum_to_rank(const ExprPtr& op_or_op_product,
                              const unsigned long k) {
     assert(op_or_op_product.is<op_t>() || op_or_op_product.is<Product>());
     return can_change_qns(op_or_op_product, excitation_type_qns(k));
   }
   
   bool lowers_rank_to_vacuum(const ExprPtr& op_or_op_product,
                              const unsigned long k) {
     assert(op_or_op_product.is<op_t>() || op_or_op_product.is<Product>());
     return can_change_qns(op_or_op_product, qns_t{}, excitation_type_qns(k));
   }
   
   #include <SeQuant/domain/mbpt/vac_av.ipp>
   
   namespace tensor {
   
   ExprPtr vac_av(ExprPtr expr, std::vector<std::pair<int, int>> nop_connections,
                  bool use_top) {
     simplify(expr);
     auto isr = get_default_context().index_space_registry();
     const auto spinorbital =
         get_default_context().spbasis() == SPBasis::spinorbital;
     // convention is to use different label for spin-orbital and spin-free RDM
     const auto rdm_label = spinorbital ? optype2label.at(OpType::RDM) : L"Γ";
   
     // only need full contractions if don't have any density outside of
     // the orbitals occupied in the vacuum
     bool full_contractions =
         (isr->reference_occupied_space() == isr->vacuum_occupied_space()) ? true
                                                                           : false;
     // N.B. reference < vacuum is not yet supported
     if (isr->reference_occupied_space().intersection(
             isr->vacuum_occupied_space()) != isr->vacuum_occupied_space()) {
       throw std::invalid_argument(
           "mbpt::tensor::vac_av: vacuum occupied orbitals must be same as or "
           "subset of the reference orbital set.");
     }
   
     FWickTheorem wick{expr};
     wick.use_topology(use_top).set_nop_connections(nop_connections);
     wick.full_contractions(full_contractions);
     auto result = wick.compute();
     simplify(result);
   
     if (Logger::instance().wick_stats) {
       std::wcout << "WickTheorem stats: # of contractions attempted = "
                  << wick.stats().num_attempted_contractions
                  << " # of useful contractions = "
                  << wick.stats().num_useful_contractions << std::endl;
     }
     // only need to handle the special case where the dense(at least partially
     // occupied) states, contain additional functions to the vacuum_occupied.
     // including a density occupied partition using a "single-reference" method
     // will replace FNOPs with RDMs. i.e. "multi-reference" RDM replacement rules
     // work in the limit of one reference.
     if (isr->reference_occupied_space() == IndexSpace::Type{} ||
         isr->reference_occupied_space(Spin::any) ==
             isr->vacuum_occupied_space(Spin::any)) {
       return result;
     } else {
       const auto target_rdm_space_type =
           get_default_context().vacuum() == Vacuum::SingleProduct
               ? isr->intersection(isr->particle_space(Spin::any),
                                   isr->hole_space(Spin::any))
               : isr->reference_occupied_space(Spin::any);
   
       // STEP1. replace NOPs by RDM
       auto replace_nop_with_rdm = [&rdm_label, spinorbital](ExprPtr& exptr) {
         auto replace = [&rdm_label, spinorbital](const auto& nop) -> ExprPtr {
           using index_container = container::svector<Index>;
           auto braidxs = nop.annihilators() |
                          ranges::views::transform(
                              [](const auto& op) { return op.index(); }) |
                          ranges::to<index_container>();
           auto ketidxs = nop.creators() |
                          ranges::views::transform(
                              [](const auto& op) { return op.index(); }) |
                          ranges::to<index_container>();
           assert(braidxs.size() ==
                  ketidxs.size());  // need to handle particle # violating case?
           const auto rank = braidxs.size();
           return ex<Tensor>(
               rdm_label, bra(std::move(braidxs)), ket(std::move(ketidxs)),
               rank > 1 && spinorbital ? Symmetry::antisymm : Symmetry::nonsymm);
         };
   
         if (exptr.template is<FNOperator>()) {
           exptr = replace(exptr.template as<FNOperator>());
         } else if (exptr.template is<BNOperator>()) {
           exptr = replace(exptr.template as<BNOperator>());
         }
       };
       result->visit(replace_nop_with_rdm, true);
   
       // STEP 2: project RDM indices onto the target RDM subspace
       // since RDM indices only make sense within a single TN expand + flatten
       // first, then do the projection individually for each TN
       expand(result);
       // flatten(result);  // TODO where is flatten?
       auto project_rdm_indices_to_target = [&](ExprPtr& exptr) {
         auto impl_for_single_tn = [&](ProductPtr& product_ptr) {
           // enlist all indices and count their instances
           auto for_each_index_in_tn = [](const auto& product_ptr,
                                          const auto& op) {
             ranges::for_each(product_ptr->factors(), [&](auto& factor) {
               auto tensor_ptr = std::dynamic_pointer_cast<AbstractTensor>(factor);
               if (tensor_ptr) {
                 ranges::for_each(tensor_ptr->_braket(),
                                  [&](auto& idx) { op(idx, *tensor_ptr); });
               }
             });
           };
   
           // compute external indices
           container::map<Index, std::size_t> indices_w_counts;
           auto retrieve_indices_with_counts = [&indices_w_counts](const auto& idx,
                                                                   auto&) {
             auto found_it = indices_w_counts.find(idx);
             if (found_it != indices_w_counts.end()) {
               found_it->second++;
             } else {
               indices_w_counts.emplace(idx, 1);
             }
           };
           for_each_index_in_tn(product_ptr, retrieve_indices_with_counts);
   
           container::set<Index> external_indices =
               indices_w_counts | ranges::views::filter([](auto& idx_cnt) {
                 auto& [idx, cnt] = idx_cnt;
                 return cnt == 1;
               }) |
               ranges::views::keys | ranges::to<container::set<Index>>;
   
           // extract RDM-only and all indices
           container::set<Index> rdm_indices;
           std::set<Index, Index::LabelCompare> all_indices;
           auto retrieve_rdm_and_all_indices = [&rdm_indices, &all_indices,
                                                &rdm_label](const auto& idx,
                                                            const auto& tensor) {
             all_indices.insert(idx);
             if (tensor._label() == rdm_label) {
               rdm_indices.insert(idx);
             }
           };
           for_each_index_in_tn(product_ptr, retrieve_rdm_and_all_indices);
   
           // compute RDM->target replacement rules
           container::map<Index, Index> replacement_rules;
           ranges::for_each(rdm_indices, [&](const Index& idx) {
             const auto target_type =
                 isr->intersection(idx.space(), target_rdm_space_type);
             if (target_type) {
               Index target = Index::make_tmp_index(target_type);
               replacement_rules.emplace(idx, target);
             }
           });
   
           if (false) {
             std::wcout << "expr = " << product_ptr->to_latex()
                        << "\n  external_indices = ";
             ranges::for_each(external_indices, [](auto& index) {
               std::wcout << index.label() << " ";
             });
             std::wcout << "\n  replrules = ";
             ranges::for_each(replacement_rules, [](auto& index) {
               std::wcout << to_latex(index.first) << "\\to"
                          << to_latex(index.second) << "\\,";
             });
             std::wcout.flush();
           }
   
           if (!replacement_rules.empty()) {
             sequant::detail::apply_index_replacement_rules(
                 product_ptr, replacement_rules, external_indices, all_indices,
                 isr);
           }
         };
   
         if (exptr.template is<Product>()) {
           auto product_ptr = exptr.template as_shared_ptr<Product>();
           impl_for_single_tn(product_ptr);
           exptr = product_ptr;
         } else {
           assert(exptr.template is<Sum>());
           auto result = std::make_shared<Sum>();
           for (auto& summand : exptr.template as<Sum>().summands()) {
             assert(summand.template is<Product>());
             auto result_summand = summand.template as<Product>().clone();
             auto product_ptr = result_summand.template as_shared_ptr<Product>();
             impl_for_single_tn(product_ptr);
             result->append(product_ptr);
           }
           exptr = result;
         }
       };
       project_rdm_indices_to_target(result);
   
       // rename dummy indices that might have been generated by
       // project_rdm_indices_to_target
       // + may combine terms
   
       // TensorCanonicalizer is given a custom comparer that moves active
       // indices to the front external-vs-internal trait still takes precedence
       auto current_index_comparer =
           TensorCanonicalizer::instance()->index_comparer();
       TensorCanonicalizer::instance()->index_comparer(
           [&](const Index& idx1, const Index& idx2) -> bool {
             auto active_space = isr->intersection(isr->particle_space(Spin::any),
                                                   isr->hole_space(Spin::any));
             const auto idx1_active = idx1.space().type() == active_space.type();
             const auto idx2_active = idx2.space().type() == active_space.type();
             if (idx1_active) {
               if (idx2_active)
                 return current_index_comparer(idx1, idx2);
               else
                 return true;
             } else {
               if (idx2_active)
                 return false;
               else
                 return current_index_comparer(idx1, idx2);
             }
           });
       simplify(result);
       TensorCanonicalizer::instance()->index_comparer(
           std::move(current_index_comparer));
   
       if (Logger::instance().wick_stats) {
         std::wcout << "WickTheorem stats: # of contractions attempted = "
                    << wick.stats().num_attempted_contractions
                    << " # of useful contractions = "
                    << wick.stats().num_useful_contractions << std::endl;
       }
       return result;
     }
   }
   
   }  // namespace tensor
   }  // namespace op
   
   bool can_change_qns(const ExprPtr& op_or_op_product, const qns_t target_qns,
                       const qns_t source_qns = {}) {
     qns_t qns = source_qns;
     if (op_or_op_product.is<Product>()) {
       const auto& op_product = op_or_op_product.as<Product>();
       for (auto& op_ptr : ranges::views::reverse(op_product.factors())) {
         assert(op_ptr->template is<op_t>());
         const auto& op = op_ptr->template as<op_t>();
         qns = op(qns);
       }
       return qns.overlaps_with(target_qns);
     } else if (op_or_op_product.is<op_t>()) {
       const auto& op = op_or_op_product.as<op_t>();
       qns = op();
       return qns.overlaps_with(target_qns);
     } else
       throw std::invalid_argument(
           "sequant::mbpt::sr::contains_rank(op_or_op_product): op_or_op_product "
           "must be mbpt::sr::op_t or Product thereof");
   }
   
   }  // namespace sequant::mbpt