Performance: short range forces low hanging fruits

- avoid a few constructions of ParticleForce
- avoid calls to npt and gay-berne if not active
- in the linked cell algorithm, first do cell-internal interactions to have better memory locality.

Improves thins by ~5% on my system (Intel)

Author: RudolfWeeber

src/core/BoxGeometry.hpp

@@ -207,8 +207,9 @@ class BoxGeometry {
    *         periodic images, i.e. <tt>a - b</tt>.
    */
   template <typename T>
-  Utils::Vector<T, 3> get_mi_vector(const Utils::Vector<T, 3> &a,
-                                    const Utils::Vector<T, 3> &b) const {
+  [[gnu::always_inline]] inline Utils::Vector<T, 3>
+  get_mi_vector(const Utils::Vector<T, 3> &a,
+                const Utils::Vector<T, 3> &b) const {
     if (type() == BoxType::LEES_EDWARDS) {
       auto const shear_plane_normal = lees_edwards_bc().shear_plane_normal;
       auto a_tmp = a;

src/core/algorithm/link_cell.hpp

@@ -40,7 +40,9 @@ void link_cell(CellIterator first, CellIterator last,
       for (auto jt = std::next(it); jt != local_particles.end(); ++jt) {
         pair_kernel(p1, *jt);
       }
-
+    }
+    for (auto it = local_particles.begin(); it != local_particles.end(); ++it) {
+      auto &p1 = *it;
       /* Pairs with neighbors */
       for (auto &neighbor : cell->neighbors().red()) {
         for (auto &p2 : neighbor->particles()) {

src/core/cell_system/CellStructure.hpp

@@ -115,8 +115,7 @@ namespace detail {
 // NOLINTNEXTLINE(bugprone-exception-escape)
 struct MinimalImageDistance {
   BoxGeometry const box;
-
-  Distance operator()(Particle const &p1, Particle const &p2) const {
+  inline Distance operator()(Particle const &p1, Particle const &p2) const {
     return Distance(box.get_mi_vector(p1.pos(), p2.pos()));
   }
 };

src/core/constraints/ShapeBasedConstraint.cpp

@@ -29,6 +29,9 @@
 #include "errorhandling.hpp"
 #include "forces_inline.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
+#ifdef GAY_BERNE
+#include "nonbonded_interactions/gay_berne.hpp"
+#endif
 #include "system/System.hpp"
 #include "thermostat.hpp"
 
@@ -111,7 +114,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
       pf.f += thole_pair_force(p, part_rep, ia_params, dist_vec, dist,
                                *system.bonded_ias, get_ptr(coulomb_kernel));
 #endif
-      pf += calc_non_central_force(p, part_rep, ia_params, dist_vec, dist);
+#ifdef GAY_BERNE
+      if (gb_active(ia_params)) {
+        pf +=
+            gb_pair_force(p.quat(), part_rep.quat(), ia_params, dist_vec, dist);
+      }
+#endif
 
 #ifdef DPD
       if (system.thermostat->thermo_switch & THERMO_DPD) {
@@ -130,7 +138,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
         pf.f += thole_pair_force(p, part_rep, ia_params, dist_vec, -dist,
                                  *system.bonded_ias, get_ptr(coulomb_kernel));
 #endif
-        pf += calc_non_central_force(p, part_rep, ia_params, dist_vec, -dist);
+#ifdef GAY_BERNE
+        if (gb_active(ia_params)) {
+          pf += gb_pair_force(p.quat(), part_rep.quat(), ia_params, dist_vec,
+                              -dist);
+        }
+#endif
 
 #ifdef DPD
         if (system.thermostat->thermo_switch & THERMO_DPD) {

src/core/forces_inline.hpp

@@ -67,6 +67,9 @@
 #include "errorhandling.hpp"
 #include "exclusions.hpp"
 #include "thermostat.hpp"
+#ifdef NPT
+#include "npt.hpp"
+#endif
 
 #include <utils/Vector.hpp>
 
@@ -140,20 +143,6 @@ inline Utils::Vector3d calc_central_radial_force(IA_parameters const &ia_params,
   return force_factor * d;
 }
 
-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
-  return pf;
-}
-
 inline void apply_opposing_force(ParticleForce &pf, Utils::Vector3d const &d) {
 #ifdef ROTATION
   // if torque is a null vector, the opposing torque is a null vector too
@@ -203,11 +192,11 @@ inline void add_non_bonded_pair_force(
     if (do_nonbonded(p1, p2)) {
 #endif
       pf.f += calc_central_radial_force(ia_params, d, dist);
-#ifdef THOLE
-      pf.f += thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias,
-                               coulomb_kernel);
+#ifdef GAY_BERNE
+      if (gb_active(ia_params)) {
+        pf += gb_pair_force(p1.quat(), p2.quat(), ia_params, d, dist);
+      }
 #endif
-      pf += calc_non_central_force(p1, p2, ia_params, d, dist);
 #ifdef EXCLUSIONS
     }
 #endif
@@ -219,7 +208,9 @@ inline void add_non_bonded_pair_force(
   /* electrostatic is calculated by energy                             */
   /*********************************************************************/
 #ifdef NPT
-  npt_add_virial_force_contribution(pf.f, d);
+  if (npt_active()) {
+    npt_add_virial_force_contribution(pf.f, d);
+  }
 #endif
 
   /***********************************************/
@@ -231,9 +222,15 @@ inline void add_non_bonded_pair_force(
   auto const q1q2 = p1.q() * p2.q();
   if (q1q2 != 0. and coulomb_kernel != nullptr) {
     pf.f += (*coulomb_kernel)(q1q2, d, dist);
+#ifdef THOLE
+    pf.f += thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias,
+                             coulomb_kernel);
+#endif
 #ifdef NPT
-    npt_add_virial_diagonalSum_contribution(
-        (*coulomb_u_kernel)(p1, p2, q1q2, d, dist));
+    if (npt_active()) {
+      npt_add_virial_diagonalSum_contribution(
+          (*coulomb_u_kernel)(p1, p2, q1q2, d, dist));
+    }
 #endif
 #ifdef P3M
     if (elc_kernel)
@@ -345,7 +342,9 @@ inline bool add_bonded_two_body_force(
       p2.force() -= result.value();
 
 #ifdef NPT
-      npt_add_virial_force_contribution(result.value(), dx);
+      if (npt_active()) {
+        npt_add_virial_force_contribution(result.value(), dx);
+      }
 #endif
       return false;
     }

src/core/nonbonded_interactions/gay_berne.hpp

@@ -45,6 +45,11 @@
 
 #include <cmath>
 
+/** Check if Gay-Berne interaction is active */
+inline bool gb_active(IA_parameters const &ia_params) {
+  return ia_params.gay_berne.cut != INACTIVE_CUTOFF;
+}
+
 /** Calculate Gay-Berne force and torques */
 inline ParticleForce gb_pair_force(Utils::Quaternion<double> const &qi,
                                    Utils::Quaternion<double> const &qj,

src/core/npt.cpp

@@ -148,4 +148,5 @@ void System::System::npt_add_virial_contribution(Utils::Vector3d const &force,
     npt_inst_pressure->p_vir += hadamard_product(force, d);
   }
 }
+
 #endif // NPT

src/core/npt.hpp

@@ -34,6 +34,10 @@
 #include <cstddef>
 #include <vector>
 
+#include "PropagationMode.hpp"
+#include "integrators/Propagation.hpp"
+#include "system/System.hpp"
+
 namespace System {
 class System;
 } // namespace System
@@ -99,4 +103,11 @@ struct InstantaneousPressure {
   Utils::Vector3d p_vel = {0., 0., 0.};
 };
 
+/** Check if NPT integration is active. */
+inline bool npt_active() {
+  auto const &system = ::System::get_system();
+  return (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) ||
+         (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO_MTK);
+}
+
 #endif // NPT

src/core/pressure_inline.hpp

@@ -33,6 +33,9 @@
 #include "errorhandling.hpp"
 #include "exclusions.hpp"
 #include "forces_inline.hpp"
+#ifdef GAY_BERNE
+#include "nonbonded_interactions/gay_berne.hpp"
+#endif
 
 #include <utils/Vector.hpp>
 #include <utils/math/tensor_product.hpp>
@@ -67,12 +70,16 @@ inline void add_non_bonded_pair_virials(
   if (do_nonbonded(p1, p2))
 #endif
   {
-    auto const force = calc_central_radial_force(ia_params, d, dist) +
+    auto force = calc_central_radial_force(ia_params, d, dist);
 #ifdef THOLE
-                       thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias,
-                                        kernel_forces) +
+    force +=
+        thole_pair_force(p1, p2, ia_params, d, dist, bonded_ias, kernel_forces);
+#endif
+#ifdef GAY_BERNE
+    if (gb_active(ia_params)) {
+      force += gb_pair_force(p1.quat(), p2.quat(), ia_params, d, dist).f;
+    }
 #endif
-                       calc_non_central_force(p1, p2, ia_params, d, dist).f;
     auto const stress = Utils::tensor_product(d, force);
     obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(),
                                              p2.mol_id(), flatten(stress));