ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
Cluster.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
20#include "config/config.hpp"
21
22#ifdef GSL
23#include "gsl/gsl_fit.h"
24#endif
25
26#include "Cluster.hpp"
27
28#include "BoxGeometry.hpp"
29#include "Particle.hpp"
30#include "errorhandling.hpp"
31#include "particle_node.hpp"
32
33#include <utils/Vector.hpp>
34
35#include <algorithm>
36#include <cmath>
37#include <cstddef>
38#include <numeric>
39#include <stdexcept>
40#include <utility>
41#include <vector>
42
43namespace ClusterAnalysis {
44
45// Center of mass of an aggregate
49
50// Center of mass of an aggregate
52Cluster::center_of_mass_subcluster(std::vector<int> const &particle_ids) {
53 sanity_checks();
54 auto const box_geo_handle = get_box_geo();
55 auto const &box_geo = *box_geo_handle;
56 Utils::Vector3d com{};
57
58 // The distances between the particles are "folded", such that all distances
59 // are smaller than box_l/2 in a periodic system. The 1st particle
60 // of the cluster is arbitrarily chosen as reference.
61
62 auto const reference_position =
63 box_geo.folded_position(get_particle_data(particles[0]).pos());
64 double total_mass = 0.;
65 for (int pid : particle_ids) {
66 auto const folded_pos =
67 box_geo.folded_position(get_particle_data(pid).pos());
68 auto const dist_to_reference =
69 box_geo.get_mi_vector(folded_pos, reference_position);
70 com += dist_to_reference * get_particle_data(pid).mass();
71 total_mass += get_particle_data(pid).mass();
72 }
73
74 // Normalize by number of particles
75 com /= total_mass;
76
77 // Re-add reference position
78 com += reference_position;
79
80 // Fold into simulation box
81 return box_geo.folded_position(com);
82}
83
85 sanity_checks();
86 auto const box_geo_handle = get_box_geo();
87 auto const &box_geo = *box_geo_handle;
88 double ld = 0.;
89 for (auto a = particles.begin(); a != particles.end(); a++) {
90 for (auto b = a; ++b != particles.end();) {
91 auto const dist = box_geo
92 .get_mi_vector(get_particle_data(*a).pos(),
94 .norm();
95
96 // Larger than previous largest distance?
97 ld = std::max(ld, dist);
98 }
99 }
100 return ld;
101}
102
103// Radius of gyration
107
108double
109Cluster::radius_of_gyration_subcluster(std::vector<int> const &particle_ids) {
110 sanity_checks();
111 auto const box_geo_handle = get_box_geo();
112 auto const &box_geo = *box_geo_handle;
113 // Center of mass
114 Utils::Vector3d com = center_of_mass_subcluster(particle_ids);
115 double sum_sq_dist = 0.;
116 for (auto const pid : particle_ids) {
117 // calculate square length of this distance
118 sum_sq_dist +=
119 box_geo.get_mi_vector(com, get_particle_data(pid).pos()).norm2();
120 }
121
122 return sqrt(sum_sq_dist / static_cast<double>(particle_ids.size()));
123}
124
125template <typename T>
126std::vector<std::size_t> sort_indices(const std::vector<T> &v) {
127
128 // Unsorted for unsorted vector (0..n-1)
129 std::vector<std::size_t> idx(v.size());
130 std::iota(idx.begin(), idx.end(), 0);
131
132 // sort indices based on comparing values in v
133 std::sort(idx.begin(), idx.end(),
134 [&v](std::size_t i1, std::size_t i2) { return v[i1] < v[i2]; });
135 return idx;
136}
137
138std::pair<double, double> Cluster::fractal_dimension(double dr) {
139#ifdef GSL
140 sanity_checks();
141 auto const box_geo_handle = get_box_geo();
142 auto const &box_geo = *box_geo_handle;
144 // calculate Df using linear regression on the logarithms of the radii of
145 // gyration against the number of particles in sub-clusters. Particles are
146 // included step by step from the center of mass outwards
147
148 // Distances of particles from the center of mass
149 std::vector<double> distances;
150
151 for (auto const &it : particles) {
152 distances.push_back(box_geo.get_mi_vector(com, get_particle_data(it).pos())
153 .norm()); // add distance from the current particle
154 // to the com in the distances vectors
155 }
156
157 // Get particle indices in the cluster which yield distances sorted in
158 // ascending order from center of mass.
159 auto particle_indices = sort_indices(distances);
160
161 // Particle ids in the current sub-cluster
162 std::vector<int> subcluster_ids;
163
164 std::vector<double> log_pcounts; // particle count
165 std::vector<double> log_diameters; // corresponding radii of gyration
166 double last_dist = 0;
167 for (auto const idx : particle_indices) {
168 subcluster_ids.push_back(particles[idx]);
169 if (distances[idx] < last_dist + dr)
170 continue;
171
172 last_dist = distances[idx];
173 if (subcluster_ids.size() == 1)
174 continue;
175 auto const current_rg = radius_of_gyration_subcluster(subcluster_ids);
176 log_pcounts.push_back(log(static_cast<double>(subcluster_ids.size())));
177 log_diameters.push_back(log(current_rg * 2.0));
178 }
179 double c0, c1, cov00, cov01, cov11, sumsq;
180 gsl_fit_linear(log_diameters.data(), 1, log_pcounts.data(), 1,
181 log_pcounts.size(), &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
182 double mean_sq_residual = sumsq / static_cast<double>(log_diameters.size());
183 return {c1, mean_sq_residual};
184#else
185 runtimeErrorMsg() << "GSL (gnu scientific library) is required for fractal "
186 "dimension calculation.";
187 return {0, 0};
188#endif
189}
190
191void Cluster::sanity_checks() const {
192 auto const box_geo_handle = get_box_geo();
193 auto const &box_geo = *box_geo_handle;
194 if (box_geo.type() != BoxType::CUBOID) {
195 throw std::runtime_error(
196 "Cluster analysis is not compatible with non-cuboid box types");
197 }
198}
199
200} // namespace ClusterAnalysis
Vector implementation and trait types for boost qvm interoperability.
Utils::Vector3d center_of_mass()
Definition Cluster.cpp:46
std::pair< double, double > fractal_dimension(double dr)
Calculate the fractal dimension N(r) via r^d, where N(r) counts the number of particles in a sphere o...
Definition Cluster.cpp:138
Utils::Vector3d center_of_mass_subcluster(std::vector< int > const &particle_ids)
Calculate the center of mass of the cluster.
Definition Cluster.cpp:52
std::vector< int > particles
Ids of the particles in the cluster.
double longest_distance()
Longest distance between any combination of two particles.
Definition Cluster.cpp:84
double radius_of_gyration()
Calculate radius of gyration of the cluster.
Definition Cluster.cpp:104
double radius_of_gyration_subcluster(std::vector< int > const &particle_ids)
Definition Cluster.cpp:109
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
std::vector< std::size_t > sort_indices(const std::vector< T > &v)
Definition Cluster.cpp:126
const Particle & get_particle_data(int p_id)
Get particle data.
Particles creation and deletion.
auto const & mass() const
Definition Particle.hpp:452
auto const & pos() const
Definition Particle.hpp:431