Merge branch 'python' into tutorial_electrokinetics

doc/sphinx/installation.rst

@@ -8,10 +8,10 @@ This chapter will describe how to get, compile and run the software.
 |es| releases are available as source code packages from the homepage [1]_.
 This is where new users should get the code. The code within release packages
 is tested and known to run on a number of platforms.
-Alternatively, people that want to use the newest features of |es| or that
-want to start contributing to the software can instead obtain the
+Alternatively, people who want to use the newest features of |es| or
+start contributing to the software can instead obtain the
 current development code via the version control system software  [2]_
-from |es|'s project page at Github  [3]_. This code might be not as well
+from |es|'s project page at GitHub  [3]_. This code might be not as well
 tested and documented as the release code; it is recommended to use this
 code only if you have already gained some experience in using |es|.
 
@@ -93,7 +93,7 @@ To compile |es| on Ubuntu 22.04 LTS, install the following dependencies:
 .. code-block:: bash
 
     sudo apt install build-essential cmake cython3 python3-pip python3-numpy \
-      libboost-all-dev openmpi-common fftw3-dev libhdf5-dev libhdf5-openmpi-dev \
+      libboost-all-dev openmpi-common fftw3-dev libfftw3-mpi-dev libhdf5-dev libhdf5-openmpi-dev \
       python3-scipy python3-opengl libgsl-dev freeglut3
 
 Optionally the ccmake utility can be installed for easier configuration:
@@ -221,7 +221,7 @@ To use Jupyter Notebook, install the following packages:
 
 .. code-block:: bash
 
-    pip3 install --user nbformat notebook 'jupyter_contrib_nbextensions==0.5.1'
+    pip3 install --user 'nbformat==5.1.3' 'nbconvert==6.4.5' 'notebook==6.4.8' 'jupyter_contrib_nbextensions==0.5.1'
     jupyter contrib nbextension install --user
     jupyter nbextension enable rubberband/main
     jupyter nbextension enable exercise2/main

src/core/Observable_stat.cpp

@@ -58,8 +58,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions
 
-  auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
-                          n_dipolar + n_vs + n_ext_fields;
+  auto const n_elements = n_kinetic + n_bonded + 2ul * 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
@@ -78,8 +78,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
 }
 
 Utils::Span<double>
-Observable_stat::non_bonded_contribution(Utils::Span<double> base_pointer,
-                                         int type1, int type2) const {
+Observable_stat::get_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};

src/core/Observable_stat.hpp

@@ -38,8 +38,9 @@ class Observable_stat {
   std::size_t m_chunk_size;
 
   /** Get contribution from a non-bonded interaction */
-  Utils::Span<double> non_bonded_contribution(Utils::Span<double> base_pointer,
-                                              int type1, int type2) const;
+  Utils::Span<double>
+  get_non_bonded_contribution(Utils::Span<double> base_pointer, int type1,
+                              int type2) const;
 
 public:
   explicit Observable_stat(std::size_t chunk_size);
@@ -91,29 +92,30 @@ class Observable_stat {
     return {bonded.data() + offset, m_chunk_size};
   }
 
-  void add_non_bonded_contribution(int type1, int type2,
+  void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2,
                                    Utils::Span<const double> data) {
     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);
+    auto const span = (molid1 == molid2) ? non_bonded_intra : non_bonded_inter;
+    auto const dest = get_non_bonded_contribution(span, type1, type2);
 
     boost::transform(dest, data, dest.begin(), std::plus<>{});
   }
 
-  void add_non_bonded_contribution(int type1, int type2, double data) {
-    add_non_bonded_contribution(type1, type2, {&data, 1});
+  void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2,
+                                   double data) {
+    add_non_bonded_contribution(type1, type2, molid1, molid2, {&data, 1});
   }
 
   /** Get contribution from a non-bonded intramolecular interaction */
   Utils::Span<double> non_bonded_intra_contribution(int type1,
                                                     int type2) const {
-    return non_bonded_contribution(non_bonded_intra, type1, type2);
+    return get_non_bonded_contribution(non_bonded_intra, type1, type2);
   }
 
   /** Get contribution from a non-bonded intermolecular interaction */
   Utils::Span<double> non_bonded_inter_contribution(int type1,
                                                     int type2) const {
-    return non_bonded_contribution(non_bonded_inter, type1, type2);
+    return get_non_bonded_contribution(non_bonded_inter, type1, type2);
   }
 
   /** MPI reduction. */

src/core/accumulators.cpp

@@ -60,24 +60,33 @@ int auto_update_next_update() {
                            });
 }
 
+namespace detail {
+struct MatchPredicate {
+  AccumulatorBase const *m_acc;
+  template <typename T> bool operator()(T const &a) const {
+    return a.acc == m_acc;
+  }
+};
+} // namespace detail
+
 void auto_update_add(AccumulatorBase *acc) {
-  assert(acc);
-  assert(std::find_if(auto_update_accumulators.begin(),
-                      auto_update_accumulators.end(), [acc](auto const &a) {
-                        return a.acc == acc;
-                      }) == auto_update_accumulators.end());
+  assert(not auto_update_contains(acc));
   auto_update_accumulators.emplace_back(acc);
 }
