ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
thermostats/brownian_inline.hpp
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
22#pragma once
23
24#include "config/config.hpp"
25
26#include "Particle.hpp"
27#include "random.hpp"
28#include "rotation.hpp"
29#include "thermostat.hpp"
30
31#include <utils/Vector.hpp>
32
33#include <cmath>
34
35/** Determine position: viscous drag driven by conservative forces.
36 * From eq. (14.39) in @cite schlick10a.
37 * @param[in] brownian_gamma Brownian translational gamma
38 * @param[in] p Particle
39 * @param[in] dt Time step
40 */
41inline Utils::Vector3d bd_drag(Thermostat::GammaType const &brownian_gamma,
42 Particle const &p, double dt) {
43 // The friction tensor Z from the Eq. (14.31) of schlick10a:
45
46#ifdef THERMOSTAT_PER_PARTICLE
47 if (p.gamma() >= Thermostat::GammaType{}) {
48 gamma = p.gamma();
49 } else
50#endif
51 {
52 gamma = brownian_gamma;
53 }
54
55#ifdef PARTICLE_ANISOTROPY
56 // Particle frictional isotropy check.
57 auto const aniso_flag = (gamma[0] != gamma[1]) || (gamma[1] != gamma[2]);
58 Utils::Vector3d delta_pos_lab;
59 if (aniso_flag) {
60 auto const force_body = convert_vector_space_to_body(p, p.force());
61 auto const delta_pos_body = hadamard_division(force_body * dt, gamma);
62 delta_pos_lab = convert_vector_body_to_space(p, delta_pos_body);
63 }
64#endif
65
66 Utils::Vector3d position = {};
67 for (unsigned int j = 0; j < 3; j++) {
68 // Second (deterministic) term of the Eq. (14.39) of schlick10a.
69 // Only a conservative part of the force is used here
70#ifdef PARTICLE_ANISOTROPY
71 if (aniso_flag) {
72 if (!p.is_fixed_along(j)) {
73 position[j] = delta_pos_lab[j];
74 }
75 } else {
76 if (!p.is_fixed_along(j)) {
77 position[j] = p.force()[j] * dt / gamma[j];
78 }
79 }
80#else
81 if (!p.is_fixed_along(j)) {
82 position[j] = p.force()[j] * dt / gamma;
83 }
84#endif // PARTICLE_ANISOTROPY
85 }
86 return position;
87}
88
89/** Set the terminal velocity driven by the conservative forces drag.
90 * From eq. (14.34) in @cite schlick10a.
91 * @param[in] brownian_gamma Brownian translational gamma
92 * @param[in] p Particle
93 */
95 Particle const &p) {
96 // The friction tensor Z from the eq. (14.31) of schlick10a:
98
99#ifdef THERMOSTAT_PER_PARTICLE
100 if (p.gamma() >= Thermostat::GammaType{}) {
101 gamma = p.gamma();
102 } else
103#endif
104 {
105 gamma = brownian_gamma;
106 }
107
108#ifdef PARTICLE_ANISOTROPY
109 // Particle frictional isotropy check.
110 auto const aniso_flag = (gamma[0] != gamma[1]) || (gamma[1] != gamma[2]);
111 Utils::Vector3d vel_lab;
112 if (aniso_flag) {
113 auto const force_body = convert_vector_space_to_body(p, p.force());
114 auto const vel_body = hadamard_division(force_body, gamma);
115 vel_lab = convert_vector_body_to_space(p, vel_body);
116 }
117#endif
118
120 for (unsigned int j = 0; j < 3; j++) {
121 // First (deterministic) term of the eq. (14.34) of schlick10a taking
122 // into account eq. (14.35). Only conservative part of the force is used
123 // here.
124#ifdef PARTICLE_ANISOTROPY
125 if (aniso_flag) {
126 if (!p.is_fixed_along(j)) {
127 velocity[j] = vel_lab[j];
128 }
129 } else {
130 if (!p.is_fixed_along(j)) {
131 velocity[j] = p.force()[j] / gamma[j];
132 }
133 }
134#else // PARTICLE_ANISOTROPY
135 if (!p.is_fixed_along(j)) {
136 velocity[j] = p.force()[j] / gamma;
137 }
138#endif // PARTICLE_ANISOTROPY
139 }
140 return velocity;
141}
142
143/** Determine the positions: random walk part.
144 * From eq. (14.37) in @cite schlick10a.
145 * @param[in] brownian Parameters
146 * @param[in] p Particle
147 * @param[in] dt Time step
148 * @param[in] kT Thermal energy
149 */
151 Particle const &p, double dt, double kT) {
152 Thermostat::GammaType sigma_pos = brownian.sigma_pos;
153#ifdef THERMOSTAT_PER_PARTICLE
154 // override default if particle-specific gamma
155 if (p.gamma() >= Thermostat::GammaType{}) {
156 if (kT > 0.0) {
157 sigma_pos = BrownianThermostat::sigma(kT, p.gamma());
158 } else {
159 sigma_pos = Thermostat::GammaType{};
160 }
161 }
162#endif // THERMOSTAT_PER_PARTICLE
163
164 // Eq. (14.37) is factored by the Gaussian noise (12.22) with its squared
165 // magnitude defined in the second eq. (14.38), schlick10a.
166 Utils::Vector3d delta_pos_body{};
167 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_WALK>(
168 brownian.rng_counter(), brownian.rng_seed(), p.id());
169 for (unsigned int j = 0; j < 3; j++) {
170 if (!p.is_fixed_along(j)) {
171#ifndef PARTICLE_ANISOTROPY
172 if (sigma_pos > 0.0) {
173 delta_pos_body[j] = sigma_pos * sqrt(dt) * noise[j];
174 } else {
175 delta_pos_body[j] = 0.0;
176 }
177#else
178 if (sigma_pos[j] > 0.0) {
179 delta_pos_body[j] = sigma_pos[j] * sqrt(dt) * noise[j];
180 } else {
181 delta_pos_body[j] = 0.0;
182 }
183#endif // PARTICLE_ANISOTROPY
184 }
185 }
186
187#ifdef PARTICLE_ANISOTROPY
188 // Particle frictional isotropy check.
189 auto const aniso_flag =
190 (sigma_pos[0] != sigma_pos[1]) || (sigma_pos[1] != sigma_pos[2]);
191 Utils::Vector3d delta_pos_lab;
192 if (aniso_flag) {
193 delta_pos_lab = convert_vector_body_to_space(p, delta_pos_body);
194 }
195#endif
196
197 Utils::Vector3d position = {};
198 for (unsigned int j = 0; j < 3; j++) {
199 if (!p.is_fixed_along(j)) {
200#ifdef PARTICLE_ANISOTROPY
201 position[j] += aniso_flag ? delta_pos_lab[j] : delta_pos_body[j];
202#else
203 position[j] += delta_pos_body[j];
204#endif
205 }
206 }
207 return position;
208}
209
210/** Determine the velocities: random walk part.
211 * From eq. (10.2.16) in @cite pottier10a.
212 * @param[in] brownian Parameters
213 * @param[in] p Particle
214 */
216 Particle const &p) {
217 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_INC>(
218 brownian.rng_counter(), brownian.rng_seed(), p.id());
220 for (unsigned int j = 0; j < 3; j++) {
221 if (!p.is_fixed_along(j)) {
222 // Random (heat) velocity. See eq. (10.2.16) taking into account eq.
223 // (10.2.18) and (10.2.29), pottier10a. Note, that the pottier10a units
224 // system (see Eq. (10.1.1) there) has been adapted towards the ESPResSo
225 // and the referenced above schlick10a one, which is defined by the eq.
226 // (14.31) of schlick10a. A difference is the mass factor to the friction
227 // tensor. The noise is Gaussian according to the convention at p. 237
228 // (last paragraph), pottier10a.
229 velocity[j] += brownian.sigma_vel * noise[j] / sqrt(p.mass());
230 }
231 }
232 return velocity;
233}
234
235#ifdef ROTATION
236
237/** Determine quaternions: viscous drag driven by conservative torques.
238 * An analogy of eq. (14.39) in @cite schlick10a.
239 * @param[in] brownian_gamma_rotation Brownian rotational gamma
240 * @param[in] p Particle
241 * @param[in] dt Time step
242 */
244bd_drag_rot(Thermostat::GammaType const &brownian_gamma_rotation, Particle &p,
245 double dt) {
247
248#ifdef THERMOSTAT_PER_PARTICLE
249 if (p.gamma_rot() >= Thermostat::GammaType{}) {
250 gamma = p.gamma_rot();
251 } else
252#endif
253 {
254 gamma = brownian_gamma_rotation;
255 }
256
257 Utils::Vector3d dphi = {};
258 for (unsigned int j = 0; j < 3; j++) {
259 if (p.can_rotate_around(j)) {
260 // only a conservative part of the torque is used here
261#ifndef PARTICLE_ANISOTROPY
262 dphi[j] = p.torque()[j] * dt / gamma;
263#else
264 dphi[j] = p.torque()[j] * dt / gamma[j];
265#endif // PARTICLE_ANISOTROPY
266 }
267 }
268 dphi = mask(p.rotation(), dphi);
269 double dphi_m = dphi.norm();
270 if (dphi_m != 0.) {
271 auto const dphi_u = dphi / dphi_m;
272 return local_rotate_particle_body(p, dphi_u, dphi_m);
273 }
274 return p.quat();
275}
276
277/** Set the terminal angular velocity driven by the conservative torques drag.
278 * An analogy of the 1st term of eq. (14.34) in @cite schlick10a.
279 * @param[in] brownian_gamma_rotation Brownian rotational gamma
280 * @param[in] p Particle
281 */
282inline Utils::Vector3d
283bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation,
284 Particle const &p) {
286
287#ifdef THERMOSTAT_PER_PARTICLE
288 if (p.gamma_rot() >= Thermostat::GammaType{}) {
289 gamma = p.gamma_rot();
290 } else
291#endif
292 {
293 gamma = brownian_gamma_rotation;
294 }
295
296 Utils::Vector3d omega = {};
297 for (unsigned int j = 0; j < 3; j++) {
298 if (p.can_rotate_around(j)) {
299#ifdef PARTICLE_ANISOTROPY
300 omega[j] = p.torque()[j] / gamma[j];
301#else
302 omega[j] = p.torque()[j] / gamma;
303#endif // PARTICLE_ANISOTROPY
304 }
305 }
306 return mask(p.rotation(), omega);
307}
308
309/** Determine the quaternions: random walk part.
310 * An analogy of eq. (14.37) in @cite schlick10a.
311 * @param[in] brownian Parameters
312 * @param[in] p Particle
313 * @param[in] dt Time step
314 * @param[in] kT Thermal energy
315 */
318 double dt, double kT) {
319
320 Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation;
321#ifdef THERMOSTAT_PER_PARTICLE
322 // override default if particle-specific gamma
323 if (p.gamma_rot() >= Thermostat::GammaType{}) {
324 if (kT > 0.) {
325 sigma_pos = BrownianThermostat::sigma(kT, p.gamma_rot());
326 } else {
327 sigma_pos = {}; // just an indication of the infinity
328 }
329 }
330#endif // THERMOSTAT_PER_PARTICLE
331
332 Utils::Vector3d dphi = {};
333 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_INC>(
334 brownian.rng_counter(), brownian.rng_seed(), p.id());
335 for (unsigned int j = 0; j < 3; j++) {
336 if (p.can_rotate_around(j)) {
337#ifndef PARTICLE_ANISOTROPY
338 if (sigma_pos > 0.0) {
339 dphi[j] = noise[j] * sigma_pos * sqrt(dt);
340 }
341#else
342 if (sigma_pos[j] > 0.0) {
343 dphi[j] = noise[j] * sigma_pos[j] * sqrt(dt);
344 }
345#endif // PARTICLE_ANISOTROPY
346 }
347 }
348 dphi = mask(p.rotation(), dphi);
349 // making the algorithm independent of the order of the rotations
350 double dphi_m = dphi.norm();
351 if (dphi_m != 0) {
352 auto const dphi_u = dphi / dphi_m;
353 return local_rotate_particle_body(p, dphi_u, dphi_m);
354 }
355 return p.quat();
356}
357
358/** Determine the angular velocities: random walk part.
359 * An analogy of eq. (10.2.16) in @cite pottier10a.
360 * @param[in] brownian Parameters
361 * @param[in] p Particle
362 */
363inline Utils::Vector3d
365 auto const sigma_vel = brownian.sigma_vel_rotation;
366
367 Utils::Vector3d domega{};
368 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_WALK>(
369 brownian.rng_counter(), brownian.rng_seed(), p.id());
370 for (unsigned int j = 0; j < 3; j++) {
371 if (p.can_rotate_around(j)) {
372 domega[j] = sigma_vel * noise[j] / sqrt(p.rinertia()[j]);
373 }
374 }
375 return mask(p.rotation(), domega);
376}
377#endif // ROTATION
Vector implementation and trait types for boost qvm interoperability.
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
Random number generation using Philox.
static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs)
Velocity of the virtual site.
Definition relative.cpp:64
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
Definition rotation.hpp:57
Utils::Vector3d convert_vector_space_to_body(const Particle &p, const Utils::Vector3d &v)
Definition rotation.hpp:61
Utils::Quaternion< double > local_rotate_particle_body(Particle const &p, const Utils::Vector3d &axis_body_frame, const double phi)
Rotate the particle p around the body-frame defined NORMALIZED axis aBodyFrame by amount phi.
Definition rotation.hpp:119
uint64_t rng_counter() const
Get current value of the RNG.
uint32_t rng_seed() const
Thermostat for Brownian dynamics.
GammaType sigma_pos
Translational noise standard deviation.
double sigma_vel_rotation
Angular velocity noise standard deviation.
double sigma_vel
Translational velocity noise standard deviation.
GammaType sigma_pos_rotation
Rotational noise standard deviation.
static GammaType sigma(double kT, GammaType const &gamma)
Calculate the noise prefactor.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & rinertia() const
Definition Particle.hpp:502
auto const & mass() const
Definition Particle.hpp:452
auto const & quat() const
Definition Particle.hpp:477
auto const & rotation() const
Definition Particle.hpp:458
auto const & gamma() const
Definition Particle.hpp:531
bool can_rotate_around(unsigned int const axis) const
Definition Particle.hpp:461
auto const & gamma_rot() const
Definition Particle.hpp:534
auto const & torque() const
Definition Particle.hpp:479
bool is_fixed_along(unsigned int const axis) const
Definition Particle.hpp:542
auto const & force() const
Definition Particle.hpp:435
auto const & id() const
Definition Particle.hpp:414
Quaternion representation.
Utils::Vector3d bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p)
Determine the angular velocities: random walk part.
Utils::Quaternion< double > bd_drag_rot(Thermostat::GammaType const &brownian_gamma_rotation, Particle &p, double dt)
Determine quaternions: viscous drag driven by conservative torques.
Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, Particle const &p, double dt, double kT)
Determine the positions: random walk part.
Utils::Vector3d bd_drag_vel(Thermostat::GammaType const &brownian_gamma, Particle const &p)
Set the terminal velocity driven by the conservative forces drag.
Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, Particle const &p)
Determine the velocities: random walk part.
Utils::Quaternion< double > bd_random_walk_rot(BrownianThermostat const &brownian, Particle const &p, double dt, double kT)
Determine the quaternions: random walk part.
Utils::Vector3d bd_drag(Thermostat::GammaType const &brownian_gamma, Particle const &p, double dt)
Determine position: viscous drag driven by conservative forces.
Utils::Vector3d bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation, Particle const &p)
Set the terminal angular velocity driven by the conservative torques drag.