2684 coulomb sr test

src/core/CMakeLists.txt

@@ -58,7 +58,7 @@ set(EspressoCore_SRC
   bonded_interactions/angle_cossquare.cpp
   bonded_interactions/angle_harmonic.cpp
   bonded_interactions/bonded_coulomb.cpp
-  bonded_interactions/bonded_coulomb_p3m_sr.cpp
+  bonded_interactions/bonded_coulomb_sr.cpp
   bonded_interactions/bonded_interaction_data.cpp
   bonded_interactions/bonded_tab.cpp
   bonded_interactions/dihedral.cpp

src/core/bonded_interactions/bonded_coulomb_p3m_sr.cpp → src/core/bonded_interactions/bonded_coulomb_sr.cpp

@@ -20,21 +20,21 @@
 */
 /** \file
  *
- *  Implementation of \ref bonded_coulomb_p3m_sr.hpp
+ *  Implementation of \ref bonded_coulomb_sr.hpp
  */
-#include "bonded_coulomb_p3m_sr.hpp"
+#include "bonded_coulomb_sr.hpp"
 #include "communication.hpp"
 
-#ifdef P3M
+#ifdef ELECTROSTATICS
 
-int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2) {
+int bonded_coulomb_sr_set_params(int bond_type, double q1q2) {
   if (bond_type < 0)
     return ES_ERROR;
 
   make_bond_type_exist(bond_type);
 
-  bonded_ia_params[bond_type].p.bonded_coulomb_p3m_sr.q1q2 = q1q2;
-  bonded_ia_params[bond_type].type = BONDED_IA_BONDED_COULOMB_P3M_SR;
+  bonded_ia_params[bond_type].p.bonded_coulomb_sr.q1q2 = q1q2;
+  bonded_ia_params[bond_type].type = BONDED_IA_BONDED_COULOMB_SR;
   bonded_ia_params[bond_type].num = 1;
 
   /* broadcast interaction parameters */

src/core/bonded_interactions/bonded_coulomb_p3m_sr.hpp → src/core/bonded_interactions/bonded_coulomb_sr.hpp

@@ -18,24 +18,24 @@
   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
-#ifndef _BONDED_COULOMB_P3M_SR_HPP
-#define _BONDED_COULOMB_P3M_SR_HPP
+#ifndef _BONDED_COULOMB_SR_HPP
+#define _BONDED_COULOMB_SR_HPP
 /** \file
- *  Routines to calculate the BONDED_COULOMB_P3M_SR Energy or/and
- * BONDED_COULOMB_P3M_SR force for a particle pair. This is only the shortrange
- * part of P3M and first used to subtract certain intramolecular interactions in
- * combination with Thole damping \ref forces.cpp
+ *  Routines to calculate the BONDED_COULOMB_SR Energy or/and
+ * BONDED_COULOMB_SR force for a particle pair. This is only the shortrange
+ * part of any coulomb interaction and first used to subtract certain
+ * intramolecular interactions in combination with Thole damping \ref forces.cpp
  */
 
 /************************************************************/
 
 #include "config.hpp"
 
-#ifdef P3M
+#ifdef ELECTROSTATICS
 
 #include "bonded_interaction_data.hpp"
 #include "debug.hpp"
-#include "electrostatics_magnetostatics/p3m.hpp"
+#include "electrostatics_magnetostatics/coulomb_inline.hpp"
 #include "particle_data.hpp"
 #include "utils.hpp"
 
@@ -44,9 +44,9 @@
  *  @retval ES_OK on success
  *  @retval ES_ERROR on error
  */
-int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2);
+int bonded_coulomb_sr_set_params(int bond_type, double q1q2);
 
-/** Computes the BONDED_COULOMB_P3M_SR pair force.
+/** Computes the BONDED_COULOMB_SR pair force.
  *  @param[in]  p1        First particle.
  *  @param[in]  p2        Second particle.
  *  @param[in]  iaparams  Interaction parameters.
@@ -55,49 +55,51 @@ int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2);
  *  @retval 0
  */
 inline int
-calc_bonded_coulomb_p3m_sr_pair_force(Particle const *p1, Particle const *p2,
-                                      Bonded_ia_parameters const *iaparams,
-                                      double const dx[3], double force[3]) {
+calc_bonded_coulomb_sr_pair_force(Particle *p1, Particle *p2,
+                                  Bonded_ia_parameters const *iaparams,
+                                  double dx[3], double force[3]) {
   double dist2 = sqrlen(dx);
   double dist = sqrt(dist2);
-  if (dist < p3m.params.r_cut) {
-    // Set to zero because p3m adds forces
-    force[0] = force[1] = force[2] = 0.;
+  if (dist < coulomb_cutoff) {
+    // TODO ugly workaround
+    Vector3d forcevec{};
 
-    p3m_add_pair_force(iaparams->p.bonded_coulomb_p3m_sr.q1q2, dx, dist2, dist,
-                       force);
+    Coulomb::calc_pair_force(p1, p2, iaparams->p.bonded_coulomb_sr.q1q2, dx,
+                             dist, dist2, forcevec);
+    force[0] = forcevec[0];
+    force[1] = forcevec[1];
+    force[2] = forcevec[2];
 
     ONEPART_TRACE(if (p1->p.identity == check_id) fprintf(
         stderr,
-        "%d: OPT: BONDED_COULOMB_P3M_SR f = (%.3e,%.3e,%.3e) with part id=%d "
+        "%d: OPT: BONDED_COULOMB_SR f = (%.3e,%.3e,%.3e) with part id=%d "
         "at dist %f\n",
         this_node, p1->f.f[0], p1->f.f[1], p1->f.f[2], p2->p.identity, dist2));
     ONEPART_TRACE(if (p2->p.identity == check_id) fprintf(
         stderr,
-        "%d: OPT: BONDED_COULOMB_P3M_SR f = (%.3e,%.3e,%.3e) with part id=%d "
+        "%d: OPT: BONDED_COULOMB_SR f = (%.3e,%.3e,%.3e) with part id=%d "
         "at dist %f\n",
         this_node, p2->f.f[0], p2->f.f[1], p2->f.f[2], p1->p.identity, dist2));
   }
 
   return 0;
 }
 
-/** Computes the BONDED_COULOMB_P3M_SR pair energy.
+/** Computes the BONDED_COULOMB_SR pair energy.
  *  @param[in]  p1        First particle.
  *  @param[in]  p2        Second particle.
  *  @param[in]  iaparams  Interaction parameters.
  *  @param[in]  dx        %Distance between the particles.
  *  @param[out] _energy   Energy.
  *  @retval 0
  */
-inline int
-bonded_coulomb_p3m_sr_pair_energy(Particle const *p1, Particle const *p2,
-                                  Bonded_ia_parameters const *iaparams,
-                                  double const dx[3], double *_energy) {
+inline int bonded_coulomb_sr_pair_energy(Particle *p1, Particle *p2,
+                                         Bonded_ia_parameters const *iaparams,
+                                         double dx[3], double *_energy) {
   double dist2 = sqrlen(dx);
   double dist = sqrt(dist2);
-  *_energy = p3m_pair_energy(iaparams->p.bonded_coulomb_p3m_sr.q1q2, dist);
-
+  *_energy = Coulomb::add_pair_energy(
+      p1, p2, iaparams->p.bonded_coulomb_sr.q1q2, dx, dist, dist2);
   return 0;
 }
 

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -23,6 +23,8 @@ enum BondedInteraction {
   BONDED_IA_QUARTIC,
   /** Type of bonded interaction is a BONDED_COULOMB */
   BONDED_IA_BONDED_COULOMB,
+  /** Type of bonded interaction is a BONDED_COULOMB_SR */
+  BONDED_IA_BONDED_COULOMB_SR,
   /** Type of bonded interaction is a dihedral potential. */
   BONDED_IA_DIHEDRAL,
   /** Type of tabulated bonded interaction potential,
@@ -58,8 +60,6 @@ enum BondedInteraction {
   BONDED_IA_UMBRELLA,
   /** Type of bonded interaction is thermalized distance bond. */
   BONDED_IA_THERMALIZED_DIST,
-  /** Type of bonded interaction is a BONDED_COULOMB_P3M_SR */
-  BONDED_IA_BONDED_COULOMB_P3M_SR,
 };
 
 /** Specify tabulated bonded interactions  */
@@ -156,13 +156,11 @@ struct Bonded_coulomb_bond_parameters {
   double prefactor;
 };
 
-#ifdef P3M
-/** Parameters for Coulomb bond p3m short-range Potential */
-struct Bonded_coulomb_p3m_sr_bond_parameters {
+/** Parameters for Coulomb bond short-range Potential */
+struct Bonded_coulomb_sr_bond_parameters {
   /** charge factor */
   double q1q2;
 };
-#endif
 
 /** Parameters for three-body angular potential.
  *  @note
@@ -309,6 +307,7 @@ union Bond_parameters {
 #endif
   Quartic_bond_parameters quartic;
   Bonded_coulomb_bond_parameters bonded_coulomb;
+  Bonded_coulomb_sr_bond_parameters bonded_coulomb_sr;
   Angle_bond_parameters angle;
   Angle_harmonic_bond_parameters angle_harmonic;
   Angle_cosine_bond_parameters angle_cosine;
@@ -321,9 +320,6 @@ union Bond_parameters {
   Umbrella_bond_parameters umbrella;
 #endif
   Thermalized_bond_parameters thermalized_bond;
-#ifdef P3M
-  Bonded_coulomb_p3m_sr_bond_parameters bonded_coulomb_p3m_sr;
-#endif
   Subt_lj_bond_parameters subt_lj;
   Rigid_bond_parameters rigid_bond;
   IBM_Triel_Parameters ibm_triel;

src/core/electrostatics_magnetostatics/coulomb_inline.hpp

@@ -15,9 +15,9 @@
 
 namespace Coulomb {
 // forces_inline
-inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
-                            double dist2, Vector3d &force) {
-  auto const q1q2 = p1->p.q * p2->p.q;
+inline void calc_pair_force(Particle *p1, Particle *p2, double const q1q2,
+                            double *d, double dist, double const dist2,
+                            Vector3d &force) {
 
   if (q1q2 != 0) {
     Vector3d f{};
@@ -54,10 +54,10 @@ inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
       add_mmm2d_coulomb_pair_force(q1q2, d, dist2, dist, f.data());
       break;
     case COULOMB_DH:
-      add_dh_coulomb_pair_force(p1, p2, d, dist, f.data());
+      add_dh_coulomb_pair_force(q1q2, d, dist, f.data());
       break;
     case COULOMB_RF:
-      add_rf_coulomb_pair_force(p1, p2, d, dist, f.data());
+      add_rf_coulomb_pair_force(q1q2, d, dist, f.data());
       break;
 #ifdef SCAFACOS
     case COULOMB_SCAFACOS:
@@ -73,8 +73,8 @@ inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
 }
 
 // pressure_inline.hpp
-inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
-                              double dist, double dist2,
+inline void add_pair_pressure(Particle *p1, Particle *p2, double q1q2,
+                              double *d, double dist, double dist2,
                               Observable_stat &virials,
                               Observable_stat &p_tensor) {
   switch (coulomb.method) {
@@ -88,7 +88,7 @@ inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
   case COULOMB_DH:
   case COULOMB_RF: {
     Vector3d force{};
-    calc_pair_force(p1, p2, d, dist, dist2, force);
+    calc_pair_force(p1, p2, q1q2, d, dist, dist2, force);
 
     /* Calculate the virial pressure */
     for (int k = 0; k < 3; k++) {
@@ -107,40 +107,40 @@ inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
 }
 
 // energy_inline
-inline double add_pair_energy(Particle *p1, Particle *p2, double *d,
-                              double dist, double dist2) {
+inline double add_pair_energy(Particle *p1, Particle *p2, double const q1q2,
+                              double *d, double dist, double dist2) {
   /* real space Coulomb */
   auto E = [&]() {
     switch (coulomb.method) {
 #ifdef P3M
     case COULOMB_P3M_GPU:
     case COULOMB_P3M:
-      return p3m_pair_energy(p1->p.q * p2->p.q, dist);
+      // TODO some energy functions include the prefactor, some don't
+      return p3m_pair_energy(q1q2, dist);
     case COULOMB_ELC_P3M:
       if (elc_params.dielectric_contrast_on) {
         return 0.5 * ELC_P3M_dielectric_layers_energy_contribution(p1, p2) +
-               p3m_pair_energy(p1->p.q * p2->p.q, dist);
+               p3m_pair_energy(q1q2, dist);
       } else {
-        return p3m_pair_energy(p1->p.q * p2->p.q, dist);
+        return p3m_pair_energy(q1q2, dist);
       }
 #endif
 #ifdef SCAFACOS
     case COULOMB_SCAFACOS:
       return Scafacos::pair_energy(p1, p2, dist);
 #endif
     case COULOMB_DH:
-      return dh_coulomb_pair_energy(p1, p2, dist);
+      return dh_coulomb_pair_energy(q1q2, dist);
     case COULOMB_RF:
-      return rf_coulomb_pair_energy(p1, p2, dist);
+      return rf_coulomb_pair_energy(q1q2, dist);
     case COULOMB_MMM1D:
-      return mmm1d_coulomb_pair_energy(p1, p2, d, dist2, dist);
+      return mmm1d_coulomb_pair_energy(q1q2, d, dist2, dist);
     case COULOMB_MMM2D:
-      return mmm2d_coulomb_pair_energy(p1->p.q * p2->p.q, d, dist2, dist);
+      return mmm2d_coulomb_pair_energy(q1q2, d, dist2, dist);
     default:
       return 0.;
     }
   }();
-
   return coulomb.prefactor * E;
 }
 } // namespace Coulomb

src/core/electrostatics_magnetostatics/debye_hueckel.hpp

@@ -25,7 +25,6 @@
  *  for a particle pair.
  */
 #include "config.hpp"
-#include "electrostatics_magnetostatics/coulomb.hpp"
 
 #ifdef ELECTROSTATICS
 
@@ -56,32 +55,30 @@ int dh_set_params(double kappa, double r_cut);
     @param dist      Distance between p1 and p2.
     @param force     returns the force on particle 1.
 */
-inline void add_dh_coulomb_pair_force(Particle *p1, Particle *p2,
-                                      double const d[3], double dist,
-                                      double force[3]) {
+inline void add_dh_coulomb_pair_force(double const q1q2, double const d[3],
+                                      double const dist, double force[3]) {
   if (dist < dh_params.r_cut) {
     double fac;
     if (dh_params.kappa > 0.0) {
       /* debye hueckel case: */
       double kappa_dist = dh_params.kappa * dist;
-      fac = coulomb.prefactor * p1->p.q * p2->p.q *
-            (exp(-kappa_dist) / (dist * dist * dist)) * (1.0 + kappa_dist);
+      fac =
+          q1q2 * (exp(-kappa_dist) / (dist * dist * dist)) * (1.0 + kappa_dist);
     } else {
       /* pure Coulomb case: */
-      fac = coulomb.prefactor * p1->p.q * p2->p.q / (dist * dist * dist);
+      fac = q1q2 / (dist * dist * dist);
     }
     for (int j = 0; j < 3; j++)
       force[j] += fac * d[j];
   }
 }
 
-inline double dh_coulomb_pair_energy(Particle *p1, Particle *p2, double dist) {
+inline double dh_coulomb_pair_energy(double const q1q2, double const dist) {
   if (dist < dh_params.r_cut) {
     if (dh_params.kappa > 0.0)
-      return coulomb.prefactor * p1->p.q * p2->p.q *
-             exp(-dh_params.kappa * dist) / dist;
+      return q1q2 * exp(-dh_params.kappa * dist) / dist;
 
-    return coulomb.prefactor * p1->p.q * p2->p.q / dist;
+    return q1q2 / dist;
   }
   return 0.0;
 }

src/core/electrostatics_magnetostatics/icc.cpp

@@ -70,7 +70,7 @@ inline void add_non_bonded_pair_force_iccp3m(Particle *p1, Particle *p2,
                                              double d[3], double dist,
                                              double dist2) {
   Vector3d force{};
-  Coulomb::calc_pair_force(p1, p2, d, dist, dist2, force);
+  Coulomb::calc_pair_force(p1, p2, p1->p.q * p2->p.q, d, dist, dist2, force);
 
   p1->f.f += force;
   p2->f.f -= force;