Skip to content

Commit

Permalink
Merge branch 'feature/pt2' of github.com:wavefunction91/MACIS into fe…
Browse files Browse the repository at this point in the history 
…ature/pt2
  • Loading branch information
David Williams-Young committed Oct 30, 2023
2 parents e068dc5 + 12cb760 commit ddb59ed
Show file tree
Hide file tree
Showing 12 changed files with 348 additions and 338 deletions.
125 changes: 62 additions & 63 deletions include/macis/asci/determinant_search.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ asci_contrib_container<wfn_t<N>> asci_contributions_constraint(
using wfn_traits = wavefunction_traits<wfn_t<N>>;
using spin_wfn_type = spin_wfn_t<wfn_t<N>>;
using spin_wfn_traits = wavefunction_traits<spin_wfn_type>;
using wfn_comp = typename wfn_traits::spin_comparator;
using wfn_comp = typename wfn_traits::spin_comparator;
if(!std::is_sorted(cdets_begin, cdets_end, wfn_comp{}))
throw std::runtime_error("PT2 Only Works with Sorted Wfns");

Expand Down Expand Up @@ -264,20 +264,19 @@ asci_contrib_container<wfn_t<N>> asci_contributions_constraint(
}
};


auto uniq_alpha = get_unique_alpha(cdets_begin, cdets_end);
const size_t nuniq_alpha = uniq_alpha.size();
std::vector<wfn_t<N>> uniq_alpha_wfn(nuniq_alpha);
std::transform(uniq_alpha.begin(), uniq_alpha.end(), uniq_alpha_wfn.begin(),
[](const auto& p) { return wfn_traits::from_spin(p.first,0); });

std::transform(
uniq_alpha.begin(), uniq_alpha.end(), uniq_alpha_wfn.begin(),
[](const auto& p) { return wfn_traits::from_spin(p.first, 0); });

using unique_alpha_data = std::vector<beta_coeff_data>;
std::vector<unique_alpha_data> uad(nuniq_alpha);
for(auto i = 0, iw = 0; i < nuniq_alpha; ++i) {
std::vector<uint32_t> occ_alpha, vir_alpha;
spin_wfn_traits::state_to_occ_vir(norb, uniq_alpha[i].first,
occ_alpha, vir_alpha);
spin_wfn_traits::state_to_occ_vir(norb, uniq_alpha[i].first, occ_alpha,
vir_alpha);

const auto nbeta = uniq_alpha[i].second;
uad[i].reserve(nbeta);
Expand All @@ -289,7 +288,7 @@ asci_contrib_container<wfn_t<N>> asci_contributions_constraint(

#endif

//const auto n_occ_alpha = wfn_traits::count(uniq_alpha_wfn[0]);
// const auto n_occ_alpha = wfn_traits::count(uniq_alpha_wfn[0]);
const auto n_occ_alpha = spin_wfn_traits::count(uniq_alpha[0].first);
const auto n_vir_alpha = norb - n_occ_alpha;
const auto n_sing_alpha = n_occ_alpha * n_vir_alpha;
Expand Down Expand Up @@ -327,9 +326,9 @@ asci_contrib_container<wfn_t<N>> asci_contributions_constraint(
}

auto gen_c_st = clock_type::now();
auto constraints =
dist_constraint_general<wfn_t<N>>(asci_settings.constraint_level, norb, n_sing_beta,
n_doub_beta, uniq_alpha_wfn, comm);
auto constraints = dist_constraint_general<wfn_t<N>>(
asci_settings.constraint_level, norb, n_sing_beta, n_doub_beta,
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());
Expand Down Expand Up @@ -453,65 +452,65 @@ asci_contrib_container<wfn_t<N>> asci_contributions_constraint(
} // Unique Alpha Loop
#else

for(size_t i_alpha = 0, iw = 0; i_alpha < nuniq_alpha; ++i_alpha) {
const auto& alpha_det = uniq_alpha[i_alpha].first;
const auto occ_alpha = bits_to_indices(alpha_det);
const bool alpha_satisfies_con = satisfies_constraint(alpha_det, con);

const auto& bcd = uad[i_alpha];
const size_t nbeta = bcd.size();
for(size_t j_beta = 0; j_beta < nbeta; ++j_beta, ++iw) {
const auto w = *(cdets_begin + iw);
const auto c = C[iw];
const auto& beta_det = bcd[j_beta].beta_string;
const auto h_diag = bcd[j_beta].h_diag;
const auto& occ_beta = bcd[j_beta].occ_beta;
const auto& vir_beta = bcd[j_beta].vir_beta;
const auto& orb_ens_alpha = bcd[j_beta].orb_ens_alpha;
const auto& orb_ens_beta = bcd[j_beta].orb_ens_beta;

// AA excitations
generate_constraint_singles_contributions_ss(
c, w, 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);
for(size_t i_alpha = 0, iw = 0; i_alpha < nuniq_alpha; ++i_alpha) {
const auto& alpha_det = uniq_alpha[i_alpha].first;
const auto occ_alpha = bits_to_indices(alpha_det);
const bool alpha_satisfies_con = satisfies_constraint(alpha_det, con);

const auto& bcd = uad[i_alpha];
const size_t nbeta = bcd.size();
for(size_t j_beta = 0; j_beta < nbeta; ++j_beta, ++iw) {
const auto w = *(cdets_begin + iw);
const auto c = C[iw];
const auto& beta_det = bcd[j_beta].beta_string;
const auto h_diag = bcd[j_beta].h_diag;
const auto& occ_beta = bcd[j_beta].occ_beta;
const auto& vir_beta = bcd[j_beta].vir_beta;
const auto& orb_ens_alpha = bcd[j_beta].orb_ens_alpha;
const auto& orb_ens_beta = bcd[j_beta].orb_ens_beta;

// AA excitations
generate_constraint_singles_contributions_ss(
c, w, 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);

// AAAA excitations
generate_constraint_doubles_contributions_ss(
c, w, con, occ_alpha, occ_beta, orb_ens_alpha.data(),
G_pqrs, norb, h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs);
// AAAA excitations
generate_constraint_doubles_contributions_ss(
c, w, con, occ_alpha, occ_beta, orb_ens_alpha.data(), G_pqrs, norb,
h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs);

// AABB excitations
generate_constraint_doubles_contributions_os(
c, w, 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);
// AABB excitations
generate_constraint_doubles_contributions_os(
c, w, 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(alpha_satisfies_con) {
// BB excitations
append_singles_asci_contributions<Spin::Beta>(
c, w, beta_det, occ_beta, vir_beta, occ_alpha,
orb_ens_beta.data(), T_pq, norb, G_red, norb, V_red, norb,
h_el_tol, h_diag, E_ASCI, ham_gen, asci_pairs);

if(alpha_satisfies_con) {
// BB excitations
append_singles_asci_contributions<Spin::Beta>(
c, w, beta_det, occ_beta, vir_beta, occ_alpha,
orb_ens_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<Spin::Beta>(
c, w, beta_det, alpha_det, occ_beta, vir_beta,
occ_alpha, orb_ens_beta.data(), G_pqrs, norb, h_el_tol, h_diag,
E_ASCI, ham_gen, asci_pairs);

// No excitation (push inf to remove from list)
asci_pairs.push_back(
{w, std::numeric_limits<double>::infinity(), 1.0});
}
// BBBB excitations
append_ss_doubles_asci_contributions<Spin::Beta>(
c, w, beta_det, alpha_det, occ_beta, vir_beta, occ_alpha,
orb_ens_beta.data(), G_pqrs, norb, h_el_tol, h_diag, E_ASCI,
ham_gen, asci_pairs);

// No excitation (push inf to remove from list)
asci_pairs.push_back(
{w, std::numeric_limits<double>::infinity(), 1.0});
}
}

// Prune Down Contributions
if(asci_pairs.size() > asci_settings.pair_size_max) {
throw std::runtime_error("DIE DIE DIE");
}

} // Unique Alpha Loop
} // Unique Alpha Loop

#endif

Expand Down Expand Up @@ -700,9 +699,9 @@ std::vector<wfn_t<N>> asci_search(
// Remove core dets
// XXX: This assumes the constraint search for now
{
auto inf_ptr = std::partition(asci_pairs.begin(), asci_pairs.end(),
[](auto& p){ return !std::isinf(p.rv()); });
asci_pairs.erase(inf_ptr, asci_pairs.end());
auto inf_ptr = std::partition(asci_pairs.begin(), asci_pairs.end(),
[](auto& p) { return !std::isinf(p.rv()); });
asci_pairs.erase(inf_ptr, asci_pairs.end());
}
#endif

Expand Down
10 changes: 5 additions & 5 deletions include/macis/asci/determinant_sort.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,14 +40,15 @@ template <typename WfnIterator>
void reorder_ci_on_alpha(WfnIterator begin, WfnIterator end, double* C) {
using wfn_type = typename WfnIterator::value_type;
using wfn_traits = wavefunction_traits<wfn_type>;
using cmp_type = typename wfn_traits::spin_comparator;
const size_t ndets = std::distance(begin,end);
using cmp_type = typename wfn_traits::spin_comparator;
const size_t ndets = std::distance(begin, end);

cmp_type comparator{};
std::vector<size_t> idx(ndets);
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(),
[&](auto i, auto j) { return comparator(*(begin+i),*(begin+j)); });
std::sort(idx.begin(), idx.end(), [&](auto i, auto j) {
return comparator(*(begin + i), *(begin + j));
});

std::vector<double> reorder_C(ndets);
std::vector<wfn_type> reorder_dets(ndets);
Expand All @@ -58,7 +59,6 @@ void reorder_ci_on_alpha(WfnIterator begin, WfnIterator end, double* C) {

std::copy(reorder_dets.begin(), reorder_dets.end(), begin);
std::copy(reorder_C.begin(), reorder_C.end(), C);

}

template <typename PairIterator>
Expand Down
6 changes: 3 additions & 3 deletions include/macis/asci/iteration.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,17 @@ auto asci_iter(ASCISettings asci_settings, MCSCFSettings mcscf_settings,
size_t nkeep = std::min(asci_settings.ncdets_max, wfn.size());

// Sort kept dets on alpha string
if(wfn.size() > 1)
if(wfn.size() > 1)
reorder_ci_on_alpha(wfn.begin(), wfn.begin() + nkeep, X.data());

// Perform the ASCI search
wfn = asci_search(asci_settings, ndets_max, wfn.begin(), wfn.begin() + nkeep,
E0, X, norb, ham_gen.T(), ham_gen.G_red(), ham_gen.V_red(),
ham_gen.G(), ham_gen.V(), ham_gen MACIS_MPI_CODE(, comm));

//std::sort(wfn.begin(), wfn.end(), bitset_less_comparator<N>{});
// std::sort(wfn.begin(), wfn.end(), bitset_less_comparator<N>{});
using wfn_traits = wavefunction_traits<wfn_t<N>>;
using wfn_comp = typename wfn_traits::spin_comparator;
using wfn_comp = typename wfn_traits::spin_comparator;
std::sort(wfn.begin(), wfn.end(), wfn_comp{});

// Rediagonalize
Expand Down
32 changes: 13 additions & 19 deletions include/macis/asci/mask_constraints.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -629,18 +629,15 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr,
#else
template <typename WfnType, typename ContainerType>
auto gen_constraints_general(size_t nlevels, size_t norb, size_t ns_othr,
size_t nd_othr,
const ContainerType& unique_alpha,
size_t nd_othr, const ContainerType& unique_alpha,
int world_size) {

using wfn_traits = wavefunction_traits<WfnType>;
using constraint_type = alpha_constraint<wfn_traits>;
using string_type = typename constraint_type::constraint_type;

constexpr bool flat_container = std::is_same_v<
std::decay_t<WfnType>,
std::decay_t<typename ContainerType::value_type>
>;
constexpr bool flat_container =
std::is_same_v<std::decay_t<WfnType>,
std::decay_t<typename ContainerType::value_type>>;

// Generate triplets + heuristic
std::vector<std::pair<constraint_type, size_t>> constraint_sizes;
Expand All @@ -653,12 +650,12 @@ auto gen_constraints_general(size_t nlevels, size_t norb, size_t ns_othr,

size_t nw = 0;
for(const auto& alpha : unique_alpha) {
if constexpr (flat_container)
if constexpr(flat_container)
nw += constraint_histogram(wfn_traits::alpha_string(alpha), ns_othr,
nd_othr, constraint);
else
nw += alpha.second *
constraint_histogram(alpha.first, ns_othr, nd_othr, constraint);
nw += alpha.second * constraint_histogram(alpha.first, ns_othr,
nd_othr, constraint);
}
if(nw) constraint_sizes.emplace_back(constraint, nw);
total_work += nw;
Expand Down Expand Up @@ -697,12 +694,12 @@ auto gen_constraints_general(size_t nlevels, size_t norb, size_t ns_othr,
size_t nw = 0;

for(const auto& alpha : unique_alpha) {
if constexpr (flat_container)
if constexpr(flat_container)
nw += constraint_histogram(wfn_traits::alpha_string(alpha), ns_othr,
nd_othr, c_next);
else
nw += alpha.second *
constraint_histogram(alpha.first, ns_othr, nd_othr, c_next);
nw += alpha.second *
constraint_histogram(alpha.first, ns_othr, nd_othr, c_next);
}
if(nw) constraint_sizes.emplace_back(c_next, nw);
total_work += nw;
Expand All @@ -719,19 +716,17 @@ auto gen_constraints_general(size_t nlevels, size_t norb, size_t ns_othr,

template <typename WfnType, typename ContainerType>
auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr,
size_t nd_othr,
const ContainerType& unique_alpha,
size_t nd_othr, const ContainerType& unique_alpha,
MPI_Comm comm) {

using wfn_traits = wavefunction_traits<WfnType>;
using constraint_type = alpha_constraint<wfn_traits>;

auto world_rank = comm_rank(comm);
auto world_size = comm_size(comm);

// Generate constraints subject to expected workload
auto constraint_sizes = gen_constraints_general<WfnType>(nlevels, norb, ns_othr,
nd_othr, unique_alpha, world_size);
auto constraint_sizes = gen_constraints_general<WfnType>(
nlevels, norb, ns_othr, nd_othr, unique_alpha, world_size);

// Global workloads
std::vector<size_t> workloads(world_size, 0);
Expand All @@ -757,7 +752,6 @@ auto dist_constraint_general(size_t nlevels, size_t norb, size_t ns_othr,
// workloads[world_rank], total_work);

return constraints;

}
#endif
#endif
Expand Down
Loading

0 comments on commit ddb59ed

Please sign in to comment.