Skip to content
Permalink
Browse files

Merge #3341

3341: core: thermostat: Removed heat_up/cool_down mechanism for Langevin ty… r=KaiSzuttor a=fweik

…pe thermostats

The Langevin-type thermostats (excluding LB), had broken initial forces.
They did still employ the heat_up/cool_down mechanism which is no longer
correct with the philox rng. This is potentially a serious bug, because it makes
the thermostat silently act with a different temperature than set for at least
parts of the time steps.


Co-authored-by: Florian Weik <fweik@icp.uni-stuttgart.de>
  • Loading branch information
bors and fweik committed Dec 3, 2019
2 parents e47dbe4 + 9c4b4e9 commit 993cb7697e84bc7085f197740e8308e22d75598b
@@ -106,16 +106,6 @@ int thermalized_bond_set_params(int bond_type, double temp_com,
return ES_OK;
}

void thermalized_bond_heat_up() {
double pref_scale = sqrt(3);
thermalized_bond_update_params(pref_scale);
}

void thermalized_bond_cool_down() {
double pref_scale = 1.0 / sqrt(3);
thermalized_bond_update_params(pref_scale);
}

void thermalized_bond_init() {

for (auto &bonded_ia_param : bonded_ia_params) {
@@ -59,9 +59,6 @@ bool thermalized_bond_is_seed_required();
uint64_t thermalized_bond_get_rng_state();
void thermalized_bond_set_rng_state(uint64_t counter);

/* Setup */
void thermalized_bond_heat_up();
void thermalized_bond_cool_down();
void thermalized_bond_update_params(double pref_scale);
void thermalized_bond_init();

@@ -72,7 +69,6 @@ void thermalized_bond_init();
2. Salt (decorrelates different counter)
3. Particle ID, Partner ID
*/

inline Utils::Vector3d v_noise(int particle_id, int partner_id) {

using rng_type = r123::Philox4x64;
@@ -103,16 +103,6 @@ void dpd_set_rng_state(const uint64_t counter) {

uint64_t dpd_get_rng_state() { return dpd_rng_counter->value(); }

void dpd_heat_up() {
double pref_scale = sqrt(3);
dpd_update_params(pref_scale);
}

void dpd_cool_down() {
double pref_scale = 1.0 / sqrt(3);
dpd_update_params(pref_scale);
}

int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
int wf, double tgamma, double tr_c, int twf) {
IA_parameters &ia_params = *get_ia_param_safe(part_type_a, part_type_b);
@@ -43,8 +43,6 @@ struct DPDParameters {
double pref = 0.0;
};

void dpd_heat_up();
void dpd_cool_down();
int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
int wf, double tgamma, double tr_c, int twf);
void dpd_init();
@@ -178,8 +178,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
*/
if (reuse_forces == -1 || (recalc_forces && reuse_forces != 1)) {
ESPRESSO_PROFILER_MARK_BEGIN("Initial Force Calculation");
thermo_heat_up();

lb_lbcoupling_deactivate();

#ifdef VIRTUAL_SITES
@@ -197,8 +195,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
#endif
}

thermo_cool_down();

ESPRESSO_PROFILER_MARK_END("Initial Force Calculation");
}

@@ -161,44 +161,3 @@ void thermo_init() {
thermo_init_npt_isotropic();
#endif
}

void langevin_heat_up() {
langevin_pref2_buffer = langevin_pref2;
langevin_pref2 *= sqrt(3);

langevin_pref2_rotation_buffer = langevin_pref2_rotation;
langevin_pref2_rotation *= sqrt(3);
}

void thermo_heat_up() {
if (thermo_switch & THERMO_LANGEVIN) {
langevin_heat_up();
}
#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_heat_up();
}
#endif
if (n_thermalized_bonds) {
thermalized_bond_heat_up();
}
}

void langevin_cool_down() {
langevin_pref2 = langevin_pref2_buffer;
langevin_pref2_rotation = langevin_pref2_rotation_buffer;
}

void thermo_cool_down() {
if (thermo_switch & THERMO_LANGEVIN) {
langevin_cool_down();
}
#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_cool_down();
}
#endif
if (n_thermalized_bonds) {
thermalized_bond_cool_down();
}
}
@@ -108,18 +108,6 @@ uint64_t langevin_get_rng_state();
start of integration */
void thermo_init();

/** very nasty: if we recalculate force when leaving/reentering the integrator,
a(t) and a((t-dt)+dt) are NOT equal in the vv algorithm. The random numbers
are drawn twice, resulting in a different variance of the random force.
This is corrected by additional heat when restarting the integrator here.
Currently only works for the Langevin thermostat, although probably also
others are affected.
*/
void thermo_heat_up();

/** pendant to \ref thermo_heat_up */
void thermo_cool_down();

#ifdef NPT
/** add velocity-dependent noise and friction for NpT-sims to the particle's
velocity
@@ -186,8 +186,16 @@ def test_03__friction_rot(self):
np.copy(system.part[0].omega_body),
o0 * np.exp(-gamma_r_i / rinertia * system.time), atol=5E-4)

def test_04__global_langevin(self):
"""Test for global Langevin parameters."""
def check_global_langevin(self, recalc_forces, loops):
"""Test velocity distribution for global Langevin parameters.
Parameters
----------
recalc_forces : :obj:`bool`
True if the forces should be recalculated after every step.
loops : :obj:`int`
Number of sampling loops
"""
N = 200
system = self.system
system.part.clear()
@@ -208,11 +216,10 @@ def test_04__global_langevin(self):
system.integrator.run(20)

# Sampling
loops = 150
v_stored = np.zeros((N * loops, 3))
omega_stored = np.zeros((N * loops, 3))
for i in range(loops):
system.integrator.run(1)
system.integrator.run(1, recalc_forces=recalc_forces)
v_stored[i * N:(i + 1) * N, :] = system.part[:].v
if espressomd.has_features("ROTATION"):
omega_stored[i * N:(i + 1) * N, :] = system.part[:].omega_body
@@ -226,6 +233,16 @@ def test_04__global_langevin(self):
self.check_velocity_distribution(
omega_stored, v_minmax, bins, error_tol, kT)

def test_global_langevin(self):
"""Test velocity distribution for global Langevin parameters."""
self.check_global_langevin(False, loops=150)

def test_global_langevin_initial_forces(self):
"""Test velocity distribution for global Langevin parameters,
when using the initial force calculation.
"""
self.check_global_langevin(True, loops=170)

@utx.skipIfMissingFeatures("LANGEVIN_PER_PARTICLE")
def test_05__langevin_per_particle(self):
"""Test for Langevin particle. Covers all combinations of

0 comments on commit 993cb76

Please sign in to comment.
You can’t perform that action at this time.