ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
CellSystem.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 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 "CellSystem.hpp"
21
25
26#include "core/BoxGeometry.hpp"
31#include "core/cells.hpp"
36#include "core/tuning.hpp"
37
38#include <utils/Vector.hpp>
40
41#include <boost/mpi/collectives/gather.hpp>
42#include <boost/variant.hpp>
43
44#include <algorithm>
45#include <cassert>
46#include <iterator>
47#include <optional>
48#include <set>
49#include <sstream>
50#include <stdexcept>
51#include <string>
52#include <unordered_map>
53#include <utility>
54#include <vector>
55
56static int coord(std::string const &s) {
57 if (s == "x")
58 return 0;
59 if (s == "y")
60 return 1;
61 if (s == "z")
62 return 2;
63 throw std::invalid_argument("Invalid Cartesian coordinate: '" + s + "'");
64}
65
66static std::string coord_letter(int c) {
67 if (c == 0)
68 return "x";
69 if (c == 1)
70 return "y";
71 assert(c == 2);
72 return "z";
73}
74
75namespace ScriptInterface {
76namespace CellSystem {
77
80 {"use_verlet_lists",
81 [this](Variant const &v) {
82 get_cell_structure().use_verlet_list = get_value<bool>(v);
83 },
84 [this]() { return get_cell_structure().use_verlet_list; }},
85 {"node_grid",
86 [this](Variant const &v) {
87 context()->parallel_try_catch([this, &v]() {
88 auto const error_msg = std::string("Parameter 'node_grid'");
93 if (n_nodes_new != n_nodes_old) {
94 std::stringstream reason;
95 reason << ": MPI world size " << n_nodes_old << " incompatible "
96 << "with new node grid [" << new_node_grid << "]";
97 throw std::invalid_argument(error_msg + reason.str());
98 }
99 try {
101 get_system().on_node_grid_change();
102 } catch (...) {
104 get_system().on_node_grid_change();
105 throw;
106 }
107 });
108 },
109 []() { return ::communicator.node_grid; }},
110 {"skin",
111 [this](Variant const &v) {
112 auto const new_skin = get_value<double>(v);
113 if (new_skin < 0.) {
114 if (context()->is_head_node()) {
115 throw std::domain_error("Parameter 'skin' must be >= 0");
116 }
117 throw Exception("");
118 }
119 get_cell_structure().set_verlet_skin(new_skin);
120 },
121 [this]() { return get_cell_structure().get_verlet_skin(); }},
122 {"decomposition_type", AutoParameter::read_only,
123 [this]() {
124 return cs_type_to_name.at(get_cell_structure().decomposition_type());
125 }},
126 {"n_square_types", AutoParameter::read_only,
127 [this]() {
128 if (get_cell_structure().decomposition_type() !=
130 return Variant{none};
131 }
132 auto const hd = get_hybrid_decomposition();
133 auto const ns_types = hd.get_n_square_types();
134 return Variant{std::vector<int>(ns_types.begin(), ns_types.end())};
135 }},
136 {"fully_connected_boundary", AutoParameter::read_only,
137 [this]() {
138 if (get_cell_structure().decomposition_type() !=
140 return Variant{none};
141 }
142 auto const rd = get_regular_decomposition();
143 auto const fcb = rd.fully_connected_boundary();
144 if (not fcb)
145 return Variant{none};
146 return Variant{std::unordered_map<std::string, Variant>{
147 {{"boundary", Variant{coord_letter((*fcb).first)}},
148 {"direction", Variant{coord_letter((*fcb).second)}}}}};
149 }},
150 {"cutoff_regular", AutoParameter::read_only,
151 [this]() {
152 if (get_cell_structure().decomposition_type() !=
154 return Variant{none};
155 }
156 auto const hd = get_hybrid_decomposition();
157 return Variant{hd.get_cutoff_regular()};
158 }},
159 {"max_cut_nonbonded", AutoParameter::read_only,
160 [this]() { return get_system().nonbonded_ias->maximal_cutoff(); }},
161 {"max_cut_bonded", AutoParameter::read_only,
162 [this]() { return get_system().bonded_ias->maximal_cutoff(); }},
163 {"interaction_range", AutoParameter::read_only,
164 [this]() { return get_system().get_interaction_range(); }},
165 });
166}
167
168Variant CellSystem::do_call_method(std::string const &name,
169 VariantMap const &params) {
170 if (name == "initialize") {
171 auto const cs_name = get_value<std::string>(params, "name");
172 auto const cs_type = cs_name_to_type.at(cs_name);
173 initialize(cs_type, params);
174 return {};
175 }
176 if (name == "resort") {
177 auto const global_flag = get_value_or<bool>(params, "global_flag", true);
178 return mpi_resort_particles(global_flag);
179 }
180 if (name == "get_state") {
181 auto state = get_parameters();
182 auto const cs_type = get_cell_structure().decomposition_type();
184 auto const rd = get_regular_decomposition();
185 state["cell_grid"] = Variant{rd.cell_grid};
186 state["cell_size"] = Variant{rd.cell_size};
187 } else if (cs_type == CellStructureType::HYBRID) {
188 auto const hd = get_hybrid_decomposition();
189 state["cell_grid"] = Variant{hd.get_cell_grid()};
190 state["cell_size"] = Variant{hd.get_cell_size()};
191 mpi_resort_particles(true); // needed to get correct particle counts
192 state["parts_per_decomposition"] =
193 Variant{std::unordered_map<std::string, Variant>{
194 {"regular", hd.count_particles_in_regular()},
195 {"n_square", hd.count_particles_in_n_square()}}};
196 }
197 state["verlet_reuse"] = get_cell_structure().get_verlet_reuse();
198 state["n_nodes"] = context()->get_comm().size();
199 return state;
200 }
201 if (name == "get_pairs") {
202 std::vector<Variant> out;
203 context()->parallel_try_catch([this, &params, &out]() {
204 auto &system = get_system();
205 system.on_observable_calc();
206 std::vector<std::pair<int, int>> pair_list;
207 auto const distance = get_value<double>(params, "distance");
208 if (boost::get<std::string>(&params.at("types")) != nullptr) {
209 auto const key = get_value<std::string>(params, "types");
210 if (key != "all") {
211 throw std::invalid_argument("Unknown argument types='" + key + "'");
212 }
213 pair_list = get_pairs(system, distance);
214 } else {
215 auto const types = get_value<std::vector<int>>(params, "types");
217 }
219 std::transform(pair_list.begin(), pair_list.end(),
220 std::back_inserter(out),
221 [](std::pair<int, int> const &pair) {
222 return std::vector<int>{pair.first, pair.second};
223 });
224 });
225 return out;
226 }
227 if (name == "get_neighbors") {
228 std::vector<std::vector<int>> neighbors_global;
229 context()->parallel_try_catch([this, &neighbors_global, &params]() {
230 auto &system = get_system();
231 system.on_observable_calc();
232 auto const dist = get_value<double>(params, "distance");
233 auto const pid = get_value<int>(params, "pid");
234 auto const ret = get_short_range_neighbors(system, pid, dist);
235 std::vector<int> neighbors_local;
236 if (ret) {
238 }
239 boost::mpi::gather(context()->get_comm(), neighbors_local,
241 });
242 std::vector<int> neighbors;
243 for (auto const &neighbors_local : neighbors_global) {
244 if (not neighbors_local.empty()) {
245 neighbors = neighbors_local;
246 break;
247 }
248 }
249 return neighbors;
250 }
251 if (name == "non_bonded_loop_trace") {
252 auto &system = get_system();
253 system.on_observable_calc();
254 std::vector<Variant> out;
255 auto pair_list =
256 non_bonded_loop_trace(system, context()->get_comm().rank());
257 Utils::Mpi::gather_buffer(pair_list, context()->get_comm());
258 std::transform(pair_list.begin(), pair_list.end(), std::back_inserter(out),
259 [](PairInfo const &pair) {
260 return std::vector<Variant>{pair.id1, pair.id2,
261 pair.pos1, pair.pos2,
262 pair.vec21, pair.node};
263 });
264 return out;
265 }
266 if (name == "tune_skin") {
267 auto &system = get_system();
268 system.tune_verlet_skin(
269 get_value<double>(params, "min_skin"),
270 get_value<double>(params, "max_skin"), get_value<double>(params, "tol"),
271 get_value<int>(params, "int_steps"),
272 get_value_or<bool>(params, "adjust_max_skin", false));
273 return get_cell_structure().get_verlet_skin();
274 }
275 if (name == "get_max_range") {
276 return get_cell_structure().max_range();
277 }
278 return {};
279}
280
281std::vector<int> CellSystem::mpi_resort_particles(bool global_flag) const {
282 auto &cell_structure = get_cell_structure();
283 cell_structure.resort_particles(global_flag);
285 auto const size = static_cast<int>(cell_structure.local_particles().size());
286 std::vector<int> n_part_per_node;
287 boost::mpi::gather(context()->get_comm(), size, n_part_per_node, 0);
288 return n_part_per_node;
289}
290
291void CellSystem::initialize(CellStructureType const &cs_type,
292 VariantMap const &params) {
293 auto const verlet = get_value_or<bool>(params, "use_verlet_lists", true);
294 auto &system = get_system();
295 m_cell_structure->use_verlet_list = verlet;
296 if (cs_type == CellStructureType::HYBRID) {
297 auto const cutoff_regular = get_value<double>(params, "cutoff_regular");
298 auto const ns_types =
299 get_value_or<std::vector<int>>(params, "n_square_types", {});
300 auto n_square_types = std::set<int>{ns_types.begin(), ns_types.end()};
301 m_cell_structure->set_hybrid_decomposition(cutoff_regular, n_square_types);
302 } else if (cs_type == CellStructureType::REGULAR) {
303 std::optional<std::pair<int, int>> fcb_pair = std::nullopt;
304 if (params.contains("fully_connected_boundary") and
305 not is_none(params.at("fully_connected_boundary"))) {
306 auto const variant =
307 get_value<VariantMap>(params, "fully_connected_boundary");
308 context()->parallel_try_catch([&fcb_pair, &variant]() {
309 fcb_pair = {{coord(boost::get<std::string>(variant.at("boundary"))),
310 coord(boost::get<std::string>(variant.at("direction")))}};
311 });
312 }
313 context()->parallel_try_catch([this, &fcb_pair]() {
314 m_cell_structure->set_regular_decomposition(
315 get_system().get_interaction_range(), fcb_pair);
316 });
317 } else {
318 system.set_cell_structure_topology(cs_type);
319 }
320}
321
322void CellSystem::configure(Particles::ParticleHandle &particle) {
323 particle.attach(m_system);
324}
325
326void CellSystem::configure(Particles::ParticleSlice &slice) {
327 slice.attach(m_system);
328}
329
330} // namespace CellSystem
331} // namespace ScriptInterface
CellStructureType
Cell structure topology.
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
static int coord(std::string const &s)
static std::string coord_letter(int c)
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
std::optional< std::vector< int > > get_short_range_neighbors(System::System const &system, int const pid, double const distance)
Get ids of particles that are within a certain distance of another particle.
Definition cells.cpp:119
std::vector< PairInfo > non_bonded_loop_trace(System::System const &system, int const rank)
Returns pairs of particle ids, positions and distance as seen by the non-bonded loop.
Definition cells.cpp:177
std::vector< std::pair< int, int > > get_pairs_of_types(System::System const &system, double const distance, std::vector< int > const &types)
Get pairs closer than distance if both their types are in types.
Definition cells.cpp:167
std::vector< std::pair< int, int > > get_pairs(System::System const &system, double const distance)
Get pairs closer than distance from the cells.
Definition cells.cpp:159
This file contains everything related to the global cell structure / cell system.
void add_parameters(std::vector< AutoParameter > &&params)
Variant do_call_method(std::string const &name, VariantMap const &params) override
virtual void parallel_try_catch(std::function< void()> const &cb) const =0
virtual boost::mpi::communicator const & get_comm() const =0
VariantMap get_parameters() const
Get current parameters.
boost::string_ref name() const
Context * context() const
Responsible context.
void attach(std::weak_ptr<::System::System > system)
void attach(std::weak_ptr<::System::System > system)
Communicator communicator
This file contains the asynchronous MPI communication.
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
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
bool is_none(Variant const &v)
Definition Variant.hpp:115
constexpr const None none
None-"literal".
Definition Variant.hpp:50
System & get_system()
void gather_buffer(std::vector< T, Allocator > &buffer, boost::mpi::communicator const &comm, int root=0)
Gather buffer with different size on each node.
T product(Vector< T, N > const &v)
Definition Vector.hpp:374
Various procedures concerning interactions between particles.
static auto & get_cell_structure()
void clear_particle_node()
Invalidate particle_node.
Particles creation and deletion.
static SteepestDescentParameters params
Currently active steepest descent instance.
Utils::Vector3i node_grid
void set_node_grid(Utils::Vector3i const &value)
Set new Cartesian topology.
static constexpr const ReadOnly read_only