Azam/feature/hardcoded biorthogonal coefficients

CMakeLists.txt

@@ -453,6 +453,7 @@ if (Boost_IS_MODULARIZED)
 endif()
 if (TARGET tiledarray)
     target_link_libraries(SeQuant PUBLIC tiledarray)
+    target_compile_definitions(SeQuant PUBLIC SEQUANT_HAS_TILEDARRAY=1 SEQUANT_HAS_BTAS=1)
 endif ()
 if (SEQUANT_HAS_EXECUTION_HEADER_STANDALONE OR SEQUANT_HAS_EXECUTION_HEADER_WITH_TBB)
     target_compile_definitions(SeQuant PUBLIC SEQUANT_HAS_EXECUTION_HEADER)

SeQuant/core/eval/eval.hpp

@@ -1,9 +1,10 @@
 #ifndef SEQUANT_EVAL_EVAL_HPP
 #define SEQUANT_EVAL_EVAL_HPP
 
+#include <SeQuant/core/eval/eval_fwd.hpp>
+
 #include <SeQuant/core/container.hpp>
 #include <SeQuant/core/eval/cache_manager.hpp>
-#include <SeQuant/core/eval/eval_fwd.hpp>
 #include <SeQuant/core/eval/result.hpp>
 #include <SeQuant/core/eval_node.hpp>
 #include <SeQuant/core/expr.hpp>
@@ -60,9 +61,6 @@ enum struct EvalMode {
   SumInplace,
   Symmetrize,
   Antisymmetrize,
-  /// NNS projection of Wang-Knizia biorthogonalization
-  /// @sa ResultPtr::biorthogonal_nns_project
-  BiorthogonalNNSProject,
   Unknown
 };
 
@@ -80,18 +78,17 @@ enum struct EvalMode {
 }
 
 [[nodiscard]] constexpr auto to_string(EvalMode mode) noexcept {
-  return (mode == EvalMode::Constant)                 ? "Constant"
-         : (mode == EvalMode::Variable)               ? "Variable"
-         : (mode == EvalMode::Tensor)                 ? "Tensor"
-         : (mode == EvalMode::Permute)                ? "Permute"
-         : (mode == EvalMode::Product)                ? "Product"
-         : (mode == EvalMode::MultByPhase)            ? "MultByPhase"
-         : (mode == EvalMode::Sum)                    ? "Sum"
-         : (mode == EvalMode::SumInplace)             ? "SumInplace"
-         : (mode == EvalMode::Symmetrize)             ? "Symmetrize"
-         : (mode == EvalMode::Antisymmetrize)         ? "Antisymmetrize"
-         : (mode == EvalMode::BiorthogonalNNSProject) ? "BiorthogonalNNSProject"
-                                                      : "??";
+  return (mode == EvalMode::Constant)         ? "Constant"
+         : (mode == EvalMode::Variable)       ? "Variable"
+         : (mode == EvalMode::Tensor)         ? "Tensor"
+         : (mode == EvalMode::Permute)        ? "Permute"
+         : (mode == EvalMode::Product)        ? "Product"
+         : (mode == EvalMode::MultByPhase)    ? "MultByPhase"
+         : (mode == EvalMode::Sum)            ? "Sum"
+         : (mode == EvalMode::SumInplace)     ? "SumInplace"
+         : (mode == EvalMode::Symmetrize)     ? "Symmetrize"
+         : (mode == EvalMode::Antisymmetrize) ? "Antisymmetrize"
+                                              : "??";
 }
 
 enum struct CacheMode { Store, Access, Release };
@@ -543,30 +540,6 @@ ResultPtr evaluate_antisymm(Args&&... args) {
   return result;
 }
 
-/// \brief Calls sequant::evaluate followed by
-/// ResultPtr::biorthogonal_nns_project \return Evaluated result as ResultPtr.
-/// \sa ResultPtr::biorthogonal_nns_project
-template <Trace EvalTrace = Trace::Default, typename... Args>
-ResultPtr evaluate_biorthogonal_nns_project(Args&&... args) {
-  ResultPtr pre = evaluate<EvalTrace>(std::forward<Args>(args)...);
-  SEQUANT_ASSERT(pre);
-
-  auto const& n0 = node0(arg0(std::forward<Args>(args)...));
-
-  ResultPtr result;
-  auto time = timed_eval_inplace([&]() {
-    result = pre->biorthogonal_nns_project(n0->as_tensor().bra_rank());
-  });
-
-  // logging
-  if constexpr (trace(EvalTrace)) {
-    auto stat = log::EvalStat{.mode = log::EvalMode::BiorthogonalNNSProject,
-                              .time = time,
-                              .memory = log::bytes(pre, result)};
-    log::eval(stat, n0->label());
-  }
-  return result;
-}
 }  // namespace sequant
 
 #endif  // SEQUANT_EVAL_EVAL_HPP

SeQuant/core/eval/result.hpp

