diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index 6a17cad57d8..a1a5921ea39 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -134,31 +134,6 @@ Particle *CellStructure::add_particle(Particle &&p) { append_indexed_particle(cell->particles(), std::move(p))); } -template bool CellStructure::run_on_particle_short_range_neighbors(Particle const &p, Kernel kernel) { - auto const cell = particle_to_cell(p); - if (cell == nullptr) { - return false; - } - // Iterate over particles inside cell - for (auto const &part : cell->particles()) { - if (part == p) { - continue; - } - kernel(part); - } - // Iterate over all neighbors - for (auto &neighbor : cell->neighbors().all()) { - // Iterate over particles in neighbors - for (auto const &part : neighbor->particles()) { - // TODO: WIP, remove later - assert(part != p); - - kernel(part); - } - } - return true; -} - int CellStructure::get_max_local_particle_id() const { auto it = std::find_if(m_particle_index.rbegin(), m_particle_index.rend(), [](const Particle *p) { return p != nullptr; }); diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index e1ddf259ccf..85671c729b0 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -642,8 +642,30 @@ struct CellStructure { * @param kernel Function to apply to all particles inside cell and neighbors * @return false if cell is not found via particle_to_cell, otherwise true */ - template - bool run_on_particle_short_range_neighbors(Particle const &p, Kernel kernel); - + template + bool run_on_particle_short_range_neighbors(const Particle &p, Kernel kernel) { + auto const cell = particle_to_cell(p); + if (cell == nullptr) { + return false; + } + // Iterate over particles inside cell + for (auto const &part : cell->particles()) { + if (part == p) { + continue; + } + kernel(part); + } + // Iterate over all neighbors + for (auto &neighbor : cell->neighbors().all()) { + // Iterate over particles in neighbors + for (auto const &part : neighbor->particles()) { + // TODO: WIP, remove later + assert(part != p); + + kernel(part); + } + } + return true; + } }; #endif // ESPRESSO_CELLSTRUCTURE_HPP diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 8c612ce47b8..ceacb64f381 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -28,12 +28,12 @@ #include "Particle.hpp" #include "communication.hpp" +#include "energy_inline.hpp" #include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" #include "particle_data.hpp" -#include "energy_inline.hpp" #include "DomainDecomposition.hpp" #include "ParticleDecomposition.hpp" @@ -85,29 +85,32 @@ std::vector> get_pairs_filtered(double const distance, return ret; } -std::vector get_short_range_neighbors(const Particle& p, const double& distance) { - std::vector ret; - auto const cutoff2 = distance * distance; - auto kernel = [&ret, &p, &cutoff2, df = detail::MinimalImageDistance{}](const Particle &p1) { - if (df(p1, p).dist2 < cutoff2) { - ret.emplace_back(p1.p.identity); - }; +template +std::vector get_short_range_neighbors(const Particle &p, + const double &distance, + Filter filter) { + std::vector ret; + auto const cutoff2 = distance * distance; + auto kernel = [&ret, &p, &cutoff2, &filter](const Particle &p1) { + if (filter(p1, p).dist2 < cutoff2) { + ret.emplace_back(p1.p.identity); }; - cell_structure.run_on_particle_short_range_neighbors(p, kernel); - return ret; + }; + cell_structure.run_on_particle_short_range_neighbors(p, kernel); + return ret; } -double particle_short_range_energy_contribution(const Particle& p) { - double ret = 0.0; - auto kernel = [&ret, &p, df = detail::MinimalImageDistance{}](const Particle &p1) { - const Utils::Vector3d vec = df(p,p1).vec21; - const auto dist2 = df(p,p1).dist2; - const auto dist = sqrt(dist2); - // Add energy for current particle pair to result - ret += compute_non_bonded_pair_energy(p, p1, vec, dist, dist2); - - }; - return ret; +double particle_short_range_energy_contribution(const Particle &p) { + double ret = 0.0; + auto kernel = [&ret, &p, + df = detail::MinimalImageDistance{}](const Particle &p1) { + const Utils::Vector3d vec = df(p, p1).vec21; + const auto dist2 = df(p, p1).dist2; + const auto dist = sqrt(dist2); + // Add energy for current particle pair to result + ret += compute_non_bonded_pair_energy(p, p1, vec, dist, dist2); + }; + return ret; } namespace boost { @@ -137,6 +140,22 @@ std::vector non_bonded_loop_trace() { return ret; } +boost::optional> +mpi_get_short_range_neighbors_local(const Particle &p, const double distance) { + auto neighbors_list = + get_short_range_neighbors(p, distance, detail::MinimalImageDistance{}); + Utils::Mpi::gather_buffer(neighbors_list, comm_cart); + return neighbors_list; +} + +REGISTER_CALLBACK_ONE_RANK(mpi_get_short_range_neighbors_local) + +std::vector mpi_get_short_range_neighbors(const Particle &p, + const double distance) { + return mpi_call(::Communication::Result::one_rank, + mpi_get_short_range_neighbors_local, p, distance); +} + static auto mpi_get_pairs_local(double const distance) { auto pairs = get_pairs_filtered(distance, [](Particle const &) { return true; }); diff --git a/src/core/cells.hpp b/src/core/cells.hpp index 81e575e6039..d08e1e4f38e 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -88,6 +88,14 @@ void cells_update_ghosts(unsigned data_parts); */ std::vector> mpi_get_pairs(double distance); +/** + * @brief Get ids of clos-erange neighbours of a particle P, that is, particles + * within the cell of particle P and its neighbouring cells than @p distance + * + */ +std::vector mpi_get_short_range_neighbors(const Particle &p, + const double distance); + /** * @brief Get pairs closer than @p distance if both their types are in @p types * @@ -110,26 +118,17 @@ void check_resort_particles(); */ std::vector mpi_resort_particles(int global_flag); -/** - * @brief Get short range neighbors of particle p - * - * This function iterates over all particles in the cell of p and its neighboring cells and adds particles close than the specified distance to the return list. - * - * @param p Particle to find short range neighbors for - * @param distance Short range cutoff distance - * @return Vector containing the particle ids of the neighbors - */ -std::vector get_short_range_neighbors(const Particle& p, const double& distance); - /** * @brief Compute short range energy of particle p * - * Iterates through particles inside cell and neighboring cells and computes energy contribution for particle p + * Iterates through particles inside cell and neighboring cells and computes + * energy contribution for particle p * * @param p Particle to calculate energy for - * @return Double containing short-range Coulomb and non-bonded energy of the particle + * @return Double containing short-range Coulomb and non-bonded energy of the + * particle */ -double particle_short_range_energy_contribution(const Particle& p); +double particle_short_range_energy_contribution(const Particle &p); /** * @brief Find the cell in which a particle is stored. diff --git a/src/python/espressomd/cellsystem.pxd b/src/python/espressomd/cellsystem.pxd index ecf56f006c8..93dee961058 100644 --- a/src/python/espressomd/cellsystem.pxd +++ b/src/python/espressomd/cellsystem.pxd @@ -35,6 +35,10 @@ cdef extern from "cells.hpp": cdef extern from "communication.hpp": int n_nodes +cdef extern from "particle_data.hpp": + ctypedef struct particle "Particle" + const particle & get_particle_data(int pid) except + + cdef extern from "cells.hpp": int CELL_STRUCTURE_DOMDEC int CELL_STRUCTURE_NSQUARE @@ -48,6 +52,7 @@ cdef extern from "cells.hpp": const DomainDecomposition * get_domain_decomposition() vector[pair[int, int]] mpi_get_pairs(double distance) except + + vector[int] mpi_get_short_range_neighbors(const particle &p, const double distance) except + vector[pair[int, int]] mpi_get_pairs_of_types(double distance, vector[int] types) except + vector[PairInfo] mpi_non_bonded_loop_trace() vector[int] mpi_resort_particles(int global_flag) @@ -74,4 +79,4 @@ cdef extern from "bonded_interactions/bonded_interaction_data.hpp": double maximal_cutoff_bonded() cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp": - double maximal_cutoff_nonbonded() + double maximal_cutoff_nonbonded() \ No newline at end of file diff --git a/src/python/espressomd/cellsystem.pyx b/src/python/espressomd/cellsystem.pyx index 9291fb0acfd..90c18aae460 100644 --- a/src/python/espressomd/cellsystem.pyx +++ b/src/python/espressomd/cellsystem.pyx @@ -28,6 +28,7 @@ from .cellsystem cimport mpi_set_skin, skin from .utils import handle_errors from .utils cimport Vector3i from .utils cimport check_type_or_throw_except, make_array_locked +from .particle_data cimport ParticleHandle cdef class CellSystem: def set_domain_decomposition(self, use_verlet_lists=True): @@ -128,6 +129,12 @@ cdef class CellSystem: pairs = self._get_pairs_of_types(distance, types) handle_errors("") return pairs + + def get_short_range_neighbors(self, ParticleHandle particle, distance): + cdef vector[int] list_of_neighbours + pass_part = get_particle_data(particle.id) + list_of_neighbours = mpi_get_short_range_neighbors(pass_part, distance) + return list_of_neighbours def _get_pairs_of_types(self, distance, types): """