ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
statistics.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21/** \file
22 * Statistical tools to analyze simulations.
23 *
24 * The corresponding header file is statistics.hpp.
25 */
26
28
29#include "BoxGeometry.hpp"
30#include "Particle.hpp"
32#include "communication.hpp"
33#include "errorhandling.hpp"
34#include "system/System.hpp"
35
36#include <utils/Vector.hpp>
37#include <utils/contains.hpp>
38#include <utils/math/sqr.hpp>
40
41#include <boost/mpi/collectives/broadcast.hpp>
42#include <boost/mpi/collectives/reduce.hpp>
43
44#include <cassert>
45#include <cmath>
46#include <cstdlib>
47#include <functional>
48#include <limits>
49#include <numbers>
50#include <stdexcept>
51#include <tuple>
52#include <utility>
53#include <vector>
54
56/** @brief Serialize @c std::tuple. */
57template <typename Archive, typename... T>
58void serialize(Archive &ar, std::tuple<T...> &pack, unsigned int const) {
59 std::apply([&](auto &...item) { ((ar & item), ...); }, pack);
60}
61} // namespace boost::serialization
62
63/** @brief Decay a tuple of only 1 type to that type. */
64template <typename... F> struct DecayTupleResult {
65 using type = std::tuple<std::invoke_result_t<F, Particle const &>...>;
66};
67
68template <typename F> struct DecayTupleResult<F> {
69 using type = std::invoke_result_t<F, Particle const &>;
70};
71
72/**
73 * @brief Gather particle traits to MPI rank 0.
74 * When only one trait is requested, return a vector of type T.
75 * When multiple traits are requested, return a vector of tuples.
76 */
77template <class... Trait>
78static auto gather_traits_for_types(System::System const &system,
79 std::vector<int> const &p_types,
80 Trait &&...trait) {
81 std::vector<typename DecayTupleResult<Trait...>::type> buffer{};
82
83 for (auto const &p : system.cell_structure->local_particles()) {
84 if (Utils::contains(p_types, p.type())) {
85 buffer.emplace_back(trait(p)...);
86 }
87 }
88
90 if (::comm_cart.rank() != 0) {
91 buffer.clear();
92 }
93 return buffer;
94}
95
96double mindist(System::System const &system, std::vector<int> const &set1,
97 std::vector<int> const &set2) {
98 using Utils::contains;
99
100 std::vector<int> buf_type{};
101 std::vector<Utils::Vector3d> buf_pos{};
102
103 auto const &box_geo = *system.box_geo;
104 auto const accept_all = set1.empty() or set2.empty();
105 for (auto const &p : system.cell_structure->local_particles()) {
106 if (accept_all or contains(set1, p.type()) or contains(set2, p.type())) {
107 buf_type.emplace_back(p.type());
108 buf_pos.emplace_back(box_geo.unfolded_position(p.pos(), p.image_box()));
109 }
110 }
111
114 if (::comm_cart.rank() != 0) {
115 buf_type.clear();
116 buf_pos.clear();
117 }
118
119 auto mindist_sq = std::numeric_limits<double>::infinity();
120 for (std::size_t j = 0ul; j < buf_type.size(); ++j) {
121 /* check which sets particle j belongs to (bit 0: set1, bit1: set2) */
122 auto in_set = 0;
123 if (set1.empty() || contains(set1, buf_type[j]))
124 in_set = 1;
125 if (set2.empty() || contains(set2, buf_type[j]))
126 in_set |= 2;
127 assert(in_set != 0);
128
129 for (auto i = j + 1ul; i != buf_type.size(); ++i) {
130 /* accept a pair if particle j is in set1 and particle i in set2 or vice
131 * versa. */
132 if (((in_set & 1) and (set2.empty() or contains(set2, buf_type[i]))) or
133 ((in_set & 2) and (set1.empty() or contains(set1, buf_type[i])))) {
134 mindist_sq = std::min(
135 mindist_sq, box_geo.get_mi_vector(buf_pos[j], buf_pos[i]).norm2());
136 }
137 }
138 }
139
140 return std::sqrt(mindist_sq);
141}
142
144 bool include_particles,
145 bool include_lbfluid) {
146 Utils::Vector3d momentum{};
147 if (include_particles) {
148 auto const particles = system.cell_structure->local_particles();
149 momentum =
150 std::accumulate(particles.begin(), particles.end(), Utils::Vector3d{},
151 [](Utils::Vector3d const &m, Particle const &p) {
152 return m + p.mass() * p.v();
153 });
154 }
155 if (include_lbfluid and system.lb.is_solver_set()) {
156 momentum += system.lb.get_momentum() * system.lb.get_lattice_speed();
157 }
158 return momentum;
159}
160
162 auto const &box_geo = *system.box_geo;
163 auto const &cell_structure = *system.cell_structure;
164 Utils::Vector3d local_com{};
165 double local_mass = 0.;
166
167 for (auto const &p : cell_structure.local_particles()) {
168 if ((p.type() == p_type or p_type == -1) and not p.is_virtual()) {
169 local_com += box_geo.unfolded_position(p.pos(), p.image_box()) * p.mass();
170 local_mass += p.mass();
171 }
172 }
173 Utils::Vector3d com{};
174 double mass = 1.; // placeholder value to avoid division by zero
175 boost::mpi::reduce(::comm_cart, local_com, com, std::plus<>(), 0);
176 boost::mpi::reduce(::comm_cart, local_mass, mass, std::plus<>(), 0);
177 return com / mass;
178}
179
181 auto const &box_geo = *system.box_geo;
182 auto const &cell_structure = *system.cell_structure;
183 Utils::Vector3d am{};
184
185 for (auto const &p : cell_structure.local_particles()) {
186 if ((p.type() == p_type or p_type == -1) and not p.is_virtual()) {
187 auto const pos = box_geo.unfolded_position(p.pos(), p.image_box());
188 am += p.mass() * vector_product(pos, p.v());
189 }
190 }
191 return am;
192}
193
195 std::vector<int> const &p_types) {
196 auto const &box_geo = *system.box_geo;
197 auto const trait_pos = [&box_geo](Particle const &p) {
198 return box_geo.unfolded_position(p.pos(), p.image_box());
199 };
200 auto const buf_pos = gather_traits_for_types(system, p_types, trait_pos);
201
202 Utils::Vector9d mat{};
203 if (::comm_cart.rank() == 0) {
204 auto const center =
205 std::accumulate(buf_pos.begin(), buf_pos.end(), Utils::Vector3d{}) /
206 static_cast<double>(buf_pos.size());
207 // compute covariance matrix
208 for (unsigned int i = 0u; i < 3u; ++i) {
209 for (unsigned int j = 0u; j < 3u; ++j) {
210 if (i > j) {
211 mat[i * 3u + j] = mat[j * 3u + i];
212 } else {
213 mat[i * 3u + j] = std::accumulate(
214 buf_pos.begin(), buf_pos.end(), 0.,
215 [i, j, &center](double acc, Utils::Vector3d const &pos) {
216 return acc + (pos[i] - center[i]) * (pos[j] - center[j]);
217 });
218 }
219 }
220 }
221 mat /= static_cast<double>(buf_pos.size());
222 }
223 return mat;
224}
225
227 int p_type) {
228 auto const &box_geo = *system.box_geo;
229 auto const &cell_structure = *system.cell_structure;
230 Utils::Vector9d mat{};
231 auto com = center_of_mass(system, p_type);
232 boost::mpi::broadcast(::comm_cart, com, 0);
233
234 for (auto const &p : cell_structure.local_particles()) {
235 if (p.type() == p_type and not p.is_virtual()) {
236 auto const pos = box_geo.unfolded_position(p.pos(), p.image_box()) - com;
237 auto const mass = p.mass();
238 mat[0] += mass * (pos[1] * pos[1] + pos[2] * pos[2]);
239 mat[4] += mass * (pos[0] * pos[0] + pos[2] * pos[2]);
240 mat[8] += mass * (pos[0] * pos[0] + pos[1] * pos[1]);
241 mat[1] -= mass * (pos[0] * pos[1]);
242 mat[2] -= mass * (pos[0] * pos[2]);
243 mat[5] -= mass * (pos[1] * pos[2]);
244 }
245 }
246 /* use symmetry */
247 mat[3] = mat[1];
248 mat[6] = mat[2];
249 mat[7] = mat[5];
250 return mat;
251}
252
253std::vector<int> nbhood(System::System const &system,
254 Utils::Vector3d const &pos, double dist) {
255 std::vector<int> buf_pid{};
256 auto const dist_sq = dist * dist;
257 auto const &box_geo = *system.box_geo;
258
259 for (auto const &p : system.cell_structure->local_particles()) {
260 auto const r_sq = box_geo.get_mi_vector(pos, p.pos()).norm2();
261 if (r_sq < dist_sq) {
262 buf_pid.push_back(p.id());
263 }
264 }
265
267 if (::comm_cart.rank() != 0) {
268 buf_pid.clear();
269 }
270
271 return buf_pid;
272}
273
274std::vector<std::vector<double>>
276 std::vector<int> const &p1_types,
277 std::vector<int> const &p2_types, double r_min,
278 double r_max, int r_bins, bool log_flag, bool int_flag) {
279
280 auto const &box_geo = *system.box_geo;
281 auto const trait_id = [](Particle const &p) { return p.id(); };
282 auto const trait_pos = [&box_geo](Particle const &p) {
283 return box_geo.unfolded_position(p.pos(), p.image_box());
284 };
285 auto const buf1 =
286 gather_traits_for_types(system, p1_types, trait_id, trait_pos);
287 auto const buf2 =
288 gather_traits_for_types(system, p2_types, trait_id, trait_pos);
289 auto const r_max2 = Utils::sqr(r_max);
290 auto const r_min2 = Utils::sqr(r_min);
291 auto const start_dist2 = Utils::sqr(r_max + 1.);
292 auto const inv_bin_width =
293 (log_flag) ? static_cast<double>(r_bins) / std::log(r_max / r_min)
294 : static_cast<double>(r_bins) / (r_max - r_min);
295
296 long cnt = 0;
297 double low = 0.0;
298 std::vector<double> distribution(r_bins);
299
300 for (auto const &[pid1, pos1] : buf1) {
301 auto min_dist2 = start_dist2;
302 /* particle loop: p2_types */
303 for (auto const &[pid2, pos2] : buf2) {
304 if (pid1 != pid2) {
305 auto const act_dist2 = box_geo.get_mi_vector(pos1, pos2).norm2();
306 if (act_dist2 < min_dist2) {
307 min_dist2 = act_dist2;
308 }
309 }
310 }
311 if (min_dist2 <= r_max2) {
312 if (min_dist2 >= r_min2) {
313 auto const min_dist = std::sqrt(min_dist2);
314 /* calculate bin index */
315 auto const ind = static_cast<int>(
316 ((log_flag) ? std::log(min_dist / r_min) : (min_dist - r_min)) *
317 inv_bin_width);
318 if (ind >= 0 and ind < r_bins) {
319 distribution[ind] += 1.0;
320 }
321 } else {
322 low += 1.0;
323 }
324 }
325 cnt++;
326 }
327
328 if (cnt != 0) {
329 // normalization
330 low /= static_cast<double>(cnt);
331 for (int i = 0; i < r_bins; i++) {
332 distribution[i] /= static_cast<double>(cnt);
333 }
334
335 // integration
336 if (int_flag) {
337 distribution[0] += low;
338 for (int i = 0; i < r_bins - 1; i++)
339 distribution[i + 1] += distribution[i];
340 }
341 }
342
343 std::vector<double> radii(r_bins);
344 if (log_flag) {
345 auto const log_fac = std::pow(r_max / r_min, 1. / r_bins);
346 radii[0] = r_min * std::sqrt(log_fac);
347 for (int i = 1; i < r_bins; ++i) {
348 radii[i] = radii[i - 1] * log_fac;
349 }
350 } else {
351 auto const bin_width = (r_max - r_min) / static_cast<double>(r_bins);
352 for (int i = 0; i < r_bins; ++i) {
353 radii[i] = r_min + bin_width / 2. + static_cast<double>(i) * bin_width;
354 }
355 }
356
357 return {radii, distribution};
358}
359
360std::vector<std::vector<double>>
361structure_factor(System::System const &system, std::vector<int> const &p_types,
362 int order) {
363 auto const &box_geo = *system.box_geo;
364 auto const trait_pos = [&box_geo](Particle const &p) {
365 return box_geo.unfolded_position(p.pos(), p.image_box());
366 };
367 auto const buf_pos = gather_traits_for_types(system, p_types, trait_pos);
368 auto const order_sq = Utils::sqr(static_cast<std::size_t>(order));
369 auto const twoPI_L = 2. * std::numbers::pi * system.box_geo->length_inv()[0];
370 std::vector<double> ff(2ul * order_sq + 1ul);
371 std::vector<double> wavevectors;
372 std::vector<double> intensities;
373
374 if (::comm_cart.rank() == 0) {
375 for (int i = 0; i <= order; i++) {
376 for (int j = -order; j <= order; j++) {
377 for (int k = -order; k <= order; k++) {
378 auto const n = i * i + j * j + k * k;
379 if ((static_cast<std::size_t>(n) <= order_sq) && (n >= 1)) {
380 double C_sum = 0.0, S_sum = 0.0;
381 for (auto const &pos : buf_pos) {
382 auto const qr = twoPI_L * (Utils::Vector3i{{i, j, k}} * pos);
383 C_sum += cos(qr);
384 S_sum += sin(qr);
385 }
386 ff[2 * n - 2] += C_sum * C_sum + S_sum * S_sum;
387 ff[2 * n - 1]++;
388 }
389 }
390 }
391 }
392
393 std::size_t length = 0;
394 for (std::size_t qi = 0; qi < order_sq; qi++) {
395 if (ff[2 * qi + 1] != 0) {
396 ff[2 * qi] /= static_cast<double>(buf_pos.size()) * ff[2 * qi + 1];
397 ++length;
398 }
399 }
400
401 wavevectors.resize(length);
402 intensities.resize(length);
403
404 int cnt = 0;
405 for (std::size_t i = 0; i < order_sq; i++) {
406 if (ff[2 * i + 1] != 0) {
407 wavevectors[cnt] = twoPI_L * std::sqrt(static_cast<double>(i + 1));
408 intensities[cnt] = ff[2 * i];
409 cnt++;
410 }
411 }
412 }
413
414 return {std::move(wavevectors), std::move(intensities)};
415}
Vector implementation and trait types for boost qvm interoperability.
Main system class.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
__device__ void vector_product(float const *a, float const *b, float *out)
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
bool contains(InputIt first, InputIt last, T const &value)
Check whether an iterator range contains a value.
Definition contains.hpp:36
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
void serialize(Archive &ar, std::tuple< T... > &pack, unsigned int const)
Serialize std::tuple.
Utils::Vector3d center_of_mass(System::System const &system, int p_type)
Calculate the center of mass of particles of a certain type.
Utils::Vector3d angular_momentum(System::System const &system, int p_type)
Calculate the angular momentum of particles of a certain type.
std::vector< int > nbhood(System::System const &system, Utils::Vector3d const &pos, double dist)
Find all particles within a given radius dist around a position pos.
Utils::Vector9d moment_of_inertia_matrix(System::System const &system, int p_type)
Calculate the moment of inertia of particles of a certain type.
Utils::Vector9d gyration_tensor(System::System const &system, std::vector< int > const &p_types)
Calculate the gyration tensor of particles of certain types.
Utils::Vector3d calc_linear_momentum(System::System const &system, bool include_particles, bool include_lbfluid)
Calculate total momentum of the system (particles & LB fluid).
std::vector< std::vector< double > > structure_factor(System::System const &system, std::vector< int > const &p_types, int order)
Calculate the spherically averaged structure factor.
double mindist(System::System const &system, std::vector< int > const &set1, std::vector< int > const &set2)
Calculate the minimal distance of two particles with types in set1 and set2, respectively.
static auto gather_traits_for_types(System::System const &system, std::vector< int > const &p_types, Trait &&...trait)
Gather particle traits to MPI rank 0.
std::vector< std::vector< double > > calc_part_distribution(System::System const &system, std::vector< int > const &p1_types, std::vector< int > const &p2_types, double r_min, double r_max, int r_bins, bool log_flag, bool int_flag)
Calculate the distribution of particles around others.
Statistical tools to analyze simulations.
std::invoke_result_t< F, Particle const & > type
Decay a tuple of only 1 type to that type.
std::tuple< std::invoke_result_t< F, Particle const & >... > type
bool is_solver_set() const
Return true if a LB solver is active.
Definition lb/Solver.cpp:64
auto get_lattice_speed() const
Get the lattice speed (agrid/tau).
Utils::Vector3d get_momentum() const
Struct holding all information for one particle.
Definition Particle.hpp:395