ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
RDF.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19#include "RDF.hpp"
20
21#include "BoxGeometry.hpp"
22#include "fetch_particles.hpp"
23#include "system/System.hpp"
24
25#include <utils/Vector.hpp>
28
29#include <boost/mpi/communicator.hpp>
30#include <boost/range/combine.hpp>
31
32#include <cmath>
33#include <cstddef>
34#include <numbers>
35#include <vector>
36
37namespace Observables {
38std::vector<double>
39RDF::operator()(boost::mpi::communicator const &comm) const {
40 auto const &local_particles_1 = fetch_particles(ids1());
41
42 if (ids2().empty()) {
43 return this->evaluate(comm, local_particles_1, {}, {});
44 }
45
46 auto const &local_particles_2 = fetch_particles(ids2());
47
48 return this->evaluate(comm, local_particles_1, local_particles_2, {});
49}
50
51std::vector<double>
52RDF::evaluate(boost::mpi::communicator const &comm,
53 ParticleReferenceRange const &local_particles_1,
54 ParticleReferenceRange const &local_particles_2,
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);
60
61 if (comm.rank() != 0) {
62 return {};
63 }
64
65 auto const &box_geo = *System::get_system().box_geo;
66 auto const bin_width = (max_r - min_r) / static_cast<double>(n_r_bins);
67 auto const inv_bin_width = 1.0 / bin_width;
68 std::vector<double> res(n_r_bins, 0.0);
69 long int cnt = 0;
70 auto op = [this, inv_bin_width, &cnt, &res, &box_geo](auto const &pos1,
71 auto const &pos2) {
72 auto const dist = box_geo.get_mi_vector(pos1, pos2).norm();
73 if (dist > min_r && dist < max_r) {
74 auto const ind =
75 static_cast<int>(std::floor((dist - min_r) * inv_bin_width));
76 res[ind]++;
77 }
78 cnt++;
79 };
80
81 if (local_particles_2.empty()) {
82 Utils::for_each_pair(positions_1, op);
83 } else {
84 auto const combine_1 = boost::combine(ids1(), positions_1);
85 auto const combine_2 = boost::combine(ids2(), positions_2);
86
87 auto op2 = [&op](auto const &it1, auto const &it2) {
88 auto const &[id1, pos1] = it1;
89 auto const &[id2, pos2] = it2;
90
91 op(pos1, pos2);
92 };
93
94 auto cmp = [](auto const &it1, auto const &it2) {
95 auto const &[id1, pos1] = it1;
96 auto const &[id2, pos2] = it2;
97
98 return id1 != id2;
99 };
100 Utils::for_each_cartesian_pair_if(combine_1, combine_2, op2, cmp);
101 }
102 if (cnt == 0)
103 return res;
104 // normalization
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));
113 }
114
115 return res;
116}
117} // namespace Observables
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
Definition RDF.cpp:39
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
System & get_system()
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...