+
 void auto_update_remove(AccumulatorBase *acc) {
-  assert(std::find_if(auto_update_accumulators.begin(),
-                      auto_update_accumulators.end(), [acc](auto const &a) {
-                        return a.acc == acc;
-                      }) != auto_update_accumulators.end());
+  assert(auto_update_contains(acc));
+  auto const beg = auto_update_accumulators.begin();
+  auto const end = auto_update_accumulators.end();
   auto_update_accumulators.erase(
-      boost::remove_if(
-          auto_update_accumulators,
-          [acc](AutoUpdateAccumulator const &au) { return au.acc == acc; }),
-      auto_update_accumulators.end());
+      std::remove_if(beg, end, detail::MatchPredicate{acc}), end);
+}
+
+bool auto_update_contains(AccumulatorBase const *acc) noexcept {
+  assert(acc);
+  auto const beg = auto_update_accumulators.begin();
+  auto const end = auto_update_accumulators.end();
+  return std::find_if(beg, end, detail::MatchPredicate{acc}) != end;
 }
 
 } // namespace Accumulators

src/core/accumulators.hpp

@@ -31,6 +31,7 @@ namespace Accumulators {
  */
 void auto_update(boost::mpi::communicator const &comm, int steps);
 int auto_update_next_update();
+bool auto_update_contains(AccumulatorBase const *) noexcept;
 void auto_update_add(AccumulatorBase *);
 void auto_update_remove(AccumulatorBase *);
 

src/core/constraints/Constraints.hpp

@@ -49,22 +49,20 @@ template <class ParticleRange, class Constraint> class Constraints {
   container_type m_constraints;
 
 public:
+  bool contains(std::shared_ptr<Constraint> const &constraint) const noexcept {
+    return std::find(begin(), end(), constraint) != end();
+  }
   void add(std::shared_ptr<Constraint> const &constraint) {
     if (not constraint->fits_in_box(box_geo.length())) {
       throw std::runtime_error("Constraint not compatible with box size.");
     }
-    assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) ==
-           m_constraints.end());
-
+    assert(not contains(constraint));
     m_constraints.emplace_back(constraint);
     on_constraint_change();
   }
   void remove(std::shared_ptr<Constraint> const &constraint) {
-    assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) !=
-           m_constraints.end());
-    m_constraints.erase(
-        std::remove(m_constraints.begin(), m_constraints.end(), constraint),
-        m_constraints.end());
+    assert(contains(constraint));
+    m_constraints.erase(std::remove(begin(), end(), constraint), end());
     on_constraint_change();
   }
 
@@ -102,7 +100,7 @@ template <class ParticleRange, class Constraint> class Constraints {
   }
 
   void on_boxl_change() const {
-    if (not this->empty()) {
+    if (not m_constraints.empty()) {
       throw std::runtime_error("The box size can not be changed because there "
                                "are active constraints.");
     }

src/core/constraints/ShapeBasedConstraint.cpp

@@ -164,7 +164,9 @@ void ShapeBasedConstraint::add_energy(const Particle &p,
     }
   }
   // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
-  if (part_rep.type() >= 0)
-    obs_energy.add_non_bonded_contribution(p.type(), part_rep.type(), energy);
+  if (part_rep.type() >= 0) {
+    obs_energy.add_non_bonded_contribution(
+        p.type(), part_rep.type(), p.mol_id(), part_rep.mol_id(), energy);
+  }
 }
 } // namespace Constraints

src/core/energy_inline.hpp

@@ -183,7 +183,7 @@ inline void add_non_bonded_pair_energy(
   if (do_nonbonded(p1, p2))
 #endif
     obs_energy.add_non_bonded_contribution(
-        p1.type(), p2.type(),
+        p1.type(), p2.type(), p1.mol_id(), p2.mol_id(),
         calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist,
                                     coulomb_kernel));
 

src/core/event.cpp

@@ -31,6 +31,7 @@
 #include "collision.hpp"
 #include "communication.hpp"
 #include "config/config.hpp"
+#include "constraints.hpp"
 #include "cuda/init.hpp"
 #include "cuda/utils.hpp"
 #include "electrostatics/coulomb.hpp"
@@ -159,6 +160,12 @@ void on_particle_charge_change() {
   partCfg().invalidate();
 }
 
