diff --git a/src/COMPRESS/dump_cfg_gz.cpp b/src/COMPRESS/dump_cfg_gz.cpp
index dcf85e1fc5..3b2f10ddb0 100644
--- a/src/COMPRESS/dump_cfg_gz.cpp
+++ b/src/COMPRESS/dump_cfg_gz.cpp
@@ -91,22 +91,20 @@ void DumpCFGGZ::write_header(bigint n)
   //   so molecules are not split across periodic box boundaries
 
   double scale = 1.0;
-  if (atom->peri_flag)
-    scale = atom->pdscale;
-  else if (unwrapflag == 1)
-    scale = UNWRAPEXPAND;
+  if (atom->peri_flag) scale = atom->pdscale;
+  else if (unwrapflag == 1) scale = UNWRAPEXPAND;
 
   std::string header = fmt::format("Number of particles = {}\n", n);
-  header += fmt::format("A = {0:g} Angstrom (basic length-scale)\n", scale);
-  header += fmt::format("H0(1,1) = {0:g} A\n", domain->xprd);
-  header += fmt::format("H0(1,2) = 0 A \n");
-  header += fmt::format("H0(1,3) = 0 A \n");
-  header += fmt::format("H0(2,1) = {0:g} A \n", domain->xy);
-  header += fmt::format("H0(2,2) = {0:g} A\n", domain->yprd);
-  header += fmt::format("H0(2,3) = 0 A \n");
-  header += fmt::format("H0(3,1) = {0:g} A \n", domain->xz);
-  header += fmt::format("H0(3,2) = {0:g} A \n", domain->yz);
-  header += fmt::format("H0(3,3) = {0:g} A\n", domain->zprd);
+  header += fmt::format("A = {:g} Angstrom (basic length-scale)\n", scale);
+  header += fmt::format("H0(1,1) = {:g} A\n", domain->xprd);
+  header += fmt::format("H0(1,2) = 0 A\n");
+  header += fmt::format("H0(1,3) = 0 A\n");
+  header += fmt::format("H0(2,1) = {:g} A\n", domain->xy);
+  header += fmt::format("H0(2,2) = {:g} A\n", domain->yprd);
+  header += fmt::format("H0(2,3) = 0 A\n");
+  header += fmt::format("H0(3,1) = {:g} A\n", domain->xz);
+  header += fmt::format("H0(3,2) = {:g} A\n", domain->yz);
+  header += fmt::format("H0(3,3) = {:g} A\n", domain->zprd);
   header += fmt::format(".NO_VELOCITY.\n");
   header += fmt::format("entry_count = {}\n", nfield - 2);
   for (int i = 0; i < nfield - 5; i++) header += fmt::format("auxiliary[{}] = {}\n", i, auxname[i]);
diff --git a/src/COMPRESS/dump_cfg_zstd.cpp b/src/COMPRESS/dump_cfg_zstd.cpp
index 8a3dc32e8d..865a8ccf87 100644
--- a/src/COMPRESS/dump_cfg_zstd.cpp
+++ b/src/COMPRESS/dump_cfg_zstd.cpp
@@ -105,16 +105,16 @@ void DumpCFGZstd::write_header(bigint n)
     scale = UNWRAPEXPAND;
 
   std::string header = fmt::format("Number of particles = {}\n", n);
-  header += fmt::format("A = {0:g} Angstrom (basic length-scale)\n", scale);
-  header += fmt::format("H0(1,1) = {0:g} A\n", domain->xprd);
-  header += fmt::format("H0(1,2) = 0 A \n");
-  header += fmt::format("H0(1,3) = 0 A \n");
-  header += fmt::format("H0(2,1) = {0:g} A \n", domain->xy);
-  header += fmt::format("H0(2,2) = {0:g} A\n", domain->yprd);
-  header += fmt::format("H0(2,3) = 0 A \n");
-  header += fmt::format("H0(3,1) = {0:g} A \n", domain->xz);
-  header += fmt::format("H0(3,2) = {0:g} A \n", domain->yz);
-  header += fmt::format("H0(3,3) = {0:g} A\n", domain->zprd);
+  header += fmt::format("A = {:g} Angstrom (basic length-scale)\n", scale);
+  header += fmt::format("H0(1,1) = {:g} A\n", domain->xprd);
+  header += fmt::format("H0(1,2) = 0 A\n");
+  header += fmt::format("H0(1,3) = 0 A\n");
+  header += fmt::format("H0(2,1) = {:g} A\n", domain->xy);
+  header += fmt::format("H0(2,2) = {:g} A\n", domain->yprd);
+  header += fmt::format("H0(2,3) = 0 A\n");
+  header += fmt::format("H0(3,1) = {:g} A\n", domain->xz);
+  header += fmt::format("H0(3,2) = {:g} A\n", domain->yz);
+  header += fmt::format("H0(3,3) = {:g} A\n", domain->zprd);
   header += fmt::format(".NO_VELOCITY.\n");
   header += fmt::format("entry_count = {}\n", nfield - 2);
   for (int i = 0; i < nfield - 5; i++) header += fmt::format("auxiliary[{}] = {}\n", i, auxname[i]);
