Merge branch 'python' into ekin_immutable

samples/grand_canonical.py

@@ -97,6 +97,9 @@
 print(RE.get_status())
 system.setup_type_map([0, 1, 2])
 
+# Set the hidden particle type to the lowest possible number to speed
+# up the simulation
+RE.set_non_interacting_type(max(types) + 1)
 
 RE.reaction(10000)
 

samples/reaction_ensemble.py

@@ -57,9 +57,16 @@
 
 # Particle setup
 #############################################################
-# type 0 = HA
-# type 1 = A-
-# type 2 = H+
+types = {
+    "HA": 0,
+    "A-": 1,
+    "H+": 2,
+}
+charge_dict = {
+    0: 0,
+    1: -1,
+    2: +1,
+}
 
 N0 = 50  # number of titratable units
 K_diss = 0.0088
@@ -75,28 +82,36 @@
         kT=1,
         exclusion_radius=1,
         seed=77)
-    RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1],
-                    product_types=[1, 2], product_coefficients=[1, 1],
-                    default_charges={0: 0, 1: -1, 2: +1})
+    RE.add_reaction(gamma=K_diss,
+                    reactant_types=[types["HA"]],
+                    reactant_coefficients=[1],
+                    product_types=[types["A-"], types["H+"]],
+                    product_coefficients=[1, 1],
+                    default_charges=charge_dict)
 elif args.mode == "constant_pH_ensemble":
     RE = espressomd.reaction_ensemble.ConstantpHEnsemble(
         kT=1, exclusion_radius=1, seed=77)
     RE.constant_pH = 2
-    RE.add_reaction(gamma=K_diss, reactant_types=[0],
-                    product_types=[1, 2],
-                    default_charges={0: 0, 1: -1, 2: +1})
+    RE.add_reaction(gamma=K_diss, reactant_types=[types["HA"]],
+                    product_types=[types["A-"], types["H+"]],
+                    default_charges=charge_dict)
 else:
     raise RuntimeError(
         "Please provide either --reaction_ensemble or --constant_pH_ensemble as argument ")
 
 print(RE.get_status())
-system.setup_type_map([0, 1, 2])
+system.setup_type_map(list(types.values()))
+
+
+# Set the hidden particle type to the lowest possible number_of_particles
+# to speed up the simulation
+RE.set_non_interacting_type(max(types.values()) + 1)
 
 for i in range(10000):
     RE.reaction()
     if i % 100 == 0:
-        print("HA", system.number_of_particles(type=0), "A-",
-              system.number_of_particles(type=1), "H+", system.number_of_particles(type=2))
+        print("HA", system.number_of_particles(type=types["HA"]), "A-",
+              system.number_of_particles(type=types["A-"]), "H+", system.number_of_particles(type=types["H+"]))
 
 print("reaction 0 has acceptance rate: ", RE.get_acceptance_rate_reaction(0))
 print("reaction 1 has acceptance rate: ", RE.get_acceptance_rate_reaction(1))

samples/reaction_ensemble_complex_reaction.py

@@ -95,6 +95,10 @@
 
 numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []}
 
+# Set the hidden particle type to the lowest possible number to speed 
+# up the simulation
+RE.set_non_interacting_type(max(types) + 1)
+
 # warmup
 RE.reaction(200)
 

samples/widom_insertion.py

@@ -113,6 +113,11 @@
 print(widom.get_status())
 system.setup_type_map([0, 1, 2])
 
+
+# Set the hidden particle type to the lowest possible number to speed 
+# up the simulation
+widom.set_non_interacting_type(max(types) + 1)
+
 n_iterations = 100
 for i in range(n_iterations):
     for j in range(50):
@@ -125,5 +130,6 @@
               f"A- {system.number_of_particles(type=1)}",
               f"H+ {system.number_of_particles(type=2)}")
 
+
 print("excess chemical potential for an ion pair ",
       widom.measure_excess_chemical_potential(insertion_reaction_id))

src/core/CMakeLists.txt

@@ -41,7 +41,6 @@ set(EspressoCore_SRC
     CellStructure.cpp
     PartCfg.cpp
     AtomDecomposition.cpp
-    reduce_observable_stat.cpp
     EspressoSystemStandAlone.cpp
     DomainDecomposition.cpp)
 

src/core/Observable_stat.cpp

@@ -24,32 +24,42 @@
 #include "config.hpp"
 
 #include "bonded_interactions/bonded_interaction_data.hpp"
+#include "communication.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include <utils/Span.hpp>
+#include <utils/index.hpp>
+
+#include <boost/mpi/collectives/reduce.hpp>
 
 #include <cassert>
 #include <cstddef>
+#include <functional>
 #include <vector>
 
