22#ifdef STOKESIAN_DYNAMICS
25#include "stokesian_dynamics/sd_cpu.hpp"
29#include "communication.hpp"
30#include "system/System.hpp"
31#include "thermostat.hpp"
37#include <boost/serialization/is_bitwise_serializable.hpp>
62 template <
class Archive>
void serialize(Archive &ar,
long int ) {
75static std::vector<double>
v_sd{};
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)");
87template <
typename ParticleIterable>
96 for (
auto &p : parts) {
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];
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];
111 double viscosity, std::unordered_map<int, double> radii,
int flags)
112 : viscosity{viscosity}, radii{radii}, flags{flags} {
114 throw std::domain_error(
"Viscosity has an invalid value: " +
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));
129 double const time_step,
double const kT) {
131 static std::vector<SD_particle_data> parts_buffer{};
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); });
142 std::size_t n_part = parts_buffer.size();
144 static std::vector<double> x_host{};
145 static std::vector<double> f_host{};
146 static std::vector<double> a_host{};
148 x_host.resize(6 * n_part);
149 f_host.resize(6 * n_part);
150 a_host.resize(n_part);
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];
158 x_host[6 * i + 3] = 1;
159 x_host[6 * i + 4] = 0;
160 x_host[6 * i + 5] = 0;
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];
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];
176 std::sqrt(kT / time_step),
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.
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.
Struct holding all information for one particle.
SD_particle_data()=default
SD_particle_data(Particle const &p)
void serialize(Archive &ar, long int)
std::unordered_map< int, double > radii
StokesianDynamicsParameters(double viscosity, std::unordered_map< int, double > radii, int flags)
Thermostat for Stokesian dynamics.