Refactor thermostats and add unit tests (#3438)

- Gather thermostat global variables in structs (#2628)
- Make thermostat functions pure functions to eliminate side effects
- Add Boost unit tests for thermostats
- Reduce thermostat code complexity

src/core/forces_inline.hpp

@@ -38,6 +38,7 @@
 #include "forces.hpp"
 #include "immersed_boundary/ibm_tribend.hpp"
 #include "immersed_boundary/ibm_triel.hpp"
+#include "integrators/langevin_inline.hpp"
 #include "metadynamics.hpp"
 #include "nonbonded_interactions/bmhtf-nacl.hpp"
 #include "nonbonded_interactions/buckingham.hpp"
@@ -101,16 +102,17 @@ inline ParticleForce external_force(Particle const &p) {
 }
 
 inline ParticleForce thermostat_force(Particle const &p) {
+  extern LangevinThermostat langevin;
   if (!(thermo_switch & THERMO_LANGEVIN)) {
     return {};
   }
 
 #ifdef ROTATION
-  return {
-      friction_thermo_langevin(p),
-      convert_vector_body_to_space(p, friction_thermo_langevin_rotation(p))};
+  return {friction_thermo_langevin(langevin, p),
+          convert_vector_body_to_space(
+              p, friction_thermo_langevin_rotation(langevin, p))};
 #else
-  return friction_thermo_langevin(p);
+  return friction_thermo_langevin(langevin, p);
 #endif
 }
 

src/core/global.cpp

@@ -47,6 +47,9 @@
 #include <unordered_map>
 
 extern double force_cap;
+extern LangevinThermostat langevin;
+extern BrownianThermostat brownian;
+extern IsotropicNptThermostat npt_iso;
 
 namespace {
 
@@ -78,12 +81,12 @@ const std::unordered_map<int, Datafield> fields{
       "cell_grid"}}, /* 1  from cells.cpp */
 #ifndef PARTICLE_ANISOTROPY
     {FIELD_LANGEVIN_GAMMA,
-     {&langevin_gamma, Datafield::Type::DOUBLE, 1,
-      "gamma"}}, /* 5  from thermostat.cpp */
+     {&langevin.gamma, Datafield::Type::DOUBLE, 1,
+      "langevin.gamma"}}, /* 5  from thermostat.cpp */
 #else
     {FIELD_LANGEVIN_GAMMA,
-     {langevin_gamma.data(), Datafield::Type::DOUBLE, 3,
-      "gamma"}}, /* 5  from thermostat.cpp */
+     {langevin.gamma.data(), Datafield::Type::DOUBLE, 3,
+      "langevin.gamma"}}, /* 5  from thermostat.cpp */
 #endif // PARTICLE_ANISOTROPY
     {FIELD_INTEG_SWITCH,
      {&integ_switch, Datafield::Type::INT, 1,
@@ -107,11 +110,11 @@ const std::unordered_map<int, Datafield> fields{
      {node_grid.data(), Datafield::Type::INT, 3,
       "node_grid"}}, /* 20 from grid.cpp */
     {FIELD_NPTISO_G0,
-     {&nptiso_gamma0, Datafield::Type::DOUBLE, 1,
-      "nptiso_gamma0"}}, /* 21 from thermostat.cpp */
+     {&npt_iso.gamma0, Datafield::Type::DOUBLE, 1,
+      "npt_iso.gamma0"}}, /* 21 from thermostat.cpp */
     {FIELD_NPTISO_GV,
-     {&nptiso_gammav, Datafield::Type::DOUBLE, 1,
-      "nptiso_gammav"}}, /* 22 from thermostat.cpp */
+     {&npt_iso.gammav, Datafield::Type::DOUBLE, 1,
+      "npt_iso.gammav"}}, /* 22 from thermostat.cpp */
     {FIELD_NPTISO_PEXT,
      {&nptiso.p_ext, Datafield::Type::DOUBLE, 1,
       "npt_p_ext"}}, /* 23 from pressure.cpp */
