PT2

CMakeLists.txt

@@ -11,6 +11,7 @@ project(MACIS VERSION 0.1 LANGUAGES C CXX)
 option( MACIS_ENABLE_MPI    "Enable MPI Bindings"    ON )
 option( MACIS_ENABLE_OPENMP "Enable OpenMP Bindings" ON )
 option( MACIS_ENABLE_BOOST  "Enable Boost" ON )
+option( MACIS_ENABLE_TREXIO "Enable TrexIO" ON )
 
 # CMake Modules
 list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake )

include/macis/asci/alpha_constraint.hpp

@@ -0,0 +1,82 @@
+/*
+ * MACIS Copyright (c) 2023, The Regents of the University of California,
+ * through Lawrence Berkeley National Laboratory (subject to receipt of
+ * any required approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * See LICENSE.txt for details
+ */
+
+#pragma once
+
+#include <macis/wfn/raw_bitset.hpp>
+namespace macis {
+
+template <typename WfnTraits>
+class alpha_constraint {
+ public:
+  using wfn_traits = WfnTraits;
+  using wfn_type = typename WfnTraits::wfn_type;
+  using spin_wfn_type = spin_wfn_t<wfn_type>;
+
+  using constraint_type = spin_wfn_type;
+  using constraint_traits = wavefunction_traits<spin_wfn_type>;
+
+ private:
+  constraint_type C_;
+  constraint_type B_;
+  uint32_t C_min_;
+  uint32_t count_;
+
+ public:
+  alpha_constraint(constraint_type C, constraint_type B, uint32_t C_min)
+      : C_(C), B_(B), C_min_(C_min), count_(constraint_traits::count(C)) {}
+
+  alpha_constraint(const alpha_constraint&) = default;
+  alpha_constraint& operator=(const alpha_constraint&) = default;
+
+  alpha_constraint(alpha_constraint&& other) noexcept = default;
+  alpha_constraint& operator=(alpha_constraint&&) noexcept = default;
+
+  inline auto C() const { return C_; }
+  inline auto B() const { return B_; }
+  inline auto C_min() const { return C_min_; }
+  inline auto count() const { return count_; }
+
+  inline spin_wfn_type c_mask_union(spin_wfn_type state) const {
+    return state & C_;
+  }
+  inline spin_wfn_type b_mask_union(spin_wfn_type state) const {
+    return state & B_;
+  }
+
+  inline spin_wfn_type symmetric_difference(spin_wfn_type state) const {
+    return state ^ C_;
+  }
+  // inline spin_wfn_type symmetric_difference(wfn_type state) const {
+  //   return symmetric_difference(wfn_traits::alpha_string(state));
+  // }
+
+  template <typename WfnType>
+  inline auto overlap(WfnType state) const {
+    return constraint_traits::count(c_mask_union(state));
+  }
+
+  template <typename WfnType>
+  inline bool satisfies_constraint(WfnType state) const {
+    return overlap(state) == count_ and
+           constraint_traits::count(symmetric_difference(state) >> C_min_) == 0;
+  }
+
+  static alpha_constraint make_triplet(unsigned i, unsigned j, unsigned k) {
+    constraint_type C = 0;
+    C.flip(i).flip(j).flip(k);
+    // constraint_type B = 1;
+    // static_assert(B.size() <= 64, "ULLONG NOT POSSIBLE HERE");
+    // B <<= k;
+    // B = B.to_ullong() - 1;
+    constraint_type B = full_mask<B.size()>(k);
+    return alpha_constraint(C, B, k);
+  }
+};
+
+}  // namespace macis

include/macis/asci/determinant_contributions.hpp

