Unverified Commit 9e1d560f authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

refactor using fmtlib and utils::

parent 284b1618
Loading
Loading
Loading
Loading
+26 −52
Original line number Diff line number Diff line
@@ -36,6 +36,7 @@
#include "error.h"
#include "rigid_const.h"
#include "utils.h"
#include "fmt/format.h"

using namespace LAMMPS_NS;
using namespace FixConst;
@@ -414,9 +415,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"temp") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 &&
          strcmp(style,"rigid/nvt/omp") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0)
      if (!utils::strmatch(style,"^rigid/n.t"))
        error->all(FLERR,"Illegal fix rigid command");
      tstat_flag = 1;
      t_start = force->numeric(FLERR,arg[iarg+1]);
@@ -426,9 +425,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"iso") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      pcouple = XYZ;
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
@@ -444,9 +441,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"aniso") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
@@ -461,9 +456,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"x") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      p_start[0] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
@@ -473,9 +466,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"y") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      p_start[1] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
@@ -485,9 +476,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"z") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      p_start[2] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
@@ -524,9 +513,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"tparam") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0 &&
          strcmp(style,"rigid/nvt/omp") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0)
      if (!utils::strmatch(style,"^rigid/n.t"))
        error->all(FLERR,"Illegal fix rigid command");
      t_chain = force->inumeric(FLERR,arg[iarg+1]);
      t_iter = force->inumeric(FLERR,arg[iarg+2]);
@@ -535,9 +522,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"pchain") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
      if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0 &&
          strcmp(style,"rigid/npt/omp") != 0 &&
          strcmp(style,"rigid/nph/omp") != 0)
      if (!utils::strmatch(style,"^rigid/np."))
        error->all(FLERR,"Illegal fix rigid command");
      p_chain = force->inumeric(FLERR,arg[iarg+1]);
      iarg += 2;
@@ -623,10 +608,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
  int nsum = 0;
  for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];

  if (me == 0) {
    if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
    if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
  }
  if (me == 0)
    utils::logmesg(lmp,fmt::format("  {} rigid bodies with {} atoms\n",
                                   nbody,nsum));
}

/* ---------------------------------------------------------------------- */
@@ -720,12 +704,9 @@ void FixRigid::init()
    for (i = 0; i < modify->nfix; i++) {
      if (modify->fix[i]->rigid_flag) rflag = 1;
      if (rflag && (modify->fmask[i] & POST_FORCE) &&
          !modify->fix[i]->rigid_flag) {
        char str[128];
        snprintf(str,128,"Fix %s alters forces after fix rigid",
                 modify->fix[i]->id);
        error->warning(FLERR,str);
      }
          !modify->fix[i]->rigid_flag)
        error->warning(FLERR,fmt::format("Fix {} alters forces after fix "
                                         "rigid",modify->fix[i]->id));
    }
  }

