From 354f540ba0094f0c9c46f083ba55b9fbe5c17ed0 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Wed, 28 Apr 2021 20:51:45 +0200 Subject: [PATCH 1/3] core: rollout of length_inv and length_half --- .../electrostatics_magnetostatics/elc.cpp | 10 ++--- .../mdlc_correction.cpp | 14 ++++--- .../p3m-dipolar.cpp | 24 ++++++------ .../electrostatics_magnetostatics/p3m.cpp | 38 +++++++++---------- src/core/grid.cpp | 4 +- src/core/integrators/velocity_verlet_npt.cpp | 2 +- src/core/statistics.cpp | 2 +- 7 files changed, 48 insertions(+), 46 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index 98358dc5f7a..da1067ca528 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/src/core/electrostatics_magnetostatics/elc.cpp @@ -222,7 +222,7 @@ static void add_dipole_force(const ParticleRange &particles) { /* for nonneutral systems, this shift gives the background contribution (rsp. for this shift, the DM of the background is zero) */ - double const shift = 0.5 * box_geo.length()[2]; + double const shift = box_geo.length_half()[2]; // collect moments @@ -284,7 +284,7 @@ static double dipole_energy(const ParticleRange &particles) { constexpr std::size_t size = 7; /* for nonneutral systems, this shift gives the background contribution (rsp. for this shift, the DM of the background is zero) */ - double const shift = 0.5 * box_geo.length()[2]; + double const shift = box_geo.length_half()[2]; // collect moments @@ -355,7 +355,7 @@ static double dipole_energy(const ParticleRange &particles) { /*****************************************************************/ inline double image_sum_b(double q, double z) { - double const shift = 0.5 * box_geo.length()[2]; + double const shift = box_geo.length_half()[2]; double const fac = elc_params.delta_mid_top * elc_params.delta_mid_bot; double const image_sum = (q / (1.0 - fac) * (z - 2.0 * fac * box_geo.length()[2] / (1.0 - fac))) - @@ -364,7 +364,7 @@ inline double image_sum_b(double q, double z) { } inline double image_sum_t(double q, double z) { - double const shift = 0.5 * box_geo.length()[2]; + double const shift = box_geo.length_half()[2]; double const fac = elc_params.delta_mid_top * elc_params.delta_mid_bot; double const image_sum = (q / (1.0 - fac) * (z + 2.0 * fac * box_geo.length()[2] / (1.0 - fac))) - @@ -379,7 +379,7 @@ static double z_energy(const ParticleRange &particles) { /* for nonneutral systems, this shift gives the background contribution (rsp. for this shift, the DM of the background is zero) */ - double const shift = 0.5 * box_geo.length()[2]; + double const shift = box_geo.length_half()[2]; if (elc_params.dielectric_contrast_on) { if (elc_params.const_pot) { diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp index 6e7f4bde1b3..50b7ba5bad4 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp @@ -126,8 +126,8 @@ double get_DLC_dipolar(int kcut, std::vector &fs, double s1z, s2z, s3z, s4z; double ss; - auto const facux = 2.0 * Utils::pi() / box_geo.length()[0]; - auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1]; + auto const facux = 2.0 * Utils::pi() * box_geo.length_inv()[0]; + auto const facuy = 2.0 * Utils::pi() * box_geo.length_inv()[1]; double energy = 0.0; for (int ix = -kcut; ix <= +kcut; ix++) { @@ -243,7 +243,8 @@ double get_DLC_dipolar(int kcut, std::vector &fs, // Multiply by the factors we have left during the loops - auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]); + auto const piarea = + Utils::pi() * box_geo.length_inv()[0] * box_geo.length_inv()[1]; for (int j = 0; j < n_local_particles; j++) { fs[j] *= piarea; @@ -259,8 +260,8 @@ double get_DLC_dipolar(int kcut, std::vector &fs, * %Algorithm implemented accordingly to @cite brodka04a. */ double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) { - auto const facux = 2.0 * Utils::pi() / box_geo.length()[0]; - auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1]; + auto const facux = 2.0 * Utils::pi() * box_geo.length_inv()[0]; + auto const facuy = 2.0 * Utils::pi() * box_geo.length_inv()[1]; double energy = 0.0; for (int ix = -kcut; ix <= +kcut; ix++) { @@ -314,7 +315,8 @@ double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) { // Multiply by the factors we have left during the loops - auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]); + auto const piarea = + Utils::pi() * box_geo.length_inv()[0] * box_geo.length_inv()[1]; energy *= (-piarea); return (this_node == 0) ? energy : 0.0; } diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index d2fe0cef479..79b9cdee07d 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -229,7 +229,7 @@ void dp3m_init() { void dp3m_set_tune_params(double r_cut, int mesh, int cao, double accuracy) { if (r_cut >= 0) { dp3m.params.r_cut = r_cut; - dp3m.params.r_cut_iL = r_cut / box_geo.length()[0]; + dp3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0]; } if (mesh >= 0) @@ -268,7 +268,7 @@ void dp3m_set_params(double r_cut, int mesh, int cao, double alpha, Dipole::set_method_local(DIPOLAR_P3M); dp3m.params.r_cut = r_cut; - dp3m.params.r_cut_iL = r_cut / box_geo.length()[0]; + dp3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0]; dp3m.params.mesh[2] = dp3m.params.mesh[1] = dp3m.params.mesh[0] = mesh; dp3m.params.cao = cao; dp3m.params.alpha = alpha; @@ -443,7 +443,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, } } node_k_space_energy_dip *= - dipole_prefac * Utils::pi() / box_geo.length()[0]; + dipole_prefac * Utils::pi() * box_geo.length_inv()[0]; boost::mpi::reduce(comm_cart, node_k_space_energy_dip, k_space_energy_dip, std::plus<>(), 0); @@ -454,7 +454,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, k_space_energy_dip -= dipole.prefactor * (dp3m.sum_mu2 * 2 * - pow(dp3m.params.alpha_L / box_geo.length()[0], 3) * + pow(dp3m.params.alpha_L * box_geo.length_inv()[0], 3) * Utils::sqrt_pi_i() / 3.0); auto const volume = box_geo.volume(); @@ -534,7 +534,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, /* Assign force component from mesh to particle */ Utils::integral_parameter( dp3m.params.cao, - dipole_prefac * (2 * Utils::pi() / box_geo.length()[0]), d_rs, + dipole_prefac * (2 * Utils::pi() * box_geo.length_inv()[0]), d_rs, particles); } @@ -623,8 +623,8 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, /* Assign force component from mesh to particle */ Utils::integral_parameter( dp3m.params.cao, - dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), d_rs, - particles); + dipole_prefac * pow(2 * Utils::pi() * box_geo.length_inv()[0], 2), + d_rs, particles); } } /* if (dp3m.sum_mu2 > 0) */ } /* if (force_flag) */ @@ -1143,8 +1143,8 @@ int dp3m_adaptive_tune(bool verbose) { r_cut_iL_min = 0; r_cut_iL_max = std::min(min_local_box_l, min_box_l / 2) - skin; - r_cut_iL_min *= 1. / box_geo.length()[0]; - r_cut_iL_max *= 1. / box_geo.length()[0]; + r_cut_iL_min *= box_geo.length_inv()[0]; + r_cut_iL_max *= box_geo.length_inv()[0]; } else { r_cut_iL_min = r_cut_iL_max = dp3m.params.r_cut_iL; @@ -1382,7 +1382,7 @@ double dp3m_rtbisection(double box_size, double prefac, double r_cut_iL, void dp3m_init_a_ai_cao_cut() { for (int i = 0; i < 3; i++) { - dp3m.params.ai[i] = (double)dp3m.params.mesh[i] / box_geo.length()[i]; + dp3m.params.ai[i] = dp3m.params.mesh[i] * box_geo.length_inv()[i]; dp3m.params.a[i] = 1.0 / dp3m.params.ai[i]; dp3m.params.cao_cut[i] = 0.5 * dp3m.params.a[i] * dp3m.params.cao; } @@ -1394,7 +1394,7 @@ bool dp3m_sanity_checks_boxl() { bool ret = false; for (int i = 0; i < 3; i++) { /* check k-space cutoff */ - if (dp3m.params.cao_cut[i] >= 0.5 * box_geo.length()[i]) { + if (dp3m.params.cao_cut[i] >= box_geo.length_half()[i]) { runtimeErrorMsg() << "dipolar P3M_init: k-space cutoff " << dp3m.params.cao_cut[i] << " is larger than half of box dimension " @@ -1469,7 +1469,7 @@ void dp3m_scaleby_box_l() { } dp3m.params.r_cut = dp3m.params.r_cut_iL * box_geo.length()[0]; - dp3m.params.alpha = dp3m.params.alpha_L / box_geo.length()[0]; + dp3m.params.alpha = dp3m.params.alpha_L * box_geo.length_inv()[0]; dp3m_init_a_ai_cao_cut(); p3m_calc_lm_ld_pos(dp3m.local_mesh, dp3m.params); dp3m_sanity_checks_boxl(); diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 58eb290e686..25758ec84c9 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -216,7 +216,7 @@ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double accuracy) { if (r_cut >= 0) { p3m.params.r_cut = r_cut; - p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]); + p3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0]; } else if (r_cut == -1.0) { p3m.params.r_cut = 0; p3m.params.r_cut_iL = 0; @@ -260,7 +260,7 @@ void p3m_set_params(double r_cut, const int *mesh, int cao, double alpha, coulomb.method = COULOMB_P3M; p3m.params.r_cut = r_cut; - p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]); + p3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0]; p3m.params.mesh[2] = mesh[2]; p3m.params.mesh[1] = mesh[1]; p3m.params.mesh[0] = mesh[0]; @@ -434,14 +434,14 @@ Utils::Vector9d p3m_calc_kspace_pressure_tensor() { for (j[1] = 0; j[1] < p3m.fft.plan[3].new_mesh[RY]; j[1]++) { for (j[2] = 0; j[2] < p3m.fft.plan[3].new_mesh[RZ]; j[2]++) { auto const kx = 2.0 * Utils::pi() * - p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] / - box_geo.length()[RX]; + p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] * + box_geo.length_inv()[RX]; auto const ky = 2.0 * Utils::pi() * - p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] / - box_geo.length()[RY]; + p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] * + box_geo.length_inv()[RY]; auto const kz = 2.0 * Utils::pi() * - p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] / - box_geo.length()[RZ]; + p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] * + box_geo.length_inv()[RZ]; auto const sqk = Utils::sqr(kx) + Utils::sqr(ky) + Utils::sqr(kz); auto const node_k_space_energy = @@ -509,8 +509,8 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag, int d_rs = (d + p3m.ks_pnum) % 3; /* directions */ auto const k = 2.0 * Utils::pi() * - p3m.d_op[d_rs][j[d] + p3m.fft.plan[3].start[d]] / - box_geo.length()[d_rs]; + p3m.d_op[d_rs][j[d] + p3m.fft.plan[3].start[d]] * + box_geo.length_inv()[d_rs]; /* i*k*(Re+i*Im) = - Im*k + i*Re*k (i=sqrt(-1)) */ p3m.E_mesh[d_rs][2 * ind + 0] = -k * phi_hat.imag(); @@ -675,7 +675,7 @@ static double p3m_mcr_time(const int mesh[3], int cao, double r_cut_iL, p3m.params.mesh[2] = mesh[2]; p3m.params.cao = cao; p3m.params.alpha_L = alpha_L; - p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]); + p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0]; /* initialize p3m structures */ mpi_bcast_coulomb_params(); @@ -1004,7 +1004,7 @@ int p3m_adaptive_tune(bool verbose) { */ } else if (p3m.params.mesh[1] == -1 && p3m.params.mesh[2] == -1) { double mesh_density = mesh_density_min = mesh_density_max = - p3m.params.mesh[0] / box_geo.length()[0]; + p3m.params.mesh[0] * box_geo.length_inv()[0]; p3m.params.mesh[1] = static_cast(std::round(mesh_density * box_geo.length()[1])); p3m.params.mesh[2] = @@ -1020,7 +1020,7 @@ int p3m_adaptive_tune(bool verbose) { } } else { mesh_density_min = mesh_density_max = - p3m.params.mesh[0] / box_geo.length()[0]; + p3m.params.mesh[0] * box_geo.length_inv()[0]; if (verbose) { std::printf("fixed mesh %d %d %d\n", p3m.params.mesh[0], @@ -1034,8 +1034,8 @@ int p3m_adaptive_tune(bool verbose) { r_cut_iL_min = 0; r_cut_iL_max = std::min(min_local_box_l, min_box_l / 2.0) - skin; - r_cut_iL_min *= (1. / box_geo.length()[0]); - r_cut_iL_max *= (1. / box_geo.length()[0]); + r_cut_iL_min *= box_geo.length_inv()[0]; + r_cut_iL_max *= box_geo.length_inv()[0]; } else { r_cut_iL_min = r_cut_iL_max = p3m.params.r_cut_iL; @@ -1136,7 +1136,7 @@ int p3m_adaptive_tune(bool verbose) { p3m.params.mesh[2] = mesh[2]; p3m.params.cao = cao; p3m.params.alpha_L = alpha_L; - p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]); + p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0]; p3m.params.accuracy = accuracy; /* broadcast tuned p3m parameters */ mpi_bcast_coulomb_params(); @@ -1244,7 +1244,7 @@ void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3], void p3m_init_a_ai_cao_cut() { for (int i = 0; i < 3; i++) { - p3m.params.ai[i] = (double)p3m.params.mesh[i] / box_geo.length()[i]; + p3m.params.ai[i] = p3m.params.mesh[i] * box_geo.length_inv()[i]; p3m.params.a[i] = 1.0 / p3m.params.ai[i]; p3m.params.cao_cut[i] = 0.5 * p3m.params.a[i] * p3m.params.cao; } @@ -1254,7 +1254,7 @@ bool p3m_sanity_checks_boxl() { bool ret = false; for (int i = 0; i < 3; i++) { /* check k-space cutoff */ - if (p3m.params.cao_cut[i] >= 0.5 * box_geo.length()[i]) { + if (p3m.params.cao_cut[i] >= box_geo.length_half()[i]) { runtimeErrorMsg() << "P3M_init: k-space cutoff " << p3m.params.cao_cut[i] << " is larger than half of box dimension " << box_geo.length()[i]; @@ -1346,7 +1346,7 @@ void p3m_scaleby_box_l() { } p3m.params.r_cut = p3m.params.r_cut_iL * box_geo.length()[0]; - p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]); + p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0]; p3m_init_a_ai_cao_cut(); p3m_calc_lm_ld_pos(p3m.local_mesh, p3m.params); p3m_sanity_checks_boxl(); diff --git a/src/core/grid.cpp b/src/core/grid.cpp index 843e909fc95..593fba1f562 100644 --- a/src/core/grid.cpp +++ b/src/core/grid.cpp @@ -111,8 +111,8 @@ void grid_changed_n_nodes() { } void rescale_boxl(int dir, double d_new) { - double scale = - (dir - 3) ? d_new / box_geo.length()[dir] : d_new / box_geo.length()[0]; + double scale = (dir - 3) ? d_new * box_geo.length_inv()[dir] + : d_new * box_geo.length_inv()[0]; /* If shrinking, rescale the particles first. */ if (scale <= 1.) { diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 34790cb375c..f73b0ea123c 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -115,7 +115,7 @@ void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) { L_new = pow(nptiso.volume, 1.0 / nptiso.dimension); - scal[1] = L_new / box_geo.length()[nptiso.non_const_dim]; + scal[1] = L_new * box_geo.length_inv()[nptiso.non_const_dim]; scal[0] = 1 / scal[1]; } MPI_Bcast(scal, 3, MPI_DOUBLE, 0, comm_cart); diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index d4f8571e818..e34ac40b269 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -274,7 +274,7 @@ void calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, auto const order_sq = order * order; std::vector ff(2 * order_sq + 1); - auto const twoPI_L = 2 * Utils::pi() / box_geo.length()[0]; + auto const twoPI_L = 2 * Utils::pi() * box_geo.length_inv()[0]; for (int i = 0; i <= order; i++) { for (int j = -order; j <= order; j++) { From 05f5a066e652e2fbd2ef3eb12b70980fff70bce2 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Wed, 28 Apr 2021 21:23:47 +0200 Subject: [PATCH 2/3] core: elc: removed custom inverse box dimensions --- .../electrostatics_magnetostatics/elc.cpp | 129 ++++++++++-------- 1 file changed, 73 insertions(+), 56 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index da1067ca528..36a44813e53 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/src/core/electrostatics_magnetostatics/elc.cpp @@ -50,11 +50,6 @@ #include #include -/** \name Inverse box dimensions */ -/**@{*/ -static double ux, uy, uz; -/**@}*/ - ELC_struct elc_params = {1e100, 10, 1, 0, true, true, false, 1, 1, false, 0, 0, 0, 0, 0.0}; @@ -107,12 +102,6 @@ static double dipole_energy(const ParticleRange &particles); static double z_energy(const ParticleRange &particles); static void add_z_force(const ParticleRange &particles); -void ELC_setup_constants() { - ux = 1 / box_geo.length()[0]; - uy = 1 / box_geo.length()[1]; - uz = 1 / box_geo.length()[2]; -} - /** * @brief Calculate cached sin/cos values for one direction. * @@ -215,7 +204,9 @@ inline void check_gap_elc(const Particle &p) { * See @cite yeh99a. */ static void add_dipole_force(const ParticleRange &particles) { - double const pref = coulomb.prefactor * 4 * Utils::pi() * ux * uy * uz; + double const pref = coulomb.prefactor * 4 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1] * + box_geo.length_inv()[2]; constexpr std::size_t size = 3; auto local_particles = particles; @@ -251,7 +242,7 @@ static void add_dipole_force(const ParticleRange &particles) { } gblcblk[0] *= pref; - gblcblk[1] *= pref / elc_params.h / uz; + gblcblk[1] *= pref / elc_params.h * box_geo.length()[2]; gblcblk[2] *= pref; distribute(size); @@ -280,7 +271,9 @@ static void add_dipole_force(const ParticleRange &particles) { * See @cite yeh99a. */ static double dipole_energy(const ParticleRange &particles) { - double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy * uz; + double const pref = coulomb.prefactor * 2 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1] * + box_geo.length_inv()[2]; constexpr std::size_t size = 7; /* for nonneutral systems, this shift gives the background contribution (rsp. for this shift, the DM of the background is zero) */ @@ -336,7 +329,8 @@ static double dipole_energy(const ParticleRange &particles) { if (elc_params.dielectric_contrast_on) { if (elc_params.const_pot) { // zero potential difference contribution - energy += pref / elc_params.h / uz * Utils::sqr(gblcblk[6]); + energy += + pref / elc_params.h * box_geo.length()[2] * Utils::sqr(gblcblk[6]); // external potential shift contribution energy -= 2 * elc_params.pot_diff / elc_params.h * gblcblk[6]; } @@ -374,7 +368,8 @@ inline double image_sum_t(double q, double z) { /*****************************************************************/ static double z_energy(const ParticleRange &particles) { - double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy; + double const pref = coulomb.prefactor * 2 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1]; constexpr std::size_t size = 4; /* for nonneutral systems, this shift gives the background contribution @@ -452,7 +447,8 @@ static double z_energy(const ParticleRange &particles) { /*****************************************************************/ static void add_z_force(const ParticleRange &particles) { - double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy; + double const pref = coulomb.prefactor * 2 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1]; constexpr std::size_t size = 1; if (elc_params.dielectric_contrast_on) { @@ -512,7 +508,8 @@ template void setup_PoQ(std::size_t index, double omega, const ParticleRange &particles) { assert(index >= 1); - double const pref_di = coulomb.prefactor * 4 * Utils::pi() * ux * uy; + double const pref_di = coulomb.prefactor * 4 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1]; double const pref = -pref_di / expm1(omega * box_geo.length()[2]); constexpr std::size_t size = 4; double lclimgebot[4], lclimgetop[4], lclimge[4]; @@ -659,7 +656,8 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega, const ParticleRange &particles) { assert(index_p >= 1); assert(index_q >= 1); - double const pref_di = coulomb.prefactor * 8 * Utils::pi() * ux * uy; + double const pref_di = coulomb.prefactor * 8 * Utils::pi() * + box_geo.length_inv()[0] * box_geo.length_inv()[1]; double const pref = -pref_di / expm1(omega * box_geo.length()[2]); constexpr std::size_t size = 8; double lclimgebot[8], lclimgetop[8], lclimge[8]; @@ -789,8 +787,10 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega, static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega, const ParticleRange &particles) { constexpr double c_2pi = 2 * Utils::pi(); - double const pref_x = c_2pi * ux * static_cast(index_p) / omega; - double const pref_y = c_2pi * uy * static_cast(index_q) / omega; + double const pref_x = + c_2pi * box_geo.length_inv()[0] * static_cast(index_p) / omega; + double const pref_y = + c_2pi * box_geo.length_inv()[1] * static_cast(index_q) / omega; constexpr std::size_t size = 8; std::size_t ic = 0; @@ -847,45 +847,54 @@ static double PQ_energy(double omega, std::size_t n_part) { void ELC_add_force(const ParticleRange &particles) { constexpr double c_2pi = 2 * Utils::pi(); - auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1); - auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1); + auto const n_scxcache = + std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1); + auto const n_scycache = + std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1); - prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy); + prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache, + box_geo.length_inv()[1]); partblk.resize(particles.size() * 8); add_dipole_force(particles); add_z_force(particles); /* the second condition is just for the case of numerical accident */ - for (std::size_t p = 1; - ux * static_cast(p - 1) < elc_params.far_cut && p <= n_scxcache; + for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast(p - 1) < + elc_params.far_cut && + p <= n_scxcache; p++) { - auto const omega = c_2pi * ux * static_cast(p); + auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast(p); setup_PoQ(p, omega, particles); distribute(4); add_PoQ_force(particles); } - for (std::size_t q = 1; - uy * static_cast(q - 1) < elc_params.far_cut && q <= n_scycache; + for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast(q - 1) < + elc_params.far_cut && + q <= n_scycache; q++) { - auto const omega = c_2pi * uy * static_cast(q); + auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast(q); setup_PoQ(q, omega, particles); distribute(4); add_PoQ_force(particles); } - for (std::size_t p = 1; - ux * static_cast(p - 1) < elc_params.far_cut && p <= n_scxcache; + for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast(p - 1) < + elc_params.far_cut && + p <= n_scxcache; p++) { for (std::size_t q = 1; - Utils::sqr(ux * static_cast(p - 1)) + - Utils::sqr(uy * static_cast(q - 1)) < + Utils::sqr(box_geo.length_inv()[0] * static_cast(p - 1)) + + Utils::sqr(box_geo.length_inv()[1] * + static_cast(q - 1)) < elc_params.far_cut2 && q <= n_scycache; q++) { - auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast(p)) + - Utils::sqr(uy * static_cast(q))); + auto const omega = + c_2pi * + sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast(p)) + + Utils::sqr(box_geo.length_inv()[1] * static_cast(q))); setup_PQ(p, q, omega, particles); distribute(8); add_PQ_force(p, q, omega, particles); @@ -898,43 +907,52 @@ double ELC_energy(const ParticleRange &particles) { auto energy = dipole_energy(particles); energy += z_energy(particles); - auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1); - auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1); - prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy); + auto const n_scxcache = + std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1); + auto const n_scycache = + std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1); + prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache, + box_geo.length_inv()[1]); auto const n_localpart = particles.size(); partblk.resize(n_localpart * 8); /* the second condition is just for the case of numerical accident */ - for (std::size_t p = 1; - ux * static_cast(p - 1) < elc_params.far_cut && p <= n_scxcache; + for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast(p - 1) < + elc_params.far_cut && + p <= n_scxcache; p++) { - auto const omega = c_2pi * ux * static_cast(p); + auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast(p); setup_PoQ(p, omega, particles); distribute(4); energy += PoQ_energy(omega, n_localpart); } - for (std::size_t q = 1; - uy * static_cast(q - 1) < elc_params.far_cut && q <= n_scycache; + for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast(q - 1) < + elc_params.far_cut && + q <= n_scycache; q++) { - auto const omega = c_2pi * uy * static_cast(q); + auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast(q); setup_PoQ(q, omega, particles); distribute(4); energy += PoQ_energy(omega, n_localpart); } - for (std::size_t p = 1; - ux * static_cast(p - 1) < elc_params.far_cut && p <= n_scxcache; + for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast(p - 1) < + elc_params.far_cut && + p <= n_scxcache; p++) { for (std::size_t q = 1; - Utils::sqr(ux * static_cast(p - 1)) + - Utils::sqr(uy * static_cast(q - 1)) < + Utils::sqr(box_geo.length_inv()[0] * static_cast(p - 1)) + + Utils::sqr(box_geo.length_inv()[1] * + static_cast(q - 1)) < elc_params.far_cut2 && q <= n_scycache; q++) { - auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast(p)) + - Utils::sqr(uy * static_cast(q))); + auto const omega = + c_2pi * + sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast(p)) + + Utils::sqr(box_geo.length_inv()[1] * static_cast(q))); setup_PQ(p, q, omega, particles); distribute(8); energy += PQ_energy(omega, n_localpart); @@ -949,7 +967,8 @@ double ELC_tune_far_cut(ELC_struct const ¶ms) { constexpr auto maximal_far_cut = 50.; double const h = params.h; double lz = box_geo.length()[2]; - double const min_inv_boxl = std::min(ux, uy); + double const min_inv_boxl = + std::min(box_geo.length_inv()[0], box_geo.length_inv()[1]); if (params.dielectric_contrast_on) { // adjust lz according to dielectric layer method @@ -965,7 +984,8 @@ double ELC_tune_far_cut(ELC_struct const ¶ms) { do { const auto prefactor = 2 * Utils::pi() * far_cut; - const auto sum = prefactor + 2 * (ux + uy); + const auto sum = + prefactor + 2 * (box_geo.length_inv()[0] + box_geo.length_inv()[1]); const auto den = -expm1(-prefactor * lz); const auto num1 = exp(prefactor * (h - lz)); const auto num2 = exp(-prefactor * (h + lz)); @@ -1018,8 +1038,6 @@ void ELC_sanity_checks(ELC_struct const ¶ms) { } void ELC_init() { - - ELC_setup_constants(); elc_params.h = box_geo.length()[2] - elc_params.gap_size; if (elc_params.dielectric_contrast_on) { @@ -1107,7 +1125,6 @@ void ELC_set_params(double maxPWerror, double gap_size, double far_cut, ELC_sanity_checks(new_elc_params); - ELC_setup_constants(); if (new_elc_params.far_calculated) { new_elc_params.far_cut = ELC_tune_far_cut(new_elc_params); } From 54ea3b5573c458dcbea581e0364c6e3eebf6a297 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Tue, 18 May 2021 16:03:11 +0200 Subject: [PATCH 3/3] core: mmm1d: removed private inverse box_l implementation --- .../electrostatics_magnetostatics/mmm1d.cpp | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/mmm1d.cpp b/src/core/electrostatics_magnetostatics/mmm1d.cpp index 3204db5018d..1453fcec6d5 100644 --- a/src/core/electrostatics_magnetostatics/mmm1d.cpp +++ b/src/core/electrostatics_magnetostatics/mmm1d.cpp @@ -75,7 +75,7 @@ /** @name inverse box dimensions and other constants */ /**@{*/ -static double uz, L2, uz2, prefuz2, prefL3_i; +static double uz2, prefuz2, prefL3_i; /**@}*/ MMM1D_struct mmm1d_params = {0.05, 1e-5, 0}; @@ -85,8 +85,9 @@ static std::vector bessel_radii; static double far_error(int P, double minrad) { // this uses an upper bound to all force components and the potential - auto const rhores = 2 * Utils::pi() * uz * minrad; - auto const pref = 4 * uz * std::max(1.0, 2 * Utils::pi() * uz); + auto const rhores = 2 * Utils::pi() * box_geo.length_inv()[2] * minrad; + auto const pref = 4 * box_geo.length_inv()[2] * + std::max(1.0, 2 * Utils::pi() * box_geo.length_inv()[2]); return pref * K1(rhores * P) * exp(rhores) / rhores * (P - 1 + 1 / rhores); } @@ -171,11 +172,9 @@ int MMM1D_init() { if (mmm1d_params.far_switch_radius_2 >= Utils::sqr(box_geo.length()[2])) mmm1d_params.far_switch_radius_2 = 0.8 * Utils::sqr(box_geo.length()[2]); - uz = 1 / box_geo.length()[2]; - L2 = box_geo.length()[2] * box_geo.length()[2]; - uz2 = uz * uz; + uz2 = Utils::sqr(box_geo.length_inv()[2]); prefuz2 = coulomb.prefactor * uz2; - prefL3_i = prefuz2 * uz; + prefL3_i = prefuz2 * box_geo.length_inv()[2]; determine_bessel_radii(mmm1d_params.maxPWerror, MAXIMAL_B_CUT); prepare_polygamma_series(mmm1d_params.maxPWerror, @@ -189,7 +188,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d, auto const n_modPsi = static_cast(modPsi.size() >> 1); auto const rxy2 = d[0] * d[0] + d[1] * d[1]; auto const rxy2_d = rxy2 * uz2; - auto const z_d = d[2] * uz; + auto const z_d = d[2] * box_geo.length_inv()[2]; Utils::Vector3d F; if (rxy2 <= mmm1d_params.far_switch_radius_2) { @@ -245,7 +244,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d, } else { /* far range formula */ auto const rxy = sqrt(rxy2); - auto const rxy_d = rxy * uz; + auto const rxy_d = rxy * box_geo.length_inv()[2]; double sr = 0, sz = 0; for (int bp = 1; bp < MAXIMAL_B_CUT; bp++) { @@ -266,7 +265,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d, sr *= uz2 * 4 * c_2pi; sz *= uz2 * 4 * c_2pi; - auto const pref = sr / rxy + 2 * uz / rxy2; + auto const pref = sr / rxy + 2 * box_geo.length_inv()[2] / rxy2; F = {pref * d[0], pref * d[1], sz}; } @@ -283,7 +282,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d, auto const n_modPsi = static_cast(modPsi.size() >> 1); auto const rxy2 = d[0] * d[0] + d[1] * d[1]; auto const rxy2_d = rxy2 * uz2; - auto const z_d = d[2] * uz; + auto const z_d = d[2] * box_geo.length_inv()[2]; double E; if (rxy2 <= mmm1d_params.far_switch_radius_2) { @@ -301,7 +300,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d, r2n *= rxy2_d; } - E *= uz; + E *= box_geo.length_inv()[2]; /* real space parts */ @@ -319,7 +318,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d, } else { /* far range formula */ auto const rxy = sqrt(rxy2); - auto const rxy_d = rxy * uz; + auto const rxy_d = rxy * box_geo.length_inv()[2]; /* The first Bessel term will compensate a little bit the log term, so add them close together */ E = -0.25 * log(rxy2_d) + 0.5 * (Utils::ln_2() - Utils::gamma()); @@ -330,7 +329,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d, auto const fq = c_2pi * bp; E += K0(fq * rxy_d) * cos(fq * z_d); } - E *= 4 * uz; + E *= 4 * box_geo.length_inv()[2]; } return chpref * E;