Separate short-range forces by signature (#4768)

Fixes #4755 Description of changes: - split short-range force kernel into a central radial kernel, a central radial with charge kernel, and a non-central kernel
espressomd · Aug 8, 2023 · 4aa06e7 · 4aa06e7
2 parents c0a3b74 + 15337d2
commit 4aa06e7
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 18 deletions.
diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp
@@ -91,8 +91,11 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
 
     if (dist > 0) {
       outer_normal_vec = -dist_vec / dist;
-      pf = calc_non_bonded_pair_force(p, part_rep, ia_params, dist_vec, dist,
-                                      get_ptr(coulomb_kernel));
+      pf = calc_central_radial_force(p, part_rep, ia_params, dist_vec, dist) +
+           calc_central_radial_charge_force(p, part_rep, ia_params, dist_vec,
+                                            dist, get_ptr(coulomb_kernel)) +
+           calc_non_central_force(p, part_rep, ia_params, dist_vec, dist);
+
 #ifdef DPD
       if (thermo_switch & THERMO_DPD) {
         dpd_force =
@@ -103,8 +106,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
 #endif
     } else if (m_penetrable && (dist <= 0)) {
       if ((!m_only_positive) && (dist < 0)) {
-        pf = calc_non_bonded_pair_force(p, part_rep, ia_params, dist_vec, -dist,
-                                        get_ptr(coulomb_kernel));
+        pf =
+            calc_central_radial_force(p, part_rep, ia_params, dist_vec, -dist) +
+            calc_central_radial_charge_force(p, part_rep, ia_params, dist_vec,
+                                             -dist, get_ptr(coulomb_kernel)) +
+            calc_non_central_force(p, part_rep, ia_params, dist_vec, -dist);
+
 #ifdef DPD
         if (thermo_switch & THERMO_DPD) {
           dpd_force = dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,

diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp
@@ -74,10 +74,11 @@
 
 #include <tuple>
 
-inline ParticleForce calc_non_bonded_pair_force(
-    Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
-    Utils::Vector3d const &d, double const dist,
-    Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel) {
+inline ParticleForce calc_central_radial_force(Particle const &p1,
+                                               Particle const &p2,
+                                               IA_parameters const &ia_params,
+                                               Utils::Vector3d const &d,
+                                               double const dist) {
 
   ParticleForce pf{};
   double force_factor = 0;
@@ -133,19 +134,39 @@ inline ParticleForce calc_non_bonded_pair_force(
 #ifdef LJCOS2
   force_factor += ljcos2_pair_force_factor(ia_params, dist);
 #endif
-/* Thole damping */
-#ifdef THOLE
-  pf.f += thole_pair_force(p1, p2, ia_params, d, dist, coulomb_kernel);
-#endif
 /* tabulated */
 #ifdef TABULATED
   force_factor += tabulated_pair_force_factor(ia_params, dist);
 #endif
+  pf.f += force_factor * d;
+  return pf;
+}
+
+inline ParticleForce calc_central_radial_charge_force(
+    Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
+    Utils::Vector3d const &d, double const dist,
+    Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel) {
+
+  ParticleForce pf{};
+/* Thole damping */
+#ifdef THOLE
+  pf.f += thole_pair_force(p1, p2, ia_params, d, dist, coulomb_kernel);
+#endif
+
+  return pf;
+}
+
+inline ParticleForce calc_non_central_force(Particle const &p1,
+                                            Particle const &p2,
+                                            IA_parameters const &ia_params,
+                                            Utils::Vector3d const &d,
+                                            double const dist) {
+
+  ParticleForce pf{};
 /* Gay-Berne */
 #ifdef GAY_BERNE
   pf += gb_pair_force(p1.quat(), p2.quat(), ia_params, d, dist);
 #endif
-  pf.f += force_factor * d;
   return pf;
 }
 
@@ -189,10 +210,15 @@ inline void add_non_bonded_pair_force(
 
   if (dist < ia_params.max_cut) {
 #ifdef EXCLUSIONS
-    if (do_nonbonded(p1, p2))
+    if (do_nonbonded(p1, p2)) {
+#endif
+      pf += calc_central_radial_force(p1, p2, ia_params, d, dist);
+      pf += calc_central_radial_charge_force(p1, p2, ia_params, d, dist,
+                                             coulomb_kernel);
+      pf += calc_non_central_force(p1, p2, ia_params, d, dist);
+#ifdef EXCLUSIONS
+    }
 #endif
-      pf += calc_non_bonded_pair_force(p1, p2, ia_params, d, dist,
-                                       coulomb_kernel);
   }
 
   /***********************************************/

diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp
@@ -64,8 +64,11 @@ inline void add_non_bonded_pair_virials(
 #endif
   {
     auto const &ia_params = get_ia_param(p1.type(), p2.type());
-    auto const force =
-        calc_non_bonded_pair_force(p1, p2, ia_params, d, dist, kernel_forces).f;
+    auto const force = calc_central_radial_force(p1, p2, ia_params, d, dist).f +
+                       calc_central_radial_charge_force(p1, p2, ia_params, d,
+                                                        dist, kernel_forces)
+                           .f +
+                       calc_non_central_force(p1, p2, ia_params, d, dist).f;
     auto const stress = Utils::tensor_product(d, force);
 
     auto const type1 = p1.mol_id();