From 354f540ba0094f0c9c46f083ba55b9fbe5c17ed0 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Wed, 28 Apr 2021 20:51:45 +0200 Subject: [PATCH] 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++) {