Skip to content

Commit

Permalink
Refactor thermostats noise sigma calculation
Browse files Browse the repository at this point in the history 
Simplify and centralize the logic to calculate the noise standard
deviation in Brownian/Langevin thermostats. Replace all the inverse
sigmas by regular sigmas in the Brownian prefactors; those were no
longer needed since the Brownian code is now separate from the
Langevin code.
  • Loading branch information
@jngrad
jngrad committed Jan 25, 2020
1 parent 53d253f commit dc47afc
Show file tree
Hide file tree
Showing 3 changed files with 70 additions and 74 deletions.
69 changes: 31 additions & 38 deletions src/core/integrators/brownian_inline.hpp
View file
Expand Up @@ -206,51 +206,48 @@ void bd_random_walk(BrownianThermostat const &brownian, Particle &p,
if (p.p.is_virtual && !thermo_virtual)
return;
// first, set defaults
Thermostat::GammaType sigma_pos_inv = brownian.sigma_pos_inv;
Thermostat::GammaType sigma_pos = brownian.sigma_pos;

// Override defaults if per-particle values for T and gamma are given
#ifdef BROWNIAN_PER_PARTICLE
auto const constexpr temp_coeff = 2.0;

if (p.p.gamma >= Thermostat::GammaType{}) {
// Is a particle-specific temperature also specified?
if (p.p.T >= 0.) {
if (p.p.T > 0.0) {
sigma_pos_inv = sqrt(p.p.gamma / (temp_coeff * p.p.T));
sigma_pos = BrownianThermostat::sigma(p.p.T, p.p.gamma);
} else {
// just an indication of the infinity
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Thermostat::GammaType{};
}
} else
// default temperature but particle-specific gamma
if (temperature > 0.0) {
sigma_pos_inv = sqrt(p.p.gamma / (temp_coeff * temperature));
sigma_pos = BrownianThermostat::sigma(temperature, p.p.gamma);
} else {
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Thermostat::GammaType{};
}
} // particle-specific gamma
else {
// No particle-specific gamma, but is there particle-specific temperature
if (p.p.T >= 0.) {
if (p.p.T > 0.0) {
sigma_pos_inv = sqrt(brownian.gamma / (temp_coeff * p.p.T));
sigma_pos = BrownianThermostat::sigma(p.p.T, brownian.gamma);
} else {
// just an indication of the infinity
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Thermostat::GammaType{};
}
} else {
// default values for both
sigma_pos_inv = brownian.sigma_pos_inv;
sigma_pos = brownian.sigma_pos;
}
}
#endif /* BROWNIAN_PER_PARTICLE */
#endif // BROWNIAN_PER_PARTICLE

bool aniso_flag; // particle anisotropy flag

#ifdef PARTICLE_ANISOTROPY
// Particle frictional isotropy check.
aniso_flag = (sigma_pos_inv[0] != sigma_pos_inv[1]) ||
(sigma_pos_inv[1] != sigma_pos_inv[2]);
aniso_flag = (sigma_pos[0] != sigma_pos[1]) || (sigma_pos[1] != sigma_pos[2]);
#else
aniso_flag = false;
#endif // PARTICLE_ANISOTROPY
Expand All @@ -267,14 +264,14 @@ void bd_random_walk(BrownianThermostat const &brownian, Particle &p,
#endif
{
#ifndef PARTICLE_ANISOTROPY
if (sigma_pos_inv > 0.0) {
delta_pos_body[j] = (1.0 / sigma_pos_inv) * sqrt(dt) * noise[j];
if (sigma_pos > 0.0) {
delta_pos_body[j] = sigma_pos * sqrt(dt) * noise[j];
} else {
delta_pos_body[j] = 0.0;
}
#else
if (sigma_pos_inv[j] > 0.0) {
delta_pos_body[j] = (1.0 / sigma_pos_inv[j]) * sqrt(dt) * noise[j];
if (sigma_pos[j] > 0.0) {
delta_pos_body[j] = sigma_pos[j] * sqrt(dt) * noise[j];
} else {
delta_pos_body[j] = 0.0;
}
Expand Down Expand Up @@ -315,17 +312,16 @@ inline void bd_random_walk_vel(BrownianThermostat const &brownian,

// Override defaults if per-particle values for T and gamma are given
#ifdef BROWNIAN_PER_PARTICLE
auto const constexpr temp_coeff = 1.0;
// Is a particle-specific temperature specified?
if (p.p.T >= 0.) {
sigma_vel = sqrt(temp_coeff * p.p.T);
sigma_vel = BrownianThermostat::sigma(p.p.T);
} else {
sigma_vel = brownian.sigma_vel;
}
#else
// defaults
sigma_vel = brownian.sigma_vel;
#endif /* BROWNIAN_PER_PARTICLE */
#endif // BROWNIAN_PER_PARTICLE

Utils::Vector3d noise = Random::v_noise_g<RNGSalt::BROWNIAN_INC>(
brownian.rng_counter->value(), p.identity());
Expand Down Expand Up @@ -438,44 +434,42 @@ void bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation,
void bd_random_walk_rot(BrownianThermostat const &brownian, Particle &p,
double dt) {
// first, set defaults
Thermostat::GammaType sigma_pos_inv = brownian.sigma_pos_rotation_inv;
Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation;

// Override defaults if per-particle values for T and gamma are given
#ifdef BROWNIAN_PER_PARTICLE
auto const constexpr temp_coeff = 2.0;

if (p.p.gamma_rot >= Thermostat::GammaType{}) {
// Is a particle-specific temperature also specified?
if (p.p.T >= 0.) {
if (p.p.T > 0.0) {
sigma_pos_inv = sqrt(p.p.gamma_rot / (temp_coeff * p.p.T));
sigma_pos = BrownianThermostat::sigma(p.p.T, p.p.gamma_rot);
} else {
// just an indication of the infinity
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Utils::Vector3d{};
}
} else if (temperature > 0.) {
// Default temperature but particle-specific gamma
sigma_pos_inv = sqrt(p.p.gamma_rot / (temp_coeff * temperature));
sigma_pos = BrownianThermostat::sigma(temperature, p.p.gamma_rot);
} else {
// just an indication of the infinity
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Utils::Vector3d{};
}
} // particle-specific gamma
else {
// No particle-specific gamma, but is there particle-specific temperature
if (p.p.T >= 0.) {
if (p.p.T > 0.0) {
sigma_pos_inv = sqrt(brownian.gamma_rotation / (temp_coeff * p.p.T));
sigma_pos = BrownianThermostat::sigma(p.p.T, brownian.gamma_rotation);
} else {
// just an indication of the infinity
sigma_pos_inv = brownian.gammatype_nan;
sigma_pos = Utils::Vector3d{};
}
} else {
// Defaut values for both
sigma_pos_inv = brownian.sigma_pos_rotation_inv;
sigma_pos = brownian.sigma_pos_rotation;
}
}
#endif /* BROWNIAN_PER_PARTICLE */
#endif // BROWNIAN_PER_PARTICLE

Utils::Vector3d dphi = {0.0, 0.0, 0.0};
Utils::Vector3d noise = Random::v_noise_g<RNGSalt::BROWNIAN_ROT_INC>(
Expand All @@ -486,14 +480,14 @@ void bd_random_walk_rot(BrownianThermostat const &brownian, Particle &p,
#endif
{
#ifndef PARTICLE_ANISOTROPY
if (sigma_pos_inv > 0.0) {
dphi[j] = noise[j] * (1.0 / sigma_pos_inv) * sqrt(dt);
if (sigma_pos > 0.0) {
dphi[j] = noise[j] * sigma_pos * sqrt(dt);
} else {
dphi[j] = 0.0;
}
#else
if (sigma_pos_inv[j] > 0.0) {
dphi[j] = noise[j] * (1.0 / sigma_pos_inv[j]) * sqrt(dt);
if (sigma_pos[j] > 0.0) {
dphi[j] = noise[j] * sigma_pos[j] * sqrt(dt);
} else {
dphi[j] = 0.0;
}
Expand All @@ -520,17 +514,16 @@ void bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle &p) {

// Override defaults if per-particle values for T and gamma are given
#ifdef BROWNIAN_PER_PARTICLE
auto const constexpr temp_coeff = 1.0;
// Is a particle-specific temperature specified?
if (p.p.T >= 0.) {
sigma_vel = sqrt(temp_coeff * p.p.T);
sigma_vel = BrownianThermostat::sigma(p.p.T);
} else {
sigma_vel = brownian.sigma_vel_rotation;
}
#else
// set defaults
sigma_vel = brownian.sigma_vel_rotation;
#endif /* BROWNIAN_PER_PARTICLE */
#endif // BROWNIAN_PER_PARTICLE

Utils::Vector3d domega;
Utils::Vector3d noise = Random::v_noise_g<RNGSalt::BROWNIAN_ROT_WALK>(
Expand Down
10 changes: 4 additions & 6 deletions src/core/integrators/langevin_inline.hpp
View file
Expand Up @@ -52,14 +52,13 @@ friction_thermo_langevin(LangevinThermostat const &langevin,
// Override defaults if per-particle values for T and gamma are given
#ifdef LANGEVIN_PER_PARTICLE
if (p.p.gamma >= Thermostat::GammaType{} or p.p.T >= 0.) {
auto const constexpr temp_coeff = 24.0;
auto const kT = p.p.T >= 0. ? p.p.T : temperature;
auto const gamma =
p.p.gamma >= Thermostat::GammaType{} ? p.p.gamma : langevin.gamma;
pref_friction = -gamma;
pref_noise = sqrt(temp_coeff * kT * gamma / time_step);
pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
}
#endif /* LANGEVIN_PER_PARTICLE */
#endif // LANGEVIN_PER_PARTICLE

// Get effective velocity in the thermostatting
#ifdef ENGINE
Expand Down Expand Up @@ -109,15 +108,14 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin,
// Override defaults if per-particle values for T and gamma are given
#ifdef LANGEVIN_PER_PARTICLE
if (p.p.gamma_rot >= Thermostat::GammaType{} or p.p.T >= 0.) {
auto const constexpr temp_coeff = 24.0;
auto const kT = p.p.T >= 0. ? p.p.T : temperature;
auto const gamma = p.p.gamma_rot >= Thermostat::GammaType{}
? p.p.gamma_rot
: langevin.gamma_rotation;
pref_friction = -gamma;
pref_noise = sqrt(temp_coeff * kT * gamma / time_step);
pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
}
#endif /* LANGEVIN_PER_PARTICLE */
#endif // LANGEVIN_PER_PARTICLE

using Random::v_noise;

Expand Down
65 changes: 35 additions & 30 deletions src/core/thermostat.hpp
View file
Expand Up @@ -69,8 +69,6 @@ constexpr double sentinel(double) { return -1.0; }
constexpr Utils::Vector3d sentinel(Utils::Vector3d) {
return {-1.0, -1.0, -1.0};
}
constexpr double set_nan(double) { return NAN; }
constexpr Utils::Vector3d set_nan(Utils::Vector3d) { return {NAN, NAN, NAN}; }
/*@}*/
} // namespace

Expand Down Expand Up @@ -107,12 +105,17 @@ struct LangevinThermostat {
*/
void recalc_prefactors() {
pref_friction = -gamma;
pref_noise = sqrt(24.0 * temperature / time_step * gamma);
pref_noise = sigma(temperature, time_step, gamma);
// If gamma_rotation is not set explicitly, use the translational one.
if (gamma_rotation < GammaType{}) {
gamma_rotation = gamma;
}
pref_noise_rotation = sqrt(24.0 * temperature / time_step * gamma_rotation);
pref_noise_rotation = sigma(temperature, time_step, gamma_rotation);
}
/** Calculate the noise standard deviation. */
static GammaType sigma(double kT, double time_step, GammaType const &gamma) {
constexpr auto const temp_coeff = 24.0;
return sqrt((temp_coeff * kT / time_step) * gamma);
}
/** @name Parameters */
/*@{*/
Expand Down Expand Up @@ -150,17 +153,12 @@ struct BrownianThermostat {
* which is only valid for the BD. Just a square root of kT, see (10.2.17)
* and comments in 2 paragraphs afterwards, @cite Pottier2010.
*/
sigma_vel = sqrt(temperature);
sigma_vel = sigma(temperature);
/** The random walk position dispersion is defined by the second eq. (14.38)
* of @cite Schlick2010. Its time interval factor will be added in the
* Brownian Dynamics functions. Its square root is the standard deviation.
*/
if (temperature > 0.0) {
sigma_pos_inv = sqrt(gamma / (2.0 * temperature));
} else {
// just an indication of the infinity
sigma_pos_inv = gammatype_nan;
}
sigma_pos = sigma(temperature, gamma);
#ifdef ROTATION
/** Note: the BD thermostat assigns the brownian viscous parameters as well.
* They correspond to the friction tensor Z from the eq. (14.31) of
Expand All @@ -170,41 +168,48 @@ struct BrownianThermostat {
if (gamma_rotation < GammaType{}) {
gamma_rotation = gamma;
}
sigma_vel_rotation = sqrt(temperature);
if (temperature > 0.0) {
sigma_pos_rotation_inv = sqrt(gamma_rotation / (2.0 * temperature));
} else {
// just an indication of the infinity
sigma_pos_rotation_inv = gammatype_nan;
}
sigma_vel_rotation = sigma(temperature);
sigma_pos_rotation = sigma(temperature, gamma_rotation);
#endif // ROTATION
}
/** Calculate the noise standard deviation. */
static GammaType sigma(double kT, GammaType const &gamma) {
constexpr auto const temp_coeff = 2.0;
return sqrt(Utils::hadamard_division(temp_coeff * kT, gamma));
}
/** Calculate the noise standard deviation. */
static double sigma(double kT) {
constexpr auto const temp_coeff = 1.0;
return sqrt(temp_coeff * kT);
}
/** @name Parameters */
/*@{*/
/** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
GammaType gamma = sentinel(GammaType{});
/** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
GammaType gamma_rotation = sentinel(GammaType{});
/*@}*/
/** Sentinel value for divisions by zero. */
GammaType const gammatype_nan = set_nan(GammaType{});
/** @name Prefactors */
/*@{*/
/** Inverse of the translational noise standard deviation.
* Stores @f$ \left(\sqrt{2D_{\text{trans}}}\right)^{-1} @f$ with
/** Translational noise standard deviation.
* Stores @f$ \sqrt{2D_{\text{trans}}} @f$ with
* @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$
* the translational diffusion coefficient
* the translational diffusion coefficient.
*/
GammaType sigma_pos_inv = sentinel(GammaType{});
/** Inverse of the rotational noise standard deviation.
* Stores @f$ \left(\sqrt{2D_{\text{rot}}}\right)^{-1} @f$ with
GammaType sigma_pos = sentinel(GammaType{});
/** Rotational noise standard deviation.
* Stores @f$ \sqrt{2D_{\text{rot}}} @f$ with
* @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$
* the rotational diffusion coefficient
* the rotational diffusion coefficient.
*/
GammaType sigma_pos_rotation = sentinel(GammaType{});
/** Translational velocity noise standard deviation.
* Stores @f$ \sqrt{k_B T} @f$.
*/
GammaType sigma_pos_rotation_inv = sentinel(GammaType{});
/** Translational velocity noise standard deviation. */
double sigma_vel = 0;
/** Angular velocity noise standard deviation. */
/** Angular velocity noise standard deviation.
* Stores @f$ \sqrt{k_B T} @f$.
*/
double sigma_vel_rotation = 0;
/*@}*/
/** RNG counter, used for both translation and rotation. */
Expand Down

0 comments on commit dc47afc

Please sign in to comment.