@@ -152,12 +155,12 @@ const std::unordered_map<int, Datafield> fields{
       "swimming_particles_exist"}}, /* from particle_data.cpp */
 #ifndef PARTICLE_ANISOTROPY
     {FIELD_LANGEVIN_GAMMA_ROTATION,
-     {&langevin_gamma_rotation, Datafield::Type::DOUBLE, 1,
-      "gamma_rot"}}, /* 55 from thermostat.cpp */
+     {&langevin.gamma_rotation, Datafield::Type::DOUBLE, 1,
+      "langevin.gamma_rotation"}}, /* 55 from thermostat.cpp */
 #else
     {FIELD_LANGEVIN_GAMMA_ROTATION,
-     {langevin_gamma_rotation.data(), Datafield::Type::DOUBLE, 3,
-      "gamma_rot"}}, /* 55 from thermostat.cpp */
+     {langevin.gamma_rotation.data(), Datafield::Type::DOUBLE, 3,
+      "langevin.gamma_rotation"}}, /* 55 from thermostat.cpp */
 #endif
 #ifdef OIF_GLOBAL_FORCES
     {FIELD_MAX_OIF_OBJECTS,

src/core/global.hpp

@@ -48,7 +48,7 @@ enum Fields {
   FIELD_BOXL = 0,
   /** index of \ref DomainDecomposition::cell_grid */
   FIELD_CELLGRID,
-  /** index of \ref langevin_gamma */
+  /** index of \ref LangevinThermostat::gamma */
   FIELD_LANGEVIN_GAMMA,
   /** index of \ref integ_switch */
   FIELD_INTEG_SWITCH,
@@ -64,9 +64,9 @@ enum Fields {
   FIELD_RIGIDBONDS,
   /** index of \ref node_grid */
   FIELD_NODEGRID,
-  /** index of \ref nptiso_gamma0 */
+  /** index of \ref IsotropicNptThermostat::gamma0 */
   FIELD_NPTISO_G0,
-  /** index of \ref nptiso_gammav */
+  /** index of \ref IsotropicNptThermostat::gammav */
   FIELD_NPTISO_GV,
   /** index of \ref nptiso_struct::p_ext */
   FIELD_NPTISO_PEXT,
@@ -91,7 +91,7 @@ enum Fields {
   FIELD_LATTICE_SWITCH,
   /** index of \ref min_global_cut */
   FIELD_MIN_GLOBAL_CUT,
-  /** index of \ref langevin_gamma_rotation */
+  /** index of \ref LangevinThermostat::gamma_rotation */
   FIELD_LANGEVIN_GAMMA_ROTATION,
   FIELD_MAX_OIF_OBJECTS, // soft objects as per the object-in-fluid method
   /** index of \ref n_thermalized_bonds */

src/core/integrate.cpp

@@ -58,6 +58,7 @@
 #include "virtual_sites.hpp"
 
 #include "integrators/brownian_inline.hpp"
+#include "integrators/langevin_inline.hpp"
 #include "integrators/steepest_descent.hpp"
 #include "integrators/velocity_verlet_inline.hpp"
 #include "integrators/velocity_verlet_npt.hpp"
@@ -140,6 +141,7 @@ bool integrator_step_1(ParticleRange &particles) {
 
 /** Calls the hook of the propagation kernels after force calculation */
 void integrator_step_2(ParticleRange &particles) {
+  extern BrownianThermostat brownian;
   switch (integ_switch) {
   case INTEG_METHOD_STEEPEST_DESCENT:
     // Nothing
@@ -154,7 +156,7 @@ void integrator_step_2(ParticleRange &particles) {
 #endif
   case INTEG_METHOD_BD:
     // the Ermak-McCammon's Brownian Dynamics requires a single step
-    brownian_dynamics_propagator(particles);
+    brownian_dynamics_propagator(brownian, particles);
     break;
   default:
     throw std::runtime_error("Unknown value for INTEG_SWITCH");