Program Listing for File spin.hpp

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

//
// Created by Eduard Valeyev on 2019-02-27.
//

#ifndef SEQUANT_DOMAIN_MBPT_SPIN_HPP
#define SEQUANT_DOMAIN_MBPT_SPIN_HPP

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

#include <SeQuant/core/container.hpp>
#include <SeQuant/core/expr.hpp>
#include <SeQuant/core/index.hpp>
#include <SeQuant/core/tensor.hpp>

#include <cassert>
#include <cstddef>
#include <initializer_list>
#include <string>
#include <vector>

namespace sequant {

namespace mbpt {

enum class Spin : bitset_t {
  alpha = 0b000001,
  beta = 0b000010,
  any = 0b000011,  // both bits set so that overlap and union work as expected
                   // (any & alpha = alpha, alpha | beta = any)
  free = any,      // syntax sugar
  spinmask = 0b000011  // coincides with any
};

// Spin is a scoped enum, hence not implicitly convertible to
// QuantumNumbersAttr::bitset_t
static_assert(!std::is_convertible_v<sequant::mbpt::Spin, bitset_t>);
// QuantumNumbersAttr::bitset_t cannot be constructed from Spin
static_assert(!std::is_constructible_v<bitset_t, sequant::mbpt::Spin>);
// but Spin can be cast to QuantumNumbersAttr::bitset_t
static_assert(meta::is_statically_castable_v<sequant::mbpt::Spin, bitset_t>);
// Spin cannot be cast to nonsense ...
static_assert(
    !meta::is_statically_castable_v<sequant::mbpt::Spin, std::string>);

inline Spin operator~(Spin s) {
  return static_cast<Spin>(~(static_cast<bitset_t>(s)));
}
inline Spin operator|(Spin s1, Spin s2) {
  return static_cast<Spin>(static_cast<bitset_t>(s1) |
                           static_cast<bitset_t>(s2));
}
inline Spin operator&(Spin s1, Spin s2) {
  return static_cast<Spin>(static_cast<bitset_t>(s1) &
                           static_cast<bitset_t>(s2));
}

inline Spin to_spin(const QuantumNumbersAttr& t) {
  assert((t.to_int32() & static_cast<int>(Spin::spinmask)) != 0);
  return static_cast<Spin>(static_cast<Spin>(t.to_int32()) & Spin::spinmask);
}

inline QuantumNumbersAttr spinannotation_remove(const QuantumNumbersAttr& t) {
  return t.intersection(QuantumNumbersAttr(~Spin::spinmask));
}

template <typename WS, typename = std::enable_if_t<
                           meta::is_wstring_or_view_v<std::decay_t<WS>>>>
std::wstring spinannotation_remove(WS&& label) {
  auto [base, orbital] = Index::make_split_label(label);
  if (base.back() == L'↑' || base.back() == L'↓') {
    base.remove_suffix(1);
  }
  return Index::make_merged_label(base, orbital);
}

template <typename WS, typename = std::enable_if_t<
                           meta::is_wstring_or_view_v<std::decay_t<WS>>>>
std::wstring spinannotation_add(WS&& label, Spin s) {
  auto [base, ordinal] = Index::make_split_label(label);
  switch (s) {
    case Spin::any:
      return std::wstring(label);
    case Spin::alpha:
      return Index::make_merged_label(std::wstring(base) + L"↑", ordinal);
    case Spin::beta:
      return Index::make_merged_label(std::wstring(base) + L"↓", ordinal);
    default:
      assert(false && "invalid quantum number");
      abort();
  }
}

template <typename WS, typename = std::enable_if_t<
                           meta::is_wstring_or_view_v<std::decay_t<WS>>>>
std::wstring spinannotation_replacе(WS&& label, Spin s) {
  auto label_sf = spinannotation_remove(std::forward<WS>(label));
  return spinannotation_add(label_sf, s);
}

}  // namespace mbpt

// make alpha-spin idx
[[nodiscard]] Index make_spinalpha(const Index& idx);

// make beta-spin idx
[[nodiscard]] Index make_spinbeta(const Index& idx);

// make null-spin idx
[[nodiscard]] Index make_spinfree(const Index& idx);

ExprPtr transform_expr(const ExprPtr& expr,
                       const container::map<Index, Index>& index_replacements,
                       Constant::scalar_type scaling_factor = 1);

ExprPtr swap_bra_ket(const ExprPtr& expr);

ExprPtr append_spin(const ExprPtr& expr,
                    const container::map<Index, Index>& index_replacements);

ExprPtr remove_spin(const ExprPtr& expr);

bool spin_symm_tensor(const Tensor& tensor);

bool same_spin_tensor(const Tensor& tensor);

bool can_expand(const Tensor& tensor);

ExprPtr expand_antisymm(const Tensor& tensor, bool skip_spinsymm = false);

ExprPtr expand_antisymm(const ExprPtr& expr, bool skip_spinsymm = false);

bool has_tensor(const ExprPtr& expr, std::wstring label);

container::svector<container::map<Index, Index>> A_maps(const Tensor& A);

ExprPtr remove_tensor(const Product& product, std::wstring label);

ExprPtr remove_tensor(const ExprPtr& expr, std::wstring label);

ExprPtr expand_A_op(const Product& product);

ExprPtr symmetrize_expr(const Product& product);

ExprPtr symmetrize_expr(const ExprPtr& expr);

ExprPtr expand_A_op(const ExprPtr& expr);

container::svector<container::map<Index, Index>> P_maps(
    const Tensor& P, bool keep_canonical = true, bool pair_wise = false);

ExprPtr expand_P_op(const Product& product, bool keep_canonical = true,
                    bool pair_wise = true);

ExprPtr expand_P_op(const ExprPtr& expr, bool keep_canonical = true,
                    bool pair_wise = true);

container::svector<container::map<Index, Index>> S_replacement_maps(
    const Tensor& S);

ExprPtr S_maps(const ExprPtr& expr);


template <typename Seq0, typename Seq1>
std::size_t count_cycles(Seq0&& v0, const Seq1& v1) {
  std::remove_reference_t<Seq0> v(std::forward<Seq0>(v0));
  using T = std::decay_t<decltype(v[0])>;
  assert(ranges::is_permutation(v, v1));

  auto make_null = []() -> T {
    if constexpr (std::is_arithmetic_v<T>) {
      return -1;
    } else if constexpr (std::is_same_v<T, Index>) {
      return L"p_50";
    } else  // unreachable
      abort();
  };

  std::size_t n_cycles = 0;
  const auto null = make_null();
  assert(ranges::contains(v, null) == false);
  assert(ranges::contains(v1, null) == false);
  for (auto it = v.begin(); it != v.end(); ++it) {
    if (*it != null) {
      n_cycles++;
      auto idx = std::distance(v.begin(), it);
      auto it0 = it;
      auto it1 = std::find(v1.begin(), v1.end(), *it0);
      auto idx1 = std::distance(v1.begin(), it1);
      do {
        it0 = std::find(v.begin(), v.end(), v[idx1]);
        it1 = std::find(v1.begin(), v1.end(), *it0);
        idx1 = std::distance(v1.begin(), it1);
        *it0 = null;
      } while (idx1 != idx);
    }
  }
  return n_cycles;
};

ExprPtr closed_shell_spintrace(
    const ExprPtr& expr,
    const container::svector<container::svector<Index>>& ext_index_groups = {});

container::svector<container::svector<Index>> external_indices(Tensor const&);

ExprPtr closed_shell_CC_spintrace(ExprPtr const& expr);

ExprPtr closed_shell_CC_spintrace_rigorous(ExprPtr const& expr);

auto index_list(const ExprPtr& expr);

Tensor swap_spin(const Tensor& t);

ExprPtr swap_spin(const ExprPtr& expr);

ExprPtr merge_tensors(const Tensor& O1, const Tensor& O2);

std::vector<ExprPtr> open_shell_A_op(const Tensor& A);

std::vector<ExprPtr> open_shell_P_op_vector(const Tensor& A);

std::vector<ExprPtr> open_shell_spintrace(
    const ExprPtr& expr,
    const container::svector<container::svector<Index>>& ext_index_groups,
    const int single_spin_case = 0);

std::vector<ExprPtr> open_shell_CC_spintrace(const ExprPtr& expr);

ExprPtr spintrace(
    const ExprPtr& expr,
    container::svector<container::svector<Index>> ext_index_groups = {},
    bool spinfree_index_spaces = true);

ExprPtr factorize_S(const ExprPtr& expression,
                    std::initializer_list<IndexList> ext_index_groups,
                    bool fast_method = true);

ExprPtr biorthogonal_transform(
    const sequant::ExprPtr& expr,
    const container::svector<container::svector<sequant::Index>>&
        ext_index_groups = {},
    double threshold = 1.e-12);

}  // namespace sequant

#endif  // SEQUANT_DOMAIN_MBPT_SPIN_HPP