Program Listing for File op.hpp

↰ Return to documentation for file (SeQuant/domain/mbpt/op.hpp)

//
// Created by Eduard Valeyev on 2019-03-26.
//

#ifndef SEQUANT_DOMAIN_MBPT_OP_HPP
#define SEQUANT_DOMAIN_MBPT_OP_HPP

#include <SeQuant/domain/mbpt/fwd.hpp>

#include <SeQuant/domain/mbpt/spin.hpp>

#include <SeQuant/core/attr.hpp>
#include <SeQuant/core/container.hpp>
#include <SeQuant/core/context.hpp>
#include <SeQuant/core/expr.hpp>
#include <SeQuant/core/hash.hpp>
#include <SeQuant/core/index.hpp>
#include <SeQuant/core/interval.hpp>
#include <SeQuant/core/math.hpp>
#include <SeQuant/core/op.hpp>
#include <SeQuant/core/rational.hpp>
#include <SeQuant/core/space.hpp>
#include <SeQuant/core/utility/strong.hpp>

#include <range/v3/iterator/basic_iterator.hpp>
#include <range/v3/range/conversion.hpp>
#include <range/v3/range/primitives.hpp>
#include <range/v3/view/map.hpp>
#include <range/v3/view/transform.hpp>
#include <range/v3/view/view.hpp>
#include <range/v3/view/zip.hpp>

#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <initializer_list>
#include <iterator>
#include <map>
#include <optional>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <utility>
#include <vector>

namespace sequant {
namespace mbpt {

DEFINE_STRONG_TYPE_FOR_INTEGER(nₚ, std::int64_t);  // define nₚ
DEFINE_STRONG_TYPE_FOR_INTEGER(nₕ, std::int64_t);  // define nₕ

#ifndef DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT
#define DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(OP)                      \
  inline ExprPtr OP(std::int64_t Rank) { return OP(nₚ(Rank), nₕ(Rank)); } \
  inline ExprPtr OP(nₚ Rank) { return OP(Rank, nₕ(Rank)); }
#endif  // DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT

template <typename QuantumNumbers>
bool is_vacuum(QuantumNumbers qns);

inline IndexSpace make_space(const IndexSpace::Type& type) {
  return get_default_context().index_space_registry()->retrieve(type,
                                                                Spin::any);
}

enum class OpType {
  h,
  f,
  f̃,
  θ,
  g,
  t,
  λ,
  A,
  S,
  L,
  R,
  R12,
  GR,
  C,
  RDM,
  RDMCumulant,
  δ,
  h_1,
  t_1,
  λ_1,
  invalid
};

inline const std::map<OpType, std::wstring> optype2label{
    {OpType::h, L"h"},
    {OpType::f, L"f"},
    {OpType::f̃, L"f̃"},
    {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::δ, L"δ"},
    {OpType::h_1, L"h¹"},
    {OpType::t_1, L"t¹"},
    {OpType::λ_1, L"λ¹"}};

inline const std::map<std::wstring, OpType> label2optype =
    ranges::views::zip(ranges::views::values(optype2label),
                       ranges::views::keys(optype2label)) |
    ranges::to<std::map<std::wstring, OpType>>();

enum class OpClass { ex, deex, gen };

std::vector<std::wstring> cardinal_tensor_labels();

std::wstring to_wstring(OpType op);
OpClass to_class(OpType op);



template <typename QuantumNumbers, Statistics S = Statistics::FermiDirac>
class Operator;

template <typename QuantumNumbers>
using FOperator = Operator<QuantumNumbers, Statistics::FermiDirac>;
template <typename QuantumNumbers>
using BOperator = Operator<QuantumNumbers, Statistics::BoseEinstein>;

template <Statistics S>
class Operator<void, S> : public Expr, public Labeled {
 protected:
  Operator() = default;
  Operator(std::function<std::wstring_view()> label_generator,
           std::function<ExprPtr()> tensor_form_generator)
      : label_generator_(std::move(label_generator)),
        tensor_form_generator_(tensor_form_generator) {}

