.. _program_listing_file_SeQuant_domain_mbpt_spin.hpp:

Program Listing for File spin.hpp
=================================

|exhale_lsh| :ref:`Return to documentation for file <file_SeQuant_domain_mbpt_spin.hpp>` (``SeQuant/domain/mbpt/spin.hpp``)

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

.. code-block:: cpp

   //
   // 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