ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
dpd.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21/** \file
22 * Implementation of dpd.hpp.
23 */
24#include "config/config.hpp"
25
26#ifdef DPD
27
28#include "dpd.hpp"
29
30#include "BoxGeometry.hpp"
32#include "cells.hpp"
34#include "random.hpp"
35#include "system/System.hpp"
36#include "thermostat.hpp"
37
38#include <utils/Vector.hpp>
39#include <utils/math/sqr.hpp>
41#include <utils/matrix.hpp>
42
43#include <boost/mpi/collectives/reduce.hpp>
44
45#include <algorithm>
46#include <cmath>
47#include <cstdint>
48#include <functional>
49
50/** Return a random uniform 3D vector with the Philox thermostat.
51 * Random numbers depend on
52 * 1. dpd_rng_counter (initialized by seed) which is increased on integration
53 * 2. Salt (decorrelates different counters)
54 * 3. Two particle IDs (order-independent, decorrelates particles, gets rid of
55 * seed-per-node)
56 */
57Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2) {
58 return Random::noise_uniform<RNGSalt::SALT_DPD>(
59 dpd.rng_counter(), dpd.rng_seed(), (pid1 < pid2) ? pid2 : pid1,
60 (pid1 < pid2) ? pid1 : pid2);
61}
62
63void dpd_init(double kT, double time_step) {
64 auto &nonbonded_ias = *System::get_system().nonbonded_ias;
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);
69
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);
74 }
75 }
76}
77
78static double weight(int type, double r_cut, double k, double r) {
79 if (type == 0) {
80 return 1.;
81 }
82 return 1. - pow((r / r_cut), k);
83}
84
86 Utils::Vector3d const &v, double dist,
87 Utils::Vector3d const &noise) {
88 if (dist < params.cutoff) {
89 auto const omega = weight(params.wf, params.cutoff, params.k, dist);
90 auto const omega2 = Utils::sqr(omega);
91
92 auto const f_d = params.gamma * omega2 * v;
93 auto const f_r = params.pref * omega * noise;
94
95 return f_r - f_d;
96 }
97
98 return {};
99}
100
102dpd_pair_force(Particle const &p1, Particle const &p2, DPDThermostat const &dpd,
103 BoxGeometry const &box_geo, IA_parameters const &ia_params,
104 Utils::Vector3d const &d, double dist, double dist2) {
105 if (ia_params.dpd.radial.cutoff <= 0.0 && ia_params.dpd.trans.cutoff <= 0.0) {
106 return {};
107 }
108
109 auto const v21 =
110 box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v());
111 auto const noise_vec =
112 (ia_params.dpd.radial.pref > 0.0 || ia_params.dpd.trans.pref > 0.0)
113 ? dpd_noise(dpd, p1.id(), p2.id())
114 : Utils::Vector3d{};
115
116 auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, noise_vec);
117 auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, noise_vec);
118
119 /* Projection operator to radial direction */
120 auto const P = tensor_product(d / dist2, d);
121 /* This is equivalent to P * f_r + (1 - P) * f_t, but with
122 * doing only one matrix-vector multiplication */
123 auto const force = P * (f_r - f_t) + f_t;
124 return force;
125}
126
128 auto &system = System::get_system();
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();
133
135 cell_structure.non_bonded_loop([&stress, &box_geo, &nonbonded_ias](
136 Particle const &p1, Particle const &p2,
137 Distance const &d) {
138 auto const v21 =
139 box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v());
140
141 auto const &ia_params = nonbonded_ias.get_ia_param(p1.type(), p2.type());
142 auto const dist = std::sqrt(d.dist2);
143
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, {});
146
147 /* Projection operator to radial direction */
148 auto const P = tensor_product(d.vec21 / d.dist2, d.vec21);
149 /* This is equivalent to P * f_r + (1 - P) * f_t, but with
150 * doing only one matrix-vector multiplication */
151 auto const f = P * (f_r - f_t) + f_t;
152
153 stress += tensor_product(d.vec21, f);
154 });
155
156 return stress;
157}
158
159/**
160 * @brief Viscous stress tensor of the DPD interaction.
161 *
162 * This calculates the total viscous stress contribution of the
163 * DPD interaction. It contains only the dissipative contributions
164 * of the interaction without noise. It's calculated as the
165 * sum over all pair virials as
166 * \f[
167 * \sigma^{\nu\mu} = V^{-1}\sum_i \sum_{j < i} r_{i,j}^{\nu} (- \gamma_{i,j}
168 * v_{i,j})^{\mu} \f] where \f$\gamma_{i,j}\f$ is the (in general tensor valued)
169 * DPD friction coefficient for particles i and j, \f$v_{i,j}\f$, \f$r_{i,j}\f$
170 * are their relative velocity and distance and \f$V\f$ is the box volume.
171 *
172 * @return Stress tensor contribution.
173 */
174Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm) {
175 auto const &box_geo = *System::get_system().box_geo;
176 auto const local_stress = dpd_viscous_stress_local();
177 std::remove_const_t<decltype(local_stress)> global_stress;
178
179 boost::mpi::reduce(comm, local_stress, global_stress, std::plus<>(), 0);
180
181 return Utils::flatten(global_stress) / box_geo.volume();
182}
183
184#endif // DPD
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 &params, Utils::Vector3d const &v, double dist, Utils::Vector3d const &noise)
Definition dpd.cpp:85
static auto dpd_viscous_stress_local()
Definition dpd.cpp:127
void dpd_init(double kT, double time_step)
Definition dpd.cpp:63
static double weight(int type, double r_cut, double k, double r)
Definition dpd.cpp:78
Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm)
Viscous stress tensor of the DPD interaction.
Definition dpd.cpp:174
Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2)
Return a random uniform 3D vector with the Philox thermostat.
Definition dpd.cpp:57
Routines to use DPD as thermostat or pair force .
Matrix implementation and trait types for boost qvm interoperability.
System & get_system()
void flatten(Range const &v, OutputIterator out)
Flatten a range of ranges.
Definition flatten.hpp:64
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
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.
Utils::Vector3d vec21
Parameters for non-bonded interactions.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & v() const
Definition Particle.hpp:433
auto const & type() const
Definition Particle.hpp:418
auto const & pos() const
Definition Particle.hpp:431
auto const & id() const
Definition Particle.hpp:414
double gamma
Dampening constant.
Matrix representation with static size.
Definition matrix.hpp:65