diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp
index 037af71990..a7ca71e080 100644
--- a/src/REPLICA/fix_pimd_langevin.cpp
+++ b/src/REPLICA/fix_pimd_langevin.cpp
@@ -113,7 +113,8 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
 
   int seed = -1;
 
-  if (domain->dimension != 3) error->universe_all(FLERR, "Fix pimd/langevin requires a 3d system");
+  if (domain->dimension != 3)
+    error->universe_all(FLERR, fmt::format("Fix {} requires a 3d system", style));
 
   for (int i = 0; i < 6; i++) {
     p_flag[i] = 0;
@@ -127,7 +128,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[i + 1], "pimd") == 0)
         method = PIMD;
       else
-        error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin");
+        error->universe_all(FLERR, fmt::format("Unknown method parameter for fix {}", style));
     } else if (strcmp(arg[i], "integrator") == 0) {
       if (strcmp(arg[i + 1], "obabo") == 0)
         integrator = OBABO;
@@ -135,8 +136,9 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
         integrator = BAOAB;
       else
         error->universe_all(FLERR,
-                            "Unknown integrator parameter for fix pimd/langevin. Only obabo and "
-                            "baoab integrators are supported!");
+                            fmt::format("Unknown integrator parameter for fix {}. Only obabo and "
+                                        "baoab integrators are supported!",
+                                        style));
     } else if (strcmp(arg[i], "ensemble") == 0) {
       if (strcmp(arg[i + 1], "nve") == 0) {
         ensemble = NVE;
@@ -156,15 +158,16 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
         pstat_flag = 1;
       } else
         error->universe_all(FLERR,
-                            "Unknown ensemble parameter for fix pimd/langevin. Only nve, nvt, nph, "
-                            "and npt ensembles are supported!");
+                            fmt::format("Unknown ensemble parameter for fix {}. Only nve, nvt, "
+                                        "nph, and npt ensembles are supported!",
+                                        style));
     } else if (strcmp(arg[i], "fmass") == 0) {
       fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
       if (fmass < 0.0 || fmass > np)
-        error->universe_all(FLERR, "Invalid fmass value for fix pimd/langevin");
+        error->universe_all(FLERR, fmt::format("Invalid fmass value for fix {}", style));
     } else if (strcmp(arg[i], "sp") == 0) {
       sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
-      if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/langevin");
+      if (sp < 0.0) error->universe_all(FLERR, fmt::format("Invalid sp value for fix {}", style));
     } else if (strcmp(arg[i], "fmmode") == 0) {
       if (strcmp(arg[i + 1], "physical") == 0)
         fmmode = PHYSICAL;
@@ -172,8 +175,9 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
         fmmode = NORMAL;
       else
         error->universe_all(FLERR,
-                            "Unknown fictitious mass mode for fix pimd/langevin. Only physical "
-                            "mass and normal mode mass are supported!");
+                            fmt::format("Unknown fictitious mass mode for fix {}. Only physical "
+                                        "mass and normal mode mass are supported!",
+                                        style));
     } else if (strcmp(arg[i], "scale") == 0) {
       if (method == PIMD)
         error->universe_all(
@@ -182,10 +186,11 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
             "scale parameter if you do want to use method pimd.");
       pilescale = utils::numeric(FLERR, arg[i + 1], false, lmp);
       if (pilescale < 0.0)
-        error->universe_all(FLERR, "Invalid PILE_L scale value for fix pimd/langevin");
+        error->universe_all(FLERR, fmt::format("Invalid PILE_L scale value for fix {}", style));
     } else if (strcmp(arg[i], "temp") == 0) {
       temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
-      if (temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd/langevin");
+      if (temp < 0.0)
+        error->universe_all(FLERR, fmt::format("Invalid temp value for fix {}", style));
     } else if (strcmp(arg[i], "lj") == 0) {
       lj_epsilon = utils::numeric(FLERR, arg[i + 1], false, lmp);
       lj_sigma = utils::numeric(FLERR, arg[i + 2], false, lmp);
@@ -207,7 +212,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
       } else if (strcmp(arg[i + 1], "BZP") == 0) {
         barostat = BZP;
       } else
-        error->universe_all(FLERR, "Unknown barostat parameter for fix pimd/langevin");
+        error->universe_all(FLERR, fmt::format("Unknown barostat parameter for fix {}", style));
     } else if (strcmp(arg[i], "iso") == 0) {
       pstyle = ISO;
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
@@ -237,7 +242,8 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
       pdim++;
     } else if (strcmp(arg[i], "taup") == 0) {
       tau_p = utils::numeric(FLERR, arg[i + 1], false, lmp);
-      if (tau_p <= 0.0) error->universe_all(FLERR, "Invalid tau_p value for fix pimd/langevin");
+      if (tau_p <= 0.0)
+        error->universe_all(FLERR, fmt::format("Invalid tau_p value for fix {}", style));
     } else if (strcmp(arg[i], "fixcom") == 0) {
       if (strcmp(arg[i + 1], "yes") == 0)
         removecomflag = 1;
@@ -403,6 +409,7 @@ FixPIMDLangevin::~FixPIMDLangevin()
   memory->destroy(tagsend);
   memory->destroy(tagrecv);
   memory->destroy(bufbeads);
+  MPI_Comm_free(&rootworld);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -422,10 +429,10 @@ int FixPIMDLangevin::setmask()
 void FixPIMDLangevin::init()
 {
   if (atom->map_style == Atom::MAP_NONE)
-    error->all(FLERR, "Fix pimd/langevin requires an atom map, see atom_modify");
+    error->all(FLERR, fmt::format("Fix {} requires an atom map, see atom_modify", style));
 
   if (universe->me == 0 && universe->uscreen)
-    fprintf(universe->uscreen, "Fix pimd/langevin: initializing Path-Integral ...\n");
+    utils::print(universe->uscreen, "Fix {}: initializing Path-Integral ...\n", style);
 
   // prepare the constants
 
@@ -448,8 +455,8 @@ void FixPIMDLangevin::init()
   fbond = _fbond * force->mvv2e;
 
   if ((universe->me == 0) && (universe->uscreen))
-    fprintf(universe->uscreen,
-            "Fix pimd/langevin: -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond);
+    utils::print(universe->uscreen, "Fix {}: -P/(beta^2 * hbar^2) = {:20.7e} (kcal/mol/A^2)\n\n",
+                 style, fbond);
 
   if (integrator == OBABO) {
     dtf = 0.5 * update->dt * force->ftm2v;
@@ -462,7 +469,7 @@ void FixPIMDLangevin::init()
     dtv2 = dtv * dtv;
     dtv3 = THIRD * dtv2 * dtv * force->ftm2v;
   } else {
-    error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin");
+    error->universe_all(FLERR, fmt::format("Unknown integrator parameter for fix {}", style));
   }
 
   comm_init();
@@ -525,7 +532,8 @@ void FixPIMDLangevin::setup(int vflag)
   } else {
     error->universe_all(
         FLERR,
-        "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+        fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!",
+                    style));
   }
   collect_xc();
   compute_spring_energy();
@@ -587,7 +595,8 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/)
     } else {
       error->universe_all(
           FLERR,
-          "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+          fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!",
+                      style));
     }
   } else if (integrator == BAOAB) {
     if (pstat_flag) {
@@ -609,7 +618,8 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/)
     } else {
       error->universe_all(
           FLERR,
-          "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+          fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!",
+                      style));
     }
     if (tstat_flag) {
       o_step();
@@ -624,12 +634,14 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/)
     } else {
       error->universe_all(
           FLERR,
-          "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+          fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!",
+                      style));
     }
   } else {
     error->universe_all(FLERR,
-                        "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab "
-                        "integrators are supported!");
+                        fmt::format("Unknown integrator parameter for fix {}. Only obabo and baoab "
+                                    "integrators are supported!",
+                                    style));
   }
   collect_xc();
 
