ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
Analysis.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2013-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 "Analysis.hpp"
21#include "ObservableStat.hpp"
22
23#include "core/BoxGeometry.hpp"
27#include "core/cells.hpp"
29#include "core/dpd.hpp"
31
33
34#include <utils/Vector.hpp>
35#include <utils/contains.hpp>
37
38#include <boost/mpi/collectives/all_reduce.hpp>
39
40#include <algorithm>
41#include <cmath>
42#include <functional>
43#include <memory>
44#include <sstream>
45#include <stdexcept>
46#include <string>
47#include <vector>
48
49namespace ScriptInterface {
50namespace Analysis {
51
52/** @brief Check if a contiguous range of particle ids exists. */
53static void check_topology(CellStructure const &cell_structure, int chain_start,
54 int chain_length, int n_chains) {
55 try {
56 if (n_chains <= 0) {
57 throw std::domain_error("Chain analysis needs at least 1 chain");
58 }
59 if (chain_length <= 0) {
60 throw std::domain_error("Chain analysis needs at least 1 bead per chain");
61 }
62
63 auto n_particles_local = 0;
64 for (int i = 0; i < n_chains; ++i) {
65 for (int j = 0; j < chain_length; ++j) {
66 auto const pid = chain_start + i * chain_length + j;
67 auto ptr = cell_structure.get_local_particle(pid);
68 if (ptr != nullptr and not ptr->is_ghost()) {
70 }
71 }
72 }
73 auto const n_particles_total = boost::mpi::all_reduce(
74 ::comm_cart, n_particles_local, std::plus<int>{});
75
77 for (int i = 0; i < chain_length * n_chains; ++i) {
78 auto const pid = chain_start + i;
79 auto ptr = cell_structure.get_local_particle(pid);
80 int local_count = 0;
81 if (ptr != nullptr and not ptr->is_ghost()) {
82 local_count = 1;
83 }
84 auto const total_count =
85 boost::mpi::all_reduce(::comm_cart, local_count, std::plus<int>{});
86 if (total_count == 0) {
87 std::stringstream error_msg;
88 error_msg << "Particle with id " << pid << " does not exist; "
89 << "cannot perform analysis on the range chain_start="
90 << chain_start << ", number_of_chains=" << n_chains
91 << ", chain_length=" << chain_length << ". "
92 << "Please provide a contiguous range of particle ids.";
93 throw std::runtime_error(error_msg.str());
94 }
95 }
96 }
97 } catch (...) {
98 if (::comm_cart.rank() == 0) {
99 throw;
100 }
101 throw Exception("");
102 }
103}
104
105void Analysis::check_particle_type(int p_type) const {
106 auto const &nonbonded_ias = get_system().nonbonded_ias;
107 if (p_type < 0 or p_type > nonbonded_ias->get_max_seen_particle_type()) {
108 std::stringstream error_msg;
109 error_msg << "Particle type " << p_type << " does not exist";
110 throw std::invalid_argument(error_msg.str());
111 }
112}
113
114Variant Analysis::do_call_method(std::string const &name,
115 VariantMap const &parameters) {
116 if (name == "linear_momentum") {
117 auto const local = calc_linear_momentum(
118 get_system(), get_value_or<bool>(parameters, "include_particles", true),
119 get_value_or<bool>(parameters, "include_lbfluid", true));
120 return mpi_reduce_sum(context()->get_comm(), local).as_vector();
121 }
122 if (name == "particle_energy") {
123 auto &system = get_system();
124 auto const pid = get_value<int>(parameters, "pid");
125 auto const local = system.particle_short_range_energy_contribution(pid);
126 return mpi_reduce_sum(context()->get_comm(), local);
127 }
128#ifdef DIPOLE_FIELD_TRACKING
129 if (name == "calc_long_range_fields") {
130 get_system().calculate_long_range_fields();
131 return {};
132 }
133#endif
134 if (name == "particle_neighbor_pids") {
135 auto &system = get_system();
136 system.on_observable_calc();
137 std::unordered_map<int, std::vector<int>> dict;
138 context()->parallel_try_catch([&]() {
139 auto neighbor_pids = get_neighbor_pids(system);
140 Utils::Mpi::gather_buffer(neighbor_pids, context()->get_comm());
141 std::for_each(neighbor_pids.begin(), neighbor_pids.end(),
142 [&dict](NeighborPIDs const &neighbor_pid) {
143 dict[neighbor_pid.pid] = neighbor_pid.neighbor_pids;
144 });
145 });
147 }
148#ifdef DPD
149 if (name == "dpd_stress") {
150 auto const result = dpd_stress(context()->get_comm());
151 return result.as_vector();
152 }
153#endif // DPD
154 if (name == "min_dist") {
155 auto const p_types1 = get_value<std::vector<int>>(parameters, "p_types1");
156 auto const p_types2 = get_value<std::vector<int>>(parameters, "p_types2");
157 for (auto const p_type : p_types1) {
158 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
159 }
160 for (auto const p_type : p_types2) {
161 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
162 }
164 }
165 if (name == "center_of_mass") {
166 auto const p_type = get_value<int>(parameters, "p_type");
167 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
168 auto const local = center_of_mass(get_system(), p_type);
169 return mpi_reduce_sum(context()->get_comm(), local).as_vector();
170 }
171 if (name == "angular_momentum") {
172 auto const p_type = get_value<int>(parameters, "p_type");
173 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
174 auto const local = angular_momentum(get_system(), p_type);
175 return mpi_reduce_sum(context()->get_comm(), local).as_vector();
176 }
177 if (name == "nbhood") {
178 auto const pos = get_value<Utils::Vector3d>(parameters, "pos");
179 auto const radius = get_value<double>(parameters, "r_catch");
180 auto const result = nbhood(get_system(), pos, radius);
181 return result;
182 }
183 if (name == "calc_re") {
184 auto const &system = get_system();
185 auto const chain_start = get_value<int>(parameters, "chain_start");
186 auto const chain_length = get_value<int>(parameters, "chain_length");
187 auto const n_chains = get_value<int>(parameters, "number_of_chains");
189 auto const result = calc_re(system, chain_start, chain_length, n_chains);
190 return std::vector<double>(result.begin(), result.end());
191 }
192 if (name == "calc_rg") {
193 auto const &system = get_system();
194 auto const chain_start = get_value<int>(parameters, "chain_start");
195 auto const chain_length = get_value<int>(parameters, "chain_length");
196 auto const n_chains = get_value<int>(parameters, "number_of_chains");
199 context()->parallel_try_catch([&]() {
200 auto const result = calc_rg(system, chain_start, chain_length, n_chains);
201 output = Variant{std::vector<double>(result.begin(), result.end())};
202 });
203 return output;
204 }
205 if (name == "calc_rh") {
206 auto const &system = get_system();
207 auto const chain_start = get_value<int>(parameters, "chain_start");
208 auto const chain_length = get_value<int>(parameters, "chain_length");
209 auto const n_chains = get_value<int>(parameters, "number_of_chains");
211 auto const result = calc_rh(system, chain_start, chain_length, n_chains);
212 return std::vector<double>(result.begin(), result.end());
213 }
214 if (name == "gyration_tensor") {
215 auto const p_types = get_value<std::vector<int>>(parameters, "p_types");
216 for (auto const p_type : p_types) {
217 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
218 }
219 auto const mat = gyration_tensor(get_system(), p_types);
220 return std::vector<double>(mat.begin(), mat.end());
221 }
222 if (name == "moment_of_inertia_matrix") {
223 auto const p_type = get_value<int>(parameters, "p_type");
224 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
225 auto const local = moment_of_inertia_matrix(get_system(), p_type);
226 return mpi_reduce_sum(context()->get_comm(), local).as_vector();
227 }
228 if (name == "structure_factor") {
229 auto const order = get_value<int>(parameters, "sf_order");
230 auto const p_types = get_value<std::vector<int>>(parameters, "sf_types");
232 if (order < 1)
233 throw std::domain_error("order has to be a strictly positive number");
234 });
235 for (auto const p_type : p_types) {
236 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
237 }
238 auto const result = structure_factor(get_system(), p_types, order);
239 return make_vector_of_variants(result);
240 }
241 if (name == "distribution") {
242 auto const &box_l = get_system().box_geo->length();
243 auto const r_max_limit =
244 0.5 * std::min(std::min(box_l[0], box_l[1]), box_l[2]);
245 auto const r_min = get_value_or<double>(parameters, "r_min", 0.);
246 auto const r_max = get_value_or<double>(parameters, "r_max", r_max_limit);
247 auto const r_bins = get_value_or<int>(parameters, "r_bins", 100);
248 auto const log_flag = get_value_or<bool>(parameters, "log_flag", false);
249 auto const int_flag = get_value_or<bool>(parameters, "int_flag", false);
250 context()->parallel_try_catch([=]() {
251 if (log_flag and r_min <= 0.) {
252 throw std::domain_error("Parameter 'r_min' must be > 0");
253 }
254 if (r_min < 0.) {
255 throw std::domain_error("Parameter 'r_min' must be >= 0");
256 }
257 if (r_min >= r_max) {
258 throw std::domain_error("Parameter 'r_max' must be > 'r_min'");
259 }
260 if (r_max > r_max_limit) {
261 throw std::domain_error("Parameter 'r_max' must be <= box_l / 2");
262 }
263 if (r_bins <= 0) {
264 throw std::domain_error("Parameter 'r_bins' must be >= 1");
265 }
266 });
267 auto const p_types1 =
268 get_value<std::vector<int>>(parameters, "type_list_a");
269 auto const p_types2 =
270 get_value<std::vector<int>>(parameters, "type_list_b");
271 for (auto const p_type : p_types1) {
272 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
273 }
274 for (auto const p_type : p_types2) {
275 context()->parallel_try_catch([&]() { check_particle_type(p_type); });
276 }
280 }
281 if (name == "calculate_energy") {
282 return m_obs_stat->do_call_method("calculate_energy", {});
283 }
284 if (name == "calculate_scalar_pressure") {
285 return m_obs_stat->do_call_method("calculate_scalar_pressure", {});
286 }
287 if (name == "calculate_pressure_tensor") {
288 return m_obs_stat->do_call_method("calculate_pressure_tensor", {});
289 }
290 return {};
291}
292
293} // namespace Analysis
294} // namespace ScriptInterface
Vector implementation and trait types for boost qvm interoperability.
std::vector< NeighborPIDs > get_neighbor_pids(System::System const &system)
Returns pairs of particle ids and neighbor particle id lists.
Definition cells.cpp:189
This file contains everything related to the global cell structure / cell system.
Variant do_call_method(std::string const &name, VariantMap const &parameters) override
Definition Analysis.cpp:114
virtual void parallel_try_catch(std::function< void()> const &cb) const =0
boost::string_ref name() const
Context * context() const
Responsible context.
boost::mpi::communicator comm_cart
The communicator.
This file contains the asynchronous MPI communication.
Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm)
Viscous stress tensor of the DPD interaction.
Definition dpd.cpp:174
Routines to use DPD as thermostat or pair force .
static void check_topology(CellStructure const &cell_structure, int chain_start, int chain_length, int n_chains)
Check if a contiguous range of particle ids exists.
Definition Analysis.cpp:53
T get_value(Variant const &v)
Extract value of specific type T from a Variant.
std::unordered_map< std::string, Variant > VariantMap
Definition Variant.hpp:69
auto make_unordered_map_of_variants(std::unordered_map< K, V > const &v)
Definition Variant.hpp:80
T mpi_reduce_sum(boost::mpi::communicator const &comm, T const &result)
Reduce object by sum on the head node.
auto make_vector_of_variants(std::vector< T > const &v)
Definition Variant.hpp:88
boost::make_recursive_variant< None, bool, int, std::size_t, double, std::string, ObjectRef, Utils::Vector3b, Utils::Vector3i, Utils::Vector2d, Utils::Vector3d, Utils::Vector4d, std::vector< int >, std::vector< double >, std::vector< boost::recursive_variant_ >, std::unordered_map< int, boost::recursive_variant_ >, std::unordered_map< std::string, boost::recursive_variant_ > >::type Variant
Possible types for parameters.
Definition Variant.hpp:67
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
Various procedures concerning interactions between particles.
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.
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::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.