+void on_particle_local_change() {
+  cells_update_ghosts(global_ghost_flags());
+
+  recalc_forces = true;
+}
+
 void on_particle_change() {
   if (cell_structure.decomposition_type() ==
       CellStructureType::CELL_STRUCTURE_HYBRID) {
@@ -241,6 +248,8 @@ void on_boxl_change(bool skip_method_adaption) {
     system.dipoles.on_boxl_change();
 #endif
   }
+
+  Constraints::constraints.on_boxl_change();
 }
 
 void on_cell_structure_change() {

src/core/event.hpp

@@ -66,6 +66,9 @@ void on_particle_change();
 /** called every time the charge of a particle has changed. */
 void on_particle_charge_change();
 
+/** called every time that local properties of a particle have changed. */
+void on_particle_local_change();
+
 /** called every time the Coulomb parameters are changed.
 
 all Coulomb methods have a short range part, aka near field

src/core/pressure_inline.hpp

@@ -70,10 +70,8 @@ inline void add_non_bonded_pair_virials(
                            .f +
                        calc_non_central_force(p1, p2, ia_params, d, dist).f;
     auto const stress = Utils::tensor_product(d, force);
-
-    auto const type1 = p1.mol_id();
-    auto const type2 = p2.mol_id();
-    obs_pressure.add_non_bonded_contribution(type1, type2, flatten(stress));
+    obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(),
+                                             p2.mol_id(), flatten(stress));
   }
 
 #ifdef ELECTROSTATICS

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -149,13 +149,19 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
     auto const random_index = i_random(number_of_particles_with_type(type));
     return get_random_p_id(type, random_index);
   };
+  auto only_local_changes = true;
   for (int i = 0; i < std::min(n_product_types, n_reactant_types); i++) {
     auto const n_product_coef = reaction.product_coefficients[i];
     auto const n_reactant_coef = reaction.reactant_coefficients[i];
     // change std::min(reactant_coefficients(i),product_coefficients(i)) many
     // particles of reactant_types(i) to product_types(i)
     auto const old_type = reaction.reactant_types[i];
     auto const new_type = reaction.product_types[i];
+#ifdef ELECTROSTATICS
+    if (charges_of_types.at(new_type) != charges_of_types.at(old_type)) {
+      only_local_changes = false;
+    }
+#endif
     for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) {
       auto const p_id = get_random_p_id_of_type(old_type);
       on_particle_type_change(p_id, old_type, new_type);
@@ -167,7 +173,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
       }
       bookkeeping.changed.emplace_back(p_id, old_type);
     }
-    on_particle_change();
     // create product_coefficients(i)-reactant_coefficients(i) many product
     // particles iff product_coefficients(i)-reactant_coefficients(i)>0,
     // iff product_coefficients(i)-reactant_coefficients(i)<0, hide this number
@@ -180,7 +185,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
         check_exclusion_range(p_id, type);
         bookkeeping.created.emplace_back(p_id);
       }
-      on_particle_change();
+      only_local_changes = false;
     } else if (delta_n < 0) {
       auto const type = reaction.reactant_types[i];
       for (int j = 0; j < -delta_n; j++) {
@@ -189,7 +194,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
         check_exclusion_range(p_id, type);
         hide_particle(p_id, type);
       }
-      on_particle_change();
+      only_local_changes = false;
     }
   }
   // create or hide particles of types with noncorresponding replacement types
@@ -204,7 +209,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
         check_exclusion_range(p_id, type);
         hide_particle(p_id, type);
       }
-      on_particle_change();
     } else {
       // create additional product_types particles
       auto const type = reaction.product_types[i];
@@ -213,9 +217,14 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction,
         check_exclusion_range(p_id, type);
         bookkeeping.created.emplace_back(p_id);
       }
-      on_particle_change();
     }
   }
+  // determine which fine-grained event to trigger
+  if (n_product_types == n_reactant_types and only_local_changes) {
+    on_particle_local_change();
+  } else {
+    on_particle_change();
+  }
 }
 
 std::unordered_map<int, int>

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -114,6 +114,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
     BOOST_REQUIRE_GE(get_particle_node_parallel(pid2), rank ? -1 : 1);
     BOOST_REQUIRE_GE(get_particle_node_parallel(pid3), rank ? -1 : 1);
   }
+  set_particle_property(pid1, &Particle::mol_id, type_a);
+  set_particle_property(pid2, &Particle::mol_id, type_b);
+  set_particle_property(pid3, &Particle::mol_id, type_b);
 
   auto const reset_particle_positions = [&start_positions]() {
     for (auto const &kv : start_positions) {
@@ -225,17 +228,18 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
     on_non_bonded_ia_change();
 
     // matrix indices and reference energy value
-    auto const max_type = type_b + 1;
-    auto const n_pairs = Utils::upper_triangular(type_b, type_b, max_type) + 1;
-    auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, max_type);
-    auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, max_type);
+    auto const size = std::max(type_a, type_b) + 1;
+    auto const n_pairs = Utils::upper_triangular(type_b, type_b, size) + 1;
+    auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, size);
+    auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, size);
     auto const frac6 = Utils::int_pow<6>(sig / r_off);
     auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift);
 
     // measure energies
     auto const obs_energy = calculate_energy();
     if (rank == 0) {
       for (int i = 0; i < n_pairs; ++i) {
+        // particles were set up with type == mol_id
         auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.;
         auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.;
         BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10);