@@ -671,7 +683,7 @@ void FixPIMDLangevin::final_integrate()
   } else if (integrator == BAOAB) {
 
   } else {
-    error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin");
+    error->universe_all(FLERR, fmt::format("Unknown integrator parameter for fix {}", style));
   }
 }
 
@@ -1045,8 +1057,9 @@ void FixPIMDLangevin::langevin_init()
     c1 = exp(-gamma * update->dt);
   else
     error->universe_all(FLERR,
-                        "Unknown integrator parameter for fix pimd/langevin. Only obabo and "
-                        "baoab integrators are supported!");
+                        fmt::format("Unknown integrator parameter for fix {}. Only obabo and baoab "
+                                    "integrators are supported!",
+                                    style));
 
   c2 = sqrt(1.0 - c1 * c1);    // note that c1 and c2 here only works for the centroid mode.
 
@@ -1068,8 +1081,9 @@ void FixPIMDLangevin::langevin_init()
           c1_k[i] = exp(-1.0 * update->dt / tau_k[i]);
         else
           error->universe_all(FLERR,
-                              "Unknown integrator parameter for fix pimd/langevin. Only obabo and "
-                              "baoab integrators are supported!");
+                              fmt::format("Unknown integrator parameter for fix {}. Only obabo and "
+                                          "baoab integrators are supported!",
+                                          style));
         c2_k[i] = sqrt(1.0 - c1_k[i] * c1_k[i]);
       }
       for (int i = 0; i < np; i++) {
@@ -1420,7 +1434,9 @@ void FixPIMDLangevin::remove_com_motion()
       }
     }
   } else {
-    error->all(FLERR, "Unknown method for fix pimd/langevin. Only nmpimd and pimd are supported!");
+    error->all(
+        FLERR,
+        fmt::format("Unknown method for fix {}. Only nmpimd and pimd are supported!", style));
   }
 }
 