@@ -2300,12 +2281,9 @@ void FixRigid::readfile(int which, double *vec,

  if (me == 0) {
    fp = fopen(inpfile,"r");
    if (fp == NULL) {
      char str[128];
      snprintf(str,128,"Cannot open fix rigid inpfile %s",inpfile);
      error->one(FLERR,str);
    }

    if (fp == NULL)
      error->one(FLERR,fmt::format("Cannot open fix rigid file {}: {}",
                                   inpfile,utils::getsyserror()));
    while (1) {
      eof = fgets(line,MAXLINE,fp);
      if (eof == NULL) error->one(FLERR,"Unexpected end of fix rigid file");
@@ -2410,19 +2388,15 @@ void FixRigid::write_restart_file(const char *file)
{
  if (me) return;

  char outfile[128];
  snprintf(outfile,128,"%s.rigid",file);
  FILE *fp = fopen(outfile,"w");
  if (fp == NULL) {
    char str[192];
    snprintf(str,192,"Cannot open fix rigid restart file %s",outfile);
    error->one(FLERR,str);
  }
  auto outfile = std::string(file) + ".rigid";
  FILE *fp = fopen(outfile.c_str(),"w");
  if (fp == NULL)
    error->one(FLERR,fmt::format("Cannot open fix rigid restart file {}: {}",
                                 outfile,utils::getsyserror()));

  fprintf(fp,"# fix rigid mass, COM, inertia tensor info for "
          "%d bodies on timestep " BIGINT_FORMAT "\n\n",
          nbody,update->ntimestep);
  fprintf(fp,"%d\n",nbody);
  fmt::print(fp,"# fix rigid mass, COM, inertia tensor info for "
             "{} bodies on timestep {}\n\n",nbody,update->ntimestep);
  fmt::print(fp,"{}\n",nbody);

  // compute I tensor against xyz axes from diagonalized I and current quat
  // Ispace = P Idiag P_transpose
+43 −68
Original line number Diff line number Diff line
@@ -40,6 +40,7 @@
#include "error.h"
#include "rigid_const.h"
#include "utils.h"
#include "fmt/format.h"

#include <map>

@@ -248,8 +249,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"temp") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (strcmp(style,"rigid/nvt/small") != 0 &&
          strcmp(style,"rigid/npt/small") != 0)
      if (!utils::strmatch(style,"^rigid/n.t/small"))
        error->all(FLERR,"Illegal fix rigid command");
      tstat_flag = 1;
      t_start = force->numeric(FLERR,arg[iarg+1]);
@@ -259,8 +259,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"iso") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (strcmp(style,"rigid/npt/small") != 0 &&
          strcmp(style,"rigid/nph/small") != 0)
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      pcouple = XYZ;
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
@@ -276,8 +275,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"aniso") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (strcmp(style,"rigid/npt/small") != 0 &&
          strcmp(style,"rigid/nph/small") != 0)
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
@@ -292,6 +290,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"x") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      p_start[0] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
      p_period[0] = force->numeric(FLERR,arg[iarg+3]);
@@ -300,6 +300,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"y") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      p_start[1] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
      p_period[1] = force->numeric(FLERR,arg[iarg+3]);
@@ -308,6 +310,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"z") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      p_start[2] = force->numeric(FLERR,arg[iarg+1]);
      p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
      p_period[2] = force->numeric(FLERR,arg[iarg+3]);
