35#include "system/System.hpp"
43#include <boost/mpi/collectives/reduce.hpp>
58 return Random::noise_uniform<RNGSalt::SALT_DPD>(
60 (pid1 < pid2) ? pid1 : pid2);
65 auto const max_type = nonbonded_ias.get_max_seen_particle_type();
66 for (
int type_a = 0; type_a <= max_type; type_a++) {
67 for (
int type_b = type_a; type_b <= max_type; type_b++) {
68 auto &ia_params = nonbonded_ias.get_ia_param(type_a, type_b);
70 ia_params.dpd.radial.pref =
71 sqrt(24.0 * kT * ia_params.dpd.radial.gamma / time_step);
72 ia_params.dpd.trans.pref =
73 sqrt(24.0 * kT * ia_params.dpd.trans.gamma / time_step);
78static double weight(
int type,
double r_cut,
double k,
double r) {
82 return 1. - pow((r / r_cut), k);
88 if (dist <
params.cutoff) {
93 auto const f_r =
params.pref * omega * noise;
111 auto const noise_vec =
120 auto const P = tensor_product(d / dist2, d);
123 auto const force =
P * (f_r - f_t) + f_t;
129 auto const &box_geo = *system.box_geo;
130 auto const &nonbonded_ias = *system.nonbonded_ias;
131 auto &cell_structure = *system.cell_structure;
132 system.on_observable_calc();
135 cell_structure.non_bonded_loop([&stress, &box_geo, &nonbonded_ias](
139 box_geo.velocity_difference(p1.
pos(), p2.
pos(), p1.
v(), p2.
v());
141 auto const &ia_params = nonbonded_ias.get_ia_param(p1.
type(), p2.
type());
142 auto const dist = std::sqrt(d.
dist2);
144 auto const f_r =
dpd_pair_force(ia_params.dpd.radial, v21, dist, {});
145 auto const f_t =
dpd_pair_force(ia_params.dpd.trans, v21, dist, {});
151 auto const f =
P * (f_r - f_t) + f_t;
153 stress += tensor_product(d.
vec21, f);
177 std::remove_const_t<
decltype(local_stress)> global_stress;
179 boost::mpi::reduce(comm, local_stress, global_stress, std::plus<>(), 0);
Vector implementation and trait types for boost qvm interoperability.
This file contains everything related to the global cell structure / cell system.
Utils::Vector3d velocity_difference(Utils::Vector3d const &x, Utils::Vector3d const &y, Utils::Vector3d const &u, Utils::Vector3d const &v) const
Calculate the velocity difference including the Lees-Edwards velocity.
std::shared_ptr< BoxGeometry > box_geo
std::shared_ptr< InteractionsNonBonded > nonbonded_ias
This file contains the defaults for ESPResSo.
Utils::Vector3d dpd_pair_force(DPDParameters const ¶ms, Utils::Vector3d const &v, double dist, Utils::Vector3d const &noise)
static auto dpd_viscous_stress_local()
void dpd_init(double kT, double time_step)
static double weight(int type, double r_cut, double k, double r)
Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm)
Viscous stress tensor of the DPD interaction.
Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2)
Return a random uniform 3D vector with the Philox thermostat.
Routines to use DPD as thermostat or pair force .
Matrix implementation and trait types for boost qvm interoperability.
void flatten(Range const &v, OutputIterator out)
Flatten a range of ranges.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Various procedures concerning interactions between particles.
Random number generation using Philox.
static SteepestDescentParameters params
Currently active steepest descent instance.
uint64_t rng_counter() const
Get current value of the RNG.
uint32_t rng_seed() const
Thermostat for dissipative particle dynamics.
Distance vector and length handed to pair kernels.
Parameters for non-bonded interactions.
Struct holding all information for one particle.
auto const & type() const
double gamma
Dampening constant.
Matrix representation with static size.