WIP add cutoff by type

src/core/collision.cpp

@@ -671,3 +671,36 @@ void handle_collisions(CellStructure &cell_structure) {
 }
 
 #endif // COLLISION_DETECTION
+
+double collision_detection_cutoff() {
+#ifdef COLLISION_DETECTION
+  auto &system = System::get_system();
+  auto &nonbonded_ias = *system.nonbonded_ias;
+  std::vector<int> v(nonbonded_ias.get_max_seen_particle_type() + 1);
+  std::iota(std::begin(v), std::end(v), 0);
+  return collision_detection_cutoff(v);
+#else
+  return -1;
+#endif
+}
+
+double collision_detection_cutoff(std::vector<int> types) {
+  double ret = -1;
+#ifdef COLLISION_DETECTION
+  if (collision_params.mode == CollisionModeType::BIND_CENTERS or
+      collision_params.mode == CollisionModeType::BIND_VS or
+      collision_params.mode == CollisionModeType::BIND_THREE_PARTICLES) {
+    ret = collision_params.distance;
+  }
+  if (collision_params.mode == CollisionModeType::GLUE_TO_SURF) {
+    if ((std::find(types.begin(), types.end(),
+                   collision_params.part_type_to_be_glued) != types.end()) and
+        (std::find(types.begin(), types.end(),
+                   collision_params.part_type_to_attach_vs_to) !=
+         types.end())) {
+      ret = collision_params.distance;
+    }
+  }
+#endif
+  return ret;
+}

src/core/collision.hpp

@@ -150,11 +150,5 @@ inline void detect_collision(Particle const &p1, Particle const &p2,
 }
 
 #endif // COLLISION_DETECTION
-
-inline double collision_detection_cutoff() {
-#ifdef COLLISION_DETECTION
-  if (collision_params.mode != CollisionModeType::OFF)
-    return collision_params.distance;
-#endif
-  return -1.;
-}
+double collision_detection_cutoff();
+double collision_detection_cutoff(std::vector<int> types);

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -436,4 +436,57 @@ class InteractionsNonBonded {
 
   /** @brief Get maximal cutoff. */
   double maximal_cutoff() const;
+
+  double maximal_cutoff(std::vector<int> types) const {
+    auto max_cut_nonbonded = INACTIVE_CUTOFF;
+    auto const n_types = static_cast<int>(types.size());
+    for (int i = 0; i < n_types; i++) {
+      for (int j = i; j < n_types; j++) {
+        if (types[i] > max_seen_particle_type or
+            types[j] > max_seen_particle_type) {
+          continue;
+        }
+        auto const data = get_ia_param(types[i], types[j]);
+        max_cut_nonbonded = std::max(max_cut_nonbonded, data.max_cut);
+      }
+    }
+    return max_cut_nonbonded;
+  }
+
+  double maximal_cutoff_all_except(std::vector<int> types) const {
+    auto max_cut_nonbonded = INACTIVE_CUTOFF;
+    for (int i = 0; i < max_seen_particle_type; i++) {
+      if (std::find(types.begin(), types.end(), i) != types.end()) {
+        // i in types
+        continue;
+      }
+      for (int j = i; j < max_seen_particle_type; j++) {
+        if (std::find(types.begin(), types.end(), j) != types.end()) {
+          // j in types
+          continue;
+        }
+        auto const data = get_ia_param(i, j);
+        max_cut_nonbonded = std::max(max_cut_nonbonded, data.max_cut);
+      }
+    }
+    return max_cut_nonbonded;
+  }
+
+  bool only_central_forces(std::vector<int> types) const {
+    // gay-berne is the only non-central force
+    bool ret = true;
+#ifdef GAY_BERNE
+    auto const n_types = static_cast<int>(types.size());
+    for (int i = 0; i < n_types; i++) {
+      for (int j = i; j < n_types; j++) {
+        auto const data = get_ia_param(i, j);
+        if (data.gay_berne.cut != INACTIVE_CUTOFF) {
+          ret = false;
+          break;
+        }
+      }
+    }
+#endif
+    return ret;
+  }
 };