diff --git a/src/COMPRESS/dump_xyz_gz.cpp b/src/COMPRESS/dump_xyz_gz.cpp
index 9bfead1c5c..4b544c0d67 100644
--- a/src/COMPRESS/dump_xyz_gz.cpp
+++ b/src/COMPRESS/dump_xyz_gz.cpp
@@ -81,11 +81,9 @@ void DumpXYZGZ::openfile()
 void DumpXYZGZ::write_header(bigint ndump)
 {
   if (me == 0) {
-    auto header = fmt::format("{}\n", ndump);
-    if (time_flag) {
-      double tcurrent = update->atime + (update->ntimestep-update->atimestep) + update->dt;
-      header += fmt::format(" Atoms. Timestep: {} Time: {:.6f}\n", update->ntimestep, tcurrent);
-    } else header += fmt::format(" Atoms. Timestep: {}\n", update->ntimestep);
+    auto header = fmt::format("{}\n Atoms. Timestep: {}", ndump, update->ntimestep);
+    if (time_flag) header += fmt::format(" Time: {:.6f}", compute_time());
+    header += "\n";
     writer.write(header.c_str(), header.length());
   }
 }
diff --git a/src/COMPRESS/dump_xyz_zstd.cpp b/src/COMPRESS/dump_xyz_zstd.cpp
index bb2640b1b8..ed87bb3f2c 100644
--- a/src/COMPRESS/dump_xyz_zstd.cpp
+++ b/src/COMPRESS/dump_xyz_zstd.cpp
@@ -89,11 +89,9 @@ void DumpXYZZstd::openfile()
 void DumpXYZZstd::write_header(bigint ndump)
 {
   if (me == 0) {
-    auto header = fmt::format("{}\n", ndump);
-    if (time_flag) {
-      double tcurrent = update->atime + (update->ntimestep-update->atimestep) + update->dt;
-      header += fmt::format(" Atoms. Timestep: {} Time: {:.6f}\n", update->ntimestep, tcurrent);
-    } else header += fmt::format(" Atoms. Timestep: {}\n", update->ntimestep);
+    auto header = fmt::format("{}\n Atoms. Timestep: {}", ndump, update->ntimestep);
+    if (time_flag) header += fmt::format(" Time: {:.6f}", compute_time());
+    header += "\n";
     writer.write(header.c_str(), header.length());
   }
 }
diff --git a/src/MPIIO/dump_atom_mpiio.cpp b/src/MPIIO/dump_atom_mpiio.cpp
index 9ba779924f..aedb4005d2 100644
--- a/src/MPIIO/dump_atom_mpiio.cpp
+++ b/src/MPIIO/dump_atom_mpiio.cpp
@@ -212,20 +212,10 @@ void DumpAtomMPIIO::init_style()
 
   delete[] format;
   if (format_line_user) {
-    int n = strlen(format_line_user) + 2;
-    format = new char[n];
-    strcpy(format, format_line_user);
-    strcat(format, "\n");
+    format = utils::strdup(std::string(format_line_user) + "\n");
   } else {
-    char *str;
-    if (image_flag == 0)
-      str = (char *) TAGINT_FORMAT " %d %g %g %g";
-    else
-      str = (char *) TAGINT_FORMAT " %d %g %g %g %d %d %d";
-    int n = strlen(str) + 2;
-    format = new char[n];
-    strcpy(format, str);
-    strcat(format, "\n");
+    if (image_flag == 0) format = utils::strdup(TAGINT_FORMAT " %d %g %g %g\n");
+    else format = utils::strdup(TAGINT_FORMAT " %d %g %g %g %d %d %d\n");
   }
 
   // setup boundary string
diff --git a/src/MPIIO/dump_cfg_mpiio.cpp b/src/MPIIO/dump_cfg_mpiio.cpp
index 88072d5adc..4f601eb865 100644
--- a/src/MPIIO/dump_cfg_mpiio.cpp
+++ b/src/MPIIO/dump_cfg_mpiio.cpp
@@ -271,41 +271,32 @@ void DumpCFGMPIIO::write_header(bigint n)
   // for unwrapped coords, set to UNWRAPEXPAND (10.0)
   //   so molecules are not split across periodic box boundaries
 
