ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
sd_interface.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 STOKESIAN_DYNAMICS
23#include "sd_interface.hpp"
24
25#include "stokesian_dynamics/sd_cpu.hpp"
26
27#include "BoxGeometry.hpp"
28#include "Particle.hpp"
29#include "communication.hpp"
30#include "system/System.hpp"
31#include "thermostat.hpp"
32
33#include <utils/Vector.hpp>
36
37#include <boost/serialization/is_bitwise_serializable.hpp>
38
39#include <algorithm>
40#include <cmath>
41#include <cstddef>
42#include <iterator>
43#include <stdexcept>
44#include <string>
45#include <utility>
46#include <vector>
47
48/* type for particle data transfer between nodes */
50 SD_particle_data() = default;
51 explicit SD_particle_data(Particle const &p)
52 : type(p.type()), pos(p.pos()), ext_force(p.force_and_torque()) {}
53
54 int type = 0;
55
56 /* particle position */
57 Utils::Vector3d pos = {0., 0., 0.};
58
59 /* external force */
61
62 template <class Archive> void serialize(Archive &ar, long int /* version */) {
63 ar & type;
64 ar & pos;
65 ar & ext_force;
66 }
67};
68
69BOOST_IS_BITWISE_SERIALIZABLE(SD_particle_data)
70
72
73/** Buffer that holds the (translational and angular) velocities of the local
74 * particles on each node, used for returning results. */
75static std::vector<double> v_sd{};
76
78 auto const &box_geo = *System::get_system().box_geo;
79 if (box_geo.periodic(0) or box_geo.periodic(1) or box_geo.periodic(2)) {
80 throw std::runtime_error(
81 "Stokesian Dynamics requires periodicity (False, False, False)");
82 }
83 ::params = obj;
84}
85
86/** Update translational and rotational velocities of all particles. */
87template <typename ParticleIterable>
88void sd_update_locally(ParticleIterable const &parts) {
89 std::size_t i = 0;
90
91 // Even though on the head node, the v_sd vector is larger than
92 // the (local) parts vector, this should still work. Because the local
93 // particles correspond to the first 6*n entries in the head node's v_sd
94 // (which holds the velocities of ALL particles).
95
96 for (auto &p : parts) {
97 // Copy velocities
98 p.v()[0] = v_sd[6 * i + 0];
99 p.v()[1] = v_sd[6 * i + 1];
100 p.v()[2] = v_sd[6 * i + 2];
101
102 p.omega()[0] = v_sd[6 * i + 3];
103 p.omega()[1] = v_sd[6 * i + 4];
104 p.omega()[2] = v_sd[6 * i + 5];
105
106 ++i;
107 }
108}
109
111 double viscosity, std::unordered_map<int, double> radii, int flags)
112 : viscosity{viscosity}, radii{radii}, flags{flags} {
113 if (viscosity < 0.) {
114 throw std::domain_error("Viscosity has an invalid value: " +
115 std::to_string(viscosity));
116 }
117 /* Check that radii are positive */
118 for (auto const &kv : radii) {
119 if (kv.second < 0.) {
120 throw std::domain_error(
121 "Particle radius for type " + std::to_string(kv.first) +
122 " has an invalid value: " + std::to_string(kv.second));
123 }
124 }
125}
126
128 StokesianThermostat const &stokesian,
129 double const time_step, double const kT) {
130
131 static std::vector<SD_particle_data> parts_buffer{};
132
133 parts_buffer.clear();
134 std::transform(particles.begin(), particles.end(),
135 std::back_inserter(parts_buffer),
136 [](auto const &p) { return SD_particle_data(p); });
137 Utils::Mpi::gather_buffer(parts_buffer, ::comm_cart, 0);
138
139 /* Buffer that holds local particle data, and all particles on the head
140 * node used for sending particle data to head node. */
141 if (::comm_cart.rank() == 0) {
142 std::size_t n_part = parts_buffer.size();
143
144 static std::vector<double> x_host{};
145 static std::vector<double> f_host{};
146 static std::vector<double> a_host{};
147
148 x_host.resize(6 * n_part);
149 f_host.resize(6 * n_part);
150 a_host.resize(n_part);
151
152 std::size_t i = 0;
153 for (auto const &p : parts_buffer) {
154 x_host[6 * i + 0] = p.pos[0];
155 x_host[6 * i + 1] = p.pos[1];
156 x_host[6 * i + 2] = p.pos[2];
157 // Actual orientation is not needed, just need default.
158 x_host[6 * i + 3] = 1;
159 x_host[6 * i + 4] = 0;
160 x_host[6 * i + 5] = 0;
161
162 f_host[6 * i + 0] = p.ext_force.f[0];
163 f_host[6 * i + 1] = p.ext_force.f[1];
164 f_host[6 * i + 2] = p.ext_force.f[2];
165
166 f_host[6 * i + 3] = p.ext_force.torque[0];
167 f_host[6 * i + 4] = p.ext_force.torque[1];
168 f_host[6 * i + 5] = p.ext_force.torque[2];
169
170 a_host[i] = params.radii.at(p.type);
171
172 ++i;
173 }
174
175 v_sd = sd_cpu(x_host, f_host, a_host, n_part, params.viscosity,
176 std::sqrt(kT / time_step),
177 static_cast<std::size_t>(stokesian.rng_counter()),
178 static_cast<std::size_t>(stokesian.rng_seed()), params.flags);
179 } else { // if (this_node == 0)
180 v_sd.resize(particles.size() * 6);
181 } // if (this_node == 0) {...} else
182
184 v_sd.data(), static_cast<int>(particles.size() * 6), ::comm_cart, 0);
185 sd_update_locally(particles);
186}
187
188#endif // STOKESIAN_DYNAMICS
Vector implementation and trait types for boost qvm interoperability.
base_type::size_type size() const
std::shared_ptr< BoxGeometry > box_geo
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
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.
void scatter_buffer(T *buffer, int n_elem, boost::mpi::communicator comm, int root=0)
Scatter buffer with different size on each node.
static StokesianDynamicsParameters params
void propagate_vel_pos_sd(ParticleRangeStokesian const &particles, StokesianThermostat const &stokesian, double const time_step, double const kT)
Takes the forces and torques on all particles and computes their velocities.
void register_integrator(StokesianDynamicsParameters const &obj)
void sd_update_locally(ParticleIterable const &parts)
Update translational and rotational velocities of all particles.
static std::vector< double > v_sd
Buffer that holds the (translational and angular) velocities of the local particles on each node,...
See for the Stokesian dynamics method used here.
uint64_t rng_counter() const
Get current value of the RNG.
uint32_t rng_seed() const
Force information on a particle.
Definition Particle.hpp:290
Struct holding all information for one particle.
Definition Particle.hpp:395
SD_particle_data()=default
SD_particle_data(Particle const &p)
Utils::Vector3d pos
void serialize(Archive &ar, long int)
ParticleForce ext_force
std::unordered_map< int, double > radii
StokesianDynamicsParameters(double viscosity, std::unordered_map< int, double > radii, int flags)
Thermostat for Stokesian dynamics.