ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/thermostat.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/** \file
25 * Implementation in \ref thermostat.cpp.
26 */
27
28#include "Particle.hpp"
29#include "PropagationMode.hpp"
30#include "rotation.hpp"
31#include "system/Leaf.hpp"
32
33#include "config/config.hpp"
34
35#include <utils/Counter.hpp>
36#include <utils/Vector.hpp>
37
38#include <cassert>
39#include <cmath>
40#include <cstdint>
41
42namespace Thermostat {
43#ifdef PARTICLE_ANISOTROPY
45#else
46using GammaType = double;
47#endif
48/**
49 * @brief Value for unset friction coefficient.
50 * Sentinel value for the Langevin/Brownian parameters,
51 * indicating that they have not been set yet.
52 */
53#ifdef PARTICLE_ANISOTROPY
54constexpr GammaType gamma_sentinel{{-1.0, -1.0, -1.0}};
55#else
56constexpr GammaType gamma_sentinel{-1.0};
57#endif
58/**
59 * @brief Value for a null friction coefficient.
60 */
61#ifdef PARTICLE_ANISOTROPY
62constexpr GammaType gamma_null{{0.0, 0.0, 0.0}};
63#else
64constexpr GammaType gamma_null{0.0};
65#endif
66
67#ifdef THERMOSTAT_PER_PARTICLE
68inline auto const &handle_particle_gamma(GammaType const &particle_gamma,
69 GammaType const &default_gamma) {
70 return particle_gamma >= gamma_null ? particle_gamma : default_gamma;
71}
72#endif
73
75 GammaType const &gamma_body) {
76#ifdef PARTICLE_ANISOTROPY
77 auto const aniso_flag =
78 (gamma_body[0] != gamma_body[1]) || (gamma_body[1] != gamma_body[2]);
79 const Utils::Matrix<double, 3, 3> gamma_matrix =
80 boost::qvm::diag_mat(gamma_body);
81 auto const gamma_space =
82 aniso_flag ? convert_body_to_space(p, gamma_matrix) : gamma_matrix;
83 return gamma_space;
84#else
85 return gamma_body;
86#endif
87}
88
89/** @brief Check that two kT values are close up to a small tolerance. */
90inline bool are_kT_equal(double old_kT, double new_kT) {
91 constexpr auto relative_tolerance = 1e-6;
92 if (old_kT == 0. and new_kT == 0.) {
93 return true;
94 }
95 if ((old_kT < 0. and new_kT >= 0.) or (old_kT >= 0. and new_kT < 0.)) {
96 return false;
97 }
98 auto const large_kT = (old_kT > new_kT) ? old_kT : new_kT;
99 auto const small_kT = (old_kT > new_kT) ? new_kT : old_kT;
100 return (large_kT / small_kT - 1. < relative_tolerance);
101}
102} // namespace Thermostat
103
105public:
106 /** Initialize or re-initialize the RNG counter with a seed. */
107 void rng_initialize(uint32_t const seed) {
108 m_rng_seed = seed;
109 m_initialized = true;
110 }
111 /** Increment the RNG counter */
112 void rng_increment() { m_rng_counter.increment(); }
113 /** Get current value of the RNG */
114 uint64_t rng_counter() const { return m_rng_counter.value(); }
115 void set_rng_counter(uint64_t value) {
116 m_rng_counter = Utils::Counter<uint64_t>(uint64_t{0u}, value);
117 }
118 /** Is the RNG seed required */
119 bool is_seed_required() const { return not m_initialized; }
120 uint32_t rng_seed() const { return m_rng_seed; }
121
122private:
123 /** @brief RNG counter. */
124 Utils::Counter<uint64_t> m_rng_counter{uint64_t{0u}, uint64_t{0u}};
125 /** @brief RNG seed. */
126 uint32_t m_rng_seed{0u};
127 bool m_initialized{false};
128};
129
130/** Thermostat for Langevin dynamics. */
132private:
134
135public:
136 /** Recalculate prefactors.
137 * Needs to be called every time the parameters are changed.
138 */
139 void recalc_prefactors(double kT, double time_step) {
141 pref_noise = sigma(kT, time_step, gamma);
142#ifdef ROTATION
143 pref_noise_rotation = sigma(kT, time_step, gamma_rotation);
144#endif // ROTATION
145 }
146 /** Calculate the noise prefactor.
147 * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma / dt} / \sigma_\eta @f$
148 * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform
149 * process @f$ \eta(t) @f$.
150 */
151 static GammaType sigma(double kT, double time_step, GammaType const &gamma) {
152 // random uniform noise has variance 1/12
153 constexpr auto const temp_coeff = 2.0 * 12.0;
154 return sqrt((temp_coeff * kT / time_step) * gamma);
155 }
156 /** @name Parameters */
157 /**@{*/
158 /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
160#ifdef ROTATION
161 /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
163#endif // ROTATION
164 /**@}*/
165 /** @name Prefactors */
166 /**@{*/
167 /** Prefactor for the friction.
168 * Stores @f$ \gamma_{\text{trans}} @f$.
169 */
171 /** Prefactor for the translational velocity noise.
172 * Stores @f$ \sqrt{2 k_B T \gamma_{\text{trans}} / dt} / \sigma_\eta @f$.
173 */
175#ifdef ROTATION
176 /** Prefactor for the angular velocity noise.
177 * Stores @f$ \sqrt{2 k_B T \gamma_{\text{rot}} / dt} / \sigma_\eta @f$.
178 */
180#endif // ROTATION
181 /**@}*/
182};
183
184/** Thermostat for Brownian dynamics.
185 * Default particle mass is assumed to be unitary in these global parameters.
186 */
188private:
190
191public:
192 /** Recalculate prefactors.
193 * Needs to be called every time the parameters are changed.
194 */
195 void recalc_prefactors(double kT) {
196 /** The heat velocity dispersion corresponds to the Gaussian noise only,
197 * which is only valid for the BD. Just a square root of kT, see (10.2.17)
198 * and comments in 2 paragraphs afterwards, @cite pottier10a.
199 */
200 sigma_vel = sigma(kT);
201 /** The random walk position dispersion is defined by the second eq. (14.38)
202 * of @cite schlick10a. Its time interval factor will be added in the
203 * Brownian Dynamics functions. Its square root is the standard deviation.
204 */
205 sigma_pos = sigma(kT, gamma);
206#ifdef ROTATION
207 /** Note: the BD thermostat assigns the brownian viscous parameters as well.
208 * They correspond to the friction tensor Z from the eq. (14.31) of
209 * @cite schlick10a.
210 */
213#endif // ROTATION
214 }
215 /** Calculate the noise prefactor.
216 * Evaluates the quantity @f$ \sqrt{2 k_B T / \gamma} / \sigma_\eta @f$
217 * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian
218 * process @f$ \eta(t) @f$.
219 */
220 static GammaType sigma(double kT, GammaType const &gamma) {
221 constexpr auto const temp_coeff = 2.0;
222 return sqrt(Utils::hadamard_division(temp_coeff * kT, gamma));
223 }
224 /** Calculate the noise prefactor.
225 * Evaluates the quantity @f$ \sqrt{k_B T} / \sigma_\eta @f$
226 * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian
227 * process @f$ \eta(t) @f$.
228 */
229 static double sigma(double kT) {
230 constexpr auto const temp_coeff = 1.0;
231 return sqrt(temp_coeff * kT);
232 }
233 /** @name Parameters */
234 /**@{*/
235 /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
237 /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
239 /**@}*/
240 /** @name Prefactors */
241 /**@{*/
242 /** Translational noise standard deviation.
243 * Stores @f$ \sqrt{2D_{\text{trans}}} @f$ with
244 * @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$
245 * the translational diffusion coefficient.
246 */
248#ifdef ROTATION
249 /** Rotational noise standard deviation.
250 * Stores @f$ \sqrt{2D_{\text{rot}}} @f$ with
251 * @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$
252 * the rotational diffusion coefficient.
253 */
255#endif // ROTATION
256 /** Translational velocity noise standard deviation.
257 * Stores @f$ \sqrt{k_B T} @f$.
258 */
259 double sigma_vel = 0.;
260#ifdef ROTATION
261 /** Angular velocity noise standard deviation.
262 * Stores @f$ \sqrt{k_B T} @f$.
263 */
265#endif // ROTATION
266 /**@}*/
267};
268
269#ifdef NPT
270/** Thermostat for isotropic NPT dynamics. */
272private:
274
275public:
276 /** Recalculate prefactors.
277 * Needs to be called every time the parameters are changed.
278 */
279 void recalc_prefactors(double kT, double piston, double time_step) {
280 assert(piston > 0.0);
281 auto const half_time_step = time_step / 2.0;
282 pref_rescale_0 = -gamma0 * half_time_step;
283 pref_noise_0 = sigma(kT, gamma0, time_step);
284 pref_rescale_V = -gammav * half_time_step / piston;
285 pref_noise_V = sigma(kT, gammav, time_step);
286 }
287 /** Calculate the noise prefactor.
288 * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma dt / 2} / \sigma_\eta @f$
289 * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform
290 * process @f$ \eta(t) @f$.
291 */
292 static double sigma(double kT, double gamma, double time_step) {
293 // random uniform noise has variance 1/12; the temperature
294 // coefficient of 2 is canceled out by the half time step
295 constexpr auto const temp_coeff = 12.0;
296 return sqrt(temp_coeff * kT * gamma * time_step);
297 }
298 /** @name Parameters */
299 /**@{*/
300 /** Friction coefficient of the particles @f$ \gamma^0 @f$ */
301 double gamma0 = 0.;
302 /** Friction coefficient for the box @f$ \gamma^V @f$ */
303 double gammav = 0.;
304 /**@}*/
305 /** @name Prefactors */
306 /**@{*/
307 /** Particle velocity rescaling at half the time step.
308 * Stores @f$ \gamma^{0}\cdot\frac{dt}{2} @f$.
309 */
310 double pref_rescale_0 = 0.;
311 /** Particle velocity rescaling noise standard deviation.
312 * Stores @f$ \sqrt{k_B T \gamma^{0} dt} / \sigma_\eta @f$.
313 */
314 double pref_noise_0 = 0.;
315 /** Volume rescaling at half the time step.
316 * Stores @f$ \frac{\gamma^{V}}{Q}\cdot\frac{dt}{2} @f$.
317 */
318 double pref_rescale_V = 0.;
319 /** Volume rescaling noise standard deviation.
320 * Stores @f$ \sqrt{k_B T \gamma^{V} dt} / \sigma_\eta @f$.
321 */
322 double pref_noise_V = 0.;
323 /**@}*/
324};
325#endif
326
327/** Thermostat for lattice-Boltzmann particle coupling. */
329 /** @name Parameters */
330 /**@{*/
331 /** Friction coefficient. */
332 double gamma = -1.;
333 /** Internal flag to disable particle coupling during force recalculation. */
334 bool couple_to_md = false;
335 /**@}*/
336};
337
339
340/** Thermostat for thermalized bonds. */
342 void recalc_prefactors(double time_step, BondedInteractionsMap &bonded_ias);
343};
344
345#ifdef DPD
346/** Thermostat for dissipative particle dynamics. */
347struct DPDThermostat : public BaseThermostat {};
348#endif
349
350#ifdef STOKESIAN_DYNAMICS
351/** Thermostat for Stokesian dynamics. */
355#endif
356
357namespace Thermostat {
358class Thermostat : public System::Leaf<Thermostat> {
359public:
360 /** @brief Thermal energy of the simulated heat bath. */
361 double kT = -1.;
362 /** @brief Bitmask of currently active thermostats. */
363 int thermo_switch = THERMO_OFF;
364 std::shared_ptr<LangevinThermostat> langevin;
365 std::shared_ptr<BrownianThermostat> brownian;
366#ifdef NPT
367 std::shared_ptr<IsotropicNptThermostat> npt_iso;
368#endif
369 std::shared_ptr<LBThermostat> lb;
370#ifdef DPD
371 std::shared_ptr<DPDThermostat> dpd;
372#endif
373#ifdef STOKESIAN_DYNAMICS
374 std::shared_ptr<StokesianThermostat> stokesian;
375#endif
376 std::shared_ptr<ThermalizedBondThermostat> thermalized_bond;
377
378 /** Increment RNG counters */
379 void philox_counter_increment();
380
381 /** Initialize constants of all thermostats. */
382 void recalc_prefactors(double time_step);
383
385 if (lb) {
386 lb->couple_to_md = true;
387 }
388 }
389 void lb_coupling_deactivate();
390};
391} // namespace Thermostat
@ THERMO_OFF
Vector implementation and trait types for boost qvm interoperability.
container for bonded interactions.
Abstract class that represents a component of the system.
std::shared_ptr< BrownianThermostat > brownian
std::shared_ptr< ThermalizedBondThermostat > thermalized_bond
std::shared_ptr< LangevinThermostat > langevin
std::shared_ptr< DPDThermostat > dpd
std::shared_ptr< LBThermostat > lb
std::shared_ptr< IsotropicNptThermostat > npt_iso
std::shared_ptr< StokesianThermostat > stokesian
T value() const
Definition Counter.hpp:43
void increment()
Definition Counter.hpp:41
This file contains the defaults for ESPResSo.
bool are_kT_equal(double old_kT, double new_kT)
Check that two kT values are close up to a small tolerance.
Utils::Vector3d GammaType
constexpr GammaType gamma_sentinel
Value for unset friction coefficient.
constexpr GammaType gamma_null
Value for a null friction coefficient.
auto handle_particle_anisotropy(Particle const &p, GammaType const &gamma_body)
auto const & handle_particle_gamma(GammaType const &particle_gamma, GammaType const &default_gamma)
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
auto hadamard_division(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:422
This file contains all subroutines required to process rotational motion.
auto convert_body_to_space(const Utils::Quaternion< double > &quat, const Utils::Matrix< T, 3, 3 > &A)
Transform matrix from body- to space-fixed frame.
Definition rotation.hpp:86
uint64_t rng_counter() const
Get current value of the RNG.
uint32_t rng_seed() const
void rng_increment()
Increment the RNG counter.
bool is_seed_required() const
Is the RNG seed required.
void rng_initialize(uint32_t const seed)
Initialize or re-initialize the RNG counter with a seed.
void set_rng_counter(uint64_t value)
Thermostat for Brownian dynamics.
GammaType gamma
Translational friction coefficient .
GammaType gamma_rotation
Rotational friction coefficient .
void recalc_prefactors(double kT)
Recalculate prefactors.
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.
static double sigma(double kT)
Calculate the noise prefactor.
Thermostat for dissipative particle dynamics.
Thermostat for isotropic NPT dynamics.
double pref_rescale_V
Volume rescaling at half the time step.
double pref_rescale_0
Particle velocity rescaling at half the time step.
double gamma0
Friction coefficient of the particles .
double pref_noise_0
Particle velocity rescaling noise standard deviation.
double gammav
Friction coefficient for the box .
double pref_noise_V
Volume rescaling noise standard deviation.
static double sigma(double kT, double gamma, double time_step)
Calculate the noise prefactor.
void recalc_prefactors(double kT, double piston, double time_step)
Recalculate prefactors.
Thermostat for lattice-Boltzmann particle coupling.
bool couple_to_md
Internal flag to disable particle coupling during force recalculation.
double gamma
Friction coefficient.
Thermostat for Langevin dynamics.
GammaType pref_friction
Prefactor for the friction.
static GammaType sigma(double kT, double time_step, GammaType const &gamma)
Calculate the noise prefactor.
GammaType gamma_rotation
Rotational friction coefficient .
GammaType gamma
Translational friction coefficient .
void recalc_prefactors(double kT, double time_step)
Recalculate prefactors.
GammaType pref_noise
Prefactor for the translational velocity noise.
GammaType pref_noise_rotation
Prefactor for the angular velocity noise.
Struct holding all information for one particle.
Definition Particle.hpp:395
Thermostat for Stokesian dynamics.
Thermostat for thermalized bonds.
void recalc_prefactors(double time_step, BondedInteractionsMap &bonded_ias)
Matrix representation with static size.
Definition matrix.hpp:65