+  double scale = 1.0;
+  if (atom->peri_flag) scale = atom->pdscale;
+  else if (unwrapflag == 1) scale = UNWRAPEXPAND;
+
+  auto header = fmt::format("Number of particles = {}\n",n);
+  header += fmt::format("A = {} Angstrom (basic length-scale)\n",scale);
+  header += fmt::format("H0(1,1) = {} A\n",domain->xprd);
+  header += fmt::format("H0(1,2) = 0 A\n");
+  header += fmt::format("H0(1,3) = 0 A\n");
+  header += fmt::format("H0(2,1) = {} A\n",domain->xy);
+  header += fmt::format("H0(2,2) = {} A\n",domain->yprd);
+  header += fmt::format("H0(2,3) = 0 A\n");
+  header += fmt::format("H0(3,1) = {} A\n",domain->xz);
+  header += fmt::format("H0(3,2) = {} A\n",domain->yz);
+  header += fmt::format("H0(3,3) = {} A\n",domain->zprd);
+  header += fmt::format(".NO_VELOCITY.\n");
+  header += fmt::format("entry_count = {}\n",nfield-2);
+  for (int i = 0; i < nfield-5; i++)
+    header += fmt::format("auxiliary[{}] = {}\n",i,auxname[i]);
+
   if (performEstimate) {
-
-    headerBuffer = (char *) malloc(MAX_TEXT_HEADER_SIZE);
-
-    headerSize = 0;
-
-    double scale = 1.0;
-    if (atom->peri_flag) scale = atom->pdscale;
-    else if (unwrapflag == 1) scale = UNWRAPEXPAND;
-
-    char str[64];
-
-    sprintf(str,"Number of particles = %s\n",BIGINT_FORMAT);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),str,n);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"A = %g Angstrom (basic length-scale)\n",scale);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,1) = %g A\n",domain->xprd);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,2) = 0 A \n");
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,3) = 0 A \n");
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,1) = %g A \n",domain->xy);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,2) = %g A\n",domain->yprd);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,3) = 0 A \n");
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,1) = %g A \n",domain->xz);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,2) = %g A \n",domain->yz);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,3) = %g A\n",domain->zprd);
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),".NO_VELOCITY.\n");
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"entry_count = %d\n",nfield-2);
-    for (int i = 0; i < nfield-5; i++)
-      headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"auxiliary[%d] = %s\n",i,auxname[i]);
-  }
-  else { // write data
-
+    headerSize = header.size();
+  } else { // write data
     if (me == 0)
-      MPI_File_write_at(mpifh,mpifo,headerBuffer,headerSize,MPI_CHAR,MPI_STATUS_IGNORE);
-    mpifo += headerSize;
-    free(headerBuffer);
+      MPI_File_write_at(mpifh,mpifo,header.c_str(),header.size(),MPI_CHAR,MPI_STATUS_IGNORE);
+    mpifo += header.size();
   }
 }
 
diff --git a/src/MPIIO/dump_xyz_mpiio.cpp b/src/MPIIO/dump_xyz_mpiio.cpp
index c976932b52..e829ff7a30 100644
--- a/src/MPIIO/dump_xyz_mpiio.cpp
+++ b/src/MPIIO/dump_xyz_mpiio.cpp
@@ -100,12 +100,11 @@ void DumpXYZMPIIO::openfile()
   if (append_flag) { // append open
     int err = MPI_File_open( world, filecurrent, MPI_MODE_CREATE | MPI_MODE_APPEND |
                              MPI_MODE_WRONLY  , MPI_INFO_NULL, &mpifh);
-    if (err != MPI_SUCCESS) error->one(FLERR, "Cannot open dump file {}",filecurrent);
+    if (err != MPI_SUCCESS) error->one(FLERR, "Cannot open dump file {}", filecurrent);
 
     int myrank;
     MPI_Comm_rank(world,&myrank);
-    if (myrank == 0)
-      MPI_File_get_size(mpifh,&mpifo);
+    if (myrank == 0) MPI_File_get_size(mpifh,&mpifo);
     MPI_Bcast(&mpifo, 1, MPI_LMP_BIGINT, 0, world);
     MPI_File_set_size(mpifh,mpifo+headerSize+sumFileSize);
     currentFileSize = mpifo+headerSize+sumFileSize;
@@ -238,19 +237,16 @@ void DumpXYZMPIIO::init_style()
 
 void DumpXYZMPIIO::write_header(bigint n)
 {
+  auto header = fmt::format("{}\n Atoms. Timestep: {}", n, update->ntimestep);
+  if (time_flag) header += fmt::format(" Time: {:.6f}", compute_time());
+  header += "\n";
+
   if (performEstimate) {
-
-    headerBuffer = (char *) malloc(MAX_TEXT_HEADER_SIZE);
-
-    headerSize = 0;
-    headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),BIGINT_FORMAT "\n",n);
-    headerSize += sprintf(&((char*)headerBuffer)[headerSize],"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep);
+    headerSize = header.size();
   } else { // write data
-
     if (me == 0)
-      MPI_File_write_at(mpifh,mpifo,headerBuffer,headerSize,MPI_CHAR,MPI_STATUS_IGNORE);
-    mpifo += headerSize;
-    free(headerBuffer);
+      MPI_File_write_at(mpifh,mpifo,header.c_str(),header.size(),MPI_CHAR,MPI_STATUS_IGNORE);
+    mpifo += header.size();
   }
 }
 
diff --git a/src/OPENMP/fix_qeq_comb_omp.cpp b/src/OPENMP/fix_qeq_comb_omp.cpp
index 82d0b0665a..18f579ad3a 100644
--- a/src/OPENMP/fix_qeq_comb_omp.cpp
+++ b/src/OPENMP/fix_qeq_comb_omp.cpp
@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -36,45 +35,41 @@ using namespace FixConst;
 
 /* ---------------------------------------------------------------------- */
 
