23#include "system/System.hpp"
29#include <boost/mpi/communicator.hpp>
30#include <boost/range/combine.hpp>
43 return this->evaluate(comm, local_particles_1, {}, {});
48 return this->evaluate(comm, local_particles_1, local_particles_2, {});
52RDF::evaluate(boost::mpi::communicator
const &comm,
56 auto const positions_1 = detail::get_all_particle_positions(
57 comm, local_particles_1,
ids1(),
traits,
false);
58 auto const positions_2 = detail::get_all_particle_positions(
59 comm, local_particles_2,
ids2(),
traits,
false);
61 if (comm.rank() != 0) {
67 auto const inv_bin_width = 1.0 / bin_width;
68 std::vector<double> res(
n_r_bins, 0.0);
70 auto op = [
this, inv_bin_width, &cnt, &res, &box_geo](
auto const &pos1,
72 auto const dist = box_geo.get_mi_vector(pos1, pos2).norm();
75 static_cast<int>(std::floor((dist -
min_r) * inv_bin_width));
81 if (local_particles_2.empty()) {
84 auto const combine_1 = boost::combine(
ids1(), positions_1);
85 auto const combine_2 = boost::combine(
ids2(), positions_2);
87 auto op2 = [&op](
auto const &it1,
auto const &it2) {
88 auto const &[id1, pos1] = it1;
89 auto const &[id2, pos2] = it2;
94 auto cmp = [](
auto const &it1,
auto const &it2) {
95 auto const &[id1, pos1] = it1;
96 auto const &[id2, pos2] = it2;
105 auto const volume = box_geo.volume();
106 for (std::size_t i = 0u; i < res.size(); ++i) {
107 auto const r_in =
static_cast<double>(i) * bin_width +
min_r;
108 auto const r_out = r_in + bin_width;
109 auto const bin_volume =
110 (4. / 3.) * std::numbers::pi *
111 (Utils::int_pow<3>(r_out) - Utils::int_pow<3>(r_in));
112 res[i] *= volume / (bin_volume *
static_cast<double>(cnt));
Vector implementation and trait types for boost qvm interoperability.
std::vector< int > & ids2()
std::vector< int > & ids1()
std::vector< double > operator()(boost::mpi::communicator const &comm) const final
std::shared_ptr< BoxGeometry > box_geo
auto fetch_particles(std::vector< int > const &ids)
Fetch a group of particles.
std::vector< std::reference_wrapper< Particle const > > ParticleReferenceRange
void for_each_pair(ForwardIterator first, ForwardIterator last, BinaryOp op)
Execute op for each pair of elements in [first, last) once.
void for_each_cartesian_pair_if(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2, ForwardIterator last2, BinaryOp op, BinaryCmp cmp)
Execute op for each pair of elements between [first1, last1) and [first2, last2) if a condition is sa...