Fix regressions in auto-exclusions (espressomd#4654)

The bond length is no longer misinterpreted as a bond id
(fixes bug introduced by a48ed7c).

Fixes espressomd#4652

Description of changes:
- the feature now works with non-contiguous particle ids
- the feature no longer adds spurious exclusions to particle ids in the range `[1, distance]`
- completely rewrite the test case to check all topologies and regular decompositions

src/core/particle_data.cpp

@@ -812,11 +812,13 @@ void auto_exclusions(int distance) {
      setting the bonds, which the user apparently accepted.
   */
   for (auto &kv : partners) {
-    auto const id = kv.first;
+    auto const pid1 = kv.first;
     auto const partner_list = kv.second;
-    for (int j : partner_list)
-      if (id < j)
-        add_particle_exclusion(id, j);
+    for (int j = 0; j < partner_list.size(); j += 2) {
+      auto const pid2 = partner_list[j];
+      if (pid1 < pid2)
+        add_particle_exclusion(pid1, pid2);
+    }
   }
 }
 #endif // EXCLUSIONS

testsuite/python/CMakeLists.txt

@@ -238,7 +238,7 @@ python_test(FILE lj.py MAX_NUM_PROC 4)
 python_test(FILE pairs.py MAX_NUM_PROC 4)
 python_test(FILE polymer_linear.py MAX_NUM_PROC 4)
 python_test(FILE polymer_diamond.py MAX_NUM_PROC 4)
-python_test(FILE auto_exclusions.py MAX_NUM_PROC 1)
+python_test(FILE auto_exclusions.py MAX_NUM_PROC 4)
 python_test(FILE observable_cylindrical.py MAX_NUM_PROC 4)
 python_test(FILE observable_cylindricalLB.py MAX_NUM_PROC 2 LABELS gpu)
 python_test(FILE analyze_chains.py MAX_NUM_PROC 1)

testsuite/python/auto_exclusions.py

@@ -17,64 +17,212 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
+import numpy as np
 import unittest as ut
 import unittest_decorators as utx
 import espressomd
 
 
 @utx.skipIfMissingFeatures("EXCLUSIONS")
-class AutoExclusions(ut.TestCase):
+class Test(ut.TestCase):
+    """
+    Check the auto-exclusions feature against various polymer topologies.
+    Verify the exclusion list is still correct when the polymer spans
+    three contiguous, yet different, cells, such that particles at one
+    end of the polymer don't have access to the particles at the other
+    end of the polymer via the ghost layer.
+    """
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system.cell_system.skin = 0.1
+    node_grid_ref = list(system.cell_system.node_grid)
 
-    def setUp(self):
+    def tearDown(self):
         self.system.part.clear()
+        self.system.cell_system.node_grid = self.node_grid_ref
+
+    def set_particles_on_cube(self):
+        # place particles on a cube centered at the origin
+        length = self.system.cell_system.skin / 2.
+        self.system.part.add(pos=length * (np.array([0, 0, 0]) - 1))
+        self.system.part.add(pos=length * (np.array([2, 0, 0]) - 1))
+        self.system.part.add(pos=length * (np.array([2, 2, 0]) - 1))
+        self.system.part.add(pos=length * (np.array([0, 2, 0]) - 1))
+        self.system.part.add(pos=length * (np.array([0, 2, 2]) - 1))
+        self.system.part.add(pos=length * (np.array([2, 2, 2]) - 1))
+        self.system.part.add(pos=length * (np.array([2, 0, 2]) - 1))
+        self.system.part.add(pos=length * (np.array([0, 0, 2]) - 1))
+
+        # particles should be equally distributed on all nodes
+        n_nodes = np.prod(self.node_grid_ref)
+        if n_nodes in (1, 2, 4, 8):
+            p_nodes = sorted(list(self.system.part.all().node))
+            p_nodes_ref = list(np.repeat(np.arange(n_nodes), 8 // n_nodes))
+            assert p_nodes == p_nodes_ref
+
+    def set_particles_on_line(self, length):
+        box_l_z = self.system.box_l[2]
+        for i in range(length):
+            self.system.part.add(pos=np.array([0, 0, box_l_z]) * i / length)
+
+        # particles should be distributed on multiple nodes
+        n_nodes = np.prod(self.node_grid_ref)
+        if n_nodes > 1:
+            assert len(set(self.system.part.all().node)) > 1
 
     def test_linear(self):
         bond = espressomd.interactions.Virtual()
-        s = self.system
-        s.bonded_inter.add(bond)
+        system = self.system
+        system.bonded_inter.add(bond)
 
-        for i in range(10):
-            s.part.add(id=i, pos=[0, 0, 0])
+        def check():
+            # topology: 0-1-2-...-n
+            length = len(system.part)
+            for pid in range(length - 1):
+                system.part.by_id(pid).add_bond((bond, pid + 1))
+            system.auto_exclusions(1)
 
-        for i in range(9):
-            s.part.by_id(i).add_bond((bond, i + 1))
+            for pid in range(1, length - 1):
+                excl = sorted(list(system.part.by_id(pid).exclusions))
+                self.assertEqual(excl, [pid - 1, pid + 1])
 
-        s.auto_exclusions(1)
+            excl = list(system.part.by_id(0).exclusions)
+            self.assertEqual(excl, [1])
 
-        for p in range(1, 9):
-            excl = s.part.by_id(p).exclusions
-            self.assertEqual(len(excl), 2)
-            self.assertIn(p - 1, excl)
-            self.assertIn(p + 1, excl)
+            excl = list(system.part.by_id(length - 1).exclusions)
+            self.assertEqual(excl, [length - 2])
 
-        excl = s.part.by_id(0).exclusions
-        self.assertEqual(len(excl), 1)
-        self.assertIn(1, excl)
+        with self.subTest(msg='cube'):
+            self.set_particles_on_cube()
+            check()
 
-        excl = s.part.by_id(9).exclusions
-        self.assertEqual(len(excl), 1)
-        self.assertIn(8, excl)
+        self.tearDown()
+
+        with self.subTest(msg='line'):
+            system.cell_system.node_grid = [1, 1, np.prod(self.node_grid_ref)]
+            self.set_particles_on_line(16)
+            check()
 
     def test_ring(self):
         bond = espressomd.interactions.Virtual()
-        s = self.system
-        s.bonded_inter.add(bond)
+        system = self.system
+        system.bonded_inter.add(bond)
+
+        def check():
+            # topology: 0-1-2-...-n-0
+            length = len(system.part)
+            for pid in range(length):
+                system.part.by_id(pid).add_bond((bond, (pid + 1) % length))
+            system.auto_exclusions(2)
+
+            for pid in range(length):
+                excl = sorted(list(system.part.by_id(pid).exclusions))
+                excl_ref = np.mod([pid - 1, pid - 2, pid + 1, pid + 2], length)
+                self.assertEqual(excl, sorted(list(excl_ref)))
 
-        for i in range(10):
-            s.part.add(id=i, pos=[0, 0, 0])
+        with self.subTest(msg='cube'):
+            self.set_particles_on_cube()
+            check()
 
-        for i in range(10):
-            s.part.by_id(i).add_bond((bond, (i + 1) % 10))
+        self.tearDown()
 
-        s.auto_exclusions(2)
+        with self.subTest(msg='line'):
+            system.cell_system.node_grid = [1, 1, np.prod(self.node_grid_ref)]
+            self.set_particles_on_line(16)
+            check()
 
-        for p in range(10):
-            excl = s.part.by_id(p).exclusions
-            self.assertEqual(len(excl), 4)
-            for i in range(1, 3):
-                self.assertIn((p - i) % 10, excl)
-                self.assertIn((p + i) % 10, excl)
+    def test_branched(self):
+        bond = espressomd.interactions.Virtual()
+        system = self.system
+        system.bonded_inter.add(bond)
+
+        length = system.cell_system.skin / 2.
+        p0 = system.part.add(pos=length * (np.array([0, 0, 0]) - 1))
+        p1 = system.part.add(pos=length * (np.array([2, 0, 0]) - 1))
+        p2 = system.part.add(pos=length * (np.array([0, 2, 0]) - 1))
+        p3 = system.part.add(pos=length * (np.array([2, 2, 0]) - 1))
+
+        # topology: 0-1(-2)-3
+        p1.add_bond((bond, p0))
+        p2.add_bond((bond, p1))
+        p3.add_bond((bond, p1))
+        system.auto_exclusions(2)
+
+        self.assertEqual(sorted(list(p0.exclusions)), [1, 2, 3])
+        self.assertEqual(sorted(list(p1.exclusions)), [0, 2, 3])
+        self.assertEqual(sorted(list(p2.exclusions)), [0, 1, 3])
+        self.assertEqual(sorted(list(p3.exclusions)), [0, 1, 2])
+
+    def test_diamond(self):
+        bond = espressomd.interactions.Virtual()
+        system = self.system
+        system.bonded_inter.add(bond)
+
+        length = system.cell_system.skin / 2.
+        p0 = system.part.add(pos=length * (np.array([0, 0, 2]) - 1))
+        p1 = system.part.add(pos=length * (np.array([0, 0, 0]) - 1))
+        p2 = system.part.add(pos=length * (np.array([2, 0, 0]) - 1))
+        p3 = system.part.add(pos=length * (np.array([0, 2, 0]) - 1))
+        p4 = system.part.add(pos=length * (np.array([2, 2, 0]) - 1))
+        p5 = system.part.add(pos=length * (np.array([2, 2, 2]) - 1))
+
+        # topology: 0-1(-2-4-5)-3-4-5
+        p1.add_bond((bond, p0))
+        p2.add_bond((bond, p1))
+        p3.add_bond((bond, p1))
+        p2.add_bond((bond, p4))
+        p3.add_bond((bond, p4))
+        p5.add_bond((bond, p4))
+        system.auto_exclusions(3)
+
+        self.assertEqual(sorted(list(p0.exclusions)), [1, 2, 3, 4])
+        self.assertEqual(sorted(list(p1.exclusions)), [0, 2, 3, 4, 5])
+        self.assertEqual(sorted(list(p2.exclusions)), [0, 1, 3, 4, 5])
+        self.assertEqual(sorted(list(p3.exclusions)), [0, 1, 2, 4, 5])
+        self.assertEqual(sorted(list(p4.exclusions)), [0, 1, 2, 3, 5])
+        self.assertEqual(sorted(list(p5.exclusions)), [1, 2, 3, 4])
+
+    def test_id_gaps(self):
+        bond = espressomd.interactions.Virtual()
+        system = self.system
+        system.bonded_inter.add(bond)
+
+        p0 = system.part.add(id=0, pos=[0.1, 0.1, 0.1])
+        p4 = system.part.add(id=4, pos=[0.2, 0.1, 0.1])
+
+        # topology: 0-4
+        p0.add_bond((bond, p4))
+        system.auto_exclusions(1)
+
+        self.assertEqual(list(p0.exclusions), [4])
+        self.assertEqual(list(p4.exclusions), [0])
+
+    def test_non_pairwise(self):
+        """
+        Check that bonds with 0, 2 and 3 partners don't generate exclusions.
+        """
+        angle = espressomd.interactions.AngleHarmonic(bend=1., phi0=2.)
+        dihe = espressomd.interactions.Dihedral(bend=1., mult=2, phase=2.)
+        if espressomd.has_features(["VIRTUAL_SITES_INERTIALESS_TRACERS"]):
+            volcons = espressomd.interactions.IBM_VolCons(softID=1, kappaV=1.)
+        system = self.system
+        system.bonded_inter.add(angle)
+        system.bonded_inter.add(dihe)
+        if espressomd.has_features(["VIRTUAL_SITES_INERTIALESS_TRACERS"]):
+            system.bonded_inter.add(volcons)
+
+        self.system.part.add(pos=[0.1, 0.1, 0.1])
+        self.system.part.add(pos=[0.2, 0.1, 0.1])
+        self.system.part.add(pos=[0.2, 0.2, 0.1])
+        self.system.part.add(pos=[0.1, 0.2, 0.1])
+        system.part.by_id(0).add_bond((angle, 1, 2))
+        system.part.by_id(0).add_bond((dihe, 1, 2, 3))
+        if espressomd.has_features(["VIRTUAL_SITES_INERTIALESS_TRACERS"]):
+            system.part.by_id(0).add_bond((volcons,))
+
+        system.auto_exclusions(2)
+
+        for p in system.part.all():
+            self.assertEqual(list(p.exclusions), [])
 
 
 if __name__ == "__main__":