ELC Coverage, integer flags converted

src/core/electrostatics_magnetostatics/elc.cpp

@@ -58,8 +58,8 @@
 static double ux, ux2, uy, uy2, uz, height_inverse;
 /*@}*/
 
-ELC_struct elc_params = {1e100, 10,    1, 0, 1, 1, false, 1,
-                         1,     false, 0, 0, 0, 0, 0.0};
+ELC_struct elc_params = {1e100, 10,    1, 0, true, true, false, 1,
+                         1,     false, 0, 0, 0,    0,    0.0};
 
 /****************************************
  * LOCAL ARRAYS
@@ -1152,14 +1152,18 @@ int ELC_tune(double error) {
     return ES_ERROR;
 
   elc_params.far_cut = min_inv_boxl;
+
   do {
-    err =
-        0.5 * (exp(2 * M_PI * elc_params.far_cut * h) / (lz - h) *
-                   (C_2PI * elc_params.far_cut + 2 * (ux + uy) + 1 / (lz - h)) /
-                   (expm1(2 * M_PI * elc_params.far_cut * lz)) +
-               exp(-2 * M_PI * elc_params.far_cut * h) / (lz + h) *
-                   (C_2PI * elc_params.far_cut + 2 * (ux + uy) + 1 / (lz + h)) /
-                   (expm1(2 * M_PI * elc_params.far_cut * lz)));
+    const auto prefactor = 2 * Utils::pi() * elc_params.far_cut;
+
+    const auto sum = prefactor + 2 * (ux + uy);
+    const auto den = -expm1(-prefactor * lz);
+    const auto num1 = exp(prefactor * (h - lz));
+    const auto num2 = exp(-prefactor * (h + lz));
+
+    err = 0.5 / den *
+          (num1 * (sum + 1 / (lz - h)) / (lz - h) +
+           num2 * (sum + 1 / (lz + h)) / (lz + h));
 
     elc_params.far_cut += min_inv_boxl;
   } while (err > error && elc_params.far_cut < MAXIMAL_FAR_CUT);
@@ -1263,7 +1267,7 @@ void ELC_on_resort_particles() {
 }
 
 int ELC_set_params(double maxPWerror, double gap_size, double far_cut,
-                   int neutralize, double delta_top, double delta_bot,
+                   bool neutralize, double delta_top, double delta_bot,
                    bool const_pot, double pot_diff) {
   elc_params.maxPWerror = maxPWerror;
   elc_params.gap_size = gap_size;
@@ -1276,7 +1280,7 @@ int ELC_set_params(double maxPWerror, double gap_size, double far_cut,
     elc_params.delta_mid_bot = delta_bot;
 
     // neutralize is automatic with dielectric contrast
-    elc_params.neutralize = 0;
+    elc_params.neutralize = false;
     // initial setup of parameters, may change later when P3M is finally tuned
     // set the space_layer to be 1/3 of the gap size, so that box = layer
     elc_params.space_layer = (1. / 3.) * gap_size;
@@ -1309,9 +1313,9 @@ int ELC_set_params(double maxPWerror, double gap_size, double far_cut,
   elc_params.far_cut = far_cut;
   if (far_cut != -1) {
     elc_params.far_cut2 = Utils::sqr(far_cut);
-    elc_params.far_calculated = 0;
+    elc_params.far_calculated = false;
   } else {
-    elc_params.far_calculated = 1;
+    elc_params.far_calculated = true;
     if (ELC_tune(elc_params.maxPWerror) == ES_ERROR) {
       runtimeErrorMsg() << "ELC tuning failed, gap size too small";
     }

src/core/electrostatics_magnetostatics/elc.hpp

@@ -48,13 +48,13 @@ typedef struct {
    */
   double gap_size;
   /** @copybrief MMM2D_struct::far_calculated */
-  int far_calculated;
+  bool far_calculated;
   /** Flag whether the box is neutralized by a homogeneous background.
    *  If true, use a homogeneous neutralizing background for nonneutral
    *  systems. Unlike the 3D case, this background adds an additional
    *  force pointing towards the system center, so be careful with this.
    */
-  int neutralize;
+  bool neutralize;
 
   /// @copybrief MMM2D_struct::dielectric_contrast_on
   bool dielectric_contrast_on;
