Fix Simple Pore

src/core/reaction_ensemble.cpp

@@ -717,9 +717,9 @@ int ReactionAlgorithm::create_particle(int desired_type) {
   // for components
   double vel[3];
   // we use mass=1 for all particles, think about adapting this
-  vel[0] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
-  vel[1] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
-  vel[2] = std::pow(2 * PI * temperature, -3.0 / 2.0) * gaussian_random();
+  vel[0] = std::sqrt(temperature) * gaussian_random();
+  vel[1] = std::sqrt(temperature) * gaussian_random();
+  vel[2] = std::sqrt(temperature) * gaussian_random();
 #ifdef ELECTROSTATICS
   double charge = charges_of_types[desired_type];
 #endif
@@ -831,6 +831,12 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
     p_id = p_id_s_changed_particles[i];
     // change particle position
     new_pos = get_random_position_in_box();
+    double vel[3];
+    auto const &p = get_particle_data(p_id);
+    vel[0] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    vel[1] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    vel[2] = std::sqrt(temperature / p.p.mass) * gaussian_random();
+    set_particle_v(p_id, vel);
     // new_pos=get_random_position_in_box_enhanced_proposal_of_small_radii();
     // //enhanced proposal of small radii
     place_particle(p_id, new_pos.data());

src/core/shapes/SimplePore.cpp

@@ -34,11 +34,26 @@ std::pair<double, double> SimplePore::dist_half_pore(double r, double z) const {
   assert(z >= 0.0);
   assert(r >= 0.0);
 
-  if (z <= c_z) {
-    /* Cylinder section */
+  /*
+   *  We have to find the line that splits area 1 (r determines distance) from
+   *  area 2 (z determines distance) inside pore. In area 3 we have to consider
+   *  z and r to determine the distance.
+   *
+   *   |        x
+   *   |   2  x
+   *   |    x
+   *  _|_ x    1
+   *    \|       ^ r
+   *  3  |-------|
+   *     |   z <-
+   */
+
+  if ((z <= c_z) && (r <= (c_z + c_r - z))) {
+    /* Cylinder section, inner */
     return {m_rad - r, 0};
-  } else if ((r >= c_r) && (z >= m_half_length)) {
-    /* Wall section */
+  } else if (((z >= c_z) && (r >= c_r)) ||
+             ((z <= c_z) && (r > (c_z + c_r - z)))) {
+    /* Wall section and outer cylinder */
     return {0, m_half_length - z};
   } else {
     /* Smoothing area */
@@ -73,11 +88,27 @@ void SimplePore::calculate_dist(const Vector3d &pos, double *dist,
   double dr, dz;
   std::tie(dr, dz) = dist_half_pore(r, std::abs(z));
 
+  double side = -1;
+  if (((dz == 0) && (r <= m_rad)) ||                  // cylinder section
+      ((dr == 0) && (std::abs(z) > m_half_length))) { // ||
+    side = 1;
+  } else {
+    // smoothing area
+    if (std::abs(z) >= c_z) {
+      double angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad);
+      double dist_offset =
+          m_smoothing_rad - (std::cos(angle) * m_smoothing_rad);
+      if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) {
+        side = 1;
+      }
+    }
+  }
+
   if (z <= 0.0) {
     dz *= -1;
   }
 
-  *dist = std::sqrt(dr * dr + dz * dz);
+  *dist = std::sqrt(dr * dr + dz * dz) * side;
   for (int i = 0; i < 3; i++) {
     vec[i] = -dr * e_r[i] + -dz * e_z[i];
   }

src/python/espressomd/reaction_ensemble.pyx

@@ -16,6 +16,10 @@ cdef class ReactionAlgorithm(object):
     reaction algorithm by setting the standard pressure, temperature, and the
     exclusion radius.
 
+    Note: When creating particles the velocities of the new particles are set
+    according the Maxwell-Boltzmann distribution. In this step the mass of the
+    new particle is assumed to equal 1.
+
 
     Parameters
     ----------
@@ -226,6 +230,11 @@ cdef class ReactionAlgorithm(object):
                 "No dictionary for relation between types and default charges provided.")
         #check electroneutrality of the provided reaction
         if(self._params["check_for_electroneutrality"]):
+            charges = np.array(list(self._params["default_charges"].values()))
+            if(np.count_nonzero(charges) == 0):
+                # all partices have zero charge
+                # no need to check electroneutrality
+                return
             total_charge_change = 0.0
             for i in range(len(self._params["reactant_coefficients"])):
                 type_here = self._params["reactant_types"][i]
@@ -235,7 +244,6 @@ cdef class ReactionAlgorithm(object):
                 type_here = self._params["product_types"][j]
                 total_charge_change += self._params["product_coefficients"][
                     j] * self._params["default_charges"][type_here]
-            charges = np.array(list(self._params["default_charges"].values()))
             min_abs_nonzero_charge = np.min(
                 np.abs(charges[np.nonzero(charges)[0]]))
             if abs(total_charge_change) / min_abs_nonzero_charge > 1e-10:

testsuite/python/simple_pore.py

@@ -31,7 +31,17 @@
            "Features not available, skipping test!")
 class SimplePoreConstraint(ut.TestCase):
 
-    def test(self):
+    def test_orientation(self):
+        pore = SimplePore(
+            axis=[1., 0., 0.], radius=2., smoothing_radius=.1, length=2., center=[5., 5., 5.])
+
+        d, _ = pore.call_method("calc_distance", position=[.0, .0, .0])
+        self.assertGreater(d, 0.)
+
+        d, _ = pore.call_method("calc_distance", position=[5., 5., .0])
+        self.assertLess(d, 0.)
+
+    def test_stability(self):
         s = espressomd.System(box_l=[1.0, 1.0, 1.0])
         s.seed = s.cell_system.get_state()['n_nodes'] * [1234]
         box_yz = 15.