From 2e4f50d048721251c3563bc07f9949b49085a129 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 26 Nov 2019 18:29:34 +0100 Subject: [PATCH 1/2] core: thermostat: Removed head_up/cool_down mechanism for Langevin type thermostats --- .../bonded_interactions/thermalized_bond.cpp | 10 ----- .../bonded_interactions/thermalized_bond.hpp | 4 -- src/core/dpd.cpp | 10 ----- src/core/dpd.hpp | 2 - src/core/integrate.cpp | 4 -- src/core/thermostat.cpp | 41 ------------------- src/core/thermostat.hpp | 12 ------ 7 files changed, 83 deletions(-) diff --git a/src/core/bonded_interactions/thermalized_bond.cpp b/src/core/bonded_interactions/thermalized_bond.cpp index af05e4757b6..297b61e4e3b 100644 --- a/src/core/bonded_interactions/thermalized_bond.cpp +++ b/src/core/bonded_interactions/thermalized_bond.cpp @@ -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) { diff --git a/src/core/bonded_interactions/thermalized_bond.hpp b/src/core/bonded_interactions/thermalized_bond.hpp index 5db22853815..07a8fa4c01c 100644 --- a/src/core/bonded_interactions/thermalized_bond.hpp +++ b/src/core/bonded_interactions/thermalized_bond.hpp @@ -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; diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index ce4b6992042..10ee586dd9e 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -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); diff --git a/src/core/dpd.hpp b/src/core/dpd.hpp index e129588f77c..32445cd1de0 100644 --- a/src/core/dpd.hpp +++ b/src/core/dpd.hpp @@ -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(); diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index cd036e68074..e605f0141f8 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -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"); } diff --git a/src/core/thermostat.cpp b/src/core/thermostat.cpp index 90b6e73a1b6..a9c1f94ae99 100644 --- a/src/core/thermostat.cpp +++ b/src/core/thermostat.cpp @@ -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(); - } -} diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 9090ce54296..8e8874cfa89 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -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 From 9c4b4e9c409e8c8b515af18a93da31a917071cde Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 28 Nov 2019 12:51:57 +0100 Subject: [PATCH 2/2] testsuite: Test for initial forces thermostat application --- testsuite/python/langevin_thermostat.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/testsuite/python/langevin_thermostat.py b/testsuite/python/langevin_thermostat.py index 7132c8fbe30..d8e2ac8a13e 100644 --- a/testsuite/python/langevin_thermostat.py +++ b/testsuite/python/langevin_thermostat.py @@ -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