@@ -16,22 +16,27 @@ namespace macis {
 template <typename WfnT>
 struct asci_contrib {
   WfnT state;
-  double rv;
+  double c_times_matel;
+  double h_diag;
+
+  auto rv() const { return c_times_matel / h_diag; }
+  auto pt2() const { return rv() * c_times_matel; }
 };
 
 template <typename WfnT>
 using asci_contrib_container = std::vector<asci_contrib<WfnT>>;
 
-template <size_t N, size_t NShift>
+template <Spin Sigma, typename WfnType, typename SpinWfnType>
 void append_singles_asci_contributions(
-    double coeff, wfn_t<2 * N> state_full, wfn_t<N> state_same,
+    double coeff, WfnType state_full, SpinWfnType state_same,
     const std::vector<uint32_t>& occ_same,
     const std::vector<uint32_t>& vir_same,
     const std::vector<uint32_t>& occ_othr, const double* eps_same,
     const double* T_pq, const size_t LDT, const double* G_kpq, const size_t LDG,
     const double* V_kpq, const size_t LDV, double h_el_tol, double root_diag,
-    double E0, HamiltonianGenerator<2 * N>& ham_gen,
-    asci_contrib_container<wfn_t<2 * N>>& asci_contributions) {
+    double E0, const HamiltonianGeneratorBase<double>& ham_gen,
+    asci_contrib_container<WfnType>& asci_contributions) {
+  using wfn_traits = wavefunction_traits<WfnType>;
   const auto LDG2 = LDG * LDG;
   const auto LDV2 = LDV * LDV;
   for(auto i : occ_same)
@@ -44,11 +49,11 @@ void append_singles_asci_contributions(
       for(auto p : occ_othr) h_el += V_ov[p];
 
       // Early Exit
-      if(std::abs(h_el) < h_el_tol) continue;
+      if(std::abs(coeff * h_el) < h_el_tol) continue;
 
       // Calculate Excited Determinant
-      auto ex_det = state_full;
-      ex_det.flip(i + NShift).flip(a + NShift);
+      auto ex_det = wfn_traits::template single_excitation_no_check<Sigma>(
+          state_full, i, a);
 
       // Calculate Excitation Sign in a Canonical Way
       auto sign = single_excitation_sign(state_same, a, i);
@@ -57,22 +62,25 @@ void append_singles_asci_contributions(
       // Calculate fast diagonal matrix element
       auto h_diag =
           ham_gen.fast_diag_single(eps_same[i], eps_same[a], i, a, root_diag);
-      h_el /= (E0 - h_diag);
+      // h_el /= (E0 - h_diag);
 
       // Append to return values
-      asci_contributions.push_back({ex_det, coeff * h_el});
+      asci_contributions.push_back({ex_det, coeff * h_el, E0 - h_diag});
 
     }  // Loop over single extitations
 }
 
-template <size_t N, size_t NShift>
+template <Spin Sigma, typename WfnType, typename SpinWfnType>
 void append_ss_doubles_asci_contributions(
-    double coeff, wfn_t<2 * N> state_full, wfn_t<N> state_spin,
-    const std::vector<uint32_t>& ss_occ, const std::vector<uint32_t>& vir,
-    const std::vector<uint32_t>& os_occ, const double* eps_same,
-    const double* G, size_t LDG, double h_el_tol, double root_diag, double E0,
-    HamiltonianGenerator<2 * N>& ham_gen,
-    asci_contrib_container<wfn_t<2 * N>>& asci_contributions) {
+    double coeff, WfnType state_full, SpinWfnType state_same,
+    SpinWfnType state_other, const std::vector<uint32_t>& ss_occ,
+    const std::vector<uint32_t>& vir, const std::vector<uint32_t>& os_occ,
+    const double* eps_same, const double* G, size_t LDG, double h_el_tol,
+    double root_diag, double E0,
+    const HamiltonianGeneratorBase<double>& ham_gen,
+    asci_contrib_container<WfnType>& asci_contributions) {
+  using wfn_traits = wavefunction_traits<WfnType>;
+  using spin_wfn_traits = wavefunction_traits<SpinWfnType>;
   const size_t nocc = ss_occ.size();
   const size_t nvir = vir.size();
 
@@ -90,8 +98,9 @@ void append_ss_doubles_asci_contributions(
           const auto jb = b + j * LDG;
           const auto G_aibj = G_ai[jb];
 
-          if(std::abs(G_aibj) < h_el_tol) continue;
+          if(std::abs(coeff * G_aibj) < h_el_tol) continue;
 
+#if 0
           // Calculate excited determinant string (spin)
           const auto full_ex_spin = wfn_t<N>(0).flip(i).flip(j).flip(a).flip(b);
           auto ex_det_spin = state_spin ^ full_ex_spin;
@@ -102,6 +111,21 @@ void append_ss_doubles_asci_contributions(
           // Calculate full excited determinant
           const auto full_ex = expand_bitset<2 * N>(full_ex_spin) << NShift;
           auto ex_det = state_full ^ full_ex;
+#else
+          // TODO: Can this be made faster since the orbital indices are known
+          //       in advance?
+          // Compute excited determinant (spin)
+          const auto full_ex_spin = spin_wfn_traits::double_excitation_no_check(
+              SpinWfnType(0), i, j, a, b);
+          const auto ex_det_spin = state_same ^ full_ex_spin;
+
+          // Calculate the sign in a canonical way
+          double sign = doubles_sign(state_same, ex_det_spin, full_ex_spin);
+
+          // Calculate full excited determinant
+          auto ex_det =
+              wfn_traits::template from_spin<Sigma>(ex_det_spin, state_other);
+#endif
 
           // Update sign of matrix element
           auto h_el = sign * G_aibj;
@@ -110,25 +134,27 @@ void append_ss_doubles_asci_contributions(
           auto h_diag =
               ham_gen.fast_diag_ss_double(eps_same[i], eps_same[j], eps_same[a],
                                           eps_same[b], i, j, a, b, root_diag);
-          h_el /= (E0 - h_diag);
+          // h_el /= (E0 - h_diag);
 
           // Append {det, c*h_el}
-          asci_contributions.push_back({ex_det, coeff * h_el});
+          asci_contributions.push_back({ex_det, coeff * h_el, E0 - h_diag});
 
         }  // Restricted BJ loop
     }      // AI Loop
 }
 
-template <size_t N>
+template <typename WfnType, typename SpinWfnType>
 void append_os_doubles_asci_contributions(
-    double coeff, wfn_t<2 * N> state_full, wfn_t<N> state_alpha,
-    wfn_t<N> state_beta, const std::vector<uint32_t>& occ_alpha,
+    double coeff, WfnType state_full, SpinWfnType state_alpha,
+    SpinWfnType state_beta, const std::vector<uint32_t>& occ_alpha,
     const std::vector<uint32_t>& occ_beta,
     const std::vector<uint32_t>& vir_alpha,
     const std::vector<uint32_t>& vir_beta, const double* eps_alpha,
     const double* eps_beta, const double* V, size_t LDV, double h_el_tol,
-    double root_diag, double E0, HamiltonianGenerator<2 * N>& ham_gen,
-    asci_contrib_container<wfn_t<2 * N>>& asci_contributions) {
+    double root_diag, double E0,
+    const HamiltonianGeneratorBase<double>& ham_gen,
+    asci_contrib_container<WfnType>& asci_contributions) {
+  using wfn_traits = wavefunction_traits<WfnType>;
   const size_t LDV2 = LDV * LDV;
   for(auto i : occ_alpha)
     for(auto a : vir_alpha) {
@@ -140,21 +166,26 @@ void append_os_doubles_asci_contributions(
           const auto jb = b + j * LDV;
           const auto V_aibj = V_ai[jb * LDV2];
 
-          if(std::abs(V_aibj) < h_el_tol) continue;
+          if(std::abs(coeff * V_aibj) < h_el_tol) continue;
 
           double sign_beta = single_excitation_sign(state_beta, b, j);
           double sign = sign_alpha * sign_beta;
-          auto ex_det = state_full;
-          ex_det.flip(a).flip(i).flip(j + N).flip(b + N);
+          // auto ex_det = state_full;
+          // ex_det.flip(a).flip(i).flip(j + N).flip(b + N);
+          auto ex_det =
+              wfn_traits::template single_excitation_no_check<Spin::Alpha>(
+                  state_full, a, i);
+          ex_det = wfn_traits::template single_excitation_no_check<Spin::Beta>(
+              ex_det, b, j);
           auto h_el = sign * V_aibj;
 
           // Evaluate fast diagonal element
           auto h_diag = ham_gen.fast_diag_os_double(eps_alpha[i], eps_beta[j],
                                                     eps_alpha[a], eps_beta[b],
                                                     i, j, a, b, root_diag);
-          h_el /= (E0 - h_diag);
+          // h_el /= (E0 - h_diag);
 
-          asci_contributions.push_back({ex_det, coeff * h_el});
+          asci_contributions.push_back({ex_det, coeff * h_el, E0 - h_diag});
         }  // BJ loop
     }      // AI loop
 }