Removed some code duplication

src/core/bonded_interactions/bonded_tab.cpp

@@ -49,11 +49,11 @@ int tabulated_bonded_set_params(int bond_type,
     bonded_ia_params[bond_type].num = 1;
   } else if (tab_type == TAB_BOND_ANGLE) {
     tab_pot->minval = 0.0;
-    tab_pot->maxval = PI + ROUND_ERROR_PREC;
+    tab_pot->maxval = Utils::pi() + ROUND_ERROR_PREC;
     bonded_ia_params[bond_type].num = 2;
   } else if (tab_type == TAB_BOND_DIHEDRAL) {
     tab_pot->minval = 0.0;
-    tab_pot->maxval = 2.0 * PI + ROUND_ERROR_PREC;
+    tab_pot->maxval = 2.0 * Utils::pi() + ROUND_ERROR_PREC;
     bonded_ia_params[bond_type].num = 3;
   } else {
     runtimeError("Unsupported tabulated bond type.");

src/core/bonded_interactions/dihedral.hpp

@@ -84,7 +84,7 @@ inline void calc_dihedral_angle(Particle const *p1, Particle const *p2,
   /* Calculate dihedral angle */
   *phi = acos(*cosphi);
   if (scalar(aXb, c) < 0.0)
-    *phi = (2.0 * PI) - *phi;
+    *phi = (2.0 * Utils::pi()) - *phi;
 }
 
 /** calculate dihedral force between particles p1, p2 p3 and p4
@@ -140,7 +140,8 @@ inline int calc_dihedral_force(Particle const *p2, Particle const *p1,
   if (fabs(sin(phi)) < TINY_SIN_VALUE) {
 #ifdef OLD_DIHEDRAL
     sinmphi_sinphi = iaparams->p.dihedral.mult *
-                     cos(2.0 * PI - iaparams->p.dihedral.mult * phi) / cos(phi);
+                     cos(2.0 * Utils::pi() - iaparams->p.dihedral.mult * phi) /
+                     cos(phi);
 #else
     /*(comes from taking the first term of the MacLaurin expansion of
       sin(n*phi - phi0) and sin(phi) and then making the division).

src/core/electrostatics_magnetostatics/mdlc_correction.cpp

@@ -608,11 +608,11 @@ int mdlc_tune(double error) {
 
   flag = 0;
   for (kc = 1; kc < limitkc; kc++) {
-    gc = kc * 2.0 * PI / lx;
+    gc = kc * 2.0 * Utils::pi() / lx;
     fa0 = sqrt(9.0 * exp(+2. * gc * h) * g1_DLC_dip(gc, lz - h) +
                22.0 * g1_DLC_dip(gc, lz) +
                9.0 * exp(-2.0 * gc * h) * g1_DLC_dip(gc, lz + h));
-    fa1 = 0.5 * sqrt(PI / (2.0 * a)) * fa0;
+    fa1 = 0.5 * sqrt(Utils::pi() / (2.0 * a)) * fa0;
     fa2 = g2_DLC_dip(gc, lz);
     de = n * (mu_max * mu_max) / (4.0 * (exp(gc * lz) - 1.0)) * (fa1 + fa2);
     if (de < error) {

src/core/electrostatics_magnetostatics/p3m-common.cpp

@@ -94,7 +94,7 @@ void p3m_add_block(double const *in, double *out, int const start[3],
 
 double p3m_analytic_cotangent_sum(int n, double mesh_i, int cao) {
   double c, res = 0.0;
-  c = Utils::sqr(cos(PI * mesh_i * (double)n));
+  c = Utils::sqr(cos(Utils::pi() * mesh_i * (double)n));
 
   switch (cao) {
   case 1: {

src/core/electrostatics_magnetostatics/p3m-dipolar.cpp

@@ -222,15 +222,15 @@ static double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]);
 /*
 
 // Do the sum over k<>0 where k={kx,ky,kz} with kx,ky,kz INTEGERS, of
-// exp(-PI**2*k**2/alpha**2/L**2)
+// exp(-Utils::pi()**2*k**2/alpha**2/L**2)
 static double dp3m_sumi1(double alpha_L){
        int k2,kx,ky,kz,kx2,ky2,limit=60;
        double suma,alpha_L2;
 
        alpha_L2= alpha_L* alpha_L;
 
        //fprintf(stderr,"alpha_L=%le\n",alpha_L);
-       //fprintf(stderr,"PI=%le\n",PI);
+       //fprintf(stderr,"Utils::pi()=%le\n",Utils::pi());
 
 
        suma=0.0;
@@ -240,7 +240,7 @@ static double dp3m_sumi1(double alpha_L){
          ky2=ky*ky;
        for(kz=-limit;kz<=limit;kz++){
            k2=kx2+ky2+kz*kz;
-           suma+=exp(-PI*PI*k2/(alpha_L*alpha_L));
+           suma+=exp(-Utils::pi()*Utils::pi()*k2/(alpha_L*alpha_L));
        }}}
        suma-=1; //It's easier to subtract the term k=0 later than put an if
 inside the loops
@@ -481,7 +481,7 @@ double dp3m_average_dipolar_self_energy(double box_l, int mesh) {
 
   MPI_Reduce(&node_phi, &phi, 1, MPI_DOUBLE, MPI_SUM, 0, comm_cart);
 
-  phi *= PI / 3. / box_l / pow(mesh, 3.0);
+  phi *= Utils::pi() / 3. / box_l / pow(mesh, 3.0);
 
   return phi;
 }
@@ -968,7 +968,7 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) {
           }
         }
       }
-      node_k_space_energy_dip *= dipole_prefac * PI / box_l[0];
+      node_k_space_energy_dip *= dipole_prefac * Utils::pi() / box_l[0];
       MPI_Reduce(&node_k_space_energy_dip, &k_space_energy_dip, 1, MPI_DOUBLE,
                  MPI_SUM, 0, comm_cart);
 
@@ -987,7 +987,7 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) {
         k_space_energy_dip -=
             dipole.prefactor *
             (dp3m.sum_mu2 * 2 * pow(dp3m.params.alpha_L * box_l_i[0], 3) *
-             wupii / 3.0);
+             Utils::sqrt_pi_i() / 3.0);
 
         double volume = box_l[0] * box_l[1] * box_l[2];
         k_space_energy_dip += dipole.prefactor * dp3m.energy_correction /
@@ -1073,7 +1073,7 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) {
         /* redistribute force component mesh */
         dp3m_spread_force_grid(dp3m.rs_mesh);
         /* Assign force component from mesh to particle */
-        P3M_assign_torques(dipole_prefac * (2 * PI / box_l[0]), d_rs);
+        P3M_assign_torques(dipole_prefac * (2 * Utils::pi() / box_l[0]), d_rs);
       }
       P3M_TRACE(fprintf(stderr, "%d: done torque calculation.\n", this_node));
 #endif /*if def ROTATION */
