Pass time step and temperature as function arguments

src/core/accumulators/Correlator.hpp

@@ -155,7 +155,7 @@ class Correlator : public AccumulatorBase {
              obs_ptr obs2, Utils::Vector3d correlation_args_ = {})
       : AccumulatorBase(delta_N), finalized(false), t(0),
         m_correlation_args(correlation_args_), m_tau_lin(tau_lin),
-        m_dt(delta_N * time_step), m_tau_max(tau_max),
+        m_dt(delta_N * get_time_step()), m_tau_max(tau_max),
         compressA_name(std::move(compress1_)),
         compressB_name(std::move(compress2_)),
         corr_operation_name(std::move(corr_operation)), A_obs(std::move(obs1)),

src/core/bonded_interactions/thermalized_bond.cpp

@@ -24,8 +24,9 @@
  */
 
 #include "thermalized_bond.hpp"
+
+#include "event.hpp"
 #include "global.hpp"
-#include "integrate.hpp"
 
 #include <cmath>
 
@@ -40,11 +41,12 @@ ThermalizedBond::ThermalizedBond(double temp_com, double gamma_com,
   this->gamma_distance = gamma_distance;
   this->r_cut = r_cut;
 
-  pref1_com = gamma_com;
-  pref2_com = std::sqrt(24.0 * gamma_com / time_step * temp_com);
-  pref1_dist = gamma_distance;
-  pref2_dist = std::sqrt(24.0 * gamma_distance / time_step * temp_distance);
+  pref1_com = -1.;
+  pref2_com = -1.;
+  pref1_dist = -1.;
+  pref2_dist = -1.;
 
   n_thermalized_bonds += 1;
   mpi_bcast_parameter(FIELD_THERMALIZEDBONDS);
+  on_parameter_change(FIELD_THERMALIZEDBONDS);
 }

src/core/bonded_interactions/thermalized_bond_utils.cpp

@@ -27,7 +27,7 @@
 
 #include <boost/variant.hpp>
 
-void thermalized_bond_init() {
+void thermalized_bond_init(double time_step) {
   for (auto &bonded_ia_param : bonded_ia_params) {
     if (auto *t = boost::get<ThermalizedBond>(&bonded_ia_param)) {
       t->pref1_com = t->gamma_com;

src/core/bonded_interactions/thermalized_bond_utils.hpp

@@ -26,6 +26,6 @@
  *
  *  Implementation in \ref thermalized_bond_utils.cpp.
  */
-void thermalized_bond_init();
+void thermalized_bond_init(double time_step);
 
 #endif

src/core/dpd.cpp

@@ -68,26 +68,24 @@ int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
                    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);
 
-  ia_params.dpd_radial = DPDParameters{
-      gamma, k, r_c, wf, sqrt(24.0 * temperature * gamma / time_step)};
-  ia_params.dpd_trans = DPDParameters{
-      tgamma, k, tr_c, twf, sqrt(24.0 * temperature * tgamma / time_step)};
+  ia_params.dpd_radial = DPDParameters{gamma, k, r_c, wf, -1.};
+  ia_params.dpd_trans = DPDParameters{tgamma, k, tr_c, twf, -1.};
 
   /* broadcast interaction parameters */
   mpi_bcast_ia_params(part_type_a, part_type_b);
 
   return ES_OK;
 }
 
-void dpd_init() {
+void dpd_init(double kT, double time_step) {
   for (int type_a = 0; type_a < max_seen_particle_type; type_a++) {
     for (int type_b = 0; type_b < max_seen_particle_type; type_b++) {
       IA_parameters &ia_params = *get_ia_param(type_a, type_b);
 
       ia_params.dpd_radial.pref =
-          sqrt(24.0 * temperature * ia_params.dpd_radial.gamma / time_step);
+          sqrt(24.0 * kT * ia_params.dpd_radial.gamma / time_step);
       ia_params.dpd_trans.pref =
-          sqrt(24.0 * temperature * ia_params.dpd_trans.gamma / time_step);
+          sqrt(24.0 * kT * ia_params.dpd_trans.gamma / time_step);
     }
   }
 }

src/core/dpd.hpp

@@ -37,7 +37,7 @@ struct IA_parameters;
 
 int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
                    double r_c, int wf, double tgamma, double tr_c, int twf);
-void dpd_init();
+void dpd_init(double kT, double time_step);
 
 Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2,
                                IA_parameters const &ia_params,

src/core/energy.cpp

@@ -105,7 +105,7 @@ void energy_calc(const double time) {
   }
 }
 
-void update_energy_local(int, int) { energy_calc(sim_time); }
+void update_energy_local(int, int) { energy_calc(get_sim_time()); }
 
 REGISTER_CALLBACK(update_energy_local)
 

src/core/event.cpp

@@ -100,7 +100,7 @@ void on_program_start() {
   }
 }
 
-void on_integration_start() {
+void on_integration_start(double time_step) {
   /********************************************/
   /* sanity checks                            */
   /********************************************/
@@ -110,20 +110,22 @@ void on_integration_start() {
   integrator_npt_sanity_checks();
 #endif
   interactions_sanity_checks();
-  lb_lbfluid_on_integration_start();
+  lb_lbfluid_sanity_checks(time_step);
 
   /********************************************/
   /* end sanity checks                        */
   /********************************************/
 
+  lb_lbfluid_on_integration_start();
+
 #ifdef CUDA
   MPI_Bcast(gpu_get_global_particle_vars_pointer_host(),
             sizeof(CUDA_global_part_vars), MPI_BYTE, 0, comm_cart);
 #endif
 
   /* Prepare the thermostat */
   if (reinit_thermo) {
-    thermo_init();
+    thermo_init(time_step);
     reinit_thermo = false;
     recalc_forces = true;
   }
@@ -328,6 +330,7 @@ void on_parameter_change(int field) {
   case FIELD_NPTISO_G0:
   case FIELD_NPTISO_GV:
   case FIELD_NPTISO_PISTON:
+  case FIELD_THERMALIZEDBONDS:
     reinit_thermo = true;
     break;
   case FIELD_FORCE_CAP:
@@ -337,7 +340,6 @@ void on_parameter_change(int field) {
   case FIELD_THERMO_SWITCH:
   case FIELD_LATTICE_SWITCH:
   case FIELD_RIGIDBONDS:
-  case FIELD_THERMALIZEDBONDS:
     break;
   case FIELD_SIMTIME:
     recalc_forces = true;

src/core/event.hpp

@@ -54,7 +54,7 @@ void on_program_start();
 /** called every time the simulation is continued/started, i.e.
  *  when switching from the script interface to the simulation core.
  */
-void on_integration_start();
+void on_integration_start(double time_step);
 
 /** called before calculating observables, i.e. energy, pressure or
  *  the integrator (forces). Initialize any methods here which are not

src/core/forces.cpp

@@ -53,7 +53,7 @@
 
 ActorList forceActors;
 
-void init_forces(const ParticleRange &particles, double time_step) {
+void init_forces(const ParticleRange &particles, double time_step, double kT) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
   /* The force initialization depends on the used thermostat and the
      thermodynamic ensemble */
@@ -67,7 +67,7 @@ void init_forces(const ParticleRange &particles, double time_step) {
      set torque to zero for all and rescale quaternions
   */
   for (auto &p : particles) {
-    p.f = init_local_particle_force(p, time_step);
+    p.f = init_local_particle_force(p, time_step, kT);
   }
 
   /* initialize ghost forces with zero
@@ -84,7 +84,7 @@ void init_forces_ghosts(const ParticleRange &particles) {
   }
 }
 
-void force_calc(CellStructure &cell_structure, double time_step) {
+void force_calc(CellStructure &cell_structure, double time_step, double kT) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
 
   espressoSystemInterface.update();
@@ -98,7 +98,7 @@ void force_calc(CellStructure &cell_structure, double time_step) {
 #ifdef ELECTROSTATICS
   icc_iteration(particles, cell_structure.ghost_particles());
 #endif
-  init_forces(particles, time_step);
+  init_forces(particles, time_step, kT);
 
   for (auto &forceActor : forceActors) {
     forceActor->computeForces(espressoSystemInterface);
@@ -134,7 +134,7 @@ void force_calc(CellStructure &cell_structure, double time_step) {
       VerletCriterion{skin, interaction_range(), coulomb_cutoff, dipole_cutoff,
                       collision_detection_cutoff()});
 
-  Constraints::constraints.add_forces(particles, sim_time);
+  Constraints::constraints.add_forces(particles, get_sim_time());
 
   if (max_oif_objects) {
     // There are two global quantities that need to be evaluated:
@@ -155,7 +155,7 @@ void force_calc(CellStructure &cell_structure, double time_step) {
   immersed_boundaries.volume_conservation(cell_structure);
 
   lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual, particles,
-                                         ghost_particles);
+                                         ghost_particles, time_step);
 
 #ifdef CUDA
   copy_forces_from_GPU(particles);

src/core/forces.hpp

@@ -55,7 +55,7 @@ void init_forces_ghosts(const ParticleRange &particles);
  *  <li> Calculate long range interaction forces
  *  </ol>
  */
-void force_calc(CellStructure &cell_structure, double time_step);
+void force_calc(CellStructure &cell_structure, double time_step, double kT);
 
 /** Calculate long range forces (P3M, ...). */
 void calc_long_range_forces(const ParticleRange &particles);

src/core/forces_inline.hpp

@@ -100,27 +100,28 @@ inline ParticleForce external_force(Particle const &p) {
   return f;
 }
 
-inline ParticleForce thermostat_force(Particle const &p, double time_step) {
+inline ParticleForce thermostat_force(Particle const &p, double time_step,
+                                      double kT) {
   extern LangevinThermostat langevin;
   if (!(thermo_switch & THERMO_LANGEVIN)) {
     return {};
   }
 
 #ifdef ROTATION
-  return {friction_thermo_langevin(langevin, p, time_step),
+  return {friction_thermo_langevin(langevin, p, time_step, kT),
           p.p.rotation ? convert_vector_body_to_space(
-                             p, friction_thermo_langevin_rotation(langevin, p,
-                                                                  time_step))
+                             p, friction_thermo_langevin_rotation(
+                                    langevin, p, time_step, kT))
                        : Utils::Vector3d{}};
 #else
-  return friction_thermo_langevin(langevin, p, time_step);
+  return friction_thermo_langevin(langevin, p, time_step, kT);
 #endif
 }
 
 /** Initialize the forces for a real particle */
 inline ParticleForce init_local_particle_force(Particle const &part,
-                                               double time_step) {
-  return thermostat_force(part, time_step) + external_force(part);
+                                               double time_step, double kT) {
+  return thermostat_force(part, time_step, kT) + external_force(part);
 }
 
 inline ParticleForce calc_non_bonded_pair_force(Particle const &p1,

src/core/grid_based_algorithms/electrokinetics_cuda.cu

@@ -33,6 +33,7 @@
 #include "grid_based_algorithms/lb_particle_coupling.hpp"
 #include "grid_based_algorithms/lbgpu.cuh"
 #include "grid_based_algorithms/lbgpu.hpp"
+#include "integrate.hpp"
 
 #include <utils/math/int_pow.hpp>
 #include <utils/math/sqr.hpp>
@@ -2246,24 +2247,24 @@ int ek_init() {
     lbpar_gpu.is_TRT = true;
 
     lb_reinit_parameters_gpu();
-    lbpar_gpu.viscosity = ek_parameters.viscosity * lbpar_gpu.time_step /
-                          Utils::sqr(lbpar_gpu.agrid);
-    lbpar_gpu.bulk_viscosity = ek_parameters.bulk_viscosity *
-                               lbpar_gpu.time_step /
-                               Utils::sqr(lbpar_gpu.agrid);
+    auto const time_step = static_cast<float>(get_time_step());
+    lbpar_gpu.viscosity =
+        ek_parameters.viscosity * time_step / Utils::sqr(lbpar_gpu.agrid);
+    lbpar_gpu.bulk_viscosity =
+        ek_parameters.bulk_viscosity * time_step / Utils::sqr(lbpar_gpu.agrid);
 
     lbpar_gpu.external_force_density =
         ek_parameters.lb_ext_force_density[0] != 0 ||
         ek_parameters.lb_ext_force_density[1] != 0 ||
         ek_parameters.lb_ext_force_density[2] != 0;
     lbpar_gpu.ext_force_density =
         Utils::Vector3f(ek_parameters.lb_ext_force_density) *
-        Utils::sqr(lbpar_gpu.agrid * lbpar_gpu.time_step);
+        Utils::sqr(lbpar_gpu.agrid * time_step);
 
     lb_reinit_parameters_gpu();
     lb_init_gpu();
 
-    ek_parameters.time_step = lbpar_gpu.time_step;
+    ek_parameters.time_step = time_step;
     ek_parameters.dim_x = lbpar_gpu.dim[0];
     ek_parameters.dim_x_padded = (ek_parameters.dim_x / 2 + 1) * 2;
     ek_parameters.dim_y = lbpar_gpu.dim[1];
@@ -3619,7 +3620,6 @@ void ek_print_lbpar() {
   printf("    float gamma_even = %f;\n", lbpar_gpu.gamma_even);
   printf("    float agrid = %f;\n", lbpar_gpu.agrid);
   printf("    float tau = %f;\n", lbpar_gpu.tau);
-  printf("    float time_step = %f;\n", lbpar_gpu.time_step);
   printf("    float bulk_viscosity = %f;\n", lbpar_gpu.bulk_viscosity);
   printf("    unsigned int dim_x = %d;\n", lbpar_gpu.dim[0]);
   printf("    unsigned int dim_y = %d;\n", lbpar_gpu.dim[1]);

src/core/grid_based_algorithms/lb.cpp

@@ -35,7 +35,6 @@
 #include "grid.hpp"
 #include "grid_based_algorithms/lb_boundaries.hpp"
 #include "halo.hpp"
-#include "integrate.hpp"
 #include "lb-d3q19.hpp"
 #include "random.hpp"
 
@@ -187,11 +186,6 @@ std::vector<LB_FluidNode> lbfields;
 
 HaloCommunicator update_halo_comm = HaloCommunicator(0);
 
-/** measures the MD time since the last fluid update */
-static double fluidstep = 0.0;
-
-/********************** The Main LB Part *************************************/
-
 /**
  * @brief Initialize fluid nodes.
  * @param[out] fields         Vector containing the fluid nodes
@@ -901,7 +895,7 @@ void lb_stream(LB_Fluid &lb_fluid, const std::array<T, 19> &populations,
 }
 
 /* Collisions and streaming (push scheme) */
-void lb_collide_stream() {
+void lb_integrate() {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
   /* loop over all lattice cells (halo excluded) */
 #ifdef LB_BOUNDARIES
@@ -977,23 +971,6 @@ void lb_collide_stream() {
 #endif
 }
 
-/** Update the lattice Boltzmann fluid.
- *
- *  This function is called from the integrator. Since the time step
- *  for the lattice dynamics can be coarser than the MD time step, we
- *  monitor the time since the last lattice update.
- */
-void lattice_boltzmann_update() {
-  auto const factor = static_cast<int>(round(lbpar.tau / time_step));
-
-  fluidstep += 1;
-  if (fluidstep >= factor) {
-    fluidstep = 0;
-
-    lb_collide_stream();
-  }
-}
-
 /***********************************************************************/
 /** \name Coupling part */
 /***********************************************************************/

src/core/grid_based_algorithms/lb.hpp

@@ -24,8 +24,7 @@
  *
  *  For performance reasons it is clever to do streaming and collision at the
  *  same time because every fluid node has to be read and written only once.
- *  This increases mainly cache efficiency. This is achieved by
- *  @ref lb_collide_stream.
+ *  This increases mainly cache efficiency.
  *
  *  The hydrodynamic fields, corresponding to density, velocity and pressure,
  *  are stored in @ref LB_FluidNode in the array @ref lbfields, the populations
@@ -177,13 +176,13 @@ extern std::vector<LB_FluidNode> lbfields;
 /************************************************************/
 /**@{*/
 
-/** Update the lattice Boltzmann system for one time step.
+/** Integrate the lattice-Boltzmann system for one time step.
  *  This function performs the collision step and the streaming step.
  *  If external force densities are present, they are applied prior to the
  *  collisions. If boundaries are present, it also applies the boundary
  *  conditions.
  */
-void lattice_boltzmann_update();
+void lb_integrate();
 
 void lb_sanity_checks(const LB_Parameters &lb_parameters);