core: thermostat: Removed heat_up/cool_down mechanism for Langevin ty…

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) {

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;

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);

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();

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");
   }
 

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();
-  }
-}

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

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