 public:
  virtual ~Operator() = default;

  std::wstring_view label() const override {
    assert(label_generator_);
    return label_generator_();
  }

  virtual ExprPtr tensor_form() const {
    assert(tensor_form_generator_);
    return tensor_form_generator_();
  }

  bool is_cnumber() const override { return false; }

 private:
  std::function<std::wstring_view()> label_generator_;
  std::function<ExprPtr()> tensor_form_generator_;
};

using FOperatorBase = FOperator<void>;
using BOperatorBase = BOperator<void>;

struct default_qns_tag {};

// clang-format off

// clang-format on
template <typename QNV = std::int64_t, typename Tag = default_qns_tag>
class QuantumNumberChange
    : public container::svector<
          boost::numeric::interval<std::make_signed_t<QNV>>, 8> {
 public:
  using QNC = std::make_signed_t<QNV>;  // change in quantum numbers
  using interval_t = boost::numeric::interval<QNC>;
  using base_type =
      container::svector<boost::numeric::interval<std::make_signed_t<QNV>>, 8>;
  using this_type = QuantumNumberChange<QNV, Tag>;

  const std::size_t size() const {
    if (get_default_context().vacuum() == Vacuum::Physical) {
      return 2;
    } else if (get_default_context().vacuum() == Vacuum::SingleProduct) {
      auto isr = get_default_context().index_space_registry();
      const auto& isr_base_spaces = isr->base_spaces();
      assert(isr_base_spaces.size() > 0);
      return isr_base_spaces.size() * 2;
    } else {
      throw std::logic_error("unknown Vacuum type");
    }
  }

  QuantumNumberChange() {
    this->resize(this->size());
    assert(this->base().size() != 0);
    std::fill(this->begin(), this->end(), interval_t{});
  }

  template <typename I, typename Range,
            typename = std::enable_if_t<
                meta::is_range_v<std::remove_reference_t<Range>> &&
                std::is_convertible_v<I, interval_t>>>
  explicit QuantumNumberChange(Range&& i) : QuantumNumberChange() {
    assert(i.size() == size());
    std::copy(i.begin(), i.end(), this->begin());
  }

  template <typename I, typename = std::enable_if_t<std::is_convertible_v<
                            std::initializer_list<I>, interval_t>>>
  explicit QuantumNumberChange(
      std::initializer_list<std::initializer_list<I>> i)
      : QuantumNumberChange() {
    assert(i.size() == size());
#ifndef NDEBUG
    if (std::find_if(i.begin(), i.end(),
                     [](const auto& ii) { return ii.size() != 2; }) != i.end())
      throw std::invalid_argument(
          "QuantumNumberChange<N>(initializer_list<initializer_list> i): each "
          "element of i must contain 2 elements");
#endif
    for (std::size_t c = 0; c != size(); ++c) {
      this->operator[](c) = interval_t{*((i.begin() + c)->begin()),
                                       *((i.begin() + c)->begin() + 1)};
    }
  }

  QuantumNumberChange& operator+=(const QuantumNumberChange& other) {
    for (std::size_t c = 0; c != size(); ++c) this->operator[](c) += other[c];
    return *this;
  }

  bool operator==(const this_type& b) const {
    return std::equal(
        this->begin(), this->end(), b.begin(),
        [](const auto& ia, const auto& ib) { return equal(ia, ib); });
  }
  bool operator!=(const this_type& b) const { return !this->operator==(b); }

  // determines the number of physical vacuum creators and annihilators for the
  // active particle and hole space from the Context. for general operators this
  // is not defined. for example: O_{e_1}^{i_1 m_1} a_{i_1 m_1}^{e_1} asking for
  // the active particle annihilators in this example is nonsense and will
  // return -1.

  interval_t ncre_particles() {
    const auto& qnvec = this->base();
    auto isr = get_default_context().index_space_registry();
    const auto& base_spaces = isr->base_spaces();
    interval_t result = 0;
    for (unsigned int i = 0; i < base_spaces.size(); i++) {
      const auto& base_space = base_spaces[i];
      const auto intersect_type =
          base_space.attr()
              .intersection(isr->particle_space(base_space.qns()).attr())
              .type();
      if (IndexSpace::Type{} != intersect_type) {
        result += qnvec[2 * i];
      }
    }
    return result;
  }

  interval_t nann_particles() {
    const auto& qnvec = this->base();
    auto isr = get_default_context().index_space_registry();
    const auto& base_spaces = isr->base_spaces();
    interval_t result = 0;
    for (unsigned int i = 0; i < base_spaces.size(); i++) {
      const auto& base_space = base_spaces[i];
      const auto intersect_type =
          base_space.attr()
              .intersection(isr->particle_space(base_space.qns()).attr())
              .type();
      if (IndexSpace::Type{} != intersect_type) {
        result += qnvec[2 * i + 1];
      }
    }
    return result;
  }

  interval_t ncre_holes() {
    const auto& qnvec = this->base();
    auto isr = get_default_context().index_space_registry();
    const auto& base_spaces = isr->base_spaces();
    interval_t result = 0;
    for (unsigned int i = 0; i < base_spaces.size(); i++) {
      const auto& base_space = base_spaces[i];
      const auto intersect_type =
          base_space.attr()
              .intersection(isr->hole_space(base_space.qns()).attr())
              .type();
      if (IndexSpace::Type{} != intersect_type) {
        result += qnvec[2 * i];
      }
    }
    return result;
  }

  interval_t nann_holes() {
    const auto& qnvec = this->base();
    auto isr = get_default_context().index_space_registry();
    const auto& base_spaces = isr->base_spaces();
    interval_t result = 0;
    for (unsigned int i = 0; i < base_spaces.size(); i++) {
      const auto& base_space = base_spaces[i];
      const auto intersect_type =
          base_space.attr()
              .intersection(isr->hole_space(base_space.qns()).attr())
              .type();
      if (IndexSpace::Type{} != intersect_type) {
        result += qnvec[2 * i + 1];
      }
    }
    return result;
  }

  template <typename I>
  bool in(I i) {
    return boost::numeric::in(static_cast<int64_t>(i), this->front());
  }

  template <typename I, typename = std::enable_if_t<std::is_integral_v<I>>>
  bool in(std::initializer_list<I> i) {
    assert(i.size() == size());
    std::array<I, 4> i_arr;
    std::copy(i.begin(), i.end(), i_arr.begin());
    return this->in(i_arr);
  }

  bool overlaps_with(base_type i) {
    for (std::size_t c = 0; c != this->size(); ++c) {
      if (!boost::numeric::overlap(i[c], this->operator[](c))) {
        return false;
      }
    }
    return true;
  }

  auto hash_value() const {
    assert(size() > 0);
    auto val = sequant::hash::value(this->operator[](0));
    for (std::size_t c = 1; c != size(); ++c) {
      sequant::hash::combine(val, sequant::hash::value(this->operator[](c)));
    }
    return val;
  }

 private:
  auto& base() { return static_cast<base_type&>(*this); }
};

template <std::size_t N, typename Tag, typename QNV>
inline bool operator==(const QuantumNumberChange<Tag, QNV>& a,
                       const QuantumNumberChange<Tag, QNV>& b) {
  return a.operator==(b);
}

template <std::size_t N, typename Tag, typename QNV>
inline bool operator!=(const QuantumNumberChange<Tag, QNV>& a,
                       const QuantumNumberChange<Tag, QNV>& b) {
  return !(a == b);
}

template <std::size_t N, typename Tag, typename QNV, typename I,
          typename = std::enable_if_t<N == 1 && std::is_integral_v<I>>>
inline bool operator==(const QuantumNumberChange<Tag, QNV>& a, I b) {
  return a.operator==(b);
}

template <std::size_t N, typename Tag, typename QNV>
inline bool equal(const QuantumNumberChange<Tag, QNV>& a,
                  const QuantumNumberChange<Tag, QNV>& b) {
  return operator==(a, b);
}

template <std::size_t N, typename Tag, typename QNV, typename I,
          typename = std::enable_if_t<N == 1 && std::is_integral_v<I>>>
inline bool operator!=(const QuantumNumberChange<Tag, QNV>& a, I b) {
  return a.operator!=(b);
}

// clang-format off
// clang-format on
using qns_t = mbpt::QuantumNumberChange<>;
using qninterval_t = typename qns_t::interval_t;
using qnc_t = qns_t;
using op_t = mbpt::Operator<qnc_t>;

qns_t combine(qns_t, qns_t);

qns_t excitation_type_qns(std::size_t k,
                          IndexSpace::QuantumNumbers SQN = Spin::any);

qns_t interval_excitation_type_qns(std::size_t k,
                                   IndexSpace::QuantumNumbers SQN = Spin::any);

qns_t deexcitation_type_qns(std::size_t k,
                            IndexSpace::QuantumNumbers SQN = Spin::any);

qns_t interval_deexcitation_type_qns(
    std::size_t k, IndexSpace::QuantumNumbers SQN = Spin::any);

qns_t general_type_qns(std::size_t k);

qns_t generic_excitation_qns(std::size_t particle_rank, std::size_t hole_rank,
                             IndexSpace particle_space, IndexSpace hole_space,
                             IndexSpace::QuantumNumbers SQN = Spin::any);

qns_t generic_deexcitation_qns(std::size_t particle_rank, std::size_t hole_rank,
                               IndexSpace particle_space, IndexSpace hole_space,
                               IndexSpace::QuantumNumbers SQN = Spin::any);

inline namespace op {
namespace tensor {
ExprPtr vac_av(ExprPtr expr,
               std::vector<std::pair<int, int>> nop_connections = {},
               bool use_top = true);
}  // namespace tensor
}  // namespace op

}  // namespace mbpt

mbpt::qns_t adjoint(mbpt::qns_t qns);

namespace mbpt {

// clang-format off

// clang-format on
template <Statistics S>
class OpMaker {
 public:
  template <typename IndexSpaceTypeRange1, typename IndexSpaceTypeRange2>
  OpMaker(OpType op, const cre<IndexSpaceTypeRange1>& cre_list,
          const ann<IndexSpaceTypeRange2>& ann_list)
      : op_(op),
        cre_spaces_(cre_list.begin(), cre_list.end()),
        ann_spaces_(ann_list.begin(), ann_list.end()) {
    assert(ncreators() > 0 || nannihilators() > 0);
  }