@@ -343,20 +347,18 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :

    } else if (strcmp(arg[iarg],"tparam") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (strcmp(style,"rigid/nvt/small") != 0 &&
          strcmp(style,"rigid/npt/small") != 0)
      if (!utils::strmatch(style,"^rigid/n.t/small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      t_chain = force->numeric(FLERR,arg[iarg+1]);
      t_iter = force->numeric(FLERR,arg[iarg+2]);
      t_order = force->numeric(FLERR,arg[iarg+3]);
      t_chain = force->inumeric(FLERR,arg[iarg+1]);
      t_iter = force->inumeric(FLERR,arg[iarg+2]);
      t_order = force->inumeric(FLERR,arg[iarg+3]);
      iarg += 4;

    } else if (strcmp(arg[iarg],"pchain") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
      if (strcmp(style,"rigid/npt/small") != 0 &&
          strcmp(style,"rigid/nph/small") != 0)
      if (!utils::strmatch(style,"^rigid/np./small"))
        error->all(FLERR,"Illegal fix rigid/small command");
      p_chain = force->numeric(FLERR,arg[iarg+1]);
      p_chain = force->inumeric(FLERR,arg[iarg+1]);
      iarg += 2;

    } else if (strcmp(arg[iarg],"gravity") == 0) {
@@ -406,14 +408,9 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
  create_bodies(bodyID);
  if (customflag) delete [] bodyID;

  double time2 = MPI_Wtime();

  if (comm->me == 0) {
    if (screen)
      fprintf(screen,"  create bodies CPU = %g seconds\n",time2-time1);
    if (logfile)
      fprintf(logfile,"  create bodies CPU = %g seconds\n",time2-time1);
  }
  if (comm->me == 0)
    utils::logmesg(lmp,fmt::format("  create bodies CPU = {:.3f} seconds\n",
                                   MPI_Wtime()-time1));

  // set nlocal_body and allocate bodies I own

@@ -471,18 +468,10 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
  MPI_Allreduce(&atomone,&atomall,1,MPI_LMP_BIGINT,MPI_SUM,world);

  if (me == 0) {
    if (screen) {
      fprintf(screen,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
              nbody,atomall);
      fprintf(screen,"  %g = max distance from body owner to body atom\n",
              maxextent);
    }
    if (logfile) {
      fprintf(logfile,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
              nbody,atomall);
      fprintf(logfile,"  %g = max distance from body owner to body atom\n",
    auto msg = fmt::format("  {} rigid bodies with {} atoms\n",nbody,atomall);
    msg += fmt::format("  {:.8} = max distance from body owner to body atom\n",
                       maxextent);
    }
    utils::logmesg(lmp,msg);
  }

  // initialize Marsaglia RNG with processor-unique seed
@@ -567,12 +556,9 @@ void FixRigidSmall::init()
    for (i = 0; i < modify->nfix; i++) {
      if (modify->fix[i]->rigid_flag) rflag = 1;
      if (rflag && (modify->fmask[i] & POST_FORCE) &&
          !modify->fix[i]->rigid_flag) {
        char str[128];
        snprintf(str,128,"Fix %s alters forces after fix rigid",
                 modify->fix[i]->id);
        error->warning(FLERR,str);
      }
          !modify->fix[i]->rigid_flag)
        error->warning(FLERR,fmt::format("Fix {} alters forces after fix "
                                         "rigid",modify->fix[i]->id));
    }
  }

@@ -2459,12 +2445,9 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody)

  if (me == 0) {
    fp = fopen(inpfile,"r");
    if (fp == NULL) {
      char str[128];
      snprintf(str,128,"Cannot open fix rigid/small inpfile %s",inpfile);
      error->one(FLERR,str);
    }

    if (fp == NULL)
      error->one(FLERR,fmt::format("Cannot open fix rigid/small file {}: {}",
                                   inpfile,utils::getsyserror()));
    while (1) {
      eof = fgets(line,MAXLINE,fp);
      if (eof == NULL)
@@ -2574,19 +2557,15 @@ void FixRigidSmall::write_restart_file(const char *file)
  // proc 0 opens file and writes header

  if (me == 0) {
    char outfile[128];
    snprintf(outfile,128,"%s.rigid",file);
    fp = fopen(outfile,"w");
    if (fp == NULL) {
      char str[128];
      snprintf(str,128,"Cannot open fix rigid restart file %s",outfile);
      error->one(FLERR,str);
    }
    auto outfile = std::string(file) + ".rigid";
    fp = fopen(outfile.c_str(),"w");
    if (fp == NULL)
      error->one(FLERR,fmt::format("Cannot open fix rigid restart file {}: {}",
                                   outfile,utils::getsyserror()));

    fprintf(fp,"# fix rigid mass, COM, inertia tensor info for "
            "%d bodies on timestep " BIGINT_FORMAT "\n\n",
            nbody,update->ntimestep);
    fprintf(fp,"%d\n",nbody);
    fmt::print(fp,"# fix rigid mass, COM, inertia tensor info for "
               "{} bodies on timestep {}\n\n",nbody,update->ntimestep);
    fmt::print(fp,"{}\n",nbody);
  }

  // communication buffer for all my rigid body info
@@ -3324,15 +3303,11 @@ void FixRigidSmall::reset_atom2body()
    atom2body[i] = -1;
    if (bodytag[i]) {
      iowner = atom->map(bodytag[i]);
      if (iowner == -1) {
        char str[128];
        sprintf(str,
                "Rigid body atoms " TAGINT_FORMAT " " TAGINT_FORMAT
                " missing on proc %d at step " BIGINT_FORMAT,
                atom->tag[i],bodytag[i],comm->me,update->ntimestep);
        error->one(FLERR,str);
      if (iowner == -1)
        error->one(FLERR,fmt::format("Rigid body atoms {} {} missing on "
                                     "proc {} at step {}",atom->tag[i],
                                     bodytag[i],comm->me,update->ntimestep));

      }
      atom2body[i] = bodyown[iowner];
    }
  }
+39 −77
Original line number Diff line number Diff line
@@ -33,6 +33,7 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"

using namespace LAMMPS_NS;
using namespace FixConst;
@@ -225,14 +226,9 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :

  find_clusters();

  double time2 = MPI_Wtime();

  if (comm->me == 0) {
    if (screen)
      fprintf(screen,"  find clusters CPU = %g seconds\n",time2-time1);
    if (logfile)
      fprintf(logfile,"  find clusters CPU = %g seconds\n",time2-time1);
  }
  if (comm->me == 0)
    utils::logmesg(lmp,fmt::format("  find clusters CPU = {:.3f} seconds\n",
                                   MPI_Wtime()-time1));

  // initialize list of SHAKE clusters to constrain

@@ -525,44 +521,31 @@ void FixShake::pre_neighbor()
      if (shake_flag[i] == 2) {
        atom1 = atom->map(shake_atom[i][0]);
        atom2 = atom->map(shake_atom[i][1]);
        if (atom1 == -1 || atom2 == -1) {
          char str[128];
          sprintf(str,"Shake atoms " TAGINT_FORMAT " " TAGINT_FORMAT
                  " missing on proc %d at step " BIGINT_FORMAT,
                  shake_atom[i][0],shake_atom[i][1],me,update->ntimestep);
          error->one(FLERR,str);
        }
        if (atom1 == -1 || atom2 == -1)
          error->one(FLERR,fmt::format("Shake atoms {} {} missing on proc "
                                       "{} at step {}",shake_atom[i][0],
                                       shake_atom[i][1],me,update->ntimestep));
        if (i <= atom1 && i <= atom2) list[nlist++] = i;
      } else if (shake_flag[i] % 2 == 1) {
        atom1 = atom->map(shake_atom[i][0]);
        atom2 = atom->map(shake_atom[i][1]);
        atom3 = atom->map(shake_atom[i][2]);
        if (atom1 == -1 || atom2 == -1 || atom3 == -1) {
          char str[128];
          sprintf(str,"Shake atoms "
                  TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
                  " missing on proc %d at step " BIGINT_FORMAT,
                  shake_atom[i][0],shake_atom[i][1],shake_atom[i][2],
                  me,update->ntimestep);
          error->one(FLERR,str);
        }
        if (atom1 == -1 || atom2 == -1 || atom3 == -1)
          error->one(FLERR,fmt::format("Shake atoms {} {} {} missing on proc "
                                       "{} at step {}",shake_atom[i][0],
                                       shake_atom[i][1],shake_atom[i][2],
                                       me,update->ntimestep));
        if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
      } else {
        atom1 = atom->map(shake_atom[i][0]);
        atom2 = atom->map(shake_atom[i][1]);
        atom3 = atom->map(shake_atom[i][2]);
        atom4 = atom->map(shake_atom[i][3]);
        if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) {
          char str[128];
          sprintf(str,"Shake atoms "
                  TAGINT_FORMAT " " TAGINT_FORMAT " "
                  TAGINT_FORMAT " " TAGINT_FORMAT
                  " missing on proc %d at step " BIGINT_FORMAT,
                  shake_atom[i][0],shake_atom[i][1],
                  shake_atom[i][2],shake_atom[i][3],
                  me,update->ntimestep);
          error->one(FLERR,str);
        }
        if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
          error->one(FLERR,fmt::format("Shake atoms {} {} {} {} missing on "
                                       "proc {} at step {}",shake_atom[i][0],
                                       shake_atom[i][1],shake_atom[i][2],
                                       shake_atom[i][3],me,update->ntimestep));
        if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
          list[nlist++] = i;
      }
@@ -697,11 +680,10 @@ void FixShake::find_clusters()
  tagint tagprev;
  double massone;

  if (me == 0 && screen) {
    if (!rattle) fprintf(screen,"Finding SHAKE clusters ...\n");
    else fprintf(screen,"Finding RATTLE clusters ...\n");
  if ((me == 0) && screen) {
    if (!rattle) fputs("Finding SHAKE clusters ...\n",screen);
    else fputs("Finding RATTLE clusters ...\n",screen);
  }

  atommols = atom->avec->onemols;

  tagint *tag = atom->tag;
@@ -1024,18 +1006,11 @@ void FixShake::find_clusters()
  MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world);

  if (me == 0) {
    if (screen) {
      fprintf(screen,"  %d = # of size 2 clusters\n",count2/2);
      fprintf(screen,"  %d = # of size 3 clusters\n",count3/3);
      fprintf(screen,"  %d = # of size 4 clusters\n",count4/4);
      fprintf(screen,"  %d = # of frozen angles\n",count1/3);
    }
    if (logfile) {
      fprintf(logfile,"  %d = # of size 2 clusters\n",count2/2);
      fprintf(logfile,"  %d = # of size 3 clusters\n",count3/3);
      fprintf(logfile,"  %d = # of size 4 clusters\n",count4/4);
      fprintf(logfile,"  %d = # of frozen angles\n",count1/3);
    }
    auto mesg = fmt::format("{:>8} = # of size 2 clusters\n",count2/2);
    mesg += fmt::format("{:>8} = # of size 3 clusters\n",count3/3);
    mesg += fmt::format("{:>8} = # of size 4 clusters\n",count4/4);
    mesg += fmt::format("{:>8} = # of frozen angles\n",count1/3);
    utils::logmesg(lmp,mesg);
  }
}

@@ -2558,34 +2533,21 @@ void FixShake::stats()
  // print stats only for non-zero counts

  if (me == 0) {

    if (screen) {
      fprintf(screen,
              "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n",
    auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n",
                            update->ntimestep);
      for (i = 1; i < nb; i++)
        if (b_count_all[i])
          fprintf(screen,"  %d %g %g %d\n",i,
                  b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i],
                  b_count_all[i]);
      for (i = 1; i < na; i++)
        if (a_count_all[i])
          fprintf(screen,"  %d %g %g\n",i,
                  a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]);
    }
    if (logfile) {
      fprintf(logfile,
              "SHAKE stats (type/ave/delta) on step " BIGINT_FORMAT "\n",
              update->ntimestep);
      for (i = 0; i < nb; i++)
        if (b_count_all[i])
          fprintf(logfile,"  %d %g %g\n",i,
                  b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i]);
      for (i = 0; i < na; i++)
        if (a_count_all[i])
          fprintf(logfile,"  %d %g %g\n",i,
                  a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]);
    }
    for (i = 1; i < nb; i++) {
      const auto bcnt = b_count_all[i];
      if (bcnt)
        mesg += fmt::format("{:>6d}   {:<9.6} {:<11.6} {:>8d}\n",i,
                            b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt);
    }
    for (i = 1; i < na; i++) {
      const auto acnt = a_count_all[i];
      if (acnt)
        mesg += fmt::format("{:>6d}   {:<9.6} {:<11.6} {:>8d}\n",i,
                            a_ave_all[i]/acnt,a_max_all[i]-a_min_all[i],acnt);
    }
    utils::logmesg(lmp,mesg);
  }

  // next timestep for stats