Skip to content

Commit

Permalink
core: rollout of length_inv and length_half
Browse files Browse the repository at this point in the history
  • Loading branch information
@reinaual
reinaual committed May 18, 2021
1 parent 5398b44 commit 354f540
Show file tree
Hide file tree
Showing 7 changed files with 48 additions and 46 deletions.
10 changes: 5 additions & 5 deletions src/core/electrostatics_magnetostatics/elc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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))) -
Expand All @@ -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))) -
Expand All @@ -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) {
Expand Down
14 changes: 8 additions & 6 deletions src/core/electrostatics_magnetostatics/mdlc_correction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &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++) {
Expand Down Expand Up @@ -243,7 +243,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &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;
Expand All @@ -259,8 +260,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &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++) {
Expand Down Expand Up @@ -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;
}
Expand Down
24 changes: 12 additions & 12 deletions src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);

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

Expand Down Expand Up @@ -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<AssignForces, 1, 7>(
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) */
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;
}
Expand All @@ -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 "
Expand Down Expand Up @@ -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();
Expand Down
38 changes: 19 additions & 19 deletions src/core/electrostatics_magnetostatics/p3m.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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 =
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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<int>(std::round(mesh_density * box_geo.length()[1]));
p3m.params.mesh[2] =
Expand All @@ -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],
Expand All @@ -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;

Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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;
}
Expand All @@ -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];
Expand Down Expand Up @@ -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();
Expand Down
4 changes: 2 additions & 2 deletions src/core/grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.) {
Expand Down
2 changes: 1 addition & 1 deletion src/core/integrators/velocity_verlet_npt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/core/statistics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,

auto const order_sq = order * order;
std::vector<double> 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++) {
Expand Down

0 comments on commit 354f540

Please sign in to comment.