Merge branch 'python' of https://github.com/espressomd/espresso into …

…simplify_cph_ensemble_with_symmetric_proposal_probability

doc/sphinx/advanced_methods.rst

@@ -1283,7 +1283,7 @@ Initialization
     system.time_step = 0.0
     system.cell_system.skin = 0.4
     ek = espressomd.electrokinetics.Electrokinetics(agrid=1.0, lb_density=1.0,
-        viscosity=1.0, friction=1.0, T=1.0, prefactor=1.0,
+        viscosity=1.0, ext_force_density = [1,0,0], friction=1.0, T=1.0, prefactor=1.0,
         stencil='linkcentered', advection=True, fluid_coupling='friction')
     system.actors.add(ek)
 
@@ -1303,7 +1303,8 @@ that your computer contains a CUDA capable GPU which is sufficiently
 modern.
 
 To set up a proper LB fluid using this command one has to specify at
-least the following options: ``agrid``, ``lb_density``, ``viscosity``, ``friction``, ``T``, and ``prefactor``. The other options can be
+least the following options: ``agrid``, ``lb_density``, ``viscosity``, 
+``friction``, ``T``, and ``prefactor``. The other options can be
 used to modify the behavior of the LB fluid. Note that the command does
 not allow the user to set the time step parameter as is the case for the
 lattice-Boltzmann command, this parameter is instead taken directly from the value set for
@@ -1382,7 +1383,7 @@ Boundaries
 ^^^^^^^^^^
 ::
 
-    ek_boundary = espressomd.electrokinetics.EKBoundary(charge_density=1.0, shape=my_shape)
+    ek_boundary = espressomd.ekboundaries.EKBoundary(charge_density=1.0, shape=my_shape)
     system.ekboundaries.add(ek_boundary)
 
 .. note:: Feature ``EK_BOUNDARIES`` required

doc/sphinx/constraints.rst

@@ -440,14 +440,20 @@ HollowConicalFrustum
 
 :class:`espressomd.shapes.HollowConicalFrustum`
 
-A hollow cone with round corners. The specific parameters
+A hollow cone with round corners.
+Can include an opening in the side (see figures below). The specific parameters
 are described in the shape's class :class:`espressomd.shapes.HollowConicalFrustum`.
 
 .. figure:: figures/shape-conical_frustum.png
    :alt: Visualization of a constraint with a HollowConicalFrustum shape.
    :align: center
    :height: 6.00000cm
 
+.. figure:: figures/shape-hollowconicalfrustum_central_angle.png
+   :alt: Visualization a HollowConicalFrustum shape with central angle
+   :align: center
+   :height: 6.00000cm
+
 .. figure:: figures/conical_frustum.png
    :alt: Schematic for the HollowConicalFrustum shape with labeled geometrical parameters.
    :align: center

doc/tutorials/active_matter/exercises/rectification_geometry.py

@@ -30,7 +30,7 @@
 from espressomd import lb
 from espressomd.lbboundaries import LBBoundary
 import espressomd.shapes
-
+import espressomd.math
 
 # Setup constants
 
@@ -90,9 +90,12 @@
 ORAD = (DIAMETER - IRAD) / np.sin(ANGLE)
 SHIFT = 0.25 * ORAD * np.cos(ANGLE)
 
+ctp = espressomd.math.CylindricalTransformationParameters(
+    axis=[-1, 0, 0], center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0])
+
 hollow_cone = LBBoundary(shape=espressomd.shapes.HollowConicalFrustum(
-    center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0],
-    axis=[-1, 0, 0], r1=ORAD, r2=IRAD, thickness=2.0, length=18,
+    cyl_transform_params=ctp,
+    r1=ORAD, r2=IRAD, thickness=2.0, length=18,
     direction=1))
 system.lbboundaries.add(hollow_cone)
 

doc/tutorials/active_matter/exercises/rectification_simulation.py

@@ -29,6 +29,7 @@
 import espressomd
 from espressomd import assert_features
 import espressomd.shapes
+import espressomd.math
 
 
 assert_features(["ENGINE", "LENNARD_JONES", "ROTATION", "MASS"])
@@ -119,6 +120,7 @@ def a2quat(phi, theta):
 system.constraints.add(shape=wall, particle_type=1)
 
 # Setup cone
+ctp = espressomd.math.CylindricalTransformationParameters(...)
 hollow_cone = espressomd.shapes.HollowConicalFrustum(...)
 system.constraints.add(shape=hollow_cone, particle_type=1)
 