-FixQEQCombOMP::FixQEQCombOMP(LAMMPS *lmp, int narg, char **arg) :
-  FixQEQComb(lmp, narg, arg)
+FixQEQCombOMP::FixQEQCombOMP(LAMMPS *lmp, int narg, char **arg) : FixQEQComb(lmp, narg, arg)
 {
-  if (narg < 5) error->all(FLERR,"Illegal fix qeq/comb/omp command");
+  if (narg < 5) error->all(FLERR, "Illegal fix qeq/comb/omp command");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixQEQCombOMP::init()
 {
-  if (!atom->q_flag)
-    error->all(FLERR,"Fix qeq/comb/omp requires atom attribute q");
+  if (!atom->q_flag) error->all(FLERR, "Fix qeq/comb/omp requires atom attribute q");
 
-  if (nullptr != force->pair_match("comb3",0))
-    error->all(FLERR,"No support for comb3 currently available in OPENMP");
+  if (nullptr != force->pair_match("comb3", 0))
+    error->all(FLERR, "No support for comb3 currently available in OPENMP");
 
-  comb = (PairComb *) force->pair_match("comb/omp",1);
+  comb = (PairComb *) force->pair_match("comb/omp", 1);
+  if (comb == nullptr) comb = (PairComb *) force->pair_match("comb", 1);
   if (comb == nullptr)
-    comb = (PairComb *) force->pair_match("comb",1);
-  if (comb == nullptr)
-    error->all(FLERR,"Must use pair_style comb or "
-               "comb/omp with fix qeq/comb/omp");
+    error->all(FLERR, "Must use pair_style comb or comb/omp with fix qeq/comb/omp");
 
-  if (utils::strmatch(update->integrate_style,"^respa")) {
-    ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
-    if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+  if (utils::strmatch(update->integrate_style, "^respa")) {
+    ilevel_respa = ((Respa *) update->integrate)->nlevels - 1;
+    if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
   }
 
   ngroup = group->count(igroup);
-  if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
+  if (ngroup == 0) error->all(FLERR, "Fix qeq/comb group has no atoms");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixQEQCombOMP::post_force(int /* vflag */)
 {
-  int i,ii,iloop,loopmax,inum,*ilist;
-  double heatpq,qmass,dtq,dtq2;
-  double enegchkall,enegmaxall;
+  int i, ii, iloop, loopmax, inum, *ilist;
+  double heatpq, qmass, dtq, dtq2;
+  double enegchkall, enegmaxall;
 
   if (update->ntimestep % nevery) return;
 
@@ -88,28 +83,28 @@ void FixQEQCombOMP::post_force(int /* vflag */)
     memory->destroy(q1);
     memory->destroy(q2);
     nmax = atom->nmax;
-    memory->create(qf,nmax,"qeq:qf");
-    memory->create(q1,nmax,"qeq:q1");
-    memory->create(q2,nmax,"qeq:q2");
+    memory->create(qf, nmax, "qeq:qf");
+    memory->create(q1, nmax, "qeq:q1");
+    memory->create(q2, nmax, "qeq:q2");
     vector_atom = qf;
   }
 
   // more loops for first-time charge equilibrium
 
   iloop = 0;
-  if (firstflag) loopmax = 500;
-  else loopmax = 200;
+  if (firstflag)
+    loopmax = 500;
+  else
+    loopmax = 200;
 
   // charge-equilibration loop
 
-  if (me == 0 && fp)
-    fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
-            update->ntimestep);
+  if (me == 0 && fp) fmt::print(fp, "Charge equilibration on step {}\n", update->ntimestep);
 
   heatpq = 0.05;
-  qmass  = 0.016;
-  dtq    = 0.01;
-  dtq2   = 0.5*dtq*dtq/qmass;
+  qmass = 0.016;
+  dtq = 0.01;
+  dtq2 = 0.5 * dtq * dtq / qmass;
 
   double enegchk = 0.0;
   double enegtot = 0.0;
@@ -126,54 +121,51 @@ void FixQEQCombOMP::post_force(int /* vflag */)
     q1[i] = q2[i] = qf[i] = 0.0;
   }
 
-  for (iloop = 0; iloop < loopmax; iloop ++) {
+  for (iloop = 0; iloop < loopmax; iloop++) {
     for (ii = 0; ii < inum; ii++) {
       i = ilist[ii];
       if (mask[i] & groupbit) {
-        q1[i] += qf[i]*dtq2 - heatpq*q1[i];
-        q[i]  += q1[i];
+        q1[i] += qf[i] * dtq2 - heatpq * q1[i];
+        q[i] += q1[i];
       }
     }
     comm->forward_comm(this);
 
-    if (comb) enegtot = comb->yasu_char(qf,igroup);
+    if (comb) enegtot = comb->yasu_char(qf, igroup);
     enegtot /= ngroup;
     enegchk = enegmax = 0.0;
 
-    for (ii = 0; ii < inum ; ii++) {
+    for (ii = 0; ii < inum; ii++) {
       i = ilist[ii];
       if (mask[i] & groupbit) {
-        q2[i] = enegtot-qf[i];
-        enegmax = MAX(enegmax,fabs(q2[i]));
+        q2[i] = enegtot - qf[i];
+        enegmax = MAX(enegmax, fabs(q2[i]));
         enegchk += fabs(q2[i]);
         qf[i] = q2[i];
       }
     }
 
-    MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
-    enegchk = enegchkall/ngroup;
-    MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
+    MPI_Allreduce(&enegchk, &enegchkall, 1, MPI_DOUBLE, MPI_SUM, world);
+    enegchk = enegchkall / ngroup;
+    MPI_Allreduce(&enegmax, &enegmaxall, 1, MPI_DOUBLE, MPI_MAX, world);
     enegmax = enegmaxall;
 
-    if (enegchk <= precision && enegmax <= 100.0*precision) break;
+    if (enegchk <= precision && enegmax <= 100.0 * precision) break;
 
     if (me == 0 && fp)
-      fprintf(fp,"  iteration: %d, enegtot %.6g, "
-              "enegmax %.6g, fq deviation: %.6g\n",
-              iloop,enegtot,enegmax,enegchk);
+      fprintf(fp, "  iteration: %d, enegtot %.6g, enegmax %.6g, fq deviation: %.6g\n", iloop,
+              enegtot, enegmax, enegchk);
 
     for (ii = 0; ii < inum; ii++) {
       i = ilist[ii];
-      if (mask[i] & groupbit)
-        q1[i] += qf[i]*dtq2 - heatpq*q1[i];
+      if (mask[i] & groupbit) q1[i] += qf[i] * dtq2 - heatpq * q1[i];
     }
   }
 
   if (me == 0 && fp) {
     if (iloop == loopmax)
-      fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
+      fprintf(fp, "Charges did not converge in %d iterations\n", iloop);
     else
-      fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
-              iloop,enegchk);
+      fprintf(fp, "Charges converged in %d iterations to %.10f tolerance\n", iloop, enegchk);
   }
 }
diff --git a/src/PHONON/fix_phonon.cpp b/src/PHONON/fix_phonon.cpp
index 3da204f6b8..b85337e5db 100644
--- a/src/PHONON/fix_phonon.cpp
+++ b/src/PHONON/fix_phonon.cpp
@@ -183,23 +183,22 @@ FixPhonon::FixPhonon(LAMMPS *lmp,  int narg, char **arg) : Fix(lmp, narg, arg)
   if (me == 0) {
     flog = fopen(logfile, "w");
     if (flog == nullptr)
-      error->one(FLERR,"Can not open output file {}: {}",
-                                   logfile,utils::getsyserror());
-    fprintf(flog,"############################################################\n");
-    fprintf(flog,"# group name of the atoms under study      : %s\n", group->names[igroup]);
-    fprintf(flog,"# total number of atoms in the group       : %d\n", ngroup);
-    fprintf(flog,"# dimension of the system                  : %d D\n", sysdim);
-    fprintf(flog,"# number of atoms per unit cell            : %d\n", nucell);
-    fprintf(flog,"# dimension of the FFT mesh                : %d x %d x %d\n", nx, ny, nz);
-    fprintf(flog,"# number of wait steps before measurement  : " BIGINT_FORMAT "\n", waitsteps);
-    fprintf(flog,"# frequency of the measurement             : %d\n", nevery);
-    fprintf(flog,"# output result after this many measurement: %d\n", nfreq);
-    fprintf(flog,"# number of processors used by this run    : %d\n", nprocs);
-    fprintf(flog,"############################################################\n");
-    fprintf(flog,"# mapping information between lattice indices and atom id\n");
-    fprintf(flog,"# nx ny nz nucell\n");
-    fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell);
-    fprintf(flog,"# l1 l2 l3 k atom_id\n");
+      error->one(FLERR,"Can not open output file {}: {}", logfile,utils::getsyserror());
+    fmt::print(flog,"############################################################\n");
+    fmt::print(flog,"# group name of the atoms under study      : {}\n", group->names[igroup]);
+    fmt::print(flog,"# total number of atoms in the group       : {}\n", ngroup);
+    fmt::print(flog,"# dimension of the system                  : {} D\n", sysdim);
+    fmt::print(flog,"# number of atoms per unit cell            : {}\n", nucell);
+    fmt::print(flog,"# dimension of the FFT mesh                : {} x {} x {}\n", nx, ny, nz);
+    fmt::print(flog,"# number of wait steps before measurement  : {}\n", waitsteps);
+    fmt::print(flog,"# frequency of the measurement             : {}\n", nevery);
+    fmt::print(flog,"# output result after this many measurement: {}\n", nfreq);
+    fmt::print(flog,"# number of processors used by this run    : {}\n", nprocs);
+    fmt::print(flog,"############################################################\n");
+    fmt::print(flog,"# mapping information between lattice indices and atom id\n");
+    fmt::print(flog,"# nx ny nz nucell\n");
+    fmt::print(flog,"{} {} {} {}\n", nx, ny, nz, nucell);
+    fmt::print(flog,"# l1 l2 l3 k atom_id\n");
     int ix, iy, iz, iu;
     for (idx = 0; idx < ngroup; ++idx) {
       itag = surf2tag[idx];
@@ -207,9 +206,9 @@ FixPhonon::FixPhonon(LAMMPS *lmp,  int narg, char **arg) : Fix(lmp, narg, arg)
       iz   = (idx/nucell)%nz;
       iy   = (idx/(nucell*nz))%ny;
       ix   = (idx/(nucell*nz*ny))%nx;
-      fprintf(flog,"%d %d %d %d " TAGINT_FORMAT "\n", ix, iy, iz, iu, itag);
+      fmt::print(flog,"{} {} {} {} {}\n", ix, iy, iz, iu, itag);
     }
-    fprintf(flog,"############################################################\n");
+    fmt::print(flog,"############################################################\n");
     fflush(flog);
   }
   surf2tag.clear();