@@ -1551,7 +1567,8 @@ void FixPIMDLangevin::compute_spring_energy()
   } else {
     error->universe_all(
         FLERR,
-        "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+        fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!",
+                     style));
   }
 }
 
@@ -1613,9 +1630,7 @@ void FixPIMDLangevin::compute_p_cv()
   } else if (method == PIMD) {
     p_cv = THIRD * inv_volume * ((2.0 * totke / np - centroid_vir) * force->nktv2p + vir) / np;
   } else {
-    error->universe_all(
-        FLERR,
-        "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!");
+    error->universe_all(FLERR, fmt::format("Unknown method parameter for fix {}. Only nmpimd and pimd are supported!", style));
   }
 }
 
diff --git a/src/REPLICA/fix_pimd_nvt.cpp b/src/REPLICA/fix_pimd_nvt.cpp
index 09c2258c17..d11bc34337 100644
--- a/src/REPLICA/fix_pimd_nvt.cpp
+++ b/src/REPLICA/fix_pimd_nvt.cpp
@@ -89,26 +89,28 @@ FixPIMDNVT::FixPIMDNVT(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
         method = CMD;
       else
         error->universe_all(
-            FLERR, fmt::format("Unknown method parameter {} for fix pimd/nvt", arg[i + 1]));
+            FLERR, fmt::format("Unknown method parameter {} for fix {}", arg[i + 1], style));
     } else if (strcmp(arg[i], "fmass") == 0) {
       fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
       if ((fmass < 0.0) || (fmass > np))
-        error->universe_all(FLERR, fmt::format("Invalid fmass value {} for fix pimd/nvt", fmass));
+        error->universe_all(FLERR, fmt::format("Invalid fmass value {} for fix {}", fmass, style));
     } else if (strcmp(arg[i], "sp") == 0) {
       sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
-      if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/nvt");
+      if (sp < 0.0) error->universe_all(FLERR, fmt::format("Invalid sp value for fix {}", style));
     } else if (strcmp(arg[i], "temp") == 0) {
       nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
-      if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd/nvt");
+      if (nhc_temp < 0.0)
+        error->universe_all(FLERR, fmt::format("Invalid temp value for fix {}", style));
     } else if (strcmp(arg[i], "nhc") == 0) {
       nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp);
-      if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd/nvt");
+      if (nhc_nchain < 2)
+        error->universe_all(FLERR, fmt::format("Invalid nhc value for fix {}", style));
     } else
-      error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd/nvt", arg[i]));
+      error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix {}", arg[i], style));
   }
 
   if (strcmp(update->unit_style, "lj") == 0)
-    error->all(FLERR, "Fix pimd/nvt does not support lj units");
+    error->all(FLERR, fmt::format("Fix {} does not support lj units", style));
 
   /* Initiation */
 
@@ -187,10 +189,10 @@ int FixPIMDNVT::setmask()
 void FixPIMDNVT::init()
 {
   if (atom->map_style == Atom::MAP_NONE)
-    error->universe_all(FLERR, "Fix pimd/nvt requires an atom map, see atom_modify");
+    error->universe_all(FLERR, fmt::format("Fix {} requires an atom map, see atom_modify", style));
 
   if (universe->me == 0 && universe->uscreen)
-    fprintf(universe->uscreen, "Fix pimd/nvt initializing Path-Integral ...\n");
+    utils::print(universe->uscreen, "Fix {} initializing Path-Integral ...\n", style);
 
   // prepare the constants
 
@@ -539,13 +541,15 @@ void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector)
 
 /* ---------------------------------------------------------------------- */
 
-void FixPIMDNVT::pre_spring_force_estimators(){
+void FixPIMDNVT::pre_spring_force_estimators()
+{
   vir_estimator();
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixPIMDNVT::vir_estimator() {
+void FixPIMDNVT::vir_estimator()
+{
   double **x = atom->x;
   double **f = atom->f;
   int nlocal = atom->nlocal;