ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
ClusterStructure.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 "ClusterStructure.hpp"
20
21#include "BoxGeometry.hpp"
22#include "Cluster.hpp"
23#include "PartCfg.hpp"
24#include "errorhandling.hpp"
25#include "particle_node.hpp"
26
28
29#include <algorithm>
30#include <memory>
31#include <stdexcept>
32#include <utility>
33#include <vector>
34
35namespace ClusterAnalysis {
36
38
40 clusters.clear();
41 cluster_id.clear();
42 m_cluster_identities.clear();
43}
44
46 return cluster_id.find(p.id()) != cluster_id.end();
47}
48
49// Analyze the cluster structure of the given particles
51 // clear data structs
52 clear();
53 sanity_checks();
54
55 // Iterate over pairs
56 auto const box_geo_handle = get_box_geo();
57 auto const &box_geo = *box_geo_handle;
58 PartCfg partCfg{box_geo};
59 Utils::for_each_pair(partCfg.begin(), partCfg.end(),
60 [this](const Particle &p1, const Particle &p2) {
61 this->add_pair(p1, p2);
62 });
63 merge_clusters();
64}
65
67 clear();
68 sanity_checks();
69 auto const box_geo_handle = get_box_geo();
70 auto const &box_geo = *box_geo_handle;
71 PartCfg partCfg{box_geo};
72 for (const auto &p : partCfg) {
73 for (auto const bond : p.bonds()) {
74 if (bond.partner_ids().size() == 1) {
75 add_pair(p, get_particle_data(bond.partner_ids()[0]));
76 }
77 }
78 }
79
80 merge_clusters();
81}
82
83void ClusterStructure::add_pair(const Particle &p1, const Particle &p2) {
84 // * check, if there's a neighbor
85 // * No: Then go on to the next particle
86 // * Yes: Then if
87 // * One of them belongs to a cluster, give the other one the same cluster
88 // id.
89 // * None of them belongs to a cluster: Give them both a new cluster id
90 // * Both belong to different clusters: Mark the clusters as identical
91 // * so that they can be put together later
92 if (!m_pair_criterion) {
93 runtimeErrorMsg() << "No cluster criterion defined";
94 return;
95 }
96 // If the two particles are neighbors...
97 if (m_pair_criterion->decide(p1, p2)) {
98
99 if // None belongs to a cluster
100 ((!part_of_cluster(p1)) && (!part_of_cluster(p2))) {
101 // Both particles belong to the same, new cluster
102 const int cid = get_next_free_cluster_id();
103
104 // assign the cluster_ids
105 cluster_id[p1.id()] = cid;
106 cluster_id[p2.id()] = cid;
107 } else if // p2 belongs to a cluster but p1 doesn't
108 (part_of_cluster(p2) && !part_of_cluster(p1)) {
109 // Give p1 the same cluster id as p2
110 cluster_id[p1.id()] = find_id_for(cluster_id.at(p2.id()));
111 } else if // i belongs to a cluster but j doesn't
112 (part_of_cluster(p1) && !part_of_cluster(p2)) {
113 // give p2 the cluster id from p1
114 cluster_id[p2.id()] = find_id_for(cluster_id.at(p1.id()));
115 } else if // Both belong to different clusters
116 (part_of_cluster(p1) && part_of_cluster(p2) &&
117 cluster_id.at(p1.id()) != cluster_id.at(p2.id())) {
118 // Clusters of p1 and p2 are one and the same. Add an identity to the list
119 // The higher number must be inserted as first value of the pair
120 // because the substitutions later have to be done in descending order
121 const int cid1 = find_id_for(cluster_id.at(p1.id()));
122 const int cid2 = find_id_for(cluster_id.at(p2.id()));
123 if (cid1 > cid2) {
124 m_cluster_identities[cid1] = cid2;
125 } else if (cid1 < cid2) {
126 m_cluster_identities[cid2] = cid1;
127 }
128 // else do nothing. The clusters are already noted for merging.
129 // Connected clusters will be merged later
130 }
131 // The case for both particles being in the same cluster does not need to be
132 // treated.
133 }
134}
135
136void ClusterStructure::merge_clusters() {
137 // Relabel particles according to the cluster identities map
138 // Also create empty cluster objects for the final cluster id
139
140 // Collect needed changes in a separate map, as doing the changes on the fly
141 // would screw up the iterators
142 std::vector<std::pair<int, int>> to_be_changed;
143
144 for (auto it : cluster_id) {
145 // particle id is in it.first and cluster id in it.second
146 // We change the cluster id according to the cluster identities
147 // map
148 const int cid = find_id_for(it.second);
149 // We note the list of changes here, so we don't modify the map
150 // while iterating
151 to_be_changed.emplace_back(it.first, cid);
152 // Empty cluster object
153 if (clusters.find(cid) == clusters.end()) {
154 clusters[cid] = std::make_shared<Cluster>(m_box_geo);
155 }
156 }
157
158 // Now act on the changes marke in above iteration
159 for (auto it : to_be_changed) {
160 cluster_id[it.first] = it.second;
161 }
162
163 // Now fill the cluster objects with particle ids
164 // Iterate over particles, fill in the cluster map
165 // to each cluster particle the corresponding cluster id
166 for (auto it : cluster_id) {
167 // If this is the first particle in this cluster, instance a new cluster
168 // object
169 if (clusters.find(it.second) == clusters.end()) {
170 clusters[it.second] = std::make_shared<Cluster>(m_box_geo);
171 }
172 clusters[it.second]->particles.push_back(it.first);
173 }
174
175 // Sort particles ids in the clusters
176 for (const auto &c : clusters) {
177 std::sort(c.second->particles.begin(), c.second->particles.end());
178 }
179}
180
181int ClusterStructure::find_id_for(int x) {
182 int tmp = x;
183 while (m_cluster_identities.find(tmp) != m_cluster_identities.end()) {
184 tmp = m_cluster_identities[tmp];
185 }
186 return tmp;
187}
188
189int ClusterStructure::get_next_free_cluster_id() {
190 // iterate over cluster_id'
191 int max_seen_cluster = 0;
192 for (auto it : cluster_id) {
193 int cid = it.second;
194 if (max_seen_cluster < cid) {
195 max_seen_cluster = cid;
196 }
197 }
198 return max_seen_cluster + 1;
199}
200
201void ClusterStructure::sanity_checks() const {
202 if (get_box_geo()->type() != BoxType::CUBOID) {
203 throw std::runtime_error(
204 "Cluster analysis is not compatible with non-cuboid box types");
205 }
206}
207
208} // namespace ClusterAnalysis
std::map< int, int > cluster_id
Map between particle ids and corresponding cluster ids.
void clear()
Clear data structures.
bool part_of_cluster(const Particle &p)
Is particle p part of a cluster.
void run_for_bonded_particles()
Run cluster analysis, consider pairs of particles connected by a bonded interaction.
void run_for_all_pairs()
Run cluster analysis, consider all particle pairs.
std::map< int, std::shared_ptr< Cluster > > clusters
Map holding the individual clusters.
Particle cache on the head node.
Definition PartCfg.hpp:34
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
void for_each_pair(ForwardIterator first, ForwardIterator last, BinaryOp op)
Execute op for each pair of elements in [first, last) once.
const Particle & get_particle_data(int p_id)
Get particle data.
Particles creation and deletion.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & id() const
Definition Particle.hpp:414