+inline std::size_t max_non_bonded_pairs() {
+  auto const n = static_cast<std::size_t>(max_seen_particle_type);
+  return (n * (n + 1)) / 2;
+}
+
 Observable_stat::Observable_stat(std::size_t chunk_size)
     : m_chunk_size(chunk_size) {
   // number of chunks for different interaction types
-  auto constexpr n_coulomb = 2;
-  auto constexpr n_dipolar = 2;
+  constexpr std::size_t n_coulomb = 2;
+  constexpr std::size_t n_dipolar = 2;
 #ifdef VIRTUAL_SITES
-  auto constexpr n_vs = 1;
+  constexpr std::size_t n_vs = 1;
 #else
-  auto constexpr n_vs = 0;
+  constexpr std::size_t n_vs = 0;
 #endif
   auto const n_bonded = bonded_ia_params.size();
   auto const n_non_bonded = max_non_bonded_pairs();
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions
 
-  // resize vector
-  auto const total = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
-                     n_dipolar + n_vs + n_ext_fields;
-  m_data = std::vector<double>(m_chunk_size * total);
+  auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
+                          n_dipolar + n_vs + n_ext_fields;
+  m_data = std::vector<double>(m_chunk_size * n_elements);
 
   // spans for the different contributions
   kinetic = Utils::Span<double>(m_data.data(), m_chunk_size);
@@ -65,3 +75,20 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
       Utils::Span<double>(non_bonded_intra.end(), n_non_bonded * m_chunk_size);
   assert(non_bonded_inter.end() == (m_data.data() + m_data.size()));
 }
+
+Utils::Span<double>
+Observable_stat::non_bonded_contribution(Utils::Span<double> base_pointer,
+                                         int type1, int type2) const {
+  auto const offset = static_cast<std::size_t>(Utils::upper_triangular(
+      std::min(type1, type2), std::max(type1, type2), max_seen_particle_type));
+  return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size};
+}
+
+void Observable_stat::mpi_reduce() {
+  if (comm_cart.rank() == 0) {
+    std::vector<double> temp(m_data);
+    boost::mpi::reduce(comm_cart, temp, m_data, std::plus<>{}, 0);
+  } else {
+    boost::mpi::reduce(comm_cart, m_data, std::plus<>{}, 0);
+  }
+}

src/core/Observable_stat.hpp

@@ -23,7 +23,6 @@
 #include <boost/range/numeric.hpp>
 
 #include <utils/Span.hpp>
-#include <utils/index.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -38,32 +37,14 @@ class Observable_stat {
   /** Number of doubles per data item */
   std::size_t m_chunk_size;
 
-  /** Calculate the maximal number of non-bonded interaction pairs in the
-   *  system.
-   */
-  static std::size_t max_non_bonded_pairs() {
-    extern int max_seen_particle_type;
-    return static_cast<std::size_t>(
-        (max_seen_particle_type * (max_seen_particle_type + 1)) / 2);
-  }
-
   /** Get contribution from a non-bonded interaction */
   Utils::Span<double> non_bonded_contribution(Utils::Span<double> base_pointer,
-                                              int type1, int type2) const {
-    extern int max_seen_particle_type;
-    int offset = Utils::upper_triangular(
-        std::min(type1, type2), std::max(type1, type2), max_seen_particle_type);
-    return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size};
-  }
+                                              int type1, int type2) const;
 
 public:
   explicit Observable_stat(std::size_t chunk_size);
 
-  auto chunk_size() const { return m_chunk_size; }
-  Utils::Span<double> data_() { return {m_data.data(), m_data.size()}; }
-  Utils::Span<const double> data_() const {
-    return {m_data.data(), m_data.size()};
-  }
+  auto get_chunk_size() const { return m_chunk_size; }
 
   /** Accumulate values.
    *  @param acc    Initial value for the accumulator.
@@ -97,7 +78,7 @@ class Observable_stat {
   Utils::Span<double> coulomb;
   /** Contribution(s) from dipolar interactions. */
   Utils::Span<double> dipolar;
-  /** Contribution(s) from virtual sites (accumulated). */
+  /** Contribution from virtual sites (accumulated). */
   Utils::Span<double> virtual_sites;
   /** Contribution from external fields (accumulated). */
   Utils::Span<double> external_fields;
@@ -107,17 +88,16 @@ class Observable_stat {
   Utils::Span<double> non_bonded_inter;
 
   /** Get contribution from a bonded interaction */
-  Utils::Span<double> bonded_contribution(int bond_id) {
-    return Utils::Span<double>(bonded.data() + m_chunk_size * bond_id,
-                               m_chunk_size);
+  Utils::Span<double> bonded_contribution(int bond_id) const {
+    auto const offset = m_chunk_size * static_cast<std::size_t>(bond_id);
+    return Utils::Span<double>(bonded.data() + offset, m_chunk_size);
   }
 
   void add_non_bonded_contribution(int type1, int type2,
                                    Utils::Span<const double> data) {
-    auto const dest =
-        (type1 == type2)
-            ? non_bonded_contribution(non_bonded_intra, type1, type2)
-            : non_bonded_contribution(non_bonded_inter, type1, type2);
+    assert(data.size() == m_chunk_size);
+    auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter;
+    auto const dest = non_bonded_contribution(source, type1, type2);
 
     boost::transform(dest, data, dest.begin(), std::plus<>{});
   }
@@ -137,6 +117,9 @@ class Observable_stat {
                                                     int type2) const {
     return non_bonded_contribution(non_bonded_inter, type1, type2);
   }
+
+  /** MPI reduction. */
+  void mpi_reduce();
 };
 
 #endif // ESPRESSO_OBSERVABLE_STAT_HPP

src/core/constraints/ShapeBasedConstraint.cpp

@@ -20,6 +20,7 @@
 #include "ShapeBasedConstraint.hpp"
 
 #include "BoxGeometry.hpp"
+#include "Observable_stat.hpp"
 #include "communication.hpp"
 #include "config.hpp"
 #include "dpd.hpp"