src/core/system/System.cpp

@@ -358,6 +358,12 @@ void System::update_local_geo() {
 }
 
 double System::maximal_cutoff() const {
+  std::vector<int> v(nonbonded_ias->get_max_seen_particle_type() + 1);
+  std::iota(std::begin(v), std::end(v), 0);
+  return System::maximal_cutoff(v);
+}
+
+double System::maximal_cutoff(std::vector<int> &types) const {
   auto max_cut = INACTIVE_CUTOFF;
   max_cut = std::max(max_cut, get_min_global_cut());
 #ifdef ELECTROSTATICS
@@ -371,9 +377,9 @@ double System::maximal_cutoff() const {
     // because bond partners are always on the local node.
     max_cut = std::max(max_cut, maximal_cutoff_bonded());
   }
-  max_cut = std::max(max_cut, nonbonded_ias->maximal_cutoff());
+  max_cut = std::max(max_cut, nonbonded_ias->maximal_cutoff(types));
 
-  max_cut = std::max(max_cut, collision_detection_cutoff());
+  max_cut = std::max(max_cut, collision_detection_cutoff(types));
   return max_cut;
 }
 

src/core/system/System.hpp

@@ -119,6 +119,7 @@ class System : public std::enable_shared_from_this<System> {
   void rebuild_cell_structure();
 
   /** @brief Calculate the maximal cutoff of all interactions. */
+  double maximal_cutoff(std::vector<int> &types) const;
   double maximal_cutoff() const;
 
   /** @brief Get the interaction range. */

src/python/espressomd/system.py

@@ -439,6 +439,9 @@ def velocity_difference(self, p1, p2):
         return self.call_method("velocity_difference", pos1=p1.pos_folded,
                                 pos2=p2.pos_folded, v1=p1.v, v2=p2.v)
 
+    def cutoff_by_types(self, types):
+        return self.call_method("cutoff_by_types", types=types)
+
     def auto_exclusions(self, distance):
         """
         Add exclusions between particles that are bonded.

src/script_interface/system/System.cpp

@@ -310,6 +310,12 @@ Variant System::do_call_method(std::string const &name,
     auto const pos2 = get_value<Utils::Vector3d>(parameters, "pos2");
     return m_instance->box_geo->get_mi_vector(pos2, pos1);
   }
+
+  if (name == "cutoff_by_types") {
+    auto types = get_value<std::vector<int>>(parameters, "types");
+    return m_instance->maximal_cutoff(types);
+  }
+
   if (name == "rotate_system") {
     rotate_system(*m_instance->cell_structure,
                   get_value<double>(parameters, "phi"),

testsuite/python/cutoffs.py

@@ -20,16 +20,25 @@
 from espressomd.interactions import FeneBond
 import numpy as np
 import unittest as ut
+import unittest_decorators as utx
 
 
 class CutOff(ut.TestCase):
 
     """
-    Test interaction cutoffs
+    Test interaction cutoffs.
+    Tests must be executed in the order given below such 
+    that the leftovers from test 2 do not break test 1.
     """
+    system = espressomd.System(box_l=3 * [50])
 
-    def test(self):
-        system = espressomd.System(box_l=[15.0, 15.0, 15.0])
+    def tearDown(self) -> None:
+        self.system.non_bonded_inter.reset()
+        self.system.bonded_inter.clear()
+        self.system.part.clear()
+
+    def test_1max_cut(self):
+        system = self.system
         system.cell_system.skin = 1
 
         # Initial state. Skin does not influence cutoffs as long as there are
@@ -68,6 +77,65 @@ def test(self):
             self.assertEqual(system.cell_system.max_cut_nonbonded, -1)
             self.assertEqual(system.cell_system.interaction_range, -1)
 
+    @utx.skipIfMissingFeatures(["LENNARD_JONES",
+                               "VIRTUAL_SITES_RELATIVE", "COLLISION_DETECTION", "DP3M"])
+    def test_2cutoff_by_type(self):
+        sys = self.system
+        min_global_cut = 0.1
+        sys.min_global_cut = min_global_cut
+        sys.non_bonded_inter[0, 0].lennard_jones.set_params(
+            sigma=0.1, epsilon=1, cutoff=1, shift="auto")
+        sys.non_bonded_inter[0, 1].lennard_jones.set_params(
+            sigma=0.1, epsilon=1, cutoff=4, shift="auto")
+        sys.non_bonded_inter[0, 2].lennard_jones.set_params(
+            sigma=0.1, epsilon=1, cutoff=3, shift="auto")
+
+        self.assertEqual(sys.max_cut_nonbonded, 4)
+        self.assertEqual(sys.cutoff_by_types([0]), 1)
+        self.assertEqual(sys.cutoff_by_types([0, 1, 2]), 4)
+        self.assertEqual(sys.cutoff_by_types([0, 2, 100]), 3)
+        self.assertEqual(sys.cutoff_by_types([1, 2]), min_global_cut)
+        self.assertEqual(sys.cutoff_by_types(
+            [9, 10, 11, 23, 1456]), min_global_cut)
+
+        sys.non_bonded_inter[0, 2].lennard_jones.set_params(
+            sigma=0.1, epsilon=1, cutoff=2, shift="auto")
+        self.assertEqual(sys.cutoff_by_types([0, 2]), 2)
+        bond = espressomd.interactions.HarmonicBond(k=1000, r_0=0)
+        sys.bonded_inter.add(bond)
+        sys.collision_detection.set_params(
+            mode="glue_to_surface",
+            distance=5.5,
+            distance_glued_particle_to_vs=0.02,
+            bond_centers=bond,
+            bond_vs=bond,
+            part_type_vs=1,
+            part_type_to_attach_vs_to=2,
+            part_type_to_be_glued=7,
+            part_type_after_glueing=8)
+        self.assertEqual(sys.cutoff_by_types([2, 7]), 5.5)
+        self.assertEqual(sys.cutoff_by_types([1, 7]), min_global_cut)
+        self.assertEqual(sys.cutoff_by_types([7, 8]), min_global_cut)
+
+        sys.collision_detection.set_params(mode="bind_centers", distance=6.6,
+                                           bond_centers=bond)
+        self.assertEqual(sys.cutoff_by_types([1, 2]), 6.6)
+        self.assertEqual(sys.cutoff_by_types([1, 7]), 6.6)
+        self.assertEqual(sys.cutoff_by_types([7, 8]), 6.6)
+        import espressomd.magnetostatics as magnetostatics
+        p3m = magnetostatics.DipolarP3M(
+            prefactor=1,
+            mesh=32,
+            alpha=0.123,
+            accuracy=1E-4,
+            r_cut=7.7,
+            cao=3,
+            tune=False)
+        sys.part.add(pos=50 * np.random.random((10, 3)),
+                     dip=np.random.random((10, 3)))
+        sys.magnetostatics.solver = p3m
+        self.assertEqual(sys.cutoff_by_types([1, 2]), 7.7)
+
 
 if __name__ == '__main__':
     ut.main()

testsuite/python/regular_decomposition.py

@@ -35,6 +35,7 @@ def setUp(self):
 
     def tearDown(self):
         self.system.part.clear()
+        self.system.non_bonded_inter.reset()
 
     def check_resort(self):
         n_part = 2351
@@ -103,7 +104,7 @@ def test_position_rounding(self):
         self.system.cell_system.skin = 0.4
         self.system.min_global_cut = 12.0 / 4.25
         self.system.part.add(pos=[25, 25, 0])
-        self.assertEqual(1, len(self.system.part))
+        self.assertEqual(1, len(self.system.part))        
 
 
 if __name__ == "__main__":