Merge 0b82008 into c103c74

Author: wavefunction91

include/macis/asci/alpha_constraint.hpp

@@ -0,0 +1,70 @@
+#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 spin_wfn_traits = wavefunction_traits<spin_wfn_type>;
+  using constraint_type = 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_(spin_wfn_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 spin_wfn_traits::count(c_mask_union(state));
+  }
+
+  template <typename WfnType>
+  inline bool satisfies_constraint(WfnType state) const {
+    return overlap(state) == count_ and 
+           spin_wfn_traits::count(symmetric_difference(state) >> C_min_) == 0;
+  }
+
+
+
+};
+
+
+}

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; }
 };
 
 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)
@@ -47,8 +52,7 @@ void append_singles_asci_contributions(
       if(std::abs(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 +61,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,
+    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,
-    HamiltonianGenerator<2 * N>& ham_gen,
-    asci_contrib_container<wfn_t<2 * N>>& asci_contributions) {
+    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();
 
@@ -92,6 +99,7 @@ void append_ss_doubles_asci_contributions(
 
           if(std::abs(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 +110,19 @@ 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 +131,26 @@ 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) {
@@ -144,17 +166,19 @@ void append_os_doubles_asci_contributions(
 
           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
 }