Skip to content

Commit

Permalink
Merge #3123
Browse files Browse the repository at this point in the history 
3123: ELC Coverage, integer flags converted r=jngrad a=reinaual

Should fix #2685 for the neutral case.
Also converted integer flags to boolean in ELC and MMM2D.




Co-authored-by: Alexander Reinauer <[email protected]>
  • Loading branch information
@bors @reinaual
bors[bot] and reinaual authored Sep 6, 2019
2 parents 95673d5 + 26fe011 commit 6c7bc5a
Show file tree
Hide file tree
Showing 8 changed files with 44 additions and 40 deletions.
30 changes: 17 additions & 13 deletions src/core/electrostatics_magnetostatics/elc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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";
}
Expand Down
8 changes: 4 additions & 4 deletions src/core/electrostatics_magnetostatics/elc.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand Down
20 changes: 10 additions & 10 deletions src/core/electrostatics_magnetostatics/mmm2d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
/*@{*/
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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 &&
Expand All @@ -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();
Expand All @@ -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;
}
}

Expand Down
8 changes: 4 additions & 4 deletions src/core/electrostatics_magnetostatics/mmm2d.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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);
Expand Down
4 changes: 2 additions & 2 deletions src/python/espressomd/electrostatic_extensions.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/python/espressomd/electrostatic_extensions.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
4 changes: 2 additions & 2 deletions src/python/espressomd/electrostatics.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down
6 changes: 3 additions & 3 deletions testsuite/python/elc_vs_mmm2d_neutral.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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"])
Expand Down

0 comments on commit 6c7bc5a

Please sign in to comment.