diff --git a/src/core/cell_system/CellStructure.cpp b/src/core/cell_system/CellStructure.cpp index 773d474eb2..91c5bee4a7 100644 --- a/src/core/cell_system/CellStructure.cpp +++ b/src/core/cell_system/CellStructure.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -248,12 +249,13 @@ void CellStructure::set_atom_decomposition() { system.on_cell_structure_change(); } -void CellStructure::set_regular_decomposition(double range) { +void CellStructure::set_regular_decomposition( + double range, std::optional> fully_connected_boundary) { auto &system = get_system(); auto &local_geo = *system.local_geo; auto const &box_geo = *system.box_geo; set_particle_decomposition(std::make_unique( - ::comm_cart, range, box_geo, local_geo)); + ::comm_cart, range, box_geo, local_geo, fully_connected_boundary)); m_type = CellStructureType::REGULAR; local_geo.set_cell_structure_type(m_type); system.on_cell_structure_change(); diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp index 005b952a22..9993fffc3d 100644 --- a/src/core/cell_system/CellStructure.hpp +++ b/src/core/cell_system/CellStructure.hpp @@ -46,6 +46,7 @@ #include #include #include +#include #include #include #include @@ -557,7 +558,9 @@ struct CellStructure : public System::Leaf { * * @param range Interaction range. */ - void set_regular_decomposition(double range); + void set_regular_decomposition( + double range, + std::optional> fully_connected_boundary); /** * @brief Set the particle decomposition to @ref HybridDecomposition. diff --git a/src/core/cell_system/HybridDecomposition.cpp b/src/core/cell_system/HybridDecomposition.cpp index a54002f8df..4c21fe6475 100644 --- a/src/core/cell_system/HybridDecomposition.cpp +++ b/src/core/cell_system/HybridDecomposition.cpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -48,7 +49,7 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm, std::set n_square_types) : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular), m_regular_decomposition(RegularDecomposition( - m_comm, cutoff_regular + skin, m_box, local_box)), + m_comm, cutoff_regular + skin, m_box, local_box, std::nullopt)), m_n_square(AtomDecomposition(m_comm, m_box)), m_n_square_types(std::move(n_square_types)), m_get_global_ghost_flags(std::move(get_ghost_flags)) { diff --git a/src/core/cell_system/RegularDecomposition.cpp b/src/core/cell_system/RegularDecomposition.cpp index 51158c580b..df48efd35c 100644 --- a/src/core/cell_system/RegularDecomposition.cpp +++ b/src/core/cell_system/RegularDecomposition.cpp @@ -397,6 +397,27 @@ void RegularDecomposition::init_cell_interactions() { auto const &node_grid = ::communicator.node_grid; auto const global_halo_offset = hadamard_product(node_pos, cell_grid) - halo; auto const global_size = hadamard_product(node_grid, cell_grid); + auto const at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) { + return (cell_idx[coord] == 0 or cell_idx[coord] == global_size[coord]); + }; + + // For the fully connected feature (cells that don't share at least a corner) + // only apply if one cell is a ghost cell (i.e. connections across the + // periodic boundary. + auto const fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a, + Utils::Vector3i b) { + if (fully_connected_boundary()) { + auto const [fc_normal, fc_dir] = *fully_connected_boundary(); + auto const involves_ghost_cell = + (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or + b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]); + if (not involves_ghost_cell) { + // check if cells do not share at least a corner + return std::abs((a - b)[fc_dir]) > 1; + } + } + return false; + }; /* Translate a node local index (relative to the origin of the local grid) * to a global index. */ @@ -418,6 +439,19 @@ void RegularDecomposition::init_cell_interactions() { return (global_index - global_halo_offset); }; + // sanity checks + if (fully_connected_boundary()) { + auto const [fc_normal, fc_dir] = *fully_connected_boundary(); + if (fc_normal == fc_dir) { + throw std::domain_error("fully_connected_boundary normal and connection " + "coordinates need to differ."); + } + if (node_grid[fc_dir] != 1) { + throw std::runtime_error( + "The MPI nodegrid must be 1 in the fully connected direction."); + } + } + /* We only consider local cells (e.g. not halo cells), which * span the range [(1,1,1), cell_grid) in local coordinates. */ auto const start = global_index(Utils::Vector3i{1, 1, 1}); @@ -431,20 +465,18 @@ void RegularDecomposition::init_cell_interactions() { Utils::Vector3i lower_index = {m - 1, n - 1, o - 1}; Utils::Vector3i upper_index = {m + 1, n + 1, o + 1}; - // /* In the fully connected case, we consider all cells - // * in the direction as neighbors, not only the nearest ones. + /* In the fully connected case, we consider all cells + * in the direction as neighbors, not only the nearest ones. // */ - // for (int i = 0; i < 3; i++) { - // if (dd.fully_connected[i]) { - // // Fully connected is only needed at the box surface - // if (i==0 and (n!=start[1] or n!=end[1]-1) and (o!=start[2] - // or o!=end[2]-1)) continue; if (i==1 and (m!=start[0] or - // m!=end[0]-1) and (o!=start[2] or o!=end[2]-1)) continue; - // if (i==2 and (m!=start[0] or m!=end[0]-1) and (n!=start[1] - // or n!=end[1]-1)) continue; lower_index[i] = 0; - // upper_index[i] = global_size[i] - 1; - // } - // } + if (fully_connected_boundary()) { + auto const [fc_boundary, fc_direction] = *fully_connected_boundary(); + + // Fully connected is only needed at the box surface + if (at_boundary(fc_boundary, {m, n, o})) { + lower_index[fc_direction] = -1; + upper_index[fc_direction] = global_size[fc_boundary]; + } + } /* In non-periodic directions, the halo needs not * be considered. */ @@ -466,6 +498,12 @@ void RegularDecomposition::init_cell_interactions() { for (int p = lower_index[2]; p <= upper_index[2]; p++) for (int q = lower_index[1]; q <= upper_index[1]; q++) for (int r = lower_index[0]; r <= upper_index[0]; r++) { + if (fully_connected_boundary()) { + // Avoid fully connecting the boundary layer and the + // next INNER layer + if (fcb_is_inner_connection({m, n, o}, {r, q, p})) + continue; + } neighbors.insert(Utils::Vector3i{r, q, p}); } @@ -629,11 +667,13 @@ GhostCommunicator RegularDecomposition::prepare_comm() { return ghost_comm; } -RegularDecomposition::RegularDecomposition(boost::mpi::communicator comm, - double range, - BoxGeometry const &box_geo, - LocalBox const &local_geo) - : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo) { +RegularDecomposition::RegularDecomposition( + boost::mpi::communicator comm, double range, BoxGeometry const &box_geo, + LocalBox const &local_geo, + std::optional> fully_connected) + : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo), + m_fully_connected_boundary(std::move(fully_connected)) { + /* set up new regular decomposition cell structure */ create_cell_grid(range); diff --git a/src/core/cell_system/RegularDecomposition.hpp b/src/core/cell_system/RegularDecomposition.hpp index 2559a6f123..f79a71ff10 100644 --- a/src/core/cell_system/RegularDecomposition.hpp +++ b/src/core/cell_system/RegularDecomposition.hpp @@ -79,6 +79,7 @@ struct RegularDecomposition : public ParticleDecomposition { boost::mpi::communicator m_comm; BoxGeometry const &m_box; LocalBox m_local_box; + std::optional> m_fully_connected_boundary = {}; std::vector cells; std::vector m_local_cells; std::vector m_ghost_cells; @@ -87,7 +88,8 @@ struct RegularDecomposition : public ParticleDecomposition { public: RegularDecomposition(boost::mpi::communicator comm, double range, - BoxGeometry const &box_geo, LocalBox const &local_geo); + BoxGeometry const &box_geo, LocalBox const &local_geo, + std::optional> fully_connected); GhostCommunicator const &exchange_ghosts_comm() const override { return m_exchange_ghosts_comm; @@ -115,6 +117,8 @@ struct RegularDecomposition : public ParticleDecomposition { Utils::Vector3d max_cutoff() const override; Utils::Vector3d max_range() const override; + auto fully_connected_boundary() const { return m_fully_connected_boundary; } + std::optional minimum_image_distance() const override { return {m_box}; } diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index e9083ff749..676b824155 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -154,7 +154,17 @@ void System::set_min_global_cut(double value) { void System::set_cell_structure_topology(CellStructureType topology) { if (topology == CellStructureType::REGULAR) { - cell_structure->set_regular_decomposition(get_interaction_range()); + if (cell_structure->decomposition_type() == CellStructureType::REGULAR) { + // get fully connected info from exising regular decomposition + auto &old_regular_decomposition = + dynamic_cast( + std::as_const(*cell_structure).decomposition()); + cell_structure->set_regular_decomposition( + get_interaction_range(), + old_regular_decomposition.fully_connected_boundary()); + } else { // prev. decomposition is not a regular decomposition + cell_structure->set_regular_decomposition(get_interaction_range(), {}); + } } else if (topology == CellStructureType::NSQUARE) { cell_structure->set_atom_decomposition(); } else { diff --git a/src/python/espressomd/cell_system.py b/src/python/espressomd/cell_system.py index 279da8d6f0..2f14a1f4b9 100644 --- a/src/python/espressomd/cell_system.py +++ b/src/python/espressomd/cell_system.py @@ -107,6 +107,12 @@ def set_regular_decomposition(self, **kwargs): use_verlet_lists : :obj:`bool`, optional Activates or deactivates the usage of Verlet lists. Defaults to ``True``. + fully_connected_boundary : :obj:`dict`, optional + If set, connects all cells on a given boundary along the given direction. + Example: ``{"boundary": "z", "direction": "x"}`` connects all + cells on the boundary normal to the z-direction along the x-axis. + This corresponds to z-axis as shear plane normal and x-axis as + shear direction in Lees-Edwards boundary conditions. """ self.call_method("initialize", name="regular_decomposition", **kwargs) diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp index 7f3cb0082a..1efd2ed3b4 100644 --- a/src/script_interface/cell_system/CellSystem.cpp +++ b/src/script_interface/cell_system/CellSystem.cpp @@ -42,7 +42,9 @@ #include #include +#include #include +#include #include #include #include @@ -51,6 +53,25 @@ #include #include +static int coord(std::string const &s) { + if (s == "x") + return 0; + if (s == "y") + return 1; + if (s == "z") + return 2; + throw std::invalid_argument("Invalid Cartesian coordinate: '" + s + "'"); +} + +static std::string coord_letter(int c) { + if (c == 0) + return "x"; + if (c == 1) + return "y"; + assert(c == 2); + return "z"; +} + namespace ScriptInterface { namespace CellSystem { @@ -112,6 +133,20 @@ CellSystem::CellSystem() { auto const ns_types = hd.get_n_square_types(); return Variant{std::vector(ns_types.begin(), ns_types.end())}; }}, + {"fully_connected_boundary", AutoParameter::read_only, + [this]() { + if (get_cell_structure().decomposition_type() != + CellStructureType::REGULAR) { + return Variant{none}; + } + auto const rd = get_regular_decomposition(); + auto const fcb = rd.fully_connected_boundary(); + if (not fcb) + return Variant{none}; + return Variant{std::unordered_map{ + {{"boundary", Variant{coord_letter((*fcb).first)}}, + {"direction", Variant{coord_letter((*fcb).second)}}}}}; + }}, {"cutoff_regular", AutoParameter::read_only, [this]() { if (get_cell_structure().decomposition_type() != @@ -264,6 +299,21 @@ void CellSystem::initialize(CellStructureType const &cs_type, get_value_or>(params, "n_square_types", {}); auto n_square_types = std::set{ns_types.begin(), ns_types.end()}; m_cell_structure->set_hybrid_decomposition(cutoff_regular, n_square_types); + } else if (cs_type == CellStructureType::REGULAR) { + std::optional> fcb_pair = std::nullopt; + if (params.contains("fully_connected_boundary") and + not is_none(params.at("fully_connected_boundary"))) { + auto const variant = + get_value(params, "fully_connected_boundary"); + context()->parallel_try_catch([&fcb_pair, &variant]() { + fcb_pair = {{coord(boost::get(variant.at("boundary"))), + coord(boost::get(variant.at("direction")))}}; + }); + } + context()->parallel_try_catch([this, &fcb_pair]() { + m_cell_structure->set_regular_decomposition( + get_system().get_interaction_range(), fcb_pair); + }); } else { system.set_cell_structure_topology(cs_type); } diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index ee90f51439..4ebdb2dafb 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -42,6 +42,8 @@ def get_lin_pos_offset(time, initial_pos_offset=None, osc_protocol = espressomd.lees_edwards.OscillatoryShear(**params_osc) off_protocol = espressomd.lees_edwards.Off() +const_offset_protocol = espressomd.lees_edwards.LinearShear( + initial_pos_offset=2.2, shear_velocity=0) def axis(coord): @@ -53,21 +55,27 @@ def axis(coord): class LeesEdwards(ut.TestCase): + box_l = [5, 5, 5] + system = espressomd.System(box_l=box_l) + node_grid = np.copy(system.cell_system.node_grid) + n_nodes = np.prod(node_grid) - system = espressomd.System(box_l=[5.0, 5.0, 5.0]) - system.cell_system.skin = 0.0 - system.cell_system.set_n_square(use_verlet_lists=True) - - time_step = 0.5 - system.time_step = time_step direction_permutations = list(itertools.permutations(["x", "y", "z"], 2)) def setUp(self): - self.system.time = 0.0 + system = self.system + system.box_l = self.box_l + system.cell_system.skin = 0. + system.cell_system.set_n_square(use_verlet_lists=True) + system.time = 0.0 + system.time_step = 0.5 + system.min_global_cut = 0. + system.cell_system.node_grid = self.node_grid def tearDown(self): system = self.system system.part.clear() + system.bonded_inter.clear() system.lees_edwards.protocol = None if espressomd.has_features("COLLISION_DETECTION"): system.collision_detection.set_params(mode="off") @@ -191,6 +199,21 @@ def test_protocols(self): shear_direction=valid, shear_plane_normal=valid, protocol=lin_protocol) + with self.assertRaisesRegex(ValueError, "fully_connected_boundary normal and connection coordinates need to differ"): + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "z"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + with self.assertRaisesRegex(ValueError, "Invalid Cartesian coordinate: 't'"): + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "t"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + if self.n_nodes > 1: + with self.assertRaisesRegex(RuntimeError, "The MPI nodegrid must be 1 in the fully connected direction"): + system.cell_system.node_grid = [1, self.n_nodes, 1] + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "y"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + def test_boundary_crossing_lin(self): """ A particle crosses the upper and lower boundary with linear shear. @@ -300,13 +323,13 @@ def test_trajectory_reconstruction(self): crossing_time = system.time system.integrator.run(1) np.testing.assert_almost_equal( - p.lees_edwards_offset, + p.lees_edwards_offset, get_lin_pos_offset(crossing_time, **params_lin)) np.testing.assert_almost_equal(p.lees_edwards_flag, -1) system.integrator.run(1) # no boundary crossing np.testing.assert_almost_equal( - p.lees_edwards_offset, + p.lees_edwards_offset, get_lin_pos_offset(crossing_time, **params_lin)) np.testing.assert_almost_equal(p.lees_edwards_flag, 0) @@ -485,6 +508,7 @@ def test_virt_sites_interaction(self): shear_velocity=2.0, initial_pos_offset=0.0) system.lees_edwards.set_boundary_conditions( shear_direction="x", shear_plane_normal="y", protocol=protocol) + system.min_global_cut = 2.5 p1 = system.part.add(pos=[2.5, 2.5, 2.5], type=10, rotation=3 * (True,), v=(0.0, -0.1, -0.25)) p2 = system.part.add(pos=(2.5, 3.5, 2.5), type=11) @@ -535,6 +559,42 @@ def test_virt_sites_interaction(self): weight_function=0, gamma=0, r_cut=0, trans_weight_function=0, trans_gamma=0, trans_r_cut=0) + @utx.skipIfMissingFeatures( + ["EXTERNAL_FORCES", "VIRTUAL_SITES_RELATIVE"]) + def test__virt_sites_rotation(self): + """ + A particle with virtual sites is placed on the boundary. We check if + the forces yield the correct torque and if a rotation frequency is + transmitted back to the virtual sites. + """ + + system = self.system + system.part.clear() + system.min_global_cut = 2.5 + + system.lees_edwards.set_boundary_conditions( + shear_direction="x", shear_plane_normal="y", protocol=lin_protocol) + + p1 = system.part.add( + id=0, pos=[2.5, 5.0, 2.5], rotation=[True] * 3) + + p2 = system.part.add(pos=(2.5, 6.0, 2.5), ext_force=(1.0, 0., 0.)) + p2.vs_auto_relate_to(0) + p3 = system.part.add(pos=(2.5, 4.0, 2.5), ext_force=(-1.0, 0., 0.)) + p3.vs_auto_relate_to(0) + + system.integrator.run(0, recalc_forces=True) + + np.testing.assert_array_almost_equal( + np.copy(p1.torque_lab), [0.0, 0.0, -2.0]) + + p1.omega_lab = (0., 0., 2.5) + system.integrator.run(0, recalc_forces=True) + for vs in p2, p3: + np.testing.assert_array_almost_equal( + system.velocity_difference(p1, vs), + np.cross(p1.omega_lab, system.distance_vec(p1, vs))) + @utx.skipIfMissingFeatures( ["EXTERNAL_FORCES", "VIRTUAL_SITES_RELATIVE", "COLLISION_DETECTION"]) def test_le_colldet(self): @@ -685,64 +745,93 @@ def test_le_breaking_bonds(self): bond_list += p.bonds np.testing.assert_array_equal(len(bond_list), 0) - def setup_lj_liquid(self): - system = self.system - system.cell_system.set_n_square(use_verlet_lists=False) - # Parameters - n = 100 - phi = 0.4 - sigma = 1. - eps = 1 - cut = sigma * 2**(1 / 6) - - # box - l = (n / 6. * np.pi * sigma**3 / phi)**(1. / 3.) - - # Setup - system.box_l = [l, l, l] - system.lees_edwards.protocol = None - - system.time_step = 0.01 - system.thermostat.turn_off() - - np.random.seed(42) - system.part.add(pos=np.random.random((n, 3)) * l) + const_offset_params = { + 'shear_velocity': 0.0, + 'shear_direction': 0, + 'shear_plane_normal': 1, + 'initial_pos_offset': 17.2} - # interactions - system.non_bonded_inter[0, 0].lennard_jones.set_params( - epsilon=eps, sigma=sigma, cutoff=cut, shift="auto") - # Remove overlap - system.integrator.set_steepest_descent( - f_max=0, gamma=0.05, max_displacement=0.05) - while system.analysis.energy()["total"] > 0.5 * n: - system.integrator.run(5) - - system.integrator.set_vv() + shear_params = { + 'shear_velocity': 0.1, + 'shear_direction': 0, + 'shear_plane_normal': 2, + 'initial_pos_offset': -np.sqrt(0.1)} @utx.skipIfMissingFeatures("LENNARD_JONES") - def test_zz_lj(self): + def run_lj_pair_visibility(self, shear_direction, shear_plane_normal): """ - Simulate an LJ liquid under linear shear and verify forces. + Simulate LJ particles coming into contact under linear shear and verify forces. This is to make sure that no pairs get lost or are outdated - in the short range loop. To have deterministic forces, velocity - capping is used rather than a thermostat. + in the short range loop. """ + shear_axis, normal_axis = axis( + shear_direction), axis(shear_plane_normal) system = self.system - self.setup_lj_liquid() + system.part.clear() + system.time = 0 + system.time_step = 0.1 protocol = espressomd.lees_edwards.LinearShear( - shear_velocity=0.3, initial_pos_offset=0.01) + shear_velocity=3, initial_pos_offset=5) system.lees_edwards.set_boundary_conditions( - shear_direction="z", shear_plane_normal="x", protocol=protocol) - system.integrator.run(1, recalc_forces=True) - tests_common.check_non_bonded_loop_trace(self, system) + shear_direction=shear_direction, shear_plane_normal=shear_plane_normal, protocol=protocol) + system.cell_system.skin = 0.2 + system.non_bonded_inter[0, 0].lennard_jones.set_params( + epsilon=1E-6, sigma=1, cutoff=1.2, shift="auto") + system.part.add( + pos=(0.1 * normal_axis, -0.8 * normal_axis), + v=(1.0 * shear_axis, -0.3 * shear_axis)) + assert np.all(system.part.all().f == 0.) + tests_common.check_non_bonded_loop_trace( + self, system, cutoff=system.non_bonded_inter[0, 0].lennard_jones.get_params()["cutoff"] + system.cell_system.skin) # Rewind the clock to get back the LE offset applied during force calc system.time = system.time - system.time_step tests_common.verify_lj_forces(system, 1E-7) - - system.thermostat.set_langevin(kT=.1, gamma=5, seed=2) - system.integrator.run(50) - tests_common.check_non_bonded_loop_trace(self, system) + have_interacted = False + for _ in range(50): + system.integrator.run(3) + if np.any(np.abs(system.part.all().f) > 0): + have_interacted = True + tests_common.check_non_bonded_loop_trace( + self, system, cutoff=system.non_bonded_inter[0, 0].lennard_jones.get_params()["cutoff"] + system.cell_system.skin) + system.time = system.time - system.time_step + tests_common.verify_lj_forces(system, 1E-7) + assert have_interacted + + def test_zz_lj_pair_visibility(self): + # check that regular decomposition without fully connected doesn't + # catch the particle + system = self.system + system.box_l = [10, 10, 10] + with self.assertRaises(AssertionError): + system.cell_system.set_regular_decomposition( + fully_connected_boundary=None) + self.assertIsNone(system.cell_system.fully_connected_boundary) + system.cell_system.node_grid = [1, self.n_nodes, 1] + self.run_lj_pair_visibility("x", "y") + + for verlet in (False, True): + for shear_direction, shear_plane_normal in self.direction_permutations: + system.cell_system.set_n_square(use_verlet_lists=verlet) + self.run_lj_pair_visibility( + shear_direction, shear_plane_normal) + + for verlet in (False, True): + for shear_direction, shear_plane_normal in self.direction_permutations: + system.cell_system.set_regular_decomposition( + fully_connected_boundary=None) + normal_axis = axis(shear_plane_normal) + system.cell_system.node_grid = [ + self.n_nodes if normal_axis[i] == 1 else 1 for i in range(3)] + fully_connected_boundary = {"boundary": shear_plane_normal, + "direction": shear_direction} + system.cell_system.set_regular_decomposition( + use_verlet_lists=verlet, + fully_connected_boundary=fully_connected_boundary) + self.assertEqual(system.cell_system.fully_connected_boundary, + fully_connected_boundary) + self.run_lj_pair_visibility( + shear_direction, shear_plane_normal) if __name__ == "__main__": diff --git a/testsuite/python/regular_decomposition.py b/testsuite/python/regular_decomposition.py index b05746f542..d94720f11e 100644 --- a/testsuite/python/regular_decomposition.py +++ b/testsuite/python/regular_decomposition.py @@ -17,8 +17,10 @@ # along with this program. If not, see . # import unittest as ut +import unittest_decorators as utx import espressomd import numpy as np +import itertools np.random.seed(42) @@ -105,6 +107,109 @@ def test_position_rounding(self): self.system.part.add(pos=[25, 25, 0]) self.assertEqual(1, len(self.system.part)) + @utx.skipIfMissingFeatures("LENNARD_JONES") + def test_fully_connected_boundary(self): + system = self.system + system.part.clear() + if system.cell_system.node_grid[1] != 1: + ng = system.cell_system.node_grid + system.cell_system.node_grid = [ng[0], 1, ng[2] * ng[1]] + system.periodic = [True] * 3 + # Check that it's initially disabled + self.assertEqual(system.cell_system.get_params()[ + "fully_connected_boundary"], None) + + # check setting and getting the parameter + system.cell_system.set_regular_decomposition( + fully_connected_boundary=dict(direction="y", boundary="z")) + self.assertEqual(system.cell_system.get_params()[ + "fully_connected_boundary"], dict(direction="y", boundary="z")) + # Check that the setting survives cell system re-initialization + system.cell_system.min_global_cut = system.box_l / 4.1 + self.assertEqual(system.cell_system.get_params()[ + "fully_connected_boundary"], dict(direction="y", boundary="z")) + + # Check particle visibility. + # Place particles on a cubic lattice and use the + # non_bonded_loop_trace() to check that all pairs are seen as expected + fc_normal = np.array((0, 0, 1)) # z + fc_normal_coord = 2 # z + fc_dir = np.array((0, 1, 0)) # y + N = 10 + system.non_bonded_inter[0, 0].lennard_jones.set_params( + sigma=1, epsilon=1, cutoff=system.box_l[0] / N + 0.01, shift="auto") + indices = [np.array((i, j, k)) for i in range(N) + for j in range(N) for k in range(N)] + + def id_for_idx(idx): return ( + idx[0] % N) * N * N + (idx[1] % N) * N + idx[2] % N + + ids = [id_for_idx(idx) for idx in indices] + dx = system.box_l / N + positions = [idx * dx for idx in indices] + system.part.add(id=ids, pos=positions) + particles = {i: system.part.by_id(i) for i in ids} + + def distance(id1, id2): + return system.distance( + particles[id1], particles[id2]) + distances = {tuple(i): distance(*i) + for i in itertools.combinations(ids, 2)} + + max_range = np.amax(system.box_l) / N + two_cells = 2 * np.amax(system.cell_system.get_state()["cell_size"]) + two_cells_2d = two_cells * np.sqrt(2) + two_cells_3d = two_cells * np.sqrt(3) + assert np.all(system.box_l / 2 > two_cells) + + # next neighbors + must_find_nn = [i for i, d in distances.items() if d <= max_range] + + # Fully connected neighbors + indices_lower_boundary = [ + idx for idx in indices if idx[fc_normal_coord] == 0] + must_find_fc = [tuple(sorted((id_for_idx(idx), id_for_idx(idx + i * fc_dir - fc_normal)))) + for idx in indices_lower_boundary for i in range(-N + 1, N)] + + # all neighbors that must be found + must_find = set(must_find_nn + must_find_fc) + + def assert_can_find(pair): + # are the particles within a range that MAY be found by the + # pair loop + p1 = particles[pair[0]] + p2 = particles[pair[1]] + d = system.distance_vec(p1, p2) + # if not accross periodic boundary: particles must be in cells + # sharing at least one corner + if np.abs( + p1.pos - p2.pos)[fc_normal_coord] < system.box_l[fc_normal_coord] / 2: + self.assertLess(np.linalg.norm(d), two_cells_3d) + # If across a the fully connected boundary + # substract the distance in the fully connected direciont (all are + # valid + d_trans = d - d * fc_dir + # in the other TWO directions, cells have to share a corner + self.assertLess(np.linalg.norm(d_trans), two_cells_2d) + + # Use the cell system trace to get all pairs + # as opposed to get_pairs() this does not have a distance check + cs_pairs = system.cell_system.non_bonded_loop_trace() + found = [] + for id1, id2, _rest1, _rest2, _rest3, _rest4 in cs_pairs: + p = tuple(sorted((id1, id2))) # Make the pair unique + found.append(p) # to check for double counting + if p in must_find: + must_find.remove(p) + else: + assert_can_find(p) # close enough so that cells share a corner + + # Check for double counting of pairs + self.assertEqual(len(found), len(set(found))) + + # check that all required pairs have been seen + self.assertEqual(must_find, set([])) + if __name__ == "__main__": ut.main()