29#include "thermostat.hpp"
46#ifdef THERMOSTAT_PER_PARTICLE
52 gamma = brownian_gamma;
55#ifdef PARTICLE_ANISOTROPY
57 auto const aniso_flag = (gamma[0] != gamma[1]) || (gamma[1] != gamma[2]);
61 auto const delta_pos_body = hadamard_division(force_body * dt, gamma);
67 for (
unsigned int j = 0; j < 3; j++) {
70#ifdef PARTICLE_ANISOTROPY
73 position[j] = delta_pos_lab[j];
77 position[j] = p.
force()[j] * dt / gamma[j];
82 position[j] = p.
force()[j] * dt / gamma;
99#ifdef THERMOSTAT_PER_PARTICLE
105 gamma = brownian_gamma;
108#ifdef PARTICLE_ANISOTROPY
110 auto const aniso_flag = (gamma[0] != gamma[1]) || (gamma[1] != gamma[2]);
114 auto const vel_body = hadamard_division(force_body, gamma);
120 for (
unsigned int j = 0; j < 3; j++) {
124#ifdef PARTICLE_ANISOTROPY
151 Particle const &p,
double dt,
double kT) {
153#ifdef THERMOSTAT_PER_PARTICLE
167 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_WALK>(
169 for (
unsigned int j = 0; j < 3; j++) {
171#ifndef PARTICLE_ANISOTROPY
172 if (sigma_pos > 0.0) {
173 delta_pos_body[j] = sigma_pos * sqrt(dt) * noise[j];
175 delta_pos_body[j] = 0.0;
178 if (sigma_pos[j] > 0.0) {
179 delta_pos_body[j] = sigma_pos[j] * sqrt(dt) * noise[j];
181 delta_pos_body[j] = 0.0;
187#ifdef PARTICLE_ANISOTROPY
189 auto const aniso_flag =
190 (sigma_pos[0] != sigma_pos[1]) || (sigma_pos[1] != sigma_pos[2]);
198 for (
unsigned int j = 0; j < 3; j++) {
200#ifdef PARTICLE_ANISOTROPY
201 position[j] += aniso_flag ? delta_pos_lab[j] : delta_pos_body[j];
203 position[j] += delta_pos_body[j];
217 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_INC>(
220 for (
unsigned int j = 0; j < 3; j++) {
248#ifdef THERMOSTAT_PER_PARTICLE
254 gamma = brownian_gamma_rotation;
258 for (
unsigned int j = 0; j < 3; j++) {
261#ifndef PARTICLE_ANISOTROPY
262 dphi[j] = p.
torque()[j] * dt / gamma;
264 dphi[j] = p.
torque()[j] * dt / gamma[j];
269 double dphi_m = dphi.
norm();
271 auto const dphi_u = dphi / dphi_m;
287#ifdef THERMOSTAT_PER_PARTICLE
293 gamma = brownian_gamma_rotation;
297 for (
unsigned int j = 0; j < 3; j++) {
299#ifdef PARTICLE_ANISOTROPY
300 omega[j] = p.
torque()[j] / gamma[j];
302 omega[j] = p.
torque()[j] / gamma;
318 double dt,
double kT) {
321#ifdef THERMOSTAT_PER_PARTICLE
333 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_INC>(
335 for (
unsigned int j = 0; j < 3; j++) {
337#ifndef PARTICLE_ANISOTROPY
338 if (sigma_pos > 0.0) {
339 dphi[j] = noise[j] * sigma_pos * sqrt(dt);
342 if (sigma_pos[j] > 0.0) {
343 dphi[j] = noise[j] * sigma_pos[j] * sqrt(dt);
350 double dphi_m = dphi.
norm();
352 auto const dphi_u = dphi / dphi_m;
368 auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_WALK>(
370 for (
unsigned int j = 0; j < 3; j++) {
372 domega[j] = sigma_vel * noise[j] / sqrt(p.
rinertia()[j]);
Vector implementation and trait types for boost qvm interoperability.
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.
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
Utils::Vector3d convert_vector_space_to_body(const Particle &p, const Utils::Vector3d &v)
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.
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.
auto const & rinertia() const
auto const & mass() const
auto const & quat() const
auto const & rotation() const
auto const & gamma() const
bool can_rotate_around(unsigned int const axis) const
auto const & gamma_rot() const
auto const & torque() const
bool is_fixed_along(unsigned int const axis) const
auto const & force() const
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.