diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index ddf4e281d56..ea0498ad998 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/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"; } diff --git a/src/core/electrostatics_magnetostatics/elc.hpp b/src/core/electrostatics_magnetostatics/elc.hpp index d75d1a0e8d8..bdb316764e8 100644 --- a/src/core/electrostatics_magnetostatics/elc.hpp +++ b/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 diff --git a/src/core/electrostatics_magnetostatics/mmm2d.cpp b/src/core/electrostatics_magnetostatics/mmm2d.cpp index 9074ea48050..227e3ce4a4f 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.cpp +++ b/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; } } diff --git a/src/core/electrostatics_magnetostatics/mmm2d.hpp b/src/core/electrostatics_magnetostatics/mmm2d.hpp index 015f57ef8d5..723bd0e42c2 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.hpp +++ b/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); diff --git a/src/python/espressomd/electrostatic_extensions.pxd b/src/python/espressomd/electrostatic_extensions.pxd index 773f28e7813..a97717edb9e 100644 --- a/src/python/espressomd/electrostatic_extensions.pxd +++ b/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 diff --git a/src/python/espressomd/electrostatic_extensions.pyx b/src/python/espressomd/electrostatic_extensions.pyx index afb1558317f..7b0e6adf667 100644 --- a/src/python/espressomd/electrostatic_extensions.pyx +++ b/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") diff --git a/src/python/espressomd/electrostatics.pxd b/src/python/espressomd/electrostatics.pxd index f653506ed74..f8657eba92e 100644 --- a/src/python/espressomd/electrostatics.pxd +++ b/src/python/espressomd/electrostatics.pxd @@ -236,7 +236,7 @@ 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; @@ -244,7 +244,7 @@ IF ELECTROSTATICS: 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(); diff --git a/testsuite/python/elc_vs_mmm2d_neutral.py b/testsuite/python/elc_vs_mmm2d_neutral.py index 9afd832c2e4..6a1bf5f5f9d 100644 --- a/testsuite/python/elc_vs_mmm2d_neutral.py +++ b/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"])