From e22433841a11cdcd4ff4bc4369cabb064f15f949 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 13 Sep 2023 15:41:01 -0700 Subject: [PATCH] changed wfn_constraint -> alpha_constraint and made appropriate changes --- include/macis/asci/alpha_constraint.hpp | 70 +++++++++++++++ include/macis/asci/determinant_search.hpp | 13 ++- include/macis/asci/mask_constraints.hpp | 101 +++++++++++----------- include/macis/asci/pt2.hpp | 10 +-- tests/asci.cxx | 24 ++--- 5 files changed, 144 insertions(+), 74 deletions(-) create mode 100644 include/macis/asci/alpha_constraint.hpp diff --git a/include/macis/asci/alpha_constraint.hpp b/include/macis/asci/alpha_constraint.hpp new file mode 100644 index 00000000..4d0ca2be --- /dev/null +++ b/include/macis/asci/alpha_constraint.hpp @@ -0,0 +1,70 @@ +#pragma once + +#include +namespace macis { + +template +class alpha_constraint { + +public: + using wfn_traits = WfnTraits; + using wfn_type = typename WfnTraits::wfn_type; + using spin_wfn_type = spin_wfn_t; + using spin_wfn_traits = wavefunction_traits; + 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 + inline auto overlap(WfnType state) const { + return spin_wfn_traits::count(c_mask_union(state)); + } + + template + inline bool satisfies_constraint(WfnType state) const { + return overlap(state) == count_ and + spin_wfn_traits::count(symmetric_difference(state) >> C_min_) == 0; + } + + + +}; + + +} diff --git a/include/macis/asci/determinant_search.hpp b/include/macis/asci/determinant_search.hpp index 447d95d2..cca0bacb 100644 --- a/include/macis/asci/determinant_search.hpp +++ b/include/macis/asci/determinant_search.hpp @@ -280,9 +280,6 @@ asci_contrib_container> asci_contributions_constraint( auto size_before = asci_pairs.size(); const double h_el_tol = asci_settings.h_el_tol; - const auto& [C, B, C_min] = con; - wfn_constraint alpha_con{ wfn_traits::alpha_string(C), wfn_traits::alpha_string(B), C_min}; - //wfn_t O = full_mask(norb); // Loop over unique alpha strings for(size_t i_alpha = 0; i_alpha < nuniq_alpha; ++i_alpha) { @@ -297,7 +294,7 @@ asci_contrib_container> asci_contributions_constraint( const auto& occ_beta = bcd.occ_beta; const auto& orb_ens_alpha = bcd.orb_ens_alpha; generate_constraint_singles_contributions_ss( - coeff, det|beta, alpha_con, occ_alpha, occ_beta, + coeff, det|beta, con, 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); } @@ -310,7 +307,7 @@ asci_contrib_container> asci_contributions_constraint( const auto& occ_beta = bcd.occ_beta; const auto& orb_ens_alpha = bcd.orb_ens_alpha; generate_constraint_doubles_contributions_ss( - coeff, det|beta, alpha_con, occ_alpha, occ_beta, + coeff, det|beta, con, occ_alpha, occ_beta, orb_ens_alpha.data(), G_pqrs, norb, h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs); } @@ -325,14 +322,14 @@ asci_contrib_container> asci_contributions_constraint( 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|beta, alpha_con, occ_alpha, occ_beta, vir_beta, + coeff, det|beta, con, 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, con)) { + if(satisfies_constraint(wfn_traits::alpha_string(det), con)) { for(const auto& bcd : uad[i_alpha].bcd) { const auto& beta = bcd.beta_string; const auto& coeff = bcd.coeff; @@ -368,7 +365,7 @@ asci_contrib_container> asci_contributions_constraint( }); asci_pairs.erase(it, asci_pairs.end()); - auto c_indices = bits_to_indices(C); + auto c_indices = bits_to_indices(con.C()); std::string c_string; for(int i = 0; i < c_indices.size(); ++i) c_string += std::to_string(c_indices[i]) + " "; diff --git a/include/macis/asci/mask_constraints.hpp b/include/macis/asci/mask_constraints.hpp index b8866cf9..a184c76a 100644 --- a/include/macis/asci/mask_constraints.hpp +++ b/include/macis/asci/mask_constraints.hpp @@ -12,33 +12,29 @@ #include #include +#include + namespace macis { -template -struct wfn_constraint { - wfn_t C; - wfn_t B; - unsigned C_min; -}; -template -bool satisfies_constraint(wfn_t det, wfn_constraint constraint) { - auto C = constraint.C; auto C_min = constraint.C_min; - return (det & C).count() == C.count() and ((det ^ C) >> C_min).count() == 0; +template +bool satisfies_constraint(wfn_t det, ConType C) { + return C.satisfies_constraint(det); } -template -auto generate_constraint_single_excitations(wfn_t det, wfn_constraint constraint) { - const auto C = constraint.C; const auto B = constraint.B; - if((det & C).count() < - (C.count() - 1)) // need to have at most one different from the constraint +template +auto generate_constraint_single_excitations(wfn_t det, ConType constraint) { + using spin_wfn_traits = typename ConType::spin_wfn_traits; + const auto C = constraint.C(); const auto B = constraint.B(); + + // need to have at most one different from the constraint + if(constraint.overlap(det) < (constraint.count()-1)) return std::make_pair(wfn_t(0), wfn_t(0)); auto o = det ^ C; auto v = (~det) & B; - if((o & C).count() == - 1) { // don't have to change this necessarily, but more clear without >= + if((o & C).count() == 1) { v = o & C; o ^= v; } @@ -50,9 +46,9 @@ auto generate_constraint_single_excitations(wfn_t det, wfn_constraint cons return std::make_pair(o, v); } -template -auto generate_constraint_double_excitations(wfn_t det, wfn_constraint constraint) { - const auto C = constraint.C; const auto B = constraint.B; +template +auto generate_constraint_double_excitations(wfn_t det, ConType constraint) { + const auto C = constraint.C(); const auto B = constraint.B(); // Occ/Vir pairs to generate excitations std::vector> O, V; @@ -102,8 +98,8 @@ auto generate_constraint_double_excitations(wfn_t det, wfn_constraint cons return std::make_tuple(O, V); } -template -void generate_constraint_singles(wfn_t det, wfn_constraint constraint, +template +void generate_constraint_singles(wfn_t det, ConType constraint, std::vector>& t_singles) { auto [o, v] = generate_constraint_single_excitations(det, constraint); const auto oc = o.count(); @@ -128,8 +124,8 @@ unsigned count_constraint_singles(Args&&... args) { return o.count() * v.count(); } -template -void generate_constraint_doubles(wfn_t det, wfn_constraint constraint, +template +void generate_constraint_doubles(wfn_t det, ConType constraint, std::vector>& t_doubles) { auto [O, V] = generate_constraint_double_excitations(det, constraint); @@ -147,9 +143,9 @@ void generate_constraint_doubles(wfn_t det, wfn_constraint constraint, * @param[in] T Triplet constraint mask * @param[in] B B mask (?) */ -template -unsigned count_constraint_doubles(wfn_t det, wfn_constraint constraint) { - const auto C = constraint.C; const auto B = constraint.B; +template +unsigned count_constraint_doubles(wfn_t det, ConType constraint) { + const auto C = constraint.C(); const auto B = constraint.B(); if((det & C) == 0) return 0; auto o = det ^ C; @@ -193,9 +189,9 @@ unsigned count_constraint_doubles(wfn_t det, wfn_constraint constraint) { return no_pairs * nv_pairs; } -template +template size_t constraint_histogram(wfn_t det, size_t n_os_singles, - size_t n_os_doubles, wfn_constraint constraint){ + size_t n_os_doubles, ConType constraint){ auto ns = count_constraint_singles(det, constraint); auto nd = count_constraint_doubles(det, constraint); @@ -461,16 +457,18 @@ auto dist_triplets_histogram(size_t norb, size_t ns_othr, size_t nd_othr, template auto make_triplet(unsigned i, unsigned j, unsigned k) { - wfn_constraint con; - - con.C = 0; - con.C.flip(i).flip(j).flip(k); - con.B = 1; - con.B <<= k; - con.B = con.B.to_ullong() - 1; - con.C_min = k; - - return con; + using wfn_type = wfn_t; + using wfn_traits = wavefunction_traits; + using constraint_type = alpha_constraint; + using string_type = typename constraint_type::constraint_type; + + string_type C = 0; + C.flip(i).flip(j).flip(k); + string_type B = 1; + B <<= k; + B = B.to_ullong() - 1; + + return constraint_type(C,B,k); } #ifdef MACIS_ENABLE_MPI @@ -479,6 +477,10 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr, size_t nd_othr, const std::vector>& unique_alpha, MPI_Comm comm) { + using wfn_type = wfn_t; + using wfn_traits = wavefunction_traits; + using constraint_type = alpha_constraint; + using string_type = typename constraint_type::constraint_type; auto world_rank = comm_rank(comm); auto world_size = comm_size(comm); @@ -488,18 +490,17 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr, std::vector workloads(world_size, 0); // Generate triplets + heuristic - std::vector, size_t>> constraint_sizes; + std::vector> constraint_sizes; constraint_sizes.reserve(norb * norb * norb); size_t total_work = 0; for(int t_i = 0; t_i < norb; ++t_i) for(int t_j = 0; t_j < t_i; ++t_j) for(int t_k = 0; t_k < t_j; ++t_k) { auto constraint = make_triplet(t_i, t_j, t_k); - const auto& [T, B, _] = constraint; size_t nw = 0; for(const auto& alpha : unique_alpha) { - nw += constraint_histogram(alpha, ns_othr, nd_othr, constraint); + nw += constraint_histogram(wfn_traits::alpha_string(alpha), ns_othr, nd_othr, constraint); } if(nw) constraint_sizes.emplace_back(constraint, nw); total_work += nw; @@ -509,7 +510,7 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr, for(size_t ilevel = 0; ilevel < nlevels; ++ilevel) { // Select constraints larger than average to be broken apart - std::vector, size_t>> tps_to_next; + std::vector> tps_to_next; { auto it = std::partition( constraint_sizes.begin(), constraint_sizes.end(), @@ -525,20 +526,20 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr, // Break apart constraints for(auto [c, nw_trip] : tps_to_next) { - const auto C_min = c.C_min; + const auto C_min = c.C_min(); // Loop over possible constraints with one more element for(auto q_l = 0; q_l < C_min; ++q_l) { // Generate masks / counts - wfn_constraint c_next = c; - c_next.C.flip(q_l); - c_next.B >>= (C_min - q_l); - c_next.C_min = q_l; + string_type cn_C = c.C(); + cn_C.flip(q_l); + string_type cn_B = c.B() >> (C_min - q_l); + constraint_type c_next(cn_C, cn_B, q_l); size_t nw = 0; for(const auto& alpha : unique_alpha) { - nw += constraint_histogram(alpha, ns_othr, nd_othr, c_next); + nw += constraint_histogram(wfn_traits::alpha_string(alpha), ns_othr, nd_othr, c_next); } if(nw) constraint_sizes.emplace_back(c_next, nw); total_work += nw; @@ -569,7 +570,7 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr, [](const auto& a, const auto& b) { return a.second > b.second; }); // Assign work - std::vector> constraints; + std::vector constraints; constraints.reserve(constraint_sizes.size() / world_size); for(auto [c, nw] : constraint_sizes) { diff --git a/include/macis/asci/pt2.hpp b/include/macis/asci/pt2.hpp index b560597c..d7cf7c46 100644 --- a/include/macis/asci/pt2.hpp +++ b/include/macis/asci/pt2.hpp @@ -130,8 +130,6 @@ double asci_pt2_constraint( 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_constraint alpha_con{ wfn_traits::alpha_string(C), wfn_traits::alpha_string(B), C_min}; // Loop over unique alpha strings for(size_t i_alpha = 0; i_alpha < nuniq_alpha; ++i_alpha) { @@ -146,7 +144,7 @@ double asci_pt2_constraint( const auto& occ_beta = bcd.occ_beta; const auto& orb_ens_alpha = bcd.orb_ens_alpha; generate_constraint_singles_contributions_ss( - coeff, det|beta, alpha_con, occ_alpha, occ_beta, + coeff, det|beta, con, 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); } @@ -159,7 +157,7 @@ double asci_pt2_constraint( const auto& occ_beta = bcd.occ_beta; const auto& orb_ens_alpha = bcd.orb_ens_alpha; generate_constraint_doubles_contributions_ss( - coeff, det|beta, alpha_con, occ_alpha, occ_beta, + coeff, det|beta, con, occ_alpha, occ_beta, orb_ens_alpha.data(), G_pqrs, norb, h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs); } @@ -174,14 +172,14 @@ double asci_pt2_constraint( 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|beta, alpha_con, occ_alpha, occ_beta, vir_beta, + coeff, det|beta, con, 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, con)) { + if(satisfies_constraint(wfn_traits::alpha_string(det), con)) { for(const auto& bcd : uad[i_alpha].bcd) { const auto& beta = bcd.beta_string; const auto& coeff = bcd.coeff; diff --git a/tests/asci.cxx b/tests/asci.cxx index f088f853..05ef72c0 100644 --- a/tests/asci.cxx +++ b/tests/asci.cxx @@ -47,16 +47,18 @@ size_t top_set_ordinal(std::bitset word, size_t NSet) { template auto make_quad(unsigned i, unsigned j, unsigned k, unsigned l) { - macis::wfn_constraint con; + using wfn_type = macis::wfn_t; + using wfn_traits = macis::wavefunction_traits; + using constraint_type = macis::alpha_constraint; + using string_type = typename constraint_type::constraint_type; - con.C = 0; - con.C.flip(i).flip(j).flip(k).flip(l); - con.B = 1; - con.B <<= l; - con.B = con.B.to_ullong() - 1; - con.C_min = l; + string_type C = 0; + C.flip(i).flip(j).flip(k).flip(l); + string_type B = 1; + B <<= l; + B = B.to_ullong() - 1; - return con; + return constraint_type(C,B,l); } TEST_CASE("Triplets") { @@ -64,6 +66,8 @@ TEST_CASE("Triplets") { size_t norb = 32; using wfn_less_comparator = macis::bitset_less_comparator; + using wfn_type = macis::wfn_t; + using wfn_traits = macis::wavefunction_traits; // Generate ficticious wfns std::vector> wfn_a = { @@ -186,7 +190,7 @@ TEST_CASE("Triplets") { const size_t n_doubles = (n_singles * (n_singles - nocc - nvir + 1)) / 4; new_triplet_hist[label] += macis::constraint_histogram( - det, n_singles, n_doubles, constraint); + wfn_traits::alpha_string(det), n_singles, n_doubles, constraint); } } @@ -213,7 +217,7 @@ TEST_CASE("Triplets") { const size_t n_doubles = (n_singles * (n_singles - nocc - nvir + 1)) / 4; new_quad_hist[label] += macis::constraint_histogram( - det, n_singles, n_doubles, constraint); + wfn_traits::alpha_string(det), n_singles, n_doubles, constraint); } }