29#include "ParticleList.hpp"
37#include "system/Leaf.hpp"
41#include <boost/container/static_vector.hpp>
42#include <boost/iterator/indirect_iterator.hpp>
43#include <boost/range/algorithm/transform.hpp>
89 auto first_non_empty =
90 std::find_if(cells.
begin(), cells.
end(),
91 [](
const Cell *c) { return not c->particles().empty(); });
111struct MinimalImageDistance {
119struct EuclidianDistance {
136 std::vector<Particle *> m_particle_index;
138 std::unique_ptr<ParticleDecomposition> m_decomposition;
144 bool m_rebuild_verlet_list =
true;
145 std::vector<std::pair<Particle *, Particle *>> m_verlet_list;
146 double m_le_pos_offset_at_last_resort = 0.;
148 double m_verlet_skin = 0.;
149 bool m_verlet_skin_set =
false;
150 double m_verlet_reuse = 0.;
169 assert(not p or p->
id() ==
id);
171 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
172 m_particle_index.resize(
static_cast<unsigned int>(
id + 1));
174 m_particle_index[
static_cast<unsigned int>(id)] = p;
217 auto &new_part = pl.
insert(std::move(p));
239 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
242 return m_particle_index[
static_cast<unsigned int>(id)];
249 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
252 return m_particle_index[
static_cast<unsigned int>(id)];
255 template <
class InputRange,
class OutputIterator>
257 boost::transform(ids, out,
357 return assert(m_decomposition), *m_decomposition;
362 return assert(m_decomposition), *m_decomposition;
370 m_resort_particles |= level;
371 assert(m_resort_particles >= level);
395 auto const lim =
Utils::sqr(m_verlet_skin / 2.) - additional_offset.norm2();
397 particles.begin(), particles.end(), [lim](
const auto &p) {
398 return ((p.pos() - p.pos_at_last_verlet_update()).norm2() > lim);
403 return m_le_pos_offset_at_last_resort;
437#ifdef BOND_CONSTRAINT
462 if (n_verlet_updates > 0) {
463 m_verlet_reuse = n_steps /
static_cast<double>(n_verlet_updates);
483 boost::container::static_vector<Particle *, 4> partners;
487 if (std::any_of(partners.begin(), partners.end(),
488 [](
Particle const *
const p) { return p == nullptr; })) {
506 template <
class Handler>
507 void execute_bond_handler(
Particle &p, Handler handler) {
508 for (
const BondView bond : p.bonds()) {
509 auto const partner_ids = bond.partner_ids();
512 auto partners = resolve_bond_partners(partner_ids);
514 auto const bond_broken =
530 void invalidate_ghosts() {
539 void set_particle_decomposition(
572 std::set<int> n_square_types);
581 template <
class Kernel>
void link_cell(Kernel kernel) {
584 auto const first = boost::make_indirect_iterator(local_cells_span.begin());
585 auto const last = boost::make_indirect_iterator(local_cells_span.end());
594 throw std::runtime_error(
"Non-cuboid box type is not compatible with a "
595 "particle decomposition that relies on "
596 "EuclideanDistance for distance calculation.");
600 [&kernel, df = detail::EuclidianDistance{}](
610 template <
class PairKernel,
class VerletCriterion>
611 void verlet_list_loop(PairKernel pair_kernel,
616 if (m_rebuild_verlet_list) {
617 m_verlet_list.clear();
620 if (verlet_criterion(p1, p2, d)) {
621 m_verlet_list.emplace_back(&p1, &p2);
622 pair_kernel(p1, p2, d);
626 m_rebuild_verlet_list =
false;
631 auto const distance_function =
633 for (
auto &pair : m_verlet_list) {
634 pair_kernel(*pair.first, *pair.second,
635 distance_function(*pair.first, *pair.second));
638 auto const distance_function = detail::EuclidianDistance{};
639 for (
auto &pair : m_verlet_list) {
640 pair_kernel(*pair.first, *pair.second,
641 distance_function(*pair.first, *pair.second));
651 template <
class BondKernel>
void bond_loop(BondKernel
const &bond_kernel) {
653 execute_bond_handler(p, bond_kernel);
661 link_cell(pair_kernel);
669 template <
class PairKernel,
class VerletCriterion>
673 verlet_list_loop(pair_kernel, verlet_criterion);
676 link_cell(pair_kernel);
716 return particle_to_cell(p);
727 template <
class Kernel>
732 if (cell ==
nullptr) {
739 auto const distance_function =
741 short_range_neighbor_loop(p, cell, kernel, distance_function);
743 auto const distance_function = detail::EuclidianDistance{};
744 short_range_neighbor_loop(p, cell, kernel, distance_function);
750 template <
class Kernel,
class DistanceFunc>
751 void short_range_neighbor_loop(
Particle const &p1,
Cell *
const cell,
752 Kernel &kernel, DistanceFunc
const &df) {
754 for (
auto const &p2 : cell->particles()) {
755 if (p1.
id() != p2.
id()) {
756 auto const vec = df(p1, p2).vec21;
761 for (
auto const neighbor : cell->neighbors().all()) {
763 if (neighbor != cell) {
764 for (
auto const &p2 : neighbor->
particles()) {
765 auto const vec = df(p1, p2).vec21;
ParticleIterator< Cell *const * > CellParticleIterator
CellStructureType
Cell structure topology.
@ NSQUARE
Atom decomposition (N-square).
unsigned map_data_parts(unsigned data_parts)
Map the data parts flags from cells to those used internally by the ghost communication.
void bond_broken_error(int id, Utils::Span< const int > partner_ids)
Immutable view on a bond.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
A distributed particle decomposition.
virtual Utils::Span< Cell *const > local_cells() const =0
Get pointer to local cells.
virtual Utils::Vector3d max_cutoff() const =0
Maximum supported cutoff.
virtual Cell * particle_to_cell(Particle const &p)=0
Determine which cell a particle id belongs to.
virtual Utils::Vector3d max_range() const =0
Range in which calculations are performed.
virtual BoxGeometry const & box() const =0
virtual boost::optional< BoxGeometry > minimum_image_distance() const =0
Return the box geometry needed for distance calculation if minimum image convention should be used ne...
Abstract class that represents a component of the system.
std::size_t capacity() const
Capacity of the container.
T & insert(T const &v)
Insert an element into the container.
std::size_t size() const
Number of elements in the container.
A stripped-down version of std::span from C++17.
DEVICE_QUALIFIER constexpr iterator end() const
DEVICE_QUALIFIER constexpr iterator begin() const
Returns true if the particles are to be considered for short range interactions.
This file contains the defaults for ESPResSo.
Ghost particles and particle exchange.
void link_cell(CellIterator first, CellIterator last, PairKernel &&pair_kernel)
Iterates over all particles in the cell range, and over all pairs within the cells and with their nei...
DataPart
Flags to select particle parts for communication.
@ DATA_PART_MOMENTUM
Particle::m.
@ DATA_PART_FORCE
Particle::f.
@ DATA_PART_PROPERTIES
Particle::p.
@ DATA_PART_BONDS
Particle::bonds.
@ DATA_PART_RATTLE
Particle::rattle.
@ DATA_PART_POSITION
Particle::r.
ParticleRange particles(Utils::Span< Cell *const > cells)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
DEVICE_QUALIFIER constexpr Span< T > make_span(T *p, std::size_t N)
Exception indicating that a particle id could not be resolved.
Describes a cell structure / cell system.
ParticleRange ghost_particles() const
Particle * get_local_particle(int id)
Get a local particle by id.
void update_particle_index(ParticleList &pl)
Update local particle index.
void check_particle_sorting() const
Check that particles are in the correct cell.
void set_regular_decomposition(double range)
Set the particle decomposition to RegularDecomposition.
void clear_resort_particles()
Set the resort level to sorted.
Cell * find_current_cell(const Particle &p)
Find cell a particle is stored in.
auto is_verlet_skin_set() const
Whether the Verlet skin is set.
ParticleDecomposition const & decomposition() const
Get the underlying particle decomposition.
void clear_particle_index()
Clear the particles index.
void update_ghosts_and_resort_particle(unsigned data_parts)
Update ghost particles, with particle resort if needed.
Particle * add_local_particle(Particle &&p)
Add a particle.
void set_verlet_skin_heuristic()
Set the Verlet skin using a heuristic.
void set_verlet_skin(double value)
Set the Verlet skin.
void ghosts_update(unsigned data_parts)
Update ghost particles.
auto get_le_pos_offset_at_last_resort() const
void get_local_particles(InputRange ids, OutputIterator out)
void update_verlet_stats(int n_steps, int n_verlet_updates)
void update_particle_index(int id, Particle *p)
Update local particle index.
void ghosts_reduce_forces()
Add forces from ghost particles to real particles.
unsigned get_resort_particles() const
Get the currently scheduled resort level.
auto get_verlet_reuse() const
Average number of integration steps the Verlet list was re-used.
void non_bonded_loop(PairKernel pair_kernel)
Non-bonded pair loop.
Utils::Vector3d max_range() const
Maximal pair range supported by current cell system.
bool check_resort_required(Utils::Vector3d const &additional_offset={}) const
Check whether a particle has moved further than half the skin since the last Verlet list update,...
const Particle * get_local_particle(int id) const
This is an overloaded member function, provided for convenience. It differs from the above function o...
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
void ghosts_count()
Synchronize number of ghosts.
void set_resort_particles(Cells::Resort level)
Increase the local resort level at least to level.
void remove_particle(int id)
Remove a particle.
Particle * add_particle(Particle &&p)
Add a particle.
void resort_particles(bool global_flag)
Resort particles.
void check_particle_index() const
Check that particle index is commensurate with particles.
auto get_verlet_skin() const
Get the Verlet skin.
void set_atom_decomposition()
Set the particle decomposition to AtomDecomposition.
bool run_on_particle_short_range_neighbors(Particle const &p, Kernel &kernel)
Run kernel on all particles inside local cell and its neighbors.
void remove_all_particles()
Remove all particles from the cell system.
ParticleRange local_particles() const
void update_particle_index(Particle &p)
Update local particle index.
void ghosts_reduce_rattle_correction()
Add rattle corrections from ghost particles to real particles.
CellStructureType decomposition_type() const
void set_hybrid_decomposition(double cutoff_regular, std::set< int > n_square_types)
Set the particle decomposition to HybridDecomposition.
int get_max_local_particle_id() const
Get the maximal particle id on this node.
Utils::Vector3d max_cutoff() const
Maximal cutoff supported by current cell system.
void non_bonded_loop(PairKernel pair_kernel, const VerletCriterion &verlet_criterion)
Non-bonded pair loop with potential use of verlet lists.
Distance vector and length handed to pair kernels.
Distance(Utils::Vector3d const &vec21)
Struct holding all information for one particle.