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 if (small_kT == 0.) {
101 return false;
102 }
103 return (large_kT / small_kT - 1. < relative_tolerance);
104}
105} // namespace Thermostat
106
108public:
109 /** Initialize or re-initialize the RNG counter with a seed. */
110 void rng_initialize(uint32_t const seed) {
111 m_rng_seed = seed;
112 m_initialized = true;
113 }
114 /** Increment the RNG counter */
115 void rng_increment() { m_rng_counter.increment(); }
116 /** Get current value of the RNG */
117 uint64_t rng_counter() const { return m_rng_counter.value(); }
118 void set_rng_counter(uint64_t value) {
119 m_rng_counter = Utils::Counter<uint64_t>(uint64_t{0u}, value);
120 }
121 /** Is the RNG seed required */
122 bool is_seed_required() const { return not m_initialized; }
123 uint32_t rng_seed() const { return m_rng_seed; }
124
125private:
126 /** @brief RNG counter. */
127 Utils::Counter<uint64_t> m_rng_counter{uint64_t{0u}, uint64_t{0u}};
128 /** @brief RNG seed. */
129 uint32_t m_rng_seed{0u};
130 bool m_initialized{false};
131};
132
133/** Thermostat for Langevin dynamics. */
135private:
137
138public:
139 /** Recalculate prefactors.
140 * Needs to be called every time the parameters are changed.
141 */
142 void recalc_prefactors(double kT, double time_step) {
144 pref_noise = sigma(kT, time_step, gamma);
145#ifdef ROTATION
146 pref_noise_rotation = sigma(kT, time_step, gamma_rotation);
147#endif // ROTATION
148 }
149 /** Calculate the noise prefactor.
150 * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma / dt} / \sigma_\eta @f$
151 * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform
152 * process @f$ \eta(t) @f$.
153 */
154 static GammaType sigma(double kT, double time_step, GammaType const &gamma) {
155 // random uniform noise has variance 1/12
156 constexpr auto const temp_coeff = 2.0 * 12.0;
157 return sqrt((temp_coeff * kT / time_step) * gamma);
158 }
159 /** @name Parameters */
160 /**@{*/
161 /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
163#ifdef ROTATION
164 /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
166#endif // ROTATION
167 /**@}*/
168 /** @name Prefactors */
169 /**@{*/
170 /** Prefactor for the friction.
171 * Stores @f$ \gamma_{\text{trans}} @f$.
172 */
174 /** Prefactor for the translational velocity noise.
175 * Stores @f$ \sqrt{2 k_B T \gamma_{\text{trans}} / dt} / \sigma_\eta @f$.
176 */
178#ifdef ROTATION
179 /** Prefactor for the angular velocity noise.
180 * Stores @f$ \sqrt{2 k_B T \gamma_{\text{rot}} / dt} / \sigma_\eta @f$.
181 */
183#endif // ROTATION
184 /**@}*/
185};
186
187/** Thermostat for Brownian dynamics.
188 * Default particle mass is assumed to be unitary in these global parameters.
189 */
191private:
193
194public:
195 /** Recalculate prefactors.
196 * Needs to be called every time the parameters are changed.
197 */
198 void recalc_prefactors(double kT) {
199 /** The heat velocity dispersion corresponds to the Gaussian noise only,
200 * which is only valid for the BD. Just a square root of kT, see (10.2.17)
201 * and comments in 2 paragraphs afterwards, @cite pottier10a.
202 */
203 sigma_vel = sigma(kT);
204 /** The random walk position dispersion is defined by the second eq. (14.38)
205 * of @cite schlick10a. Its time interval factor will be added in the
206 * Brownian Dynamics functions. Its square root is the standard deviation.
207 */
208 sigma_pos = sigma(kT, gamma);
209#ifdef ROTATION
210 /** Note: the BD thermostat assigns the brownian viscous parameters as well.
211 * They correspond to the friction tensor Z from the eq. (14.31) of
212 * @cite schlick10a.
213 */
216#endif // ROTATION
217 }
218 /** Calculate the noise prefactor.
219 * Evaluates the quantity @f$ \sqrt{2 k_B T / \gamma} / \sigma_\eta @f$
220 * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian
221 * process @f$ \eta(t) @f$.
222 */
223 static GammaType sigma(double kT, GammaType const &gamma) {
224 constexpr auto const temp_coeff = 2.0;
225 return sqrt(Utils::hadamard_division(temp_coeff * kT, gamma));
226 }
227 /** Calculate the noise prefactor.
228 * Evaluates the quantity @f$ \sqrt{k_B T} / \sigma_\eta @f$
229 * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian
230 * process @f$ \eta(t) @f$.
231 */
232 static double sigma(double kT) {
233 constexpr auto const temp_coeff = 1.0;
234 return sqrt(temp_coeff * kT);
235 }
236 /** @name Parameters */
237 /**@{*/
238 /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
240 /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
242 /**@}*/
243 /** @name Prefactors */
244 /**@{*/
245 /** Translational noise standard deviation.
246 * Stores @f$ \sqrt{2D_{\text{trans}}} @f$ with
247 * @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$
248 * the translational diffusion coefficient.
249 */
251#ifdef ROTATION
252 /** Rotational noise standard deviation.
253 * Stores @f$ \sqrt{2D_{\text{rot}}} @f$ with
254 * @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$
255 * the rotational diffusion coefficient.
256 */
258#endif // ROTATION
259 /** Translational velocity noise standard deviation.
260 * Stores @f$ \sqrt{k_B T} @f$.
261 */
262 double sigma_vel = 0.;
263#ifdef ROTATION
264 /** Angular velocity noise standard deviation.
265 * Stores @f$ \sqrt{k_B T} @f$.
266 */
268#endif // ROTATION
269 /**@}*/
270};
271
272#ifdef NPT
273/** Thermostat for isotropic NPT dynamics. */
275private:
277
278public:
279 /** Recalculate prefactors.
280 * Needs to be called every time the parameters are changed.
281 */
282 void recalc_prefactors(double kT, double piston, double time_step) {
283 assert(piston > 0.0);
284 auto const half_time_step = time_step / 2.0;
285 pref_rescale_0 = -gamma0 * half_time_step;
286 pref_noise_0 = sigma(kT, gamma0, time_step);
287 pref_rescale_V = -gammav * half_time_step / piston;
288 pref_noise_V = sigma(kT, gammav, time_step);
289 }
290 /** Calculate the noise prefactor.
291 * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma dt / 2} / \sigma_\eta @f$
292 * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform
293 * process @f$ \eta(t) @f$.
294 */
295 static double sigma(double kT, double gamma, double time_step) {
296 // random uniform noise has variance 1/12; the temperature
297 // coefficient of 2 is canceled out by the half time step
298 constexpr auto const temp_coeff = 12.0;
299 return sqrt(temp_coeff * kT * gamma * time_step);
300 }
301 /** @name Parameters */
302 /**@{*/
303 /** Friction coefficient of the particles @f$ \gamma^0 @f$ */
304 double gamma0 = 0.;
305 /** Friction coefficient for the box @f$ \gamma^V @f$ */
306 double gammav = 0.;
307 /**@}*/
308 /** @name Prefactors */
309 /**@{*/
310 /** Particle velocity rescaling at half the time step.
311 * Stores @f$ \gamma^{0}\cdot\frac{dt}{2} @f$.
312 */
313 double pref_rescale_0 = 0.;
314 /** Particle velocity rescaling noise standard deviation.
315 * Stores @f$ \sqrt{k_B T \gamma^{0} dt} / \sigma_\eta @f$.
316 */
317 double pref_noise_0 = 0.;
318 /** Volume rescaling at half the time step.
319 * Stores @f$ \frac{\gamma^{V}}{Q}\cdot\frac{dt}{2} @f$.
320 */
321 double pref_rescale_V = 0.;
322 /** Volume rescaling noise standard deviation.
323 * Stores @f$ \sqrt{k_B T \gamma^{V} dt} / \sigma_\eta @f$.
324 */
325 double pref_noise_V = 0.;
326 /**@}*/
327};
328#endif
329
330/** Thermostat for lattice-Boltzmann particle coupling. */
332 /** @name Parameters */
333 /**@{*/
334 /** Friction coefficient. */
335 double gamma = -1.;
336 /** Internal flag to disable particle coupling during force recalculation. */
337 bool couple_to_md = false;
338 /**@}*/
339};
340
342
343/** Thermostat for thermalized bonds. */
345 void recalc_prefactors(double time_step, BondedInteractionsMap &bonded_ias);
346};
347
348#ifdef DPD
349/** Thermostat for dissipative particle dynamics. */
350struct DPDThermostat : public BaseThermostat {};
351#endif
352
353#ifdef STOKESIAN_DYNAMICS
354/** Thermostat for Stokesian dynamics. */
358#endif
359
360namespace Thermostat {
361class Thermostat : public System::Leaf<Thermostat> {
362public:
363 /** @brief Thermal energy of the simulated heat bath. */
364 double kT = -1.;
365 /** @brief Bitmask of currently active thermostats. */
366 int thermo_switch = THERMO_OFF;
367 std::shared_ptr<LangevinThermostat> langevin;
368 std::shared_ptr<BrownianThermostat> brownian;
369#ifdef NPT
370 std::shared_ptr<IsotropicNptThermostat> npt_iso;
371#endif
372 std::shared_ptr<LBThermostat> lb;
373#ifdef DPD
374 std::shared_ptr<DPDThermostat> dpd;
375#endif
376#ifdef STOKESIAN_DYNAMICS
377 std::shared_ptr<StokesianThermostat> stokesian;
378#endif
379 std::shared_ptr<ThermalizedBondThermostat> thermalized_bond;
380
381 /** Increment RNG counters */
382 void philox_counter_increment();
383
384 /** Initialize constants of all thermostats. */
385 void recalc_prefactors(double time_step);
386
388 if (lb) {
389 lb->couple_to_md = true;
390 }
391 }
392 void lb_coupling_deactivate();
393};
394} // 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