@@ -737,16 +736,16 @@ void FixPhonon::postprocess( )
     fclose(fp_bin);
 
     // write log file, here however, it is the dynamical matrix that is written
-    fprintf(flog,"############################################################\n");
-    fprintf(flog,"# Current time step                      : " BIGINT_FORMAT "\n", update->ntimestep);
-    fprintf(flog,"# Total number of measurements           : %d\n", neval);
-    fprintf(flog,"# Average temperature of the measurement : %lg\n", TempAve);
-    fprintf(flog,"# Boltzmann constant under current units : %lg\n", boltz);
-    fprintf(flog,"# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]);
-    fprintf(flog,"# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]);
-    fprintf(flog,"# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]);
-    fprintf(flog,"############################################################\n");
-    fprintf(flog,"# qx\t qy \t qz \t\t Phi(q)\n");
+    fmt::print(flog,"############################################################\n");
+    fmt::print(flog,"# Current time step                      : {}\n", update->ntimestep);
+    fmt::print(flog,"# Total number of measurements           : {}\n", neval);
+    fmt::print(flog,"# Average temperature of the measurement : {}\n", TempAve);
+    fmt::print(flog,"# Boltzmann constant under current units : {}\n", boltz);
+    fmt::print(flog,"# basis vector A1 = [{} {} {}]\n", basevec[0], basevec[1], basevec[2]);
+    fmt::print(flog,"# basis vector A2 = [{} {} {}]\n", basevec[3], basevec[4], basevec[5]);
+    fmt::print(flog,"# basis vector A3 = [{} {} {}]\n", basevec[6], basevec[7], basevec[8]);
+    fmt::print(flog,"############################################################\n");
+    fmt::print(flog,"# qx\t qy \t qz \t\t Phi(q)\n");
 
     EnforceASR();
 
@@ -754,7 +753,8 @@ void FixPhonon::postprocess( )
     for (idq = 0; idq < ntotal; ++idq) {
       ndim =0;
       for (idim = 0; idim < fft_dim; ++idim)
-      for (jdim = 0; jdim < fft_dim; ++jdim) Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
+        for (jdim = 0; jdim < fft_dim; ++jdim)
+          Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
     }
 
     idq =0;
@@ -764,11 +764,10 @@ void FixPhonon::postprocess( )
         double qy = double(iy)/double(ny);
         for (int iz = 0; iz < nz; ++iz) {
           double qz = double(iz)/double(nz);
-          fprintf(flog,"%lg %lg %lg", qx, qy, qz);
+          fmt::print(flog,"{} {} {}", qx, qy, qz);
           for (idim = 0; idim < fft_dim2; ++idim)
-            fprintf(flog, " %lg %lg", std::real(Phi_all[idq][idim]),
-                                      std::imag(Phi_all[idq][idim]));
-          fprintf(flog, "\n");
+            fmt::print(flog, " {} {}", std::real(Phi_all[idq][idim]), std::imag(Phi_all[idq][idim]));
+          fmt::print(flog, "\n");
           ++idq;
         }
       }
diff --git a/src/QMMM/fix_qmmm.cpp b/src/QMMM/fix_qmmm.cpp
index 8c59077c7d..d0cfe2ed95 100644
--- a/src/QMMM/fix_qmmm.cpp
+++ b/src/QMMM/fix_qmmm.cpp
@@ -443,10 +443,7 @@ void FixQMMM::exchange_positions()
       MPI_Send(isend_buf, 4, MPI_INT,    1, QMMM_TAG_SIZE, qm_comm);
       MPI_Send(celldata,  9, MPI_DOUBLE, 1, QMMM_TAG_CELL, qm_comm);
     }
-    if (verbose > 0) {
-      if (screen) fputs("QMMM: exchange positions\n",screen);
-      if (logfile) fputs("QMMM: exchange positions\n",logfile);
-    }
+    if (verbose > 0) utils::logmesg(lmp, "QMMM: exchange positions\n");
   }
 
   if (qmmm_role == QMMM_ROLE_MASTER) {
@@ -537,10 +534,7 @@ void FixQMMM::exchange_forces()
   const int nlocal = atom->nlocal;
   const int natoms = (int) atom->natoms;
 
-  if ((comm->me) == 0 && (verbose > 0)) {
-    if (screen)  fputs("QMMM: exchange forces\n",screen);
-    if (logfile) fputs("QMMM: exchange forces\n",logfile);
-  }
+  if ((comm->me) == 0 && (verbose > 0)) utils::logmesg(lmp, "QMMM: exchange forces\n");
 
   if (qmmm_role == QMMM_ROLE_MASTER) {
     struct commdata *buf = static_cast<struct commdata *>(comm_buf);
@@ -561,21 +555,14 @@ void FixQMMM::exchange_forces()
       // so we need to apply the scaling factor to get to the
       // supported internal units ("metal" or "real")
       for (int i=0; i < num_qm; ++i) {
-        if  (verbose > 1) {
-           const char fmt[] = "[" TAGINT_FORMAT "]: QM(%g %g %g) MM(%g %g %g) /\\(%g %g %g)\n";
-           if (screen) fprintf(screen, fmt, qm_remap[i],
-                qmmm_fscale*qm_force[3*i+0], qmmm_fscale*qm_force[3*i+1], qmmm_fscale*qm_force[3*i+2],
-                mm_force_on_qm_atoms[3*i+0], mm_force_on_qm_atoms[3*i+1], mm_force_on_qm_atoms[3*i+2],
-                qmmm_fscale*qm_force[3*i+0] - mm_force_on_qm_atoms[3*i+0],
-                qmmm_fscale*qm_force[3*i+1] - mm_force_on_qm_atoms[3*i+1],
-                qmmm_fscale*qm_force[3*i+2] - mm_force_on_qm_atoms[3*i+2]);
-           if (logfile) fprintf(logfile, fmt, qm_remap[i],
-                qmmm_fscale*qm_force[3*i+0], qmmm_fscale*qm_force[3*i+1], qmmm_fscale*qm_force[3*i+2],
-                mm_force_on_qm_atoms[3*i+0], mm_force_on_qm_atoms[3*i+1], mm_force_on_qm_atoms[3*i+2],
-                qmmm_fscale*qm_force[3*i+0] - mm_force_on_qm_atoms[3*i+0],
-                qmmm_fscale*qm_force[3*i+1] - mm_force_on_qm_atoms[3*i+1],
-                qmmm_fscale*qm_force[3*i+2] - mm_force_on_qm_atoms[3*i+2]);
-        }
+        if  (verbose > 1)
+          utils::logmesg(lmp, "[{}]: QM({} {} {}) MM({} {} {}) /\\({} {} {})\n", qm_remap[i],
+                         qmmm_fscale*qm_force[3*i+0], qmmm_fscale*qm_force[3*i+1], qmmm_fscale*qm_force[3*i+2],
+                         mm_force_on_qm_atoms[3*i+0], mm_force_on_qm_atoms[3*i+1], mm_force_on_qm_atoms[3*i+2],
+                         qmmm_fscale*qm_force[3*i+0] - mm_force_on_qm_atoms[3*i+0],
+                         qmmm_fscale*qm_force[3*i+1] - mm_force_on_qm_atoms[3*i+1],
+                         qmmm_fscale*qm_force[3*i+2] - mm_force_on_qm_atoms[3*i+2]);
+
         buf[i].tag = qm_remap[i];
         buf[i].x = qmmm_fscale*qm_force[3*i+0] - mm_force_on_qm_atoms[3*i+0];
         buf[i].y = qmmm_fscale*qm_force[3*i+1] - mm_force_on_qm_atoms[3*i+1];
@@ -684,21 +671,9 @@ void FixQMMM::init()
       memory->create(qm_charge,num_qm,"qmmm:qm_charge");
       memory->create(qm_force,3*num_qm,"qmmm:qm_force");
 
-      const char fmt1[] = "Initializing QM/MM master with %d QM atoms\n";
-      const char fmt2[] = "Initializing QM/MM master with %d MM atoms\n";
-      const char fmt3[] = "Electrostatic coupling with %d atoms\n";
-
-      if (screen)  {
-        fprintf(screen,fmt1,num_qm);
-        fprintf(screen,fmt2,num_mm);
-        if (qmmm_mode == QMMM_MODE_ELEC) fprintf(screen,fmt3,num_mm-num_qm);
-      }
-
-      if (logfile) {
-        fprintf(logfile,fmt1,num_qm);
-        fprintf(logfile,fmt2,num_mm);
-        if (qmmm_mode == QMMM_MODE_ELEC) fprintf(logfile,fmt3,num_mm-num_qm);
-      }
+      utils::logmesg(lmp, "Initializing QM/MM master with {} QM atoms\n", num_qm);
+      utils::logmesg(lmp, "Initializing QM/MM master with {} MM atoms\n", num_mm);
+      utils::logmesg(lmp, "Electrostatic coupling with {} atoms\n", num_mm-num_qm);
 
     } else if (qmmm_role == QMMM_ROLE_SLAVE) {
 
@@ -711,11 +686,7 @@ void FixQMMM::init()
       memory->create(qm_coord,3*num_qm,"qmmm:qm_coord");
       memory->create(qm_force,3*num_qm,"qmmm:qm_force");
 
-      const char fmt[] = "Initializing QM/MM slave with %d QM atoms\n";
-
-      if (screen)  fprintf(screen,fmt,num_qm);
-      if (logfile) fprintf(logfile,fmt,num_qm);
-
+      utils::logmesg(lmp, "Initializing QM/MM slave with {} QM atoms\n",num_qm);
     }
 
     // communication buffer
@@ -769,15 +740,10 @@ void FixQMMM::init()
       qm_remap=taginthash_keys(qm_hash);
 
       if (verbose > 1) {
-        const char fmt[] = "qm_remap[%d]=" TAGINT_FORMAT
-                           "  qm_hash[" TAGINT_FORMAT "]=" TAGINT_FORMAT "\n";
         // print hashtable and reverse mapping
-        for (i=0; i < num_qm; ++i) {
-          if (screen) fprintf(screen,fmt,i,qm_remap[i],qm_remap[i],
-                              taginthash_lookup(qm_hash, qm_remap[i]));
-          if (logfile) fprintf(logfile,fmt,i,qm_remap[i],qm_remap[i],
-                               taginthash_lookup(qm_hash, qm_remap[i]));
-        }
+        for (i=0; i < num_qm; ++i)
+          utils::logmesg(lmp, "qm_remap[{}]={} qm_hash[{}]={}\n",
+                         i,qm_remap[i],qm_remap[i],taginthash_lookup(qm_hash, qm_remap[i]));
       }
 
     } else {
diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp
index b1c9aed58b..303327e797 100644
--- a/src/dump_cfg.cpp
+++ b/src/dump_cfg.cpp
@@ -123,21 +123,22 @@ void DumpCFG::write_header(bigint n)
   if (atom->peri_flag) scale = atom->pdscale;
   else if (unwrapflag == 1) scale = UNWRAPEXPAND;
 
-  fmt::print(fp,"Number of particles = {}\n", n);
-  fprintf(fp,"A = %g Angstrom (basic length-scale)\n",scale);
-  fprintf(fp,"H0(1,1) = %g A\n",domain->xprd);
-  fprintf(fp,"H0(1,2) = 0 A \n");
-  fprintf(fp,"H0(1,3) = 0 A \n");
-  fprintf(fp,"H0(2,1) = %g A \n",domain->xy);
-  fprintf(fp,"H0(2,2) = %g A\n",domain->yprd);
-  fprintf(fp,"H0(2,3) = 0 A \n");
-  fprintf(fp,"H0(3,1) = %g A \n",domain->xz);
-  fprintf(fp,"H0(3,2) = %g A \n",domain->yz);
-  fprintf(fp,"H0(3,3) = %g A\n",domain->zprd);
-  fprintf(fp,".NO_VELOCITY.\n");
-  fprintf(fp,"entry_count = %d\n",nfield-2);
+  auto header = fmt::format("Number of particles = {}\n",n);
+  header += fmt::format("A = {:g} Angstrom (basic length-scale)\n",scale);
+  header += fmt::format("H0(1,1) = {:g} A\n",domain->xprd);
+  header += fmt::format("H0(1,2) = 0 A\n");
+  header += fmt::format("H0(1,3) = 0 A\n");
+  header += fmt::format("H0(2,1) = {:g} A\n",domain->xy);
+  header += fmt::format("H0(2,2) = {:g} A\n",domain->yprd);
+  header += fmt::format("H0(2,3) = 0 A\n");
+  header += fmt::format("H0(3,1) = {:g} A\n",domain->xz);
+  header += fmt::format("H0(3,2) = {:g} A\n",domain->yz);
+  header += fmt::format("H0(3,3) = {:g} A\n",domain->zprd);
+  header += fmt::format(".NO_VELOCITY.\n");
+  header += fmt::format("entry_count = {}\n",nfield-2);
   for (int i = 0; i < nfield-5; i++)
-    fprintf(fp,"auxiliary[%d] = %s\n",i,auxname[i]);
+    header += fmt::format("auxiliary[{}] = {}\n",i,auxname[i]);
+  fmt::print(fp, header);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp
index e90243fd1f..d8819820ec 100644
--- a/src/dump_xyz.cpp
+++ b/src/dump_xyz.cpp
@@ -130,10 +130,10 @@ int DumpXYZ::modify_param(int narg, char **arg)
 void DumpXYZ::write_header(bigint n)
 {
   if (me == 0) {
-    if (time_flag) {
-      double tcurrent = update->atime + (update->ntimestep-update->atimestep) + update->dt;
-      fmt::print(fp,"{}\n Atoms. Timestep: {} Time: {:.6f}\n", n, update->ntimestep, tcurrent);
-    } else fmt::print(fp,"{}\n Atoms. Timestep: {}\n", n, update->ntimestep);
+    auto header = fmt::format("{}\n Atoms. Timestep: {}", n, update->ntimestep);
+    if (time_flag) header += fmt::format(" Time: {:.6f}", compute_time());
+    header += "\n";
+    fmt::print(fp, header);
   }
 }