  OpMaker(OpType op, ncre nc, nann na);

  OpMaker(OpType op, std::size_t rank);

  OpMaker(OpType op, ncre nc, nann na, const cre<IndexSpace>& cre_space,
          const ann<IndexSpace>& ann_space);

  enum class UseDepIdx {
    Bra,
    Ket,
    None
  };

  // clang-format off
  // clang-format on
  ExprPtr operator()(std::optional<UseDepIdx> dep_opt = {},
                     std::optional<Symmetry> opsymm_opt = {}) const;

  template <typename TensorGenerator>
  static ExprPtr make(const container::svector<IndexSpace>& creators,
                      const container::svector<IndexSpace>& annihilators,
                      TensorGenerator&& tensor_generator,
                      UseDepIdx dep = UseDepIdx::None) {
    const bool symm =
        get_default_context().spbasis() ==
        SPBasis::spinorbital;  // antisymmetrize if spin-orbital basis
    const auto dep_bra = dep == UseDepIdx::Bra;
    const auto dep_ket = dep == UseDepIdx::Ket;

    // not sure what it means to use nonsymmetric operator if nbra != nket
    if (!symm) assert(ranges::size(creators) == ranges::size(annihilators));

    auto make_idx_vector = [](const auto& spacetypes) {
      container::svector<Index> result;
      const auto n = spacetypes.size();
      result.reserve(n);
      for (size_t i = 0; i != n; ++i) {
        auto space = spacetypes[i];
        result.push_back(Index::make_tmp_index(space));
      }
      return result;
    };

    auto make_depidx_vector = [](const auto& spacetypes, auto&& protoidxs) {
      const auto n = spacetypes.size();
      container::svector<Index> result;
      result.reserve(n);
      for (size_t i = 0; i != n; ++i) {
        auto space = spacetypes[i];
        result.push_back(Index::make_tmp_index(space, protoidxs, true));
      }
      return result;
    };

    container::svector<Index> creidxs, annidxs;
    if (dep_bra) {
      annidxs = make_idx_vector(annihilators);
      creidxs = make_depidx_vector(creators, annidxs);
    } else if (dep_ket) {
      creidxs = make_idx_vector(creators);
      annidxs = make_depidx_vector(annihilators, creidxs);
    } else {
      creidxs = make_idx_vector(creators);
      annidxs = make_idx_vector(annihilators);
    }

    using ranges::size;
    const auto mult =
        symm ? factorial(size(creators)) * factorial(size(annihilators))
             : factorial(size(creators));
    const auto opsymm = symm ? (S == Statistics::FermiDirac ? Symmetry::antisymm
                                                            : Symmetry::symm)
                             : Symmetry::nonsymm;
    return ex<Constant>(rational{1, mult}) *
           tensor_generator(creidxs, annidxs, opsymm) *
           ex<NormalOperator<S>>(cre(creidxs), ann(annidxs),
                                 get_default_context().vacuum());
  }

