ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
statistics_chain.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 * Implementation of \ref statistics_chain.hpp "statistics_chain.hpp".
23 */
24
26
27#include "BoxGeometry.hpp"
28#include "Particle.hpp"
30#include "communication.hpp"
31#include "system/System.hpp"
32
33#include <utils/Vector.hpp>
34#include <utils/math/sqr.hpp>
36
37#include <boost/mpi/collectives/all_reduce.hpp>
38
39#include <array>
40#include <cmath>
41#include <functional>
42#include <stdexcept>
43#include <unordered_map>
44#include <vector>
45
46/**
47 * @brief Gather particle properties (or any derived quantities) on MPI rank 0.
48 * The data is collected into an unordered map.
49 */
50template <typename T> struct GatherParticleTraits {
51private:
52 std::vector<int> buffer_pid{};
53 std::vector<T> buffer_obs{};
54
55public:
56 virtual T kernel(Particle const &) const = 0;
57
58 void fetch(CellStructure const &cell_structure, int pid) {
59 auto const ptr = cell_structure.get_local_particle(pid);
60 if (ptr != nullptr and not ptr->is_ghost()) {
61 buffer_pid.emplace_back(pid);
62 buffer_obs.emplace_back(kernel(*ptr));
63 }
64 }
65
66 auto join() {
67 std::unordered_map<int, T> map{};
70 if (::comm_cart.rank() == 0) {
71 for (std::size_t i = 0u; i < buffer_pid.size(); ++i) {
72 map[buffer_pid[i]] = buffer_obs[i];
73 }
74 }
75 buffer_pid.clear();
76 buffer_obs.clear();
77 return map;
78 }
79};
80
81struct GatherPos : public GatherParticleTraits<Utils::Vector3d> {
83 GatherPos(BoxGeometry const &box_geo) : m_box_geo{box_geo} {}
84 Utils::Vector3d kernel(Particle const &p) const override {
86 }
87};
88
89struct GatherCom : public GatherParticleTraits<Utils::Vector3d> {
91 GatherCom(BoxGeometry const &box_geo) : m_box_geo{box_geo} {}
92 Utils::Vector3d kernel(Particle const &p) const override {
93 return m_box_geo.unfolded_position(p.pos(), p.image_box()) * p.mass();
94 }
95};
96
97struct GatherMass : public GatherParticleTraits<double> {
98 double kernel(Particle const &p) const override { return p.mass(); }
99};
100
101std::array<double, 4> calc_re(System::System const &system, int chain_start,
102 int chain_length, int n_chains) {
103 auto const &cell_structure = *system.cell_structure;
104 GatherPos prefetch{*system.box_geo};
105 double dist = 0.0, dist2 = 0.0, dist4 = 0.0;
106 std::array<double, 4> re{};
107
108 for (int i = 0; i < n_chains; i++) {
109 auto const pid2 = chain_start + i * chain_length;
110 auto const pid1 = pid2 + chain_length - 1;
111 prefetch.fetch(cell_structure, pid1);
112 prefetch.fetch(cell_structure, pid2);
113 }
114 auto const map = prefetch.join();
115 if (::comm_cart.rank() == 0) {
116 for (int i = 0; i < n_chains; i++) {
117 auto const pid2 = chain_start + i * chain_length;
118 auto const pid1 = pid2 + chain_length - 1;
119 auto const norm2 = (map.at(pid1) - map.at(pid2)).norm2();
120 dist += sqrt(norm2);
121 dist2 += norm2;
122 dist4 += norm2 * norm2;
123 }
124 auto const tmp = static_cast<double>(n_chains);
125 re[0] = dist / tmp;
126 re[2] = dist2 / tmp;
127 re[1] = (n_chains == 1) ? 0. : std::sqrt(re[2] - Utils::sqr(re[0]));
128 re[3] = (n_chains == 1) ? 0. : std::sqrt(dist4 / tmp - Utils::sqr(re[2]));
129 }
130 return re;
131}
132
133std::array<double, 4> calc_rg(System::System const &system, int chain_start,
134 int chain_length, int n_chains) {
135 auto const &cell_structure = *system.cell_structure;
136 GatherPos prefetch_pos{*system.box_geo};
137 GatherCom prefetch_com{*system.box_geo};
138 GatherMass prefetch_mass{};
139 double r_G = 0.0, r_G2 = 0.0, r_G4 = 0.0;
140 std::array<double, 4> rg{};
141
142 auto has_virtual = false;
143 for (int i = 0; i < n_chains * chain_length; ++i) {
144 auto const pid = chain_start + i;
145 auto const ptr = cell_structure.get_local_particle(pid);
146 if (ptr != nullptr and not ptr->is_ghost() and ptr->is_virtual()) {
147 has_virtual = true;
148 break;
149 }
150 }
151 if (boost::mpi::all_reduce(::comm_cart, has_virtual, std::logical_or<>{})) {
152 throw std::runtime_error(
153 "Center of mass is not well-defined for chains including virtual "
154 "sites. Virtual sites do not have a meaningful mass.");
155 }
156
157 for (int i = 0; i < n_chains * chain_length; ++i) {
158 auto const pid = chain_start + i;
159 prefetch_com.fetch(cell_structure, pid);
160 prefetch_pos.fetch(cell_structure, pid);
161 prefetch_mass.fetch(cell_structure, pid);
162 }
163
164 auto const map_pos = prefetch_pos.join();
165 auto const map_com = prefetch_com.join();
166 auto const map_mass = prefetch_mass.join();
167 if (::comm_cart.rank() == 0) {
168 for (int i = 0; i < n_chains; i++) {
169 double M = 0.0;
170 Utils::Vector3d r_CM{};
171 for (int j = 0; j < chain_length; j++) {
172 auto const pid = chain_start + i * chain_length + j;
173 r_CM += map_com.at(pid);
174 M += map_mass.at(pid);
175 }
176 r_CM /= M;
177 double tmp = 0.0;
178 for (int j = 0; j < chain_length; ++j) {
179 auto const pid = chain_start + i * chain_length + j;
180 auto const d = map_pos.at(pid) - r_CM;
181 tmp += d.norm2();
182 }
183 tmp /= static_cast<double>(chain_length);
184 r_G += sqrt(tmp);
185 r_G2 += tmp;
186 r_G4 += tmp * tmp;
187 }
188 auto const tmp = static_cast<double>(n_chains);
189 rg[0] = r_G / tmp;
190 rg[2] = r_G2 / tmp;
191 rg[1] = (n_chains == 1) ? 0. : std::sqrt(rg[2] - Utils::sqr(rg[0]));
192 rg[3] = (n_chains == 1) ? 0. : std::sqrt(r_G4 / tmp - Utils::sqr(rg[2]));
193 }
194 return rg;
195}
196
197std::array<double, 2> calc_rh(System::System const &system, int chain_start,
198 int chain_length, int n_chains) {
199 auto const &cell_structure = *system.cell_structure;
200 GatherPos prefetch{*system.box_geo};
201 double r_H = 0.0, r_H2 = 0.0;
202 std::array<double, 2> rh{};
203
204 auto const chain_l = static_cast<double>(chain_length);
205 auto const prefac = 0.5 * chain_l * (chain_l - 1.);
206 for (int p = 0; p < n_chains; p++) {
207 for (int i = chain_start + chain_length * p;
208 i < chain_start + chain_length * (p + 1); i++) {
209 prefetch.fetch(cell_structure, i);
210 for (int j = i + 1; j < chain_start + chain_length * (p + 1); j++) {
211 prefetch.fetch(cell_structure, j);
212 }
213 }
214 }
215 auto const map = prefetch.join();
216 if (::comm_cart.rank() == 0) {
217 for (int p = 0; p < n_chains; p++) {
218 double ri = 0.0;
219 for (int i = chain_start + chain_length * p;
220 i < chain_start + chain_length * (p + 1); i++) {
221 for (int j = i + 1; j < chain_start + chain_length * (p + 1); j++) {
222 ri += 1.0 / (map.at(i) - map.at(j)).norm();
223 }
224 }
225 auto const tmp = prefac / ri;
226 r_H += tmp;
227 r_H2 += tmp * tmp;
228 }
229 auto const tmp = static_cast<double>(n_chains);
230 rh[0] = r_H / tmp;
231 rh[1] = (n_chains == 1) ? 0. : std::sqrt(r_H2 / tmp - Utils::sqr(rh[0]));
232 }
233 return rh;
234}
Vector implementation and trait types for boost qvm interoperability.
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
Main system class.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
std::array< double, 4 > calc_rg(System::System const &system, int chain_start, int chain_length, int n_chains)
Calculate the radius of gyration.
std::array< double, 2 > calc_rh(System::System const &system, int chain_start, int chain_length, int n_chains)
Calculate the hydrodynamic radius (ref.
std::array< double, 4 > calc_re(System::System const &system, int chain_start, int chain_length, int n_chains)
Calculate the end-to-end-distance.
This file contains the code for statistics on chains.
Describes a cell structure / cell system.
Particle * get_local_particle(int id)
Get a local particle by id.
GatherCom(BoxGeometry const &box_geo)
BoxGeometry const & m_box_geo
Utils::Vector3d kernel(Particle const &p) const override
double kernel(Particle const &p) const override
Gather particle properties (or any derived quantities) on MPI rank 0.
virtual T kernel(Particle const &) const =0
void fetch(CellStructure const &cell_structure, int pid)
BoxGeometry const & m_box_geo
Utils::Vector3d kernel(Particle const &p) const override
GatherPos(BoxGeometry const &box_geo)
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & mass() const
Definition Particle.hpp:452
auto const & image_box() const
Definition Particle.hpp:444
auto const & pos() const
Definition Particle.hpp:431
DEVICE_QUALIFIER constexpr reference at(size_type i)
Definition Array.hpp:86