diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 6e978cd0961..3bbbc1b26d4 100644 --- a/src/core/particle_data.cpp +++ b/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 diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 1c7b036983f..8e0810975df 100644 --- a/testsuite/python/CMakeLists.txt +++ b/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) diff --git a/testsuite/python/auto_exclusions.py b/testsuite/python/auto_exclusions.py index 826864dd4ea..d24d6e5de1b 100644 --- a/testsuite/python/auto_exclusions.py +++ b/testsuite/python/auto_exclusions.py @@ -17,64 +17,212 @@ # along with this program. If not, see . # +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__":