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>
91 auto first_non_empty = std::ranges::find_if(
92 cells, [](
const Cell *c) {
return not c->
particles().empty(); });
112struct MinimalImageDistance {
120struct EuclidianDistance {
137 std::vector<Particle *> m_particle_index;
139 std::unique_ptr<ParticleDecomposition> m_decomposition;
145 bool m_rebuild_verlet_list =
true;
146 std::vector<std::pair<Particle *, Particle *>> m_verlet_list;
147 double m_le_pos_offset_at_last_resort = 0.;
149 double m_verlet_skin = 0.;
150 bool m_verlet_skin_set =
false;
151 double m_verlet_reuse = 0.;
170 assert(not p or p->
id() ==
id);
172 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
173 m_particle_index.resize(
static_cast<unsigned int>(
id + 1));
175 m_particle_index[
static_cast<unsigned int>(id)] = p;
218 auto &new_part = pl.
insert(std::move(p));
240 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
243 return m_particle_index[
static_cast<unsigned int>(id)];
250 if (
static_cast<unsigned int>(
id) >= m_particle_index.size())
253 return m_particle_index[
static_cast<unsigned int>(id)];
256 template <
class InputRange,
class OutputIterator>
258 std::ranges::transform(ids, out,
358 return assert(m_decomposition), *m_decomposition;
363 return assert(m_decomposition), *m_decomposition;
371 m_resort_particles |= level;
372 assert(m_resort_particles >= level);
396 auto const lim =
Utils::sqr(m_verlet_skin / 2.) - additional_offset.norm2();
398 particles.begin(), particles.end(), [lim](
const auto &p) {
399 return ((p.pos() - p.pos_at_last_verlet_update()).norm2() > lim);
404 return m_le_pos_offset_at_last_resort;
438#ifdef BOND_CONSTRAINT
463 if (n_verlet_updates > 0) {
464 m_verlet_reuse = n_steps /
static_cast<double>(n_verlet_updates);
483 auto resolve_bond_partners(std::span<const int> partner_ids) {
484 boost::container::static_vector<Particle *, 4> partners;
488 if (std::ranges::find(partners,
nullptr) != partners.end()) {
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);
513 auto const partners_span = std::span(partners.data(), partners.size());
514 auto const bond_broken = handler(p, bond.bond_id(), partners_span);
528 void invalidate_ghosts() {
537 void set_particle_decomposition(
564 std::optional<std::pair<int, int>> fully_connected_boundary);
573 std::set<int> n_square_types);
582 template <
class Kernel>
void link_cell(Kernel kernel) {
585 auto const first = boost::make_indirect_iterator(local_cells_span.begin());
586 auto const last = boost::make_indirect_iterator(local_cells_span.end());
595 throw std::runtime_error(
"Non-cuboid box type is not compatible with a "
596 "particle decomposition that relies on "
597 "EuclideanDistance for distance calculation.");
601 [&kernel, df = detail::EuclidianDistance{}](
611 template <
class PairKernel,
class VerletCriterion>
612 void verlet_list_loop(PairKernel pair_kernel,
617 if (m_rebuild_verlet_list) {
618 m_verlet_list.clear();
621 if (verlet_criterion(p1, p2, d)) {
622 m_verlet_list.emplace_back(&p1, &p2);
623 pair_kernel(p1, p2, d);
627 m_rebuild_verlet_list =
false;
632 auto const distance_function =
634 for (
auto &pair : m_verlet_list) {
635 pair_kernel(*pair.first, *pair.second,
636 distance_function(*pair.first, *pair.second));
639 auto const distance_function = detail::EuclidianDistance{};
640 for (
auto &pair : m_verlet_list) {
641 pair_kernel(*pair.first, *pair.second,
642 distance_function(*pair.first, *pair.second));
652 template <
class BondKernel>
void bond_loop(BondKernel
const &bond_kernel) {
654 execute_bond_handler(p, bond_kernel);
662 link_cell(pair_kernel);
670 template <
class PairKernel,
class VerletCriterion>
674 verlet_list_loop(pair_kernel, verlet_criterion);
677 link_cell(pair_kernel);
717 return particle_to_cell(p);
728 template <
class Kernel>
733 if (cell ==
nullptr) {
740 auto const distance_function =
742 short_range_neighbor_loop(p, cell, kernel, distance_function);
744 auto const distance_function = detail::EuclidianDistance{};
745 short_range_neighbor_loop(p, cell, kernel, distance_function);
751 template <
class Kernel,
class DistanceFunc>
752 void short_range_neighbor_loop(
Particle const &p1,
Cell *
const cell,
753 Kernel &kernel, DistanceFunc
const &df) {
755 for (
auto const &p2 : cell->particles()) {
756 if (p1.
id() != p2.
id()) {
757 auto const vec = df(p1, p2).vec21;
762 for (
auto const neighbor : cell->neighbors().all()) {
764 if (neighbor != cell) {
765 for (
auto const &p2 : neighbor->
particles()) {
766 auto const vec = df(p1, p2).vec21;
ParticleIterator< std::span< Cell *const >::iterator > 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, std::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.
auto & particles()
Particles.
A distributed particle decomposition.
virtual Utils::Vector3d max_cutoff() const =0
Maximum supported cutoff.
virtual std::span< Cell *const > local_cells() const =0
Get pointer to local cells.
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 std::optional< BoxGeometry > minimum_image_distance() const =0
Return the box geometry needed for distance calculation if minimum image convention should be used ne...
virtual BoxGeometry const & box() const =0
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.
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(std::span< Cell *const > cells)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
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 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_regular_decomposition(double range, std::optional< std::pair< int, int > > fully_connected_boundary)
Set the particle decomposition to RegularDecomposition.
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.