  template <typename TensorGenerator>
  static ExprPtr make(std::initializer_list<IndexSpace::Type> creators,
                      std::initializer_list<IndexSpace::Type> annihilators,
                      TensorGenerator&& tensor_generator,
                      UseDepIdx csv = UseDepIdx::None) {
    container::svector<IndexSpace> cre_vec(creators.begin(), creators.end());
    container::svector<IndexSpace> ann_vec(annihilators.begin(),
                                           annihilators.end());
    return OpMaker::make(cre_vec, ann_vec,
                         std::forward<TensorGenerator>(tensor_generator), csv);
  }

 protected:
  OpType op_;
  container::svector<IndexSpace> cre_spaces_;
  container::svector<IndexSpace> ann_spaces_;

  OpMaker(OpType op);

  [[nodiscard]] const auto ncreators() const { return cre_spaces_.size(); };
  [[nodiscard]] const auto nannihilators() const { return ann_spaces_.size(); };
};

extern template class OpMaker<Statistics::FermiDirac>;
extern template class OpMaker<Statistics::BoseEinstein>;

template <typename QuantumNumbers, Statistics S>
class Operator : public Operator<void, S> {
  using this_type = Operator<QuantumNumbers, S>;
  using base_type = Operator<void, S>;

 protected:
  Operator();

 public:
  Operator(std::function<std::wstring_view()> label_generator,
           std::function<ExprPtr()> tensor_form_generator,
           std::function<void(QuantumNumbers&)> qn_action);

  virtual ~Operator();

  QuantumNumbers operator()(const QuantumNumbers& qns = {}) const;

  template <typename I, typename = std::enable_if_t<std::is_integral_v<I>>>
  QuantumNumbers operator()(std::initializer_list<I> qns) const {
    QuantumNumbers result(qns);
    this->apply_to(result);
    return result;
  }

  virtual QuantumNumbers& apply_to(QuantumNumbers& qns) const;

  bool static_less_than(const Expr& that) const override;

  bool commutes_with_atom(const Expr& that) const override;

  void adjoint() override;

  bool is_adjoint_ = false;

 private:
  std::function<void(QuantumNumbers&)> qn_action_;

  bool less_than_rank_of(const this_type& that) const;

  Expr::type_id_type type_id() const override;

  ExprPtr clone() const override;

  std::wstring to_latex() const override;

  Expr::hash_type memoizing_hash() const override;

};  // class Operator

extern template class Operator<qns_t, Statistics::FermiDirac>;
extern template class Operator<qns_t, Statistics::BoseEinstein>;

inline namespace op {
namespace tensor {

using mbpt::nₕ;
using mbpt::nₚ;

// clang-format off
// clang-format on
ExprPtr H_(std::size_t k);

ExprPtr H(std::size_t k = 2);

ExprPtr F(bool use_tensor = true, IndexSpace reference_occupied = {L"", 0});

ExprPtr θ(std::size_t K);

ExprPtr T_(std::size_t K);

ExprPtr T(std::size_t K, bool skip1 = false);

ExprPtr Λ_(std::size_t K);

ExprPtr Λ(std::size_t K);

ExprPtr R_(
    nann na, ncre nc,
    const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
    const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));

ExprPtr R_(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(R_);

ExprPtr L_(
    nann na, ncre nc,
    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));

ExprPtr L_(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(L_);

// clang-format off
// clang-format on
ExprPtr P(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(P);

// clang-format off
// clang-format on
ExprPtr A(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(A);

ExprPtr S(std::int64_t K);

ExprPtr H_pt(std::size_t order, std::size_t R);

ExprPtr T_pt_(std::size_t order, std::size_t K);

ExprPtr T_pt(std::size_t order, std::size_t K, bool skip1 = false);

ExprPtr Λ_pt_(std::size_t order, std::size_t K);

ExprPtr Λ_pt(std::size_t order, std::size_t K, bool skip1 = false);
}  // namespace tensor
}  // namespace op

inline namespace op {
// clang-format off
// clang-format on
ExprPtr H_(std::size_t k);

ExprPtr H(std::size_t k = 2);

ExprPtr F(bool use_tensor = true, IndexSpace reference_occupied = {L"", 0});

ExprPtr θ(std::size_t K);

ExprPtr T_(std::size_t K);

ExprPtr T(std::size_t K, bool skip1 = false);

ExprPtr Λ_(std::size_t K);

ExprPtr Λ(std::size_t K);

ExprPtr R_(
    nann na, ncre nc,
    const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
    const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));

ExprPtr R_(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(R_);

ExprPtr L_(
    nann na, ncre nc,
    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));

ExprPtr L_(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(L_);

ExprPtr R(nann na, ncre nc,
          const cre<IndexSpace>& cre_space = cre(get_particle_space(Spin::any)),
          const ann<IndexSpace>& ann_space = ann(get_hole_space(Spin::any)));

ExprPtr R(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(R);

ExprPtr L(
    nann na, ncre nc,
    const cre<IndexSpace>& cre_space = cre(get_hole_space(Spin::any)),
    const ann<IndexSpace>& ann_space = ann(get_particle_space(Spin::any)));

ExprPtr L(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(L);

// clang-format off
// clang-format on
ExprPtr P(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(P);

// clang-format off
// clang-format on
ExprPtr A(nₚ np, nₕ nh);
DEFINE_SINGLE_SIGNED_ARGUMENT_OP_VARIANT(A);

ExprPtr S(std::int64_t K);

ExprPtr H_pt(std::size_t order, std::size_t R);

ExprPtr T_pt_(std::size_t order, std::size_t K);

ExprPtr T_pt(std::size_t order, std::size_t K, bool skip1 = false);

ExprPtr Λ_pt_(std::size_t order, std::size_t K);

ExprPtr Λ_pt(std::size_t order, std::size_t K, bool skip1 = false);

bool raises_vacuum_up_to_rank(const ExprPtr& op_or_op_product,
                              const unsigned long k);

bool lowers_rank_or_lower_to_vacuum(const ExprPtr& op_or_op_product,
                                    const unsigned long k);

bool raises_vacuum_to_rank(const ExprPtr& op_or_op_product,
                           const unsigned long k);

bool lowers_rank_to_vacuum(const ExprPtr& op_or_op_product,
                           const unsigned long k);
#include <SeQuant/domain/mbpt/vac_av.hpp>
}  // namespace op
}  // namespace mbpt
}  // namespace sequant

#endif  // SEQUANT_DOMAIN_MBPT_OP_HPP