@@ -64,7 +64,7 @@ typedef struct {
   /// @copybrief MMM2D_struct::delta_mid_bot
   double delta_mid_bot;
 
-  /// @copybrief MMM2D_struct::const_pot_on
+  /// @copybrief MMM2D_struct::const_pot
   bool const_pot;
   /// @copybrief MMM2D_struct::pot_diff
   double pot_diff;
@@ -103,7 +103,7 @@ extern ELC_struct elc_params;
  *  @retval ES_OK
  */
 int ELC_set_params(double maxPWerror, double min_dist, double far_cut,
-                   int neutralize, double delta_mid_top, double delta_mid_bot,
+                   bool neutralize, double delta_mid_top, double delta_mid_bot,
                    bool const_pot, double pot_diff);
 
 /// the force calculation

src/core/electrostatics_magnetostatics/mmm2d.cpp

@@ -135,7 +135,7 @@ static double max_near, min_far;
 ///
 static double self_energy;
 
-MMM2D_struct mmm2d_params = {1e100, 10, 1, 0, false, false, 0, 1, 1, 1};
+MMM2D_struct mmm2d_params = {1e100, 10, 1, false, false, false, 0, 1, 1, 1};
 
 /** return codes for \ref MMM2D_tune_near and \ref MMM2D_tune_far */
 /*@{*/
@@ -589,7 +589,7 @@ static void add_z_force(const ParticleRange &particles) {
   double field_tot = 0;
 
   /* Const. potential: subtract global dipole moment */
-  if (mmm2d_params.const_pot_on) {
+  if (mmm2d_params.const_pot) {
     double gbl_dm_z = 0;
     double lcl_dm_z = 0;
     for (auto const &p : particles) {
@@ -667,7 +667,7 @@ static double z_energy(const ParticleRange &particles) {
   }
 
   /* total dipole moment term, for capacitor feature */
-  if (mmm2d_params.const_pot_on) {
+  if (mmm2d_params.const_pot) {
     double gbl_dm_z = 0;
     double lcl_dm_z = 0;
 
@@ -1721,7 +1721,7 @@ void MMM2D_self_energy(const ParticleRange &particles) {
  ****************************************/
 
 int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top,
-                     double delta_bot, bool const_pot_on, double pot_diff) {
+                     double delta_bot, bool const_pot, double pot_diff) {
   int err;
 
   if (cell_structure.type != CELL_STRUCTURE_NSQUARE &&
@@ -1731,25 +1731,25 @@ int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top,
 
   mmm2d_params.maxPWerror = maxPWerror;
 
-  if (const_pot_on) {
+  if (const_pot) {
     mmm2d_params.dielectric_contrast_on = true;
     mmm2d_params.delta_mid_top = -1;
     mmm2d_params.delta_mid_bot = -1;
     mmm2d_params.delta_mult = 1;
-    mmm2d_params.const_pot_on = true;
+    mmm2d_params.const_pot = true;
     mmm2d_params.pot_diff = pot_diff;
   } else if (delta_top != 0.0 || delta_bot != 0.0) {
     mmm2d_params.dielectric_contrast_on = true;
     mmm2d_params.delta_mid_top = delta_top;
     mmm2d_params.delta_mid_bot = delta_bot;
     mmm2d_params.delta_mult = delta_top * delta_bot;
-    mmm2d_params.const_pot_on = false;
+    mmm2d_params.const_pot = false;
   } else {
     mmm2d_params.dielectric_contrast_on = false;
     mmm2d_params.delta_mid_top = 0;
     mmm2d_params.delta_mid_bot = 0;
     mmm2d_params.delta_mult = 0;
-    mmm2d_params.const_pot_on = false;
+    mmm2d_params.const_pot = false;
   }
 
   MMM2D_setup_constants();
@@ -1769,11 +1769,11 @@ int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top,
     mmm2d_params.far_cut = far_cut;
     mmm2d_params.far_cut2 = Utils::sqr(far_cut);
     if (mmm2d_params.far_cut > 0)
-      mmm2d_params.far_calculated = 0;
+      mmm2d_params.far_calculated = false;
     else {
       if ((err = MMM2D_tune_far(maxPWerror)))
         return err;
-      mmm2d_params.far_calculated = 1;
+      mmm2d_params.far_calculated = true;
     }
   }
 

src/core/electrostatics_magnetostatics/mmm2d.hpp

@@ -56,11 +56,11 @@ typedef struct {
    *  In the latter case, the cutoff will be adapted if important parameters,
    *  such as the box dimensions, change.
    */
-  int far_calculated;
+  bool far_calculated;
   /// flag whether there is any dielectric contrast in the system.
   bool dielectric_contrast_on;
   /// @brief Flag whether a const. potential is applied.
-  bool const_pot_on;
+  bool const_pot;
   /// @brief Const. potential.
   double pot_diff;
   /// dielectric contrast in the upper part of the simulation cell.
@@ -82,11 +82,11 @@ extern MMM2D_struct mmm2d_params;
  *                       Manual setting is probably only good for testing.
  *  @param delta_top     @copybrief MMM2D_struct::delta_mid_top
  *  @param delta_bot     @copybrief MMM2D_struct::delta_mid_bot
- *  @param const_pot_on  @copybrief MMM2D_struct::const_pot_on
+ *  @param const_pot  @copybrief MMM2D_struct::const_pot
  *  @param pot_diff      @copybrief MMM2D_struct::pot_diff
  */
 int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top,
-                     double delta_bot, bool const_pot_on, double pot_diff);
+                     double delta_bot, bool const_pot, double pot_diff);
 
 /** the general long range force/energy calculation */
 double MMM2D_add_far(int f, int e, const ParticleRange &particles);

src/python/espressomd/electrostatic_extensions.pxd

@@ -32,14 +32,14 @@ IF ELECTROSTATICS and P3M:
             double maxPWerror
             double gap_size
             double far_cut
-            int neutralize
+            bool neutralize
             double delta_mid_top,
             double delta_mid_bot,
             bool const_pot,
             double pot_diff
 
         int ELC_set_params(double maxPWerror, double min_dist, double far_cut,
-                           int neutralize, double delta_mid_top, double delta_mid_bot, bool const_pot, double pot_diff)
+                           bool neutralize, double delta_mid_top, double delta_mid_bot, bool const_pot, double pot_diff)
 
         # links intern C-struct with python object
         ELC_struct elc_params

src/python/espressomd/electrostatic_extensions.pyx

@@ -131,10 +131,10 @@ IF ELECTROSTATICS and P3M:
                 self._params["maxPWerror"],
                 self._params["gap_size"],
                 self._params["far_cut"],
-                int(self._params["neutralize"]),
+                self._params["neutralize"],
                 self._params["delta_mid_top"],
                 self._params["delta_mid_bot"],
-                int(self._params["const_pot"]),
+                self._params["const_pot"],
                     self._params["pot_diff"]):
                 handle_errors(
                     "ELC tuning failed, ELC is not set up to work with the GPU P3M")

src/python/espressomd/electrostatics.pxd

@@ -236,15 +236,15 @@ IF ELECTROSTATICS:
             double far_cut2;
             int far_calculated;
             bool dielectric_contrast_on;
-            bool const_pot_on;
+            bool const_pot;
             double pot_diff;
             double delta_mid_top;
             double delta_mid_bot;
             double delta_mult;
 
         cdef extern MMM2D_struct mmm2d_params;
 
-        int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top, double delta_bot, bool const_pot_on, double pot_diff);
+        int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top, double delta_bot, bool const_pot, double pot_diff);
 
         void MMM2D_init();
 

testsuite/python/elc_vs_mmm2d_neutral.py

@@ -26,8 +26,8 @@
 class ELC_vs_MMM2D_neutral(ut.TestCase):
     # Handle to espresso system
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
-    acc = 1e-6
-    elc_gap = 5.0
+    acc = 1e-7
+    elc_gap = 10.0
     box_l = 10.0
     bl2 = box_l * 0.5
     system.time_step = 0.01
@@ -127,7 +127,7 @@ def test_elc_vs_mmm2d(self):
         self.system.cell_system.node_grid = buf_node_grid
         self.system.periodicity = [1, 1, 1]
         p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=self.acc,
-                                            mesh=[16, 16, 24], cao=6)
+                                            mesh=[24, 24, 32], cao=6)
         self.system.actors.add(p3m)
 
         elc = electrostatic_extensions.ELC(**elc_param_sets["inert"])