diff --git a/src/core/collision.cpp b/src/core/collision.cpp index d52bfa747c..c84e215f99 100644 --- a/src/core/collision.cpp +++ b/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 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 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; +} \ No newline at end of file diff --git a/src/core/collision.hpp b/src/core/collision.hpp index 339e47f81a..b921430fa9 100644 --- a/src/core/collision.hpp +++ b/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 types); diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index 29474c56c2..32f4d66937 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/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 types) const { + auto max_cut_nonbonded = INACTIVE_CUTOFF; + auto const n_types = static_cast(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 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 types) const { + // gay-berne is the only non-central force + bool ret = true; +#ifdef GAY_BERNE + auto const n_types = static_cast(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; + } }; diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index 0c2462961e..2f3881a5a1 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -358,6 +358,12 @@ void System::update_local_geo() { } double System::maximal_cutoff() const { + std::vector 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 &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; } diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index a1d46118ec..e236dd2edc 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -119,6 +119,7 @@ class System : public std::enable_shared_from_this { void rebuild_cell_structure(); /** @brief Calculate the maximal cutoff of all interactions. */ + double maximal_cutoff(std::vector &types) const; double maximal_cutoff() const; /** @brief Get the interaction range. */ diff --git a/src/python/espressomd/system.py b/src/python/espressomd/system.py index b4d0dd5f18..2d43b6adb5 100644 --- a/src/python/espressomd/system.py +++ b/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. diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index 2d5ce58a8f..b51cde6bdb 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -310,6 +310,12 @@ Variant System::do_call_method(std::string const &name, auto const pos2 = get_value(parameters, "pos2"); return m_instance->box_geo->get_mi_vector(pos2, pos1); } + + if (name == "cutoff_by_types") { + auto types = get_value>(parameters, "types"); + return m_instance->maximal_cutoff(types); + } + if (name == "rotate_system") { rotate_system(*m_instance->cell_structure, get_value(parameters, "phi"), diff --git a/testsuite/python/cutoffs.py b/testsuite/python/cutoffs.py index 51c4ca0522..af2dd7ba8a 100644 --- a/testsuite/python/cutoffs.py +++ b/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() diff --git a/testsuite/python/regular_decomposition.py b/testsuite/python/regular_decomposition.py index b05746f542..313b32e801 100644 --- a/testsuite/python/regular_decomposition.py +++ b/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__":