From 15337d214a396c34b3da433dfcaa3d3ee65f2cc3 Mon Sep 17 00:00:00 2001
From: Somesh Kurahatti <skurahatti@icp.uni-stuttgart.de>
Date: Tue, 8 Aug 2023 13:15:32 +0200
Subject: [PATCH] Separate short forces by signature

---
 src/core/constraints/ShapeBasedConstraint.cpp | 15 ++++--
 src/core/forces_inline.hpp                    | 50 ++++++++++++++-----
 src/core/pressure_inline.hpp                  |  7 ++-
 3 files changed, 54 insertions(+), 18 deletions(-)

diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp
index 290e813b5a..8a3a1b8554 100644
--- 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
index 3858633060..e18a18c786 100644
--- 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
index 8467c2249b..a1f6354eb8 100644
--- 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();