@@ -1158,7 +1158,8 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) {
         dp3m_spread_force_grid(dp3m.rs_mesh_dip[1]);
         dp3m_spread_force_grid(dp3m.rs_mesh_dip[2]);
         /* Assign force component from mesh to particle */
-        dp3m_assign_forces_dip(dipole_prefac * pow(2 * PI / box_l[0], 2), d_rs);
+        dp3m_assign_forces_dip(
+            dipole_prefac * pow(2 * Utils::pi() / box_l[0], 2), d_rs);
       }
 
       P3M_TRACE(fprintf(stderr,
@@ -1450,7 +1451,7 @@ double dp3m_perform_aliasing_sums_force(const int n[3], double nominator[1]) {
   nominator[0] = 0.0;
 
   f1 = 1.0 / (double)dp3m.params.mesh[0];
-  f2 = Utils::sqr(PI / (dp3m.params.alpha_L));
+  f2 = Utils::sqr(Utils::pi() / (dp3m.params.alpha_L));
 
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
     nmx = dp3m.meshift[n[0]] + dp3m.params.mesh[0] * mx;
@@ -1538,7 +1539,7 @@ double dp3m_perform_aliasing_sums_energy(const int n[3], double nominator[1]) {
   nominator[0] = 0.0;
 
   f1 = 1.0 / (double)dp3m.params.mesh[0];
-  f2 = Utils::sqr(PI / (dp3m.params.alpha_L));
+  f2 = Utils::sqr(Utils::pi() / (dp3m.params.alpha_L));
 
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
     nmx = dp3m.meshift[n[0]] + dp3m.params.mesh[0] * mx;
@@ -2183,7 +2184,8 @@ static double dp3m_k_space_error(double box_size, double prefac, int mesh,
             he_q += d;
         }
 
-  return 8. * PI * PI / 3. * sum_q2 * sqrt(he_q / (double)n_c_part) /
+  return 8. * Utils::pi() * Utils::pi() / 3. * sum_q2 *
+         sqrt(he_q / (double)n_c_part) /
          (box_size * box_size * box_size * box_size);
 }
 
@@ -2196,7 +2198,7 @@ void dp3m_tune_aliasing_sums(int nx, int ny, int nz, int mesh, double mesh_i,
 
   double ex, ex2, nm2, U2, factor1;
 
-  factor1 = Utils::sqr(PI * alpha_L_i);
+  factor1 = Utils::sqr(Utils::pi() * alpha_L_i);
 
   *alias1 = *alias2 = 0.0;
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
@@ -2598,9 +2600,10 @@ void dp3m_compute_constants_energy_dipolar() {
   P3M_TRACE(
       fprintf(stderr, "%d: Average Dipolar Energy = %lf.\n", this_node, Ukp3m));
 
-  Eself = -(2 * pow(dp3m.params.alpha_L, 3) * wupii / 3.0);
+  Eself = -(2 * pow(dp3m.params.alpha_L, 3) * Utils::sqrt_pi_i() / 3.0);
 
-  dp3m.energy_correction = -dp3m.sum_mu2 * (Ukp3m + Eself + 2. * PI / 3.);
+  dp3m.energy_correction =
+      -dp3m.sum_mu2 * (Ukp3m + Eself + 2. * Utils::pi() / 3.);
 }
 
 /*****************************************************************************/

src/core/electrostatics_magnetostatics/p3m-dipolar.hpp

@@ -249,7 +249,7 @@ inline double dp3m_add_pair_force(Particle *p1, Particle *p2, double const *d,
     double mir = dip1 * Vector3d{d[0], d[1], d[2]};
     double mjr = dip2 * Vector3d{d[0], d[1], d[2]};
 
-    coeff = 2.0 * dp3m.params.alpha * wupii;
+    coeff = 2.0 * dp3m.params.alpha * Utils::sqrt_pi_i();
     double dist2i = 1 / dist2;
     exp_adist2 = exp(-adist * adist);
 
@@ -331,7 +331,7 @@ inline double dp3m_pair_energy(Particle *p1, Particle *p2,
     mir = dip1 * Vector3d{d[0], d[1], d[2]};
     mjr = dip2 * Vector3d{d[0], d[1], d[2]};
 
-    coeff = 2.0 * dp3m.params.alpha * wupii;
+    coeff = 2.0 * dp3m.params.alpha * Utils::sqrt_pi_i();
     dist2i = 1 / dist2;
     exp_adist2 = exp(-adist * adist);
 

src/core/electrostatics_magnetostatics/p3m.cpp

@@ -845,11 +845,11 @@ double p3m_calc_kspace_forces(int force_flag, int energy_flag) {
                comm_cart);
     if (this_node == 0) {
       /* self energy correction */
-      k_space_energy -=
-          coulomb.prefactor * (p3m.sum_q2 * p3m.params.alpha * wupii);
+      k_space_energy -= coulomb.prefactor *
+                        (p3m.sum_q2 * p3m.params.alpha * Utils::sqrt_pi_i());
       /* net charge correction */
       k_space_energy -=
-          coulomb.prefactor * p3m.square_sum_q * PI /
+          coulomb.prefactor * p3m.square_sum_q * Utils::pi() /
           (2.0 * box_l[0] * box_l[1] * box_l[2] * Utils::sqr(p3m.params.alpha));
     }
 
@@ -889,12 +889,12 @@ double p3m_calc_kspace_forces(int force_flag, int energy_flag) {
         for (j[1] = 0; j[1] < p3m.fft.plan[3].new_mesh[1]; j[1]++) {
           for (j[2] = 0; j[2] < p3m.fft.plan[3].new_mesh[2]; j[2]++) {
             /* i*k*(Re+i*Im) = - Im*k + i*Re*k     (i=sqrt(-1)) */
-            p3m.rs_mesh[ind] = -2.0 * PI *
+            p3m.rs_mesh[ind] = -2.0 * Utils::pi() *
                                (p3m.ks_mesh[ind + 1] *
                                 d_operator[j[d] + p3m.fft.plan[3].start[d]]) /
                                box_l[d_rs];
             ind++;
-            p3m.rs_mesh[ind] = 2.0 * PI * p3m.ks_mesh[ind - 1] *
+            p3m.rs_mesh[ind] = 2.0 * Utils::pi() * p3m.ks_mesh[ind - 1] *
                                d_operator[j[d] + p3m.fft.plan[3].start[d]] /
                                box_l[d_rs];
             ind++;
@@ -1141,7 +1141,7 @@ inline double perform_aliasing_sums_force(int const n[3], double numerator[3]) {
   for (i = 0; i < 3; i++)
     numerator[i] = 0.0;
 
-  f1 = Utils::sqr(PI / (p3m.params.alpha));
+  f1 = Utils::sqr(Utils::pi() / (p3m.params.alpha));
 
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
     nmx = p3m.meshift_x[n[KX]] + p3m.params.mesh[RX] * mx;
@@ -1223,7 +1223,7 @@ template <int cao> void calc_influence_function_force() {
           } else
             fak3 = fak1 / (fak2 * Utils::sqr(denominator));
 
-          p3m.g_force[ind] = 2 * fak3 / (PI);
+          p3m.g_force[ind] = 2 * fak3 / (Utils::pi());
         }
       }
     }
@@ -1267,7 +1267,7 @@ template <int cao> inline double perform_aliasing_sums_energy(int const n[3]) {
   double sx, sy, sz, f1, f2, mx, my, mz, nmx, nmy, nmz, nm2, expo;
   double limit = 30;
 
-  f1 = Utils::sqr(PI / (p3m.params.alpha));
+  f1 = Utils::sqr(Utils::pi() / (p3m.params.alpha));
 
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
     nmx = p3m.meshift_x[n[KX]] + p3m.params.mesh[RX] * mx;
@@ -1330,7 +1330,8 @@ template <int cao> void calc_influence_function_energy() {
         }
 
         else
-          p3m.g_energy[ind] = perform_aliasing_sums_energy<cao>(n) / PI;
+          p3m.g_energy[ind] =
+              perform_aliasing_sums_energy<cao>(n) / Utils::pi();
       }
     }
   }
@@ -2072,7 +2073,7 @@ void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3],
 
   double ex, ex2, nm2, U2, factor1;
 
-  factor1 = Utils::sqr(PI * alpha_L_i);
+  factor1 = Utils::sqr(Utils::pi() * alpha_L_i);
 
   *alias1 = *alias2 = 0.0;
   for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
@@ -2399,12 +2400,12 @@ void p3m_calc_kspace_stress(double *stress) {
     for (j[0] = 0; j[0] < p3m.fft.plan[3].new_mesh[RX]; j[0]++) {
       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]++) {
-          kx = 2.0 * PI * p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] /
-               box_l[RX];
-          ky = 2.0 * PI * p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] /
-               box_l[RY];
-          kz = 2.0 * PI * p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] /
-               box_l[RZ];
+          kx = 2.0 * Utils::pi() *
+               p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] / box_l[RX];
+          ky = 2.0 * Utils::pi() *
+               p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] / box_l[RY];
+          kz = 2.0 * Utils::pi() *
+               p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] / box_l[RZ];
           sqk = Utils::sqr(kx) + Utils::sqr(ky) + Utils::sqr(kz);
           if (sqk == 0) {
             node_k_space_energy = 0.0;

src/core/electrostatics_magnetostatics/p3m.hpp

@@ -231,14 +231,16 @@ inline double p3m_add_pair_force(double chgfac, double const *d, double dist2,
       double erfc_part_ri = Utils::AS_erfc_part(adist) / dist;
       double fac1 = chgfac * exp(-adist * adist);
       double fac2 =
-          fac1 * (erfc_part_ri + 2.0 * p3m.params.alpha * wupii) / dist2;
+          fac1 * (erfc_part_ri + 2.0 * p3m.params.alpha * Utils::sqrt_pi_i()) /
+          dist2;
 #else
       erfc_part_ri = erfc(adist) / dist;
-      double fac1 = cchgfac;
-      double fac2 = fac1 *
-                    (erfc_part_ri +
-                     2.0 * p3m.params.alpha * wupii * exp(-adist * adist)) /
-                    dist2;
+      double fac1 = chgfac;
+      double fac2 =
+          fac1 *
+          (erfc_part_ri +
+           2.0 * p3m.params.alpha * Utils::sqrt_pi_i() * exp(-adist * adist)) /
+          dist2;
 #endif
       for (int j = 0; j < 3; j++)
         force[j] += fac2 * d[j];