@@ -1,9 +1,10 @@
 #ifndef SEQUANT_EVAL_RESULT_HPP
 #define SEQUANT_EVAL_RESULT_HPP
 
+#include <SeQuant/core/eval/eval_fwd.hpp>
+
 #include <SeQuant/core/algorithm.hpp>
 #include <SeQuant/core/container.hpp>
-#include <SeQuant/core/eval/eval_fwd.hpp>
 #include <SeQuant/core/hash.hpp>
 #include <SeQuant/core/index.hpp>
 #include <SeQuant/core/logger.hpp>
@@ -340,157 +341,6 @@ auto particle_antisymmetrize_btas(btas::Tensor<Args...> const& arr,
   return result;
 }
 
-/// \brief This function is used to implement
-/// ResultPtr::biorthogonal_nns_project for TA::DistArray
-///
-/// \param arr The array to be "cleaned up"
-/// \param bra_rank The rank of the bra indices
-///
-/// \return The cleaned TA::DistArray.
-template <typename... Args>
-auto biorthogonal_nns_project_ta(TA::DistArray<Args...> const& arr,
-                                 size_t bra_rank) {
-  using ranges::views::iota;
-  size_t const rank = arr.trange().rank();
-  SEQUANT_ASSERT(bra_rank <= rank);
-  size_t const ket_rank = rank - bra_rank;
-
-  if (rank <= 4) {
-    return arr;
-  }
-
-  using numeric_type = typename TA::DistArray<Args...>::numeric_type;
-
-  size_t factorial_ket = 1;
-  for (size_t i = 2; i <= ket_rank; ++i) {
-    factorial_ket *= i;
-  }
-  numeric_type norm_factor = numeric_type(1) / numeric_type(factorial_ket);
-
-  TA::DistArray<Args...> result;
-
-  perm_t perm = iota(size_t{0}, rank) | ranges::to<perm_t>;
-  perm_t bra_perm = iota(size_t{0}, bra_rank) | ranges::to<perm_t>;
-  perm_t ket_perm = iota(bra_rank, rank) | ranges::to<perm_t>;
-
-  const auto lannot = ords_to_annot(perm);
-
-  auto process_permutations = [&lannot](const TA::DistArray<Args...>& input_arr,
-                                        size_t range_rank, perm_t range_perm,
-                                        const std::string& other_annot,
-                                        bool is_bra) -> TA::DistArray<Args...> {
-    if (range_rank <= 1) return input_arr;
-    TA::DistArray<Args...> result;
-
-    auto callback = [&]([[maybe_unused]] int parity) {
-      const auto range_annot = ords_to_annot(range_perm);
-      const auto annot = other_annot.empty()
-                             ? range_annot
-                             : (is_bra ? range_annot + "," + other_annot
-                                       : other_annot + "," + range_annot);
-
-      // ignore parity, all permutations get same coefficient
-      numeric_type p_ = 1;
-      if (result.is_initialized()) {
-        result(lannot) += p_ * input_arr(annot);
-      } else {
-        result(lannot) = p_ * input_arr(annot);
-      }
-    };
-    antisymmetric_permutation(ParticleRange{range_perm.begin(), range_rank},
-                              callback);
-    return result;
-  };
-
-  // identity term with coefficient +1
-  result(lannot) = arr(lannot);
-
-  // process only ket permutations with coefficient norm_factor
-  if (ket_rank > 1) {
-    const auto bra_annot = bra_rank == 0 ? "" : ords_to_annot(bra_perm);
-    auto ket_result =
-        process_permutations(arr, ket_rank, ket_perm, bra_annot, false);
-
-    result(lannot) -= norm_factor * ket_result(lannot);
-  }
-
-  TA::DistArray<Args...>::wait_for_lazy_cleanup(result.world());
-  return result;
-}
-
-/// \brief This function is used to implement
-/// ResultPtr::biorthogonal_nns_project for btas::Tensor
-///
-/// \param arr The array to be "cleaned up"
-/// \param bra_rank The rank of the bra indices
-///
-/// \return The cleaned btas::Tensor.
-template <typename... Args>
-auto biorthogonal_nns_project_btas(btas::Tensor<Args...> const& arr,
-                                   size_t bra_rank) {
-  using ranges::views::concat;
-  using ranges::views::iota;
-  size_t const rank = arr.rank();
-  SEQUANT_ASSERT(bra_rank <= rank);
-  size_t const ket_rank = rank - bra_rank;
-
-  if (rank <= 4) {
-    return arr;
-  }
-
-  using numeric_type = typename btas::Tensor<Args...>::numeric_type;
-
-  size_t factorial_ket = 1;
-  for (size_t i = 2; i <= ket_rank; ++i) {
-    factorial_ket *= i;
-  }
-  numeric_type norm_factor = numeric_type(1) / numeric_type(factorial_ket);
-
-  perm_t bra_perm = iota(size_t{0}, bra_rank) | ranges::to<perm_t>;
-  perm_t ket_perm = iota(bra_rank, rank) | ranges::to<perm_t>;
-  const auto lannot = iota(size_t{0}, rank) | ranges::to<perm_t>;
-
-  auto process_permutations = [&lannot](const btas::Tensor<Args...>& input_arr,
-                                        size_t range_rank, perm_t range_perm,
-                                        const perm_t& other_perm, bool is_bra) {
-    if (range_rank <= 1) return input_arr;
-    btas::Tensor<Args...> result{input_arr.range()};
-    result.fill(0);
-
-    auto callback = [&]([[maybe_unused]] int parity) {
-      const auto annot =
-          is_bra ? concat(range_perm, other_perm) | ranges::to<perm_t>()
-                 : concat(other_perm, range_perm) | ranges::to<perm_t>();
-
-      // ignore parity, all permutations get same coefficient
-      numeric_type p_ = 1;
-      btas::Tensor<Args...> temp;
-      btas::permute(input_arr, lannot, temp, annot);
-      btas::scal(p_, temp);
-      result += temp;
-    };
-
-    antisymmetric_permutation(ParticleRange{range_perm.begin(), range_rank},
-                              callback);
-    return result;
-  };
-
-  // identity term with coefficient +1
-  auto result = arr;
-
-  // process only ket permutations with coefficient norm_factor
-  if (ket_rank > 1) {
-    const auto bra_annot = bra_rank == 0 ? perm_t{} : bra_perm;
-    auto ket_result =
-        process_permutations(arr, ket_rank, ket_perm, bra_annot, false);
-
-    btas::scal(norm_factor, ket_result);
-    result -= ket_result;
-  }
-
-  return result;
-}
-
 template <typename... Args>
 inline void log_result(Args const&... args) noexcept {
   auto& l = Logger::instance();
@@ -631,19 +481,6 @@ class Result {
   ///
   [[nodiscard]] virtual ResultPtr antisymmetrize(size_t bra_rank) const = 0;
 
-  /// \brief Implements "biorthogonal cleanup" of closed-shell
-  /// more compact spintraced equations produced via method of
-  /// <a href="https://arxiv.org/abs/1805.00565">Wang and Knizia</a>.
-  ///
-  /// For 3-body residual (`bra_rank=3`) this implements Eq. (41) of the
-  /// Wang/Knizia paper, same as the first line of Figure 1.
-  /// For 4-body residual this implements the first line of Figure 2.
-  /// The implementation is for arbitrary ranks.
-  /// @param bra_rank the particle rank of the residual tensor (i.e.
-  ///                 its order halved)
-  [[nodiscard]] virtual ResultPtr biorthogonal_nns_project(
-      size_t bra_rank) const = 0;
-
   [[nodiscard]] bool has_value() const noexcept;
 
   [[nodiscard]] virtual ResultPtr mult_by_phase(std::int8_t) const = 0;
@@ -755,11 +592,6 @@ class ResultScalar final : public Result {
     throw unimplemented_method("antisymmetrize");
   }
 
-  [[nodiscard]] ResultPtr biorthogonal_nns_project(
-      [[maybe_unused]] size_t bra_rank) const override {
-    throw unimplemented_method("biorthogonal_nns_project");
-  }
-
   [[nodiscard]] ResultPtr mult_by_phase(std::int8_t factor) const override {
     return eval_result<ResultScalar<T>>(value() * T(factor));
   }
@@ -901,12 +733,6 @@ class ResultTensorTA final : public Result {
         particle_antisymmetrize_ta(get<ArrayT>(), bra_rank));
   }
 
-  [[nodiscard]] ResultPtr biorthogonal_nns_project(
-      size_t bra_rank) const override {
-    return eval_result<this_type>(
-        biorthogonal_nns_project_ta(get<ArrayT>(), bra_rank));
-  }
-
  private:
   [[nodiscard]] std::size_t size_in_bytes() const final {
     auto& v = get<ArrayT>();
@@ -1055,13 +881,6 @@ class ResultTensorOfTensorTA final : public Result {
     return nullptr;
   }
 
-  [[nodiscard]] ResultPtr biorthogonal_nns_project(
-      [[maybe_unused]] size_t bra_rank) const override {
-    // or? throw unimplemented_method("biorthogonal_nns_project");
-    // not implemented yet, I think I need it for CSV
-    return nullptr;
-  }
-
  private:
   [[nodiscard]] std::size_t size_in_bytes() const final {
     auto& v = get<ArrayT>();
@@ -1164,12 +983,6 @@ class ResultTensorBTAS final : public Result {
         particle_antisymmetrize_btas(get<T>(), bra_rank));
   }
 
-  [[nodiscard]] ResultPtr biorthogonal_nns_project(
-      [[maybe_unused]] size_t bra_rank) const override {
-    return eval_result<ResultTensorBTAS<T>>(
-        biorthogonal_nns_project_btas(get<T>(), bra_rank));
-  }
-
  private:
   [[nodiscard]] std::size_t size_in_bytes() const final {
     static_assert(std::is_arithmetic_v<typename T::value_type>);