doc/tutorials/active_matter/solutions/rectification_geometry.py

@@ -30,6 +30,7 @@
 from espressomd import lb
 from espressomd.lbboundaries import LBBoundary
 import espressomd.shapes
+import espressomd.math
 
 
 # Setup constants
@@ -91,9 +92,12 @@
 ORAD = (DIAMETER - IRAD) / np.sin(ANGLE)
 SHIFT = 0.25 * ORAD * np.cos(ANGLE)
 
+ctp = espressomd.math.CylindricalTransformationParameters(
+    axis=[-1, 0, 0], center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0])
+
 hollow_cone = LBBoundary(shape=espressomd.shapes.HollowConicalFrustum(
-    center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0],
-    axis=[-1, 0, 0], r1=ORAD, r2=IRAD, thickness=2.0, length=18,
+    cyl_transform_params=ctp,
+    r1=ORAD, r2=IRAD, thickness=2.0, length=18,
     direction=1))
 system.lbboundaries.add(hollow_cone)
 

doc/tutorials/active_matter/solutions/rectification_simulation.py

@@ -29,6 +29,7 @@
 import espressomd
 from espressomd import assert_features
 import espressomd.shapes
+import espressomd.math
 
 
 assert_features(["ENGINE", "LENNARD_JONES", "ROTATION", "MASS"])
@@ -118,9 +119,12 @@ def a2quat(phi, theta):
 ORAD = (DIAMETER - IRAD) / np.sin(ANGLE)
 SHIFT = 0.25 * ORAD * np.cos(ANGLE)
 
+ctp = espressomd.math.CylindricalTransformationParameters(
+    axis=[-1, 0, 0], center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0])
+
 hollow_cone = espressomd.shapes.HollowConicalFrustum(
-    center=[BOX_L[0] / 2.0 - 1.3 * SHIFT, BOX_L[1] / 2.0, BOX_L[2] / 2.0],
-    axis=[-1, 0, 0], r1=ORAD, r2=IRAD, thickness=2.0, length=18,
+    cyl_transform_params=ctp,
+    r1=ORAD, r2=IRAD, thickness=2.0, length=18,
     direction=1)
 system.constraints.add(shape=hollow_cone, particle_type=1)
 

samples/visualization_constraints.py

@@ -22,6 +22,7 @@
 import argparse
 
 import espressomd
+import espressomd.math
 import espressomd.shapes
 import espressomd.visualization_opengl
 
@@ -96,9 +97,11 @@
         particle_type=0, penetrable=True)
 
 elif args.shape == "HollowConicalFrustum":
+    ctp = espressomd.math.CylindricalTransformationParameters(
+        axis=1 / np.sqrt(2) * np.array([0.0, 1.0, 1.0]), center=[25, 25, 25], orientation=[1, 0, 0])
     system.constraints.add(shape=espressomd.shapes.HollowConicalFrustum(
         r1=12, r2=8, length=15.0, thickness=3,
-        axis=[0.0, 1.0, 1.0], center=[25, 25, 25], direction=1),
+        cyl_transform_params=ctp, direction=1, central_angle=np.pi / 2),
         particle_type=0, penetrable=True)
 
 elif args.shape == "Torus":

src/core/grid_based_algorithms/electrokinetics.hpp

@@ -48,7 +48,7 @@ typedef struct {
   float friction;
   float T;
   float prefactor;
-  float lb_force_density[3];
+  float lb_ext_force_density[3];
   unsigned int number_of_species;
   int reaction_species[3];
   float rho_reactant_reservoir;
@@ -150,6 +150,9 @@ int ek_init();
 int ek_set_agrid(float agrid);
 int ek_set_lb_density(float lb_density);
 int ek_set_viscosity(float viscosity);
+int ek_set_lb_ext_force_density(float lb_ext_force_dens_x,
+                                float lb_ext_force_dens_y,
+                                float lb_ext_force_dens_z);
 int ek_set_friction(float friction);
 int ek_set_T(float T);
 int ek_set_prefactor(float prefactor);

src/core/grid_based_algorithms/electrokinetics_cuda.cu

@@ -103,7 +103,7 @@ EK_parameters ek_parameters = {
     -1.0,
     // prefactor
     -1.0,
-    // lb_force_density
+    // lb_ext_force_density
     {0.0, 0.0, 0.0},
     // number_of_species
     0,
@@ -2241,49 +2241,33 @@ int ek_init() {
     lbpar_gpu.rho =
         (ek_parameters.lb_density < 0.0)
             ? 1.0f
-            : ek_parameters.lb_density * Utils::int_pow<3>(ek_parameters.agrid);
+            : ek_parameters.lb_density * Utils::int_pow<3>(lbpar_gpu.agrid);
 
     lbpar_gpu.is_TRT = true;
 
     lb_reinit_parameters_gpu();
     lbpar_gpu.viscosity = ek_parameters.viscosity * lbpar_gpu.time_step /
-                          Utils::sqr(ek_parameters.agrid);
+                          Utils::sqr(lbpar_gpu.agrid);
     lbpar_gpu.bulk_viscosity = ek_parameters.bulk_viscosity *
                                lbpar_gpu.time_step /
-                               Utils::sqr(ek_parameters.agrid);
-    lb_reinit_parameters_gpu();
+                               Utils::sqr(lbpar_gpu.agrid);
 
-    lb_init_gpu();
+    lbpar_gpu.external_force_density =
+        ek_parameters.lb_ext_force_density[0] != 0 ||
+        ek_parameters.lb_ext_force_density[1] != 0 ||
+        ek_parameters.lb_ext_force_density[2] != 0;
+    lbpar_gpu.ext_force_density =
+        Utils::Vector3f(ek_parameters.lb_ext_force_density) *
+        Utils::sqr(lbpar_gpu.agrid * lbpar_gpu.time_step);
 
-    if (ek_parameters.lb_force_density[0] != 0 ||
-        ek_parameters.lb_force_density[1] != 0 ||
-        ek_parameters.lb_force_density[2] != 0) {
-      lbpar_gpu.external_force_density = 1;
-      lbpar_gpu.ext_force_density[0] =
-          ek_parameters.lb_force_density[0] * ek_parameters.agrid *
-          ek_parameters.agrid * ek_parameters.time_step *
-          ek_parameters.time_step;
-      lbpar_gpu.ext_force_density[1] =
-          ek_parameters.lb_force_density[1] * ek_parameters.agrid *
-          ek_parameters.agrid * ek_parameters.time_step *
-          ek_parameters.time_step;
-      lbpar_gpu.ext_force_density[2] =
-          ek_parameters.lb_force_density[2] * ek_parameters.agrid *
-          ek_parameters.agrid * ek_parameters.time_step *
-          ek_parameters.time_step;
-      lb_reinit_extern_nodeforce_GPU(&lbpar_gpu);
-    } else {
-      lbpar_gpu.external_force_density = 0;
-      lbpar_gpu.ext_force_density[0] = 0;
-      lbpar_gpu.ext_force_density[1] = 0;
-      lbpar_gpu.ext_force_density[2] = 0;
-    }
+    lb_reinit_parameters_gpu();
+    lb_init_gpu();
 
+    ek_parameters.time_step = lbpar_gpu.time_step;
     ek_parameters.dim_x = lbpar_gpu.dim[0];
     ek_parameters.dim_x_padded = (ek_parameters.dim_x / 2 + 1) * 2;
     ek_parameters.dim_y = lbpar_gpu.dim[1];
     ek_parameters.dim_z = lbpar_gpu.dim[2];
-    ek_parameters.time_step = lbpar_gpu.time_step;
     ek_parameters.number_of_nodes =
         ek_parameters.dim_x * ek_parameters.dim_y * ek_parameters.dim_z;
 
@@ -3535,9 +3519,10 @@ void ek_print_parameters() {
   printf("  float friction = %f;\n", ek_parameters.friction);
   printf("  float T = %f;\n", ek_parameters.T);
   printf("  float prefactor = %f;\n", ek_parameters.prefactor);
-  printf("  float lb_force_density[] = {%f, %f, %f};\n",
-         ek_parameters.lb_force_density[0], ek_parameters.lb_force_density[1],
-         ek_parameters.lb_force_density[2]);
+  printf("  float lb_ext_force_density[] = {%f, %f, %f};\n",
+         ek_parameters.lb_ext_force_density[0],
+         ek_parameters.lb_ext_force_density[1],
+         ek_parameters.lb_ext_force_density[2]);
   printf("  unsigned int number_of_species = %d;\n",
          ek_parameters.number_of_species);
   printf("  int reaction_species[] = {%d, %d, %d};\n",
@@ -3697,6 +3682,15 @@ int ek_set_viscosity(float viscosity) {
   return 0;
 }
 
+int ek_set_lb_ext_force_density(float lb_ext_force_dens_x,
+                                float lb_ext_force_dens_y,
+                                float lb_ext_force_dens_z) {
+  ek_parameters.lb_ext_force_density[0] = lb_ext_force_dens_x;
+  ek_parameters.lb_ext_force_density[1] = lb_ext_force_dens_y;
+  ek_parameters.lb_ext_force_density[2] = lb_ext_force_dens_z;
+  return 0;
+}
+
 int ek_set_friction(float friction) {
 
   ek_parameters.friction = friction;

src/core/observables/PidObservable.hpp

@@ -60,8 +60,6 @@ class PidObservable : virtual public Observable {
 public:
   explicit PidObservable(std::vector<int> ids) : m_ids(std::move(ids)) {}
   std::vector<double> operator()() const final;
-
-  std::vector<int> &ids() { return m_ids; }
   std::vector<int> const &ids() const { return m_ids; }
 };
 

src/core/particle_data.cpp

@@ -1142,39 +1142,49 @@ void init_type_map(int type) {
   if (type < 0)
     throw std::runtime_error("Types may not be negative");
 
-  // fill particle map
-  if (particle_type_map.count(type) == 0)
-    particle_type_map[type] = std::unordered_set<int>();
-
+  auto &map_for_type = particle_type_map[type];
+  map_for_type.clear();
   for (auto const &p : partCfg()) {
     if (p.p.type == type)
-      particle_type_map.at(type).insert(p.p.identity);
+      map_for_type.insert(p.p.identity);
   }
 }
 
 void remove_id_from_map(int part_id, int type) {
-  if (particle_type_map.find(type) != particle_type_map.end())
-    particle_type_map.at(type).erase(part_id);
+  auto it = particle_type_map.find(type);
+  if (it != particle_type_map.end())
+    it->second.erase(part_id);
 }
 
 int get_random_p_id(int type, int random_index_in_type_map) {
-  if (random_index_in_type_map + 1 > particle_type_map.at(type).size())
+  auto it = particle_type_map.find(type);
+  if (it == particle_type_map.end()) {
+    throw std::runtime_error("The provided particle type " +
+                             std::to_string(type) +
+                             " is currently not tracked by the system.");
+  }
+
+  if (random_index_in_type_map + 1 > it->second.size())
     throw std::runtime_error("The provided index exceeds the number of "
                              "particle types listed in the particle_type_map");
-  return *std::next(particle_type_map[type].begin(), random_index_in_type_map);
+  return *std::next(it->second.begin(), random_index_in_type_map);
 }
 
 void add_id_to_type_map(int part_id, int type) {
-  if (particle_type_map.find(type) != particle_type_map.end())
-    particle_type_map.at(type).insert(part_id);
+  auto it = particle_type_map.find(type);
+  if (it != particle_type_map.end())
+    it->second.insert(part_id);
 }
 
 int number_of_particles_with_type(int type) {
-  if (particle_type_map.count(type) == 0)
+  auto it = particle_type_map.find(type);
+  if (it == particle_type_map.end()) {
     throw std::runtime_error("The provided particle type " +
                              std::to_string(type) +
                              " is currently not tracked by the system.");
-  return static_cast<int>(particle_type_map.at(type).size());
+  }
+
+  return static_cast<int>(it->second.size());
 }
 
 // The following functions are used by the python interface to obtain

src/python/espressomd/electrokinetics.pxd

@@ -39,7 +39,7 @@ IF ELECTROKINETICS and CUDA:
                 float friction
                 float T
                 float prefactor
-                float lb_force_density[3]
+                float lb_ext_force_density[3]
                 unsigned int number_of_species
                 int reaction_species[3]
                 float rho_reactant_reservoir
@@ -88,7 +88,7 @@ IF ELECTROKINETICS and CUDA:
                 float friction
                 float T
                 float prefactor
-                float lb_force_density[3]
+                float lb_ext_force_density[3]
                 unsigned int number_of_species
                 int reaction_species[3]
                 float rho_reactant_reservoir
@@ -140,6 +140,7 @@ IF ELECTROKINETICS and CUDA:
         int ek_set_lb_density(float lb_density)
         int ek_set_viscosity(float viscosity)
         int ek_set_friction(float friction)
+        int ek_set_lb_ext_force_density(float lb_ext_force_dens_x, float lb_ext_force_dens_y, float lb_ext_force_dens_z)
         int ek_set_T(float T)
         int ek_set_prefactor(float prefactor)
         int ek_set_bulk_viscosity(float bulk_viscosity)