diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index c7dd4846add74f13f53f293a1cc5ed1b871340c8..363dfe5fac67cfae3455bce30fbac24a317eae71 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -59,11 +59,12 @@ Scheduler:
 
 # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
 TimeIntegration:
-  time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.    # The end time of the simulation (in internal units).
-  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
-
+  time_begin:        0.    # The starting time of the simulation (in internal units).
+  time_end:          1.    # The end time of the simulation (in internal units).
+  dt_min:            1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:            1e-2  # The maximal time-step size of the simulation (in internal units).
+  max_dt_RMS_factor: 0.25  # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
+  
 # Parameters governing the snapshots
 Snapshots:
   basename:   output      # Common part of the name of output files
diff --git a/src/engine.c b/src/engine.c
index a3e643fea695ccc33c6b44159569a9995e4fd5ef..5dafce3ce1af480895eb9fa7d7a44e35bc18b547 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5429,7 +5429,7 @@ void engine_init(struct engine *e, struct space *s,
   e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
   e->dt_max_RMS_displacement = FLT_MAX;
   e->max_RMS_displacement_factor = parser_get_opt_param_double(
-      params, "TimeIntegration:max_dt_RMS_factor", 1.);
+      params, "TimeIntegration:max_dt_RMS_factor", 0.25);
   e->a_first_statistics =
       parser_get_opt_param_double(params, "Statistics:scale_factor_first", 0.1);
   e->time_first_statistics =
@@ -6138,6 +6138,7 @@ void engine_compute_next_statistics_time(struct engine *e) {
  */
 void engine_recompute_displacement_constraint(struct engine *e) {
 
+  /* Get the cosmological information */
   const struct cosmology *cosmo = e->cosmology;
   const float Om = cosmo->Omega_m;
   const float Ob = cosmo->Omega_b;
@@ -6151,6 +6152,16 @@ void engine_recompute_displacement_constraint(struct engine *e) {
                                       FLT_MAX,
                                       e->s->min_spart_mass,
                                       FLT_MAX};
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that the minimal mass collection worked */
+  float min_part_mass_check = FLT_MAX;
+  for (size_t i = 0; i < e->s->nr_parts; ++i)
+    min_part_mass_check =
+        min(min_part_mass_check, hydro_get_mass(&e->s->parts[i]));
+  if (min_part_mass_check != min_mass[swift_type_gas])
+    error("Error collecting minimal mass of gas particles.");
+#endif
+
 #ifdef WITH_MPI
   MPI_Allreduce(MPI_IN_PLACE, min_mass, swift_type_count, MPI_FLOAT, MPI_MIN,
                 MPI_COMM_WORLD);
@@ -6174,14 +6185,6 @@ void engine_recompute_displacement_constraint(struct engine *e) {
   float count_parts[swift_type_count] = {
       e->total_nr_parts, total_nr_dm_gparts, 0.f, 0.f, e->total_nr_sparts, 0.f};
 
-  /* Minimal mass for the two species */
-  const float min_mass_dm = min_mass[1];
-  const float min_mass_b = min(min_mass[0], min_mass[4]);
-
-  /* Inter-particle sepration for the two species */
-  const float d_dm = cbrtf(min_mass_dm / ((Om - Ob) * rho_crit));
-  const float d_b = cbrtf(min_mass_b / (Ob * rho_crit));
-
   /* Count of particles for the two species */
   const float N_dm = count_parts[1];
   const float N_b = count_parts[0] + count_parts[4];
@@ -6190,13 +6193,43 @@ void engine_recompute_displacement_constraint(struct engine *e) {
   const float vel_norm_dm = vel_norm[1];
   const float vel_norm_b = vel_norm[0] + vel_norm[4];
 
-  /* RMS peculiar motion for the two species */
-  const float rms_vel_dm = (N_dm > 0.f) ? (vel_norm_dm / N_dm) : 0.f;
-  const float rms_vel_b = (N_b > 0.f) ? (vel_norm_b / N_b) : 0.f;
+  /* Mesh forces smoothing scale */
+  const float a_smooth =
+      e->gravity_properties->a_smooth * e->s->dim[0] / e->s->cdim[0];
+
+  float dt_dm = FLT_MAX, dt_b = FLT_MAX;
+
+  /* DM case */
+  if (N_dm > 0.f) {
+
+    /* Minimal mass for the DM */
+    const float min_mass_dm = min_mass[1];
+
+    /* Inter-particle sepration for the DM */
+    const float d_dm = cbrtf(min_mass_dm / ((Om - Ob) * rho_crit));
 
-  /* Time-step based on maximum displacement */
-  const float dt_dm = a * a * d_dm / (sqrtf(rms_vel_dm) + FLT_MIN);
-  const float dt_b = a * a * d_b / (sqrtf(rms_vel_b) + FLT_MIN);
+    /* RMS peculiar motion for the DM */
+    const float rms_vel_dm = vel_norm_dm / N_dm;
+
+    /* Time-step based on maximum displacement */
+    dt_dm = a * a * min(a_smooth, d_dm) / sqrtf(rms_vel_dm);
+  }
+
+  /* Baryon case */
+  if (N_b > 0.f) {
+
+    /* Minimal mass for the baryons */
+    const float min_mass_b = min(min_mass[0], min_mass[4]);
+
+    /* Inter-particle sepration for the baryons */
+    const float d_b = cbrtf(min_mass_b / (Ob * rho_crit));
+
+    /* RMS peculiar motion for the baryons */
+    const float rms_vel_b = vel_norm_b / N_b;
+
+    /* Time-step based on maximum displacement */
+    dt_b = a * a * min(a_smooth, d_b) / sqrtf(rms_vel_b);
+  }
 
   /* Use the minimum */
   const float dt = min(dt_dm, dt_b);