Program Listing for File op.ipp

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

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

#ifndef SEQUANT_DOMAIN_MBPT_OP_IPP
#define SEQUANT_DOMAIN_MBPT_OP_IPP

#include <SeQuant/core/context.hpp>
#include <SeQuant/core/tensor_canonicalizer.hpp>
#include <SeQuant/domain/mbpt/op.hpp>

#include <range/v3/algorithm/all_of.hpp>
#include <range/v3/algorithm/any_of.hpp>
#include <range/v3/algorithm/find_if.hpp>
#include <range/v3/algorithm/is_sorted.hpp>
#include <range/v3/view/iota.hpp>

namespace sequant {
namespace mbpt {

template <typename QuantumNumbers, Statistics S>
Operator<QuantumNumbers, S>::Operator() = default;

template <typename QuantumNumbers, Statistics S>
Operator<QuantumNumbers, S>::Operator(
    std::function<std::wstring_view()> label_generator,
    std::function<ExprPtr()> tensor_form_generator,
    std::function<void(QuantumNumbers&)> qn_action)
    : base_type(std::move(label_generator), std::move(tensor_form_generator)),
      qn_action_(std::move(qn_action)) {}

template <typename QuantumNumbers, Statistics S>
Operator<QuantumNumbers, S>::Operator(
    std::function<std::wstring_view()> label_generator,
    std::function<ExprPtr()> tensor_form_generator,
    std::function<void(QuantumNumbers&)> qn_action, const OpParams& params)
    : Operator(std::move(label_generator), std::move(tensor_form_generator),
               std::move(qn_action)) {
  params.validate();
  order_ = params.order;
  // decorate label with perturbation order if non-zero
  if (order_ != 0) {
    auto decorated = detail::decorate_with_pert_order(this->label_generator_(),
                                                      static_cast<int>(order_));
    this->label_generator_ =
        [label = std::move(decorated)]() -> std::wstring_view { return label; };
  }
  // 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<container::svector<std::size_t>>();
  }
}

template <typename QuantumNumbers, Statistics S>
Operator<QuantumNumbers, S>::~Operator() = default;

template <typename QuantumNumbers, Statistics S>
QuantumNumbers Operator<QuantumNumbers, S>::operator()(
    const QuantumNumbers& qns) const {
  QuantumNumbers result(qns);
  this->apply_to(result);
  return result;
}

template <typename QuantumNumbers, Statistics S>
QuantumNumbers& Operator<QuantumNumbers, S>::apply_to(
    QuantumNumbers& qns) const {
  SEQUANT_ASSERT(qn_action_);
  if (is_vacuum(qns)) {  // action on vacuum is trivial ...
    qn_action_(qns);
  } else {  // action on a {operator. product of operators} = use Wick's theorem
    QuantumNumbers qns_this;
    qn_action_(qns_this);
    qns = combine(qns_this, qns);
  }
  return qns;
}

template <typename QuantumNumbers, Statistics S>
bool Operator<QuantumNumbers, S>::static_less_than(const Expr& that) const {
  SEQUANT_ASSERT(that.is<this_type>());
  auto& that_op = that.as<this_type>();

  // compare cardinal tensor labels first, then QN ranks
  const auto ctx = get_default_context_snapshot();
  const auto& cardinal_tensor_labels = ctx.cardinal_tensor_labels();
  const auto this_label = this->label();
  const auto that_label = that_op.label();
  if (this_label == that_label) return this->less_than_rank_of(that_op);
  const auto this_is_cardinal_it = ranges::find_if(
      cardinal_tensor_labels,
      [&this_label](const std::wstring& l) { return l == this_label; });
  const auto this_is_cardinal =
      this_is_cardinal_it != ranges::end(cardinal_tensor_labels);
  const auto that_is_cardinal_it = ranges::find_if(
      cardinal_tensor_labels,
      [&that_label](const std::wstring& l) { return l == that_label; });
  const auto that_is_cardinal =
      that_is_cardinal_it != ranges::end(cardinal_tensor_labels);
  if (this_is_cardinal && that_is_cardinal) {
    if (this_is_cardinal_it != that_is_cardinal_it)
      return this_is_cardinal_it < that_is_cardinal_it;
    else
      return this->less_than_rank_of(that_op);
  } else if (this_is_cardinal && !that_is_cardinal)
    return true;
  else if (!this_is_cardinal && that_is_cardinal)
    return false;
  else {  // !this_is_cardinal && !that_is_cardinal
    if (this_label == that_label)
      return this->less_than_rank_of(that_op);
    else
      return this_label < that_label;
  }
}

template <typename QuantumNumbers, Statistics S>
bool Operator<QuantumNumbers, S>::commutes_with_atom(const Expr& that) const {
  SEQUANT_ASSERT(that.is_cnumber() || that.is<this_type>());
  if (that.is_cnumber())
    return true;
  else {
    auto& that_op = that.as<this_type>();

    // if this has annihilators/creators in same space as that has
    // creator/annihilators return false

    const auto delta_this = (*this)();
    const auto delta_that = (that_op)();

    SEQUANT_ASSERT(delta_this.size() % 2 == 0 &&
                   delta_that.size() == delta_this.size());

    const auto nspaces = delta_this.size() / 2;
    for (std::size_t s = 0; s != nspaces; ++s) {
      const auto cre = QuantumNumbers::cre_slot(s);
      const auto ann = QuantumNumbers::ann_slot(s);
      if ((delta_this[cre].upper() > 0 && delta_that[ann].upper() > 0) ||
          (delta_this[ann].upper() > 0 && delta_that[cre].upper() > 0))
        return false;
    }
    // without a contractible pair two fermionic operators of odd rank
    // anticommute
    if constexpr (S == Statistics::FermiDirac)
      return !(this->has_odd_rank_term() && that_op.has_odd_rank_term());
    return true;
  }
}

template <typename QuantumNumbers, Statistics S>
bool Operator<QuantumNumbers, S>::has_odd_rank_term() const {
  auto compute = [this]() -> bool {
    // exact op counts give the rank directly; intervals (ops on non-base
    // spaces) do not, then the tensor form decides
    const auto dN = (*this)();
    if (ranges::all_of(dN, [](const auto& n) { return singleton(n); })) {
      std::int64_t rank = 0;
      for (const auto& n : dN) rank += n.lower();
      return rank % 2 != 0;
    }

    auto term_rank = [](const ExprPtr& term) {
      std::size_t rank = 0;
      auto count = [&rank](const ExprPtr& e) {
        if (e.is<NormalOperator<S>>()) {
          const auto& nop = e.as<NormalOperator<S>>();
          rank += nop.ncreators() + nop.nannihilators();
        }
      };
      term->visit(count, /* atoms_only = */ true);
      return rank;
    };
    const ExprPtr form = this->tensor_form();
    if (!form.is<Sum>()) return term_rank(form) % 2 != 0;
    return ranges::any_of(form.as<Sum>().summands(), [&](const ExprPtr& t) {
      return term_rank(t) % 2 != 0;
    });
  };
  if (!has_odd_rank_term_) has_odd_rank_term_ = compute();
  return *has_odd_rank_term_;
}

template <typename QuantumNumbers, Statistics S>
void Operator<QuantumNumbers, S>::adjoint() {
  const auto dN = (*this)(QuantumNumbers{});
  using qnc_t = std::decay_t<decltype(dN)>;
  static_assert(std::is_same_v<QuantumNumbers, qnc_t>,
                "mbpt::Operator::adjoint: QuantumNumbers type mismatch");

  // grab label and update according to adjoint flag
  auto lbl = std::wstring(this->label());
  SEQUANT_ASSERT(sequant::is_adjoint_label(lbl) == is_adjoint_);
  sequant::toggle_adjoint_label(lbl);

  // get order and batch_ordinals to restore later
  const auto saved_order = this->order_;
  const auto saved_batch_ordinals = this->batch_ordinals_;

  // capture the current tensor-form generator so that each call regenerates
  // fresh dummy indices; otherwise using this adjoint operator more than once
  // would yield tensors that have the same indices.
  auto base_tensor_form = this->tensor_form_generator_;
  const auto adjoint_dN = sequant::adjoint(dN);
  *this = Operator{
      [=]() -> std::wstring_view { return lbl; },  // label_generator
      [=]() -> ExprPtr {
        return sequant::adjoint(base_tensor_form());  // tensor_form_generator
      },
      [=](qnc_t& qn) {
        qn += adjoint_dN;
        return qn;  // qn_action
      }};
  this->is_adjoint_ = !this->is_adjoint_;  // toggle adjoint flag

  // restore original order and batch_ordinals
  this->order_ = saved_order;
  this->batch_ordinals_ = saved_batch_ordinals;
}

template <typename QuantumNumbers, Statistics S>
bool Operator<QuantumNumbers, S>::less_than_rank_of(
    const this_type& that) const {
  return (*this)(QuantumNumbers{}) < that(QuantumNumbers{});
}

template <typename QuantumNumbers, Statistics S>
Expr::type_id_type Operator<QuantumNumbers, S>::type_id() const {
  return Expr::get_type_id<this_type>();
};

template <typename QuantumNumbers, Statistics S>
ExprPtr Operator<QuantumNumbers, S>::clone() const {
  return ex<this_type>(*this);
}

// Expresses general operators in human interpretable form. for example:
// \hat{T}_2 is a particle conserving 2-body excitation operator a non-particle
// conserving operator \hat{R}_2_1 implies that two particles are created
// followed by a single hole creation. conversely \hat{R}_1_2 implies the that
// only one particle is annihilated followed by two holes being created. The
// rule being, that for non-particle conserving operators, the first position
// indicates where the quasiparticle is going to and the second position
// indicates where it comes from. for the case of adjoint operators, the adjoint
// is represented by the symbol ⁺ and superscripting the quasi-particle numbers.
// for example: hat{R⁺}^{1,2}} For operators in which one or more
// quasi-particles has only partial coverage in the particle_space or
// hole_space, this notation is unsuitable, and we default to level printing of
// the operator.
template <typename QuantumNumbers, Statistics S>
std::wstring Operator<QuantumNumbers, S>::to_latex() const {
  return sequant::to_latex(*this);
}

template <typename QuantumNumbers, Statistics S>
Expr::hash_type Operator<QuantumNumbers, S>::memoizing_hash() const {
  auto compute_hash = [this]() {
    using std::begin;
    using std::end;
    auto qns = (*this)(QuantumNumbers{});
    auto val = sequant::hash::value(qns);
    sequant::hash::combine(val, std::wstring(this->label()));
    if (batch_ordinals()) {
      const auto ordinals = batch_ordinals().value();
      sequant::hash::range(val, begin(ordinals), end(ordinals));
    }
    return val;
  };
  if (!this->hash_value_) {
    this->hash_value_ = compute_hash();
  } else {
    SEQUANT_ASSERT(*(this->hash_value_) == compute_hash());
  }
  return *(this->hash_value_);
}

template <typename QuantumNumbers, Statistics S>
bool Operator<QuantumNumbers, S>::static_equal(const Expr& other_expr) const {
  const auto& other =
      static_cast<const Operator<QuantumNumbers, S>&>(other_expr);
  return this->label() == other.label() &&
         (*this)(QuantumNumbers{}) == other(QuantumNumbers{}) &&
         this->batch_ordinals() == other.batch_ordinals() &&
         this->order() == other.order();
}

}  // namespace mbpt
}  // namespace sequant

#endif  // SEQUANT_DOMAIN_MBPT_OP_IPP