30#include "communication.hpp"
31#include "system/System.hpp"
37#include <boost/mpi/collectives/all_reduce.hpp>
43#include <unordered_map>
52 std::vector<int> buffer_pid{};
53 std::vector<T> buffer_obs{};
60 if (ptr !=
nullptr and not ptr->is_ghost()) {
61 buffer_pid.emplace_back(pid);
62 buffer_obs.emplace_back(
kernel(*ptr));
67 std::unordered_map<int, T> map{};
71 for (std::size_t i = 0u; i < buffer_pid.size(); ++i) {
72 map[buffer_pid[i]] = buffer_obs[i];
102 int chain_length,
int n_chains) {
105 double dist = 0.0, dist2 = 0.0, dist4 = 0.0;
106 std::array<double, 4> re{};
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);
114 auto const map = prefetch.join();
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();
122 dist4 += norm2 * norm2;
124 auto const tmp =
static_cast<double>(n_chains);
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]));
134 int chain_length,
int n_chains) {
139 double r_G = 0.0, r_G2 = 0.0, r_G4 = 0.0;
140 std::array<double, 4> rg{};
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()) {
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.");
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);
164 auto const map_pos = prefetch_pos.join();
165 auto const map_com = prefetch_com.join();
166 auto const map_mass = prefetch_mass.join();
168 for (
int i = 0; i < n_chains; i++) {
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);
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;
183 tmp /=
static_cast<double>(chain_length);
188 auto const tmp =
static_cast<double>(n_chains);
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]));
198 int chain_length,
int n_chains) {
201 double r_H = 0.0, r_H2 = 0.0;
202 std::array<double, 2> rh{};
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);
215 auto const map = prefetch.join();
217 for (
int p = 0; p < n_chains; p++) {
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();
225 auto const tmp = prefac / ri;
229 auto const tmp =
static_cast<double>(n_chains);
231 rh[1] = (n_chains == 1) ? 0. : std::sqrt(r_H2 / tmp -
Utils::sqr(rh[0]));
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.
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.
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.
auto const & mass() const
auto const & image_box() const
DEVICE_QUALIFIER constexpr reference at(size_type i)