Skip to content

Commit

Permalink
implemented py interface, fixed issues with c implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
stekajack committed Dec 14, 2021
1 parent 222a7c3 commit 5b696f9
Show file tree
Hide file tree
Showing 6 changed files with 91 additions and 64 deletions.
25 changes: 0 additions & 25 deletions src/core/CellStructure.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -134,31 +134,6 @@ Particle *CellStructure::add_particle(Particle &&p) {
append_indexed_particle(cell->particles(), std::move(p)));
}

template<class Kernel> 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; });
Expand Down
28 changes: 25 additions & 3 deletions src/core/CellStructure.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<class Kernel>
bool run_on_particle_short_range_neighbors(Particle const &p, Kernel kernel);

template <class Kernel>
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
61 changes: 40 additions & 21 deletions src/core/cells.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -85,29 +85,32 @@ std::vector<std::pair<int, int>> get_pairs_filtered(double const distance,
return ret;
}

std::vector<int> get_short_range_neighbors(const Particle& p, const double& distance) {
std::vector<int> 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 <class Filter>
std::vector<int> get_short_range_neighbors(const Particle &p,
const double &distance,
Filter filter) {
std::vector<int> 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 {
Expand Down Expand Up @@ -137,6 +140,22 @@ std::vector<PairInfo> non_bonded_loop_trace() {
return ret;
}

boost::optional<std::vector<int>>
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<int> 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; });
Expand Down
27 changes: 13 additions & 14 deletions src/core/cells.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,14 @@ void cells_update_ghosts(unsigned data_parts);
*/
std::vector<std::pair<int, int>> 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<int> 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
*
Expand All @@ -110,26 +118,17 @@ void check_resort_particles();
*/
std::vector<int> 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<int> 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.
Expand Down
7 changes: 6 additions & 1 deletion src/python/espressomd/cellsystem.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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()
7 changes: 7 additions & 0 deletions src/python/espressomd/cellsystem.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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):
"""
Expand Down

3 comments on commit 5b696f9

@stekajack
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

related to issue espressomd#4399
continuation from WIP espressomd#4401
Followed description provided in the issue description for get_short_range_neighbors
Complies fine, passes check_unit_tests.
Missing:

  • Write unit_tests/short_range_neighbors_test.cpp to test both CellStructure::run_on_particle_short_range_neighbors and get_short_range_neighbors
  • didnt finish particle_short_range_energy_contribution integration
  • test particle_short_range_energy_contribution within python interface
    Hope this helps

@jhossbach
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be better to also test the get_short_range_neighbors function within the python interface, otherwise we will have to construct a system with boxgeometry within the c++ interface.

@jhossbach
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just realized that particle_to_cell as defined in DomainDecomposition also requires particle positions... Writing unit tests might be harder than I thought. There are currently no unit tests for CellStructures to copy from. I will consult with the Espresso core devs on the next move

Please sign in to comment.