clusters seem to work (if cacheing disabled) [skip ci]

CMakeLists.txt

@@ -46,6 +46,7 @@ cmaize_find_or_build_dependency(
     BUILD_TARGET simde
     FIND_TARGET nwx::simde
     CMAKE_ARGS BUILD_TESTING=OFF
+               NWX_PLUGINPLAY_VERSION="cache_access"
 )
 
 cmaize_find_or_build_dependency(
@@ -88,11 +89,6 @@ if("${BUILD_TESTING}")
         INCLUDE_DIRS "${project_src_dir}/${PROJECT_NAME}"
         DEPENDS Catch2::Catch2 ${PROJECT_NAME}
     )
-    nwx_pybind11_tests(
-        py_ghostfragment
-        ${python_tests_dir}/integration_tests/test_ghostfragment.py
-        SUBMODULES simde chemist pluginplay parallelzone
-    )
 
     if("${INTEGRATION_TESTING}")
         cmaize_find_or_build_dependency(
@@ -102,5 +98,13 @@ if("${BUILD_TESTING}")
             FIND_TARGET nwx::nwchemex
             CMAKE_ARGS BUILD_TESTING=OFF
         )
+
+        nwx_pybind11_tests(
+            py_ghostfragment
+            ${python_tests_dir}/integration_tests/test_ghostfragment.py
+            DEPENDS nwchemex
+            SUBMODULES simde chemist pluginplay parallelzone friendzone chemcache nwchemex
+        )
+
     endif()
 endif()

src/ghostfragment/drivers/drivers.hpp

@@ -40,11 +40,24 @@ inline void set_defaults(pluginplay::ModuleManager& mm) {
     mm.change_submod("Fragment Driver", "Find broken bonds", "Broken bonds");
     mm.change_submod("Fragment Driver", "Cap broken bonds",
                      "Weighted distance");
+
+    mm.copy_module("Fragment Driver", "N-mer Driver");
+    mm.change_submod("N-mer Driver", "Molecular graph to fragments",
+                     "All nmers");
+
     mm.change_submod("FragmentedChemicalSystem Driver", "Fragmenter",
                      "Fragment Driver");
+
+    mm.copy_module("FragmentedChemicalSystem Driver", "N-mer System Driver");
+    mm.change_submod("N-mer System Driver", "Fragmenter", "N-mer Driver");
+
     mm.change_submod("Fragment Based Method", "Subsystem former",
                      "FragmentedChemicalSystem Driver");
     mm.change_submod("Fragment Based Method", "Weighter", "GMBE Weights");
+
+    mm.copy_module("Fragment Based Method", "N-mer Based Method");
+    mm.change_submod("N-mer Based Method", "Subsystem former",
+                     "N-mer System Driver");
 }
 
 } // namespace ghostfragment::drivers

src/ghostfragment/drivers/fragment.cpp

@@ -58,22 +58,42 @@ MODULE_CTOR(Fragment) {
 
 MODULE_RUN(Fragment) {
     const auto& [mol] = frags_pt::unwrap_inputs(inputs);
-
-    auto& graph_mod   = submods.at("Molecular Graph");
-    const auto& graph = graph_mod.run_as<graph_pt>(mol);
-
-    auto& frags_mod   = submods.at("Molecular graph to fragments");
+    auto& runtime     = get_runtime();
+
+    // Step 1: Form the molecular graph
+    auto& graph_mod    = submods.at("Molecular Graph");
+    const auto& graph  = graph_mod.run_as<graph_pt>(mol);
+    const auto n_nodes = graph.nodes_size();
+    const auto n_edges = graph.edges_size();
+    runtime.logger().debug("Created a graph with " + std::to_string(n_nodes) +
+                           " nodes and " + std::to_string(n_edges) + " edges.");
+
+    // Step 2: Use the graph to make fragments
+    auto& frags_mod           = submods.at("Molecular graph to fragments");
     const auto& frags_no_ints = frags_mod.run_as<graph2frags_pt>(graph);
+    const auto n_frags        = frags_no_ints.size();
+    runtime.logger().debug("Created " + std::to_string(n_frags) +
+                           " fragments.");
 
+    // Step 3: Analyze the fragments for intersections
     auto& intersect_mod = submods.at("Intersection finder");
-    const auto& frags  = intersect_mod.run_as<intersections_pt>(frags_no_ints);
+    const auto& frags   = intersect_mod.run_as<intersections_pt>(frags_no_ints);
+    const auto n_ints   = frags.size() - n_frags;
+    runtime.logger().debug("Added " + std::to_string(n_ints) +
+                           " intersections.");
 
+    // Step 4: Did forming fragments (or intersections) break bonds?
     auto& bonds_mod          = submods.at("Find broken bonds");
     const auto& conns        = graph.edges();
     const auto& broken_bonds = bonds_mod.run_as<broken_bonds_pt>(frags, conns);
+    runtime.logger().debug("Found " + std::to_string(broken_bonds.size()) +
+                           " broken bonds.");
 
+    // Step 5: Fix those broken bonds!!!!
     auto& cap_mod            = submods.at("Cap broken bonds");
     const auto& capped_frags = cap_mod.run_as<cap_pt>(frags, broken_bonds);
+    const auto n_caps        = capped_frags.cap_set().size();
+    runtime.logger().debug("Added " + std::to_string(n_caps) + " caps.");
 
     auto rv = results();
     return frags_pt::wrap_results(rv, capped_frags);

src/ghostfragment/drivers/fragment_based_method.cpp

@@ -44,11 +44,12 @@ MODULE_CTOR(FragmentBasedMethod) {
 
     add_submodule<fragmenting_pt>("Subsystem former");
     add_submodule<weight_pt>("Weighter");
-    add_submodule<basis_set_pt>("Apply basis set");
-    add_submodule<ao_energy_pt>("Energy method");
+    add_submodule<my_pt>("Energy method");
 }
 
 MODULE_RUN(FragmentBasedMethod) {
+    auto& logger = get_runtime().logger();
+
     // Step 0: Unpack input
     const auto& [sys] = my_pt::unwrap_inputs(inputs);
 
@@ -60,19 +61,25 @@ MODULE_RUN(FragmentBasedMethod) {
     auto& weight_mod    = submods.at("weighter");
     const auto& weights = weight_mod.run_as<weight_pt>(subsystems);
 
-    auto& basis_mod  = submods.at("Apply basis set");
     auto& energy_mod = submods.at("Energy method");
 
     double energy = 0.0;
+
+    auto n_subsystems              = subsystems.size();
+    decltype(n_subsystems) counter = 0;
+    auto msg = [](auto counter, auto n_subsystems, auto egy) {
+        return "Energy of subsystem " + std::to_string(counter) + " of " +
+               std::to_string(n_subsystems) + " : " + std::to_string(egy);
+    };
     for(auto&& [c_i, sys_i] : iter::zip(weights, subsystems)) {
-        auto mol_i        = sys_i.molecule().as_molecule();
-        const auto& basis = basis_mod.run_as<basis_set_pt>(mol_i);
+        auto mol_i = sys_i.molecule().as_molecule();
 
         // This is a hack until views work with values
         chemical_system_type sys_i_copy(mol_i);
 
-        const auto e_i = energy_mod.run_as<ao_energy_pt>(basis, sys_i_copy);
+        const auto e_i = energy_mod.run_as<my_pt>(sys_i_copy);
         energy += c_i * e_i;
+        logger.info(msg(counter++, n_subsystems, e_i));
     }
 
     auto rv = results();

src/ghostfragment/fragmenting/fragmenting.hpp

@@ -37,6 +37,7 @@ inline void load_modules(pluginplay::ModuleManager& mm) {
 
 inline void set_defaults(pluginplay::ModuleManager& mm) {
     mm.change_submod("Heavy Atom Partition", "Connectivity", "Covalent Radius");
+    mm.change_submod("All nmers", "Monomer maker", "Bond-Based Fragmenter");
 }
 
 } // namespace ghostfragment::fragmenting

src/ghostfragment/fragmenting/heavy_atom.cpp

@@ -56,19 +56,23 @@ MODULE_CTOR(HeavyAtom) {
 MODULE_RUN(HeavyAtom) {
     using fragmented_nuclei = typename pt::FragmentedNucleiTraits::result_type;
     using size_type         = typename fragmented_nuclei::size_type;
+    auto& logger            = get_runtime().logger();
 
     const auto& [system] = frags_pt::unwrap_inputs(inputs);
     const auto& mol      = system.molecule();
 
     auto& con_mod     = submods.at("Connectivity");
     const auto& conns = con_mod.run_as<conn_pt>(mol);
+    logger.debug("Found " + std::to_string(conns.nbonds()) + " bonds.");
 
     fragmented_nuclei frags(mol.nuclei().as_nuclei());
 
     for(size_type atom_i = 0; atom_i < mol.size(); ++atom_i) {
         std::vector<size_type> fragment;
         const auto Zi     = mol[atom_i].Z();
         const auto conn_i = conns.bonded_atoms(atom_i);
+        logger.trace("Atom " + std::to_string(atom_i) +
+                     " has Z == " + std::to_string(Zi));
         if(Zi > 1) {
             fragment.push_back(atom_i);
 

src/ghostfragment/fragmenting/intersections_by_recursion.cpp

@@ -26,27 +26,27 @@ using nuclear_index_set = typename fragments_type::nucleus_index_set;
 using size_type         = typename nuclear_index_set::size_type;
 using index_set         = std::set<size_type>;
 using intersection_set  = std::set<index_set>;
+using frag_set          = std::vector<index_set>;
 
 namespace {
-template<typename BeginItr, typename EndItr>
-void compute_intersection(const index_set& curr_frag, BeginItr&& starting_frag,
-                          EndItr&& end_itr, intersection_set& ints_so_far) {
-    while(starting_frag != end_itr) {
+void compute_intersection(const index_set& curr_frag, std::size_t starting_frag,
+                          const frag_set& frag_indices,
+                          intersection_set& ints_so_far) {
+    while(starting_frag < frag_indices.size()) {
         index_set intersection;
+        const index_set& next_frag = frag_indices[starting_frag];
         auto itr = std::inserter(intersection, intersection.begin());
         std::set_intersection(curr_frag.begin(), curr_frag.end(),
-                              starting_frag->begin(), starting_frag->end(),
-                              itr);
+                              next_frag.begin(), next_frag.end(), itr);
         ++starting_frag;
         // If it's empty and/or we've seen it befor just move on
         if(intersection.empty() || ints_so_far.count(intersection)) continue;
 
         // Add the intersection
         ints_so_far.insert(intersection);
 
-        // Copy iterator so we leave starting_frag unchanged when we recurse
-        decltype(starting_frag) new_begin(starting_frag);
-        compute_intersection(intersection, new_begin, end_itr, ints_so_far);
+        compute_intersection(intersection, starting_frag + 1, frag_indices,
+                             ints_so_far);
     }
 }
 
@@ -65,30 +65,35 @@ MODULE_CTOR(IntersectionsByRecursion) {
 }
 
 MODULE_RUN(IntersectionsByRecursion) {
+    auto& logger = get_runtime().logger();
     auto [frags] = property_type::unwrap_inputs(inputs);
 
     // It's much easier to work with nuclear indices
     std::vector<index_set> frag_indices;
 
     for(size_type i = 0; i < frags.size(); ++i) {
         const auto frag_i = frags.nuclear_indices(i);
+        std::string frag_str;
+        for(auto x : frag_i) frag_str += std::to_string(x) + ",";
+        frag_str.pop_back();
+        logger.debug("Input fragment: " + frag_str);
         frag_indices.emplace_back(frag_i.begin(), frag_i.end());
     }
 
     std::set<index_set> intersections;
 
-    auto begin = frag_indices.begin();
-    auto end   = frag_indices.end();
-
-    while(begin != end) {
-        index_set frag = *begin;
-        ++begin;
-        decltype(begin) begin_copy(begin);
-        compute_intersection(frag, begin_copy, end, intersections);
+    for(decltype(frags.size()) begin = 0; begin < frags.size(); ++begin) {
+        const index_set& frag = frag_indices[begin];
+        compute_intersection(frag, begin + 1, frag_indices, intersections);
     }
 
-    for(const auto& intersection_i : intersections)
+    for(const auto& intersection_i : intersections) {
+        std::string int_str;
+        for(auto x : intersection_i) int_str += std::to_string(x) + ",";
+        int_str.pop_back();
+        logger.debug("Found intersection: " + int_str);
         frags.insert(intersection_i.begin(), intersection_i.end());
+    }
 
     auto rv = results();
     return property_type::wrap_results(rv, frags);

src/ghostfragment/fragmenting/nmers.cpp

@@ -45,14 +45,21 @@ MODULE_CTOR(NMers) {
     description(mod_desc);
     satisfies_property_type<my_pt>();
 
-    add_input<n_type>("n").set_description("The maximum n-mer size");
+    add_input<n_type>("n")
+      .set_description("The maximum n-mer size")
+      .set_default(n_type(1));
     add_submodule<my_pt>("Monomer maker");
 }
 
 MODULE_RUN(NMers) {
+    auto& logger = get_runtime().logger();
+
     const auto& [graph] = my_pt::unwrap_inputs(inputs);
     auto n              = inputs.at("n").value<n_type>();
 
+    auto nmer_str = std::to_string(n) + "-mers";
+    logger.debug("Will be making " + nmer_str + ".");
+
     auto& monomer_mod = submods.at("Monomer maker");
     const auto& frags = monomer_mod.run_as<my_pt>(graph);
 
@@ -110,7 +117,7 @@ MODULE_RUN(NMers) {
         }
         if(i_is_good[i]) nmers.insert(nmer_i.begin(), nmer_i.end());
     }
-
+    logger.debug("Made " + std::to_string(nmers.size()) + " " + nmer_str + ".");
     auto rv = results();
     return my_pt::wrap_results(rv, nmers);
 }

src/ghostfragment/topology/covradii.cpp

@@ -55,6 +55,7 @@ MODULE_CTOR(CovRadii) {
 }
 
 MODULE_RUN(CovRadii) {
+    auto& logger          = get_runtime().logger();
     const auto& [mol]     = my_pt::unwrap_inputs(inputs);
     const auto tau        = inputs.at("tau").value<double>();
     const auto tau_plus_1 = tau + 1.0;
@@ -66,13 +67,21 @@ MODULE_RUN(CovRadii) {
     for(size_type i = 0; i < natoms; ++i) {
         const auto atom_i  = mol[i];
         const auto sigma_i = covalent_radius(atom_i.Z());
+        logger.trace("Atom " + std::to_string(i) + " has covalent radius " +
+                     std::to_string(sigma_i) + " (a.u.).");
 
         for(size_type j = i + 1; j < natoms; ++j) {
             const auto atom_j  = mol[j].as_nucleus();
             const auto sigma_j = covalent_radius(atom_j.Z());
+            logger.trace("Atom " + std::to_string(j) + " has covalent radius " +
+                         std::to_string(sigma_j) + " (a.u.).");
 
             const auto rij      = (atom_i.as_nucleus() - atom_j).magnitude();
             const auto max_bond = tau_plus_1 * (sigma_i + sigma_j);
+
+            logger.trace(std::to_string(i) + "-" + std::to_string(j) +
+                         " distance is: " + std::to_string(rij));
+
             if(rij <= max_bond) ct.add_bond(i, j);
         }
     }

src/ghostfragment/topology/nuclear_graph_from_connectivity.cpp

@@ -46,13 +46,22 @@ MODULE_CTOR(NuclearGraphFromConnectivity) {
 MODULE_RUN(NuclearGraphFromConnectivity) {
     using traits_type = pt::NuclearGraphTraits;
     using result_type = traits_type::result_type;
+    auto& logger      = get_runtime().logger();
 
     const auto& [chem_sys] = my_pt::unwrap_inputs(inputs);
-    auto& pseudo_atom_mod  = submods.at("Nodes");
-    const auto& frags      = pseudo_atom_mod.run_as<pa_pt>(chem_sys);
+
+    auto& pseudo_atom_mod = submods.at("Nodes");
+    const auto& frags     = pseudo_atom_mod.run_as<pa_pt>(chem_sys);
+    const auto n_atoms    = chem_sys.molecule().size();
+    const auto n_pas      = frags.size();
+    logger.debug("The " + std::to_string(n_atoms) +
+                 " atoms of the system were converted into " +
+                 std::to_string(n_pas) + " pseudoatoms.");
 
     auto& conn_mod         = submods.at("Connectivity");
     const auto& atom_conns = conn_mod.run_as<conn_pt>(chem_sys.molecule());
+    const auto n_bonds     = atom_conns.nbonds();
+    logger.debug("System has " + std::to_string(n_bonds) + " bonds.");
 
     const auto nnodes = frags.size();
     std::decay_t<decltype(atom_conns)> edges(nnodes);