Initial implementation of EN-PT2

include/macis/asci/determinant_sort.hpp

@@ -66,7 +66,7 @@ PairIterator sort_and_accumulate_asci_pairs(PairIterator pairs_begin,
     // Accumulate
     else {
       cur_it->c_times_matel += it->c_times_matel;
-      it->c_times_matel = 0;  // Zero out to expose potential bugs
+      it->c_times_matel = NAN;  // Zero out to expose potential bugs
     }
   }
 

include/macis/asci/mask_constraints.hpp

@@ -258,11 +258,11 @@ void generate_constraint_singles_contributions_ss(
   const size_t LDV2 = LDV * LDV;
   for(int ii = 0; ii < no; ++ii) {
     const auto i = fls(o);
-    o.flip(i);
+    o.flip(i); // Disable "i"-bit so it's not used in FLS next iteration
     auto v_cpy = v;
     for(int aa = 0; aa < nv; ++aa) {
       const auto a = fls(v_cpy);
-      v_cpy.flip(a);
+      v_cpy.flip(a); // Disable "a"-bit so it's not used in FLS next iteration
 
       double h_el = T_pq[a + i * LDT];
       const double* G_ov = G_kpq + a * LDG + i * LDG2;

include/macis/asci/pt2.hpp

@@ -0,0 +1,231 @@
+/*
+ * 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/asci/determinant_contributions.hpp>
+#include <macis/asci/determinant_sort.hpp>
+
+namespace macis {
+
+
+template <size_t N>
+double asci_pt2_constraint(
+    wavefunction_iterator_t<N> cdets_begin, wavefunction_iterator_t<N> cdets_end,
+    const double E_ASCI, const std::vector<double>& C, size_t norb,
+    const double* T_pq, const double* G_red, const double* V_red, 
+    const double* G_pqrs, const double* V_pqrs, HamiltonianGenerator<N>& ham_gen,
+    MPI_Comm comm) {
+
+  using clock_type = std::chrono::high_resolution_clock;
+  using duration_type = std::chrono::duration<double, std::milli>;
+
+  auto logger = spdlog::get("asci_search");
+  const size_t ncdets = std::distance(cdets_begin, cdets_end);
+  std::cout << "NDETS PT = " << ncdets <<  " " << C.size() << std::endl;
+  std::cout << "PT E0    = " << E_ASCI << std::endl;
+
+  std::vector<uint32_t> occ_alpha, vir_alpha;
+  std::vector<uint32_t> occ_beta, vir_beta;
+
+  // Get unique alpha strings
+  std::vector<wfn_t<N>> uniq_alpha_wfn(cdets_begin, cdets_end);
+  std::transform(uniq_alpha_wfn.begin(), uniq_alpha_wfn.end(),
+                 uniq_alpha_wfn.begin(),
+                 [=](const auto& w) { return w & full_mask<N / 2, N>(); });
+  std::sort(uniq_alpha_wfn.begin(), uniq_alpha_wfn.end(),
+            bitset_less_comparator<N>{});
+  {
+    auto it = std::unique(uniq_alpha_wfn.begin(), uniq_alpha_wfn.end());
+    uniq_alpha_wfn.erase(it, uniq_alpha_wfn.end());
+  }
+  const size_t nuniq_alpha = uniq_alpha_wfn.size();
+
+  // For each unique alpha, create a list of beta string and store metadata
+  struct beta_coeff_data {
+    wfn_t<N> beta_string;
+    std::vector<uint32_t> occ_beta;
+    std::vector<uint32_t> vir_beta;
+    std::vector<double> orb_ens_alpha;
+    std::vector<double> orb_ens_beta;
+    double coeff;
+    double h_diag;
+
+    beta_coeff_data(double c, size_t norb,
+                    const std::vector<uint32_t>& occ_alpha, wfn_t<N> w,
+                    const HamiltonianGenerator<N>& ham_gen) {
+      coeff = c;
+
+      // Compute Beta string
+      const auto beta_shift = w >> N / 2;
+      // Reduce the number of times things shift in inner loop
+      beta_string = beta_shift << N / 2;
+
+      // Compute diagonal matrix element
+      h_diag = ham_gen.matrix_element(w, w);
+
+      // Compute occ/vir for beta string
+      bitset_to_occ_vir(norb, beta_shift, occ_beta, vir_beta);
+
+      // Precompute orbital energies
+      orb_ens_alpha = ham_gen.single_orbital_ens(norb, occ_alpha, occ_beta);
+      orb_ens_beta = ham_gen.single_orbital_ens(norb, occ_beta, occ_alpha);
+    }
+  };
+
+  struct unique_alpha_data {
+    std::vector<beta_coeff_data> bcd;
+  };
+
+  std::vector<unique_alpha_data> uad(nuniq_alpha);
+  for(auto i = 0; i < nuniq_alpha; ++i) {
+    const auto wfn_a = uniq_alpha_wfn[i];
+    std::vector<uint32_t> occ_alpha, vir_alpha;
+    bitset_to_occ_vir(norb, wfn_a, occ_alpha, vir_alpha);
+    for(auto j = 0; j < ncdets; ++j) {
+      const auto w = *(cdets_begin + j);
+      if((w & full_mask<N / 2, N>()) == wfn_a) {
+        uad[i].bcd.emplace_back(C[j], norb, occ_alpha, w, ham_gen);
+      }
+    }
+  }
+
+  auto world_rank = comm_rank(comm);
+  auto world_size = comm_size(comm);
+
+  const auto n_occ_alpha = uniq_alpha_wfn[0].count();
+  const auto n_vir_alpha = norb - n_occ_alpha;
+  const auto n_sing_alpha = n_occ_alpha * n_vir_alpha;
+  const auto n_doub_alpha = (n_sing_alpha * (n_sing_alpha - norb + 1)) / 4;
+
+  logger->info("  * NS = {} ND = {}", n_sing_alpha, n_doub_alpha);
+
+  auto gen_c_st = clock_type::now();
+  auto constraints =
+      dist_constraint_general(0, norb,
+                              n_sing_alpha, n_doub_alpha, uniq_alpha_wfn, comm);
+  auto gen_c_en = clock_type::now();
+  duration_type gen_c_dur = gen_c_en - gen_c_st;
+  logger->info("  * GEN_DUR = {:.2e} ms", gen_c_dur.count());
+
+  size_t max_size = std::min(100000000ul,
+                             ncdets * (2 * n_sing_alpha +  // AA + BB
+                                       2 * n_doub_alpha +  // AAAA + BBBB
+                                       n_sing_alpha * n_sing_alpha  // AABB
+                                       ));
+  double EPT2 = 0.0;
+
+  // Process ASCI pair contributions for each constraint
+  #pragma omp parallel 
+  {
+  asci_contrib_container<wfn_t<N>> asci_pairs;
+  asci_pairs.reserve(max_size);
+  #pragma omp for reduction(+:EPT2)
+  for(size_t ic = 0; ic < constraints.size(); ++ic) {
+    const auto& con = constraints[ic];
+    //std::cout << std::distance(&constraints[0], &con) << "/" << constraints.size() << std::endl;
+    const double h_el_tol = 1e-16;
+    const auto& [C, B, C_min] = con;
+    wfn_t<N> O = full_mask<N>(norb);
+
+    // Loop over unique alpha strings
+    for(size_t i_alpha = 0; i_alpha < nuniq_alpha; ++i_alpha) {
+      const auto& det = uniq_alpha_wfn[i_alpha];
+      const auto occ_alpha = bits_to_indices(det);
+
+      // AA excitations
+      for(const auto& bcd : uad[i_alpha].bcd) {
+        const auto& beta = bcd.beta_string;
+        const auto& coeff = bcd.coeff;
+        const auto& h_diag = bcd.h_diag;
+        const auto& occ_beta = bcd.occ_beta;
+        const auto& orb_ens_alpha = bcd.orb_ens_alpha;
+        generate_constraint_singles_contributions_ss(
+            coeff, det, C, O, B, beta, occ_alpha, occ_beta,
+            orb_ens_alpha.data(), T_pq, norb, G_red, norb, V_red, norb,
+            h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs);
+      }
+
+      // AAAA excitations
+      for(const auto& bcd : uad[i_alpha].bcd) {
+        const auto& beta = bcd.beta_string;
+        const auto& coeff = bcd.coeff;
+        const auto& h_diag = bcd.h_diag;
+        const auto& occ_beta = bcd.occ_beta;
+        const auto& orb_ens_alpha = bcd.orb_ens_alpha;
+        generate_constraint_doubles_contributions_ss(
+            coeff, det, C, O, B, beta, occ_alpha, occ_beta,
+            orb_ens_alpha.data(), G_pqrs, norb, h_el_tol, h_diag, E_ASCI,
+            ham_gen, asci_pairs);
+      }
+
+      // AABB excitations
+      for(const auto& bcd : uad[i_alpha].bcd) {
+        const auto& beta = bcd.beta_string;
+        const auto& coeff = bcd.coeff;
+        const auto& h_diag = bcd.h_diag;
+        const auto& occ_beta = bcd.occ_beta;
+        const auto& vir_beta = bcd.vir_beta;
+        const auto& orb_ens_alpha = bcd.orb_ens_alpha;
+        const auto& orb_ens_beta = bcd.orb_ens_beta;
+        generate_constraint_doubles_contributions_os(
+            coeff, det, C, O, B, beta, occ_alpha, occ_beta, vir_beta,
+            orb_ens_alpha.data(), orb_ens_beta.data(), V_pqrs, norb, h_el_tol,
+            h_diag, E_ASCI, ham_gen, asci_pairs);
+      }
+
+      // If the alpha determinant satisfies the constraint,
+      // append BB and BBBB excitations
+      if(satisfies_constraint(det, C, C_min)) {
+        for(const auto& bcd : uad[i_alpha].bcd) {
+          const auto& beta = bcd.beta_string;
+          const auto& coeff = bcd.coeff;
+          const auto& h_diag = bcd.h_diag;
+          const auto& occ_beta = bcd.occ_beta;
+          const auto& vir_beta = bcd.vir_beta;
+          const auto& eps_beta = bcd.orb_ens_beta;
+
+          const auto state = det | beta;
+          const auto state_beta = bitset_hi_word(beta);
+          // BB Excitations
+          append_singles_asci_contributions<(N / 2), (N / 2)>(
+              coeff, state, state_beta, occ_beta, vir_beta, occ_alpha,
+              eps_beta.data(), T_pq, norb, G_red, norb, V_red, norb, h_el_tol,
+              h_diag, E_ASCI, ham_gen, asci_pairs);
+
+          // BBBB Excitations
+          append_ss_doubles_asci_contributions<N / 2, N / 2>(
+              coeff, state, state_beta, occ_beta, vir_beta, occ_alpha,
+              eps_beta.data(), G_pqrs, norb, h_el_tol, h_diag, E_ASCI, ham_gen,
+              asci_pairs);
+
+          // No excition - to remove for PT2
+          asci_pairs.push_back({state, std::numeric_limits<double>::infinity(), 1.0});
+        }  // Beta Loop
+      }    // Triplet Check
+
+    }    // Unique Alpha Loop
+
+    double EPT2_local = 0.0;
+    // Local S&A for each quad + update EPT2
+    {
+      auto uit = sort_and_accumulate_asci_pairs(asci_pairs.begin(), asci_pairs.end());
+      for(auto it = asci_pairs.begin(); it != uit; ++it) {
+        //if(std::find(cdets_begin, cdets_end, it->state) == cdets_end)
+        if(!std::isinf(it->c_times_matel))
+        EPT2_local += (it->c_times_matel * it->c_times_matel) / it->h_diag;
+      } 
+      asci_pairs.clear();
+    }
+
+    EPT2 += EPT2_local;
+  }  // Constraint Loop
+  }
+
+  return EPT2;
+}
+}

tests/standalone_driver.cxx

@@ -16,6 +16,7 @@
 #include <iostream>
 #include <macis/asci/grow.hpp>
 #include <macis/asci/refine.hpp>
+#include <macis/asci/pt2.hpp>
 #include <macis/hamiltonian_generator/double_loop.hpp>
 #include <macis/util/cas.hpp>
 #include <macis/util/detail/rdm_files.hpp>
@@ -277,7 +278,9 @@ int main(int argc, char** argv) {
     std::vector<double> active_ordm(n_active * n_active);
     std::vector<double> active_trdm(active_ordm.size() * active_ordm.size());
 
+    bool pt2 = true;
     double E0 = 0;
+    double EPT2 = 0;
 
     // CI
     if(job == Job::CI) {
@@ -378,9 +381,22 @@ int main(int argc, char** argv) {
           sparsexx::write_dist_mm("ham.mtx", H, 1);
         }
 #endif
+        if(pt2) {
+          EPT2 = macis::asci_pt2_constraint( dets.begin(), dets.end(), E0 -(E_inactive + E_core), C, n_active,
+            ham_gen.T(), ham_gen.G_red(), ham_gen.V_red(),ham_gen.G(), ham_gen.V(), 
+            ham_gen MACIS_MPI_CODE(, MPI_COMM_WORLD));
+        }
+      }
+
+      console->info("E(CI)     = {:.12f} Eh", E0);
+
+      if(pt2) {
+        console->info("E(PT2)    = {:.12f} Eh", EPT2);
+        console->info("E(CI+PT2) = {:.12f} Eh", E0 + EPT2);
       }
 
-      // MCSCF
+
+    // MCSCF
     } else if(job == Job::MCSCF) {
       // Possibly read active RDMs
       if(rdm_fname.size()) {
@@ -410,9 +426,9 @@ int main(int argc, char** argv) {
           NumVirtual(n_virtual), E_core, T.data(), norb, V.data(), norb,
           active_ordm.data(), n_active, active_trdm.data(),
           n_active MACIS_MPI_CODE(, MPI_COMM_WORLD));
-    }
 
-    console->info("E(CI)  = {:.12f} Eh", E0);
+      console->info("E(CASSCF)  = {:.12f} Eh", E0);
+    }
 
     // Write FCIDUMP file if requested
     if(fci_out_fname.size())