diff --git a/.github/codecov.yml b/.github/codecov.yml
new file mode 100644
index 0000000000..4645f52f7f
--- /dev/null
+++ b/.github/codecov.yml
@@ -0,0 +1,29 @@
+comment: false
+coverage:
+  notify:
+    slack:
+      default:
+        url: "secret:HWZbvgtc6OD7F3v3PfrK3/rzCJvScbh69Fi1CkLwuHK0+wIBIHVR+Q5i7q6F9Ln4OChbiRGtYAEUUsT8/jmBu4qDpIi8mx746codc0z/Z3aafLd24pBrCEPLvdCfIZxqPnw3TuUgGhwmMDZf0+thg8YNUr/MbOZ7Li2L6+ZbYuA="
+        threshold: 10%
+        only_pulls: false
+        branches:
+          - "master"
+        flags:
+          - "unit"
+        paths:
+          - "src"
+  status:
+    project:
+        default:
+            branches:
+              - "master"
+            paths:
+              - "src"
+            informational: true
+    patch:
+        default:
+            branches:
+              - "master"
+            paths:
+              - "src"
+            informational: true
diff --git a/cmake/Modules/Testing.cmake b/cmake/Modules/Testing.cmake
index 29f8b61802..6d2f39b0f3 100644
--- a/cmake/Modules/Testing.cmake
+++ b/cmake/Modules/Testing.cmake
@@ -23,14 +23,17 @@ if(ENABLE_TESTING)
         OR (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang"))
       include(CheckCXXCompilerFlag)
       set(CMAKE_CUSTOM_LINKER_DEFAULT default)
-      check_cxx_compiler_flag(-fuse-ld=lld HAVE_LLD_LINKER)
-      check_cxx_compiler_flag(-fuse-ld=gold HAVE_GOLD_LINKER)
-      check_cxx_compiler_flag(-fuse-ld=bfd HAVE_BFD_LINKER)
-      if(HAVE_LLD_LINKER)
+      check_cxx_compiler_flag(-fuse-ld=lld HAVE_LLD_LINKER_FLAG)
+      check_cxx_compiler_flag(-fuse-ld=gold HAVE_GOLD_LINKER_FLAG)
+      check_cxx_compiler_flag(-fuse-ld=bfd HAVE_BFD_LINKER_FLAG)
+      find_program(HAVE_LLD_LINKER_BIN lld ld.lld)
+      find_program(HAVE_GOLD_LINKER_BIN ld.gold)
+      find_program(HAVE_BFD_LINKER_BIN ld.bfd)
+      if(HAVE_LLD_LINKER_FLAG AND HAVE_LLD_LINKER_BIN)
         set(CMAKE_CUSTOM_LINKER_DEFAULT lld)
-      elseif(HAVE_GOLD_LINKER)
+      elseif(HAVE_GOLD_LINKER_FLAG AND HAVE_GOLD_LINKER_BIN)
         set(CMAKE_CUSTOM_LINKER_DEFAULT gold)
-      elseif(HAVE_BFD_LINKER)
+      elseif(HAVE_BFD_LINKER_FLAG AND HAVE_BFD_LINKER_BIN)
         set(CMAKE_CUSTOM_LINKER_DEFAULT bfd)
       endif()
       set(CMAKE_CUSTOM_LINKER_VALUES lld gold bfd default)
diff --git a/doc/src/Errors_warnings.rst b/doc/src/Errors_warnings.rst
index ceccc0b74c..306c9b7b31 100644
--- a/doc/src/Errors_warnings.rst
+++ b/doc/src/Errors_warnings.rst
@@ -118,6 +118,24 @@ Doc page with :doc:`ERROR messages <Errors_messages>`
    incorrect periodic images of atoms in interaction lists.  To avoid, either use
    :doc:`pair style zero <pair_zero>` with a suitable cutoff or use :doc:`comm_modify cutoff <comm_modify>`.
 
+*Communication cutoff is shorter than a bond length based estimate. This may lead to errors.*
+
+   Since LAMMPS stores topology data with individual atoms, all atoms
+   comprising a bond, angle, dihedral or improper must be present on any
+   sub-domain that "owns" the atom with the information, either as a
+   local or a ghost atom. The communication cutoff is what determines up
+   to what distance from a sub-domain boundary ghost atoms are created.
+   The communication cutoff is by default the largest non-bonded cutoff
+   plus the neighbor skin distance, but for short or non-bonded cutoffs
+   and/or long bonds, this may not be sufficient. This warning indicates
+   that there is an increased risk of a simulation stopping unexpectedly
+   because of Bond/Angle/Dihedral/Improper atoms missing.  It can be
+   silenced by manually setting the communication cutoff via
+   :doc:`comm_modify cutoff <comm_modify>`.  However, since the
+   heuristic used to determine the estimate is not always accurate, it
+   is not changed automatically and the warning may be ignored depending
+   on the specific system being simulated.
+
 *Communication cutoff is too small for SNAP micro load balancing, increased to %lf*
    Self-explanatory.
 
diff --git a/src/MC/fix_widom.h b/src/MC/fix_widom.h
index 8dd77df28c..05060e76ac 100644
--- a/src/MC/fix_widom.h
+++ b/src/MC/fix_widom.h
@@ -32,7 +32,7 @@ class FixWidom : public Fix {
   int setmask();
   void init();
   void pre_exchange();
-  
+
   void attempt_atomic_insertion();
   void attempt_molecule_insertion();
 
@@ -207,7 +207,7 @@ is on.
 
 E: Cannot use fix Widom in a 2d simulation
 
-Fix Widom is set up to run in 3d only. No 2d simulations with fix Widom 
+Fix Widom is set up to run in 3d only. No 2d simulations with fix Widom
 are allowed.
 
 E: Could not find fix Widom exclusion group ID
diff --git a/src/USER-MISC/pair_lennard_mdf.cpp b/src/USER-MISC/pair_lennard_mdf.cpp
index c539bb30e2..5e7e24af67 100644
--- a/src/USER-MISC/pair_lennard_mdf.cpp
+++ b/src/USER-MISC/pair_lennard_mdf.cpp
@@ -32,11 +32,11 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairLJ_AB_MDF::PairLJ_AB_MDF(LAMMPS *lmp) : Pair(lmp) {}
+PairLennardMDF::PairLennardMDF(LAMMPS *lmp) : Pair(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
-PairLJ_AB_MDF::~PairLJ_AB_MDF()
+PairLennardMDF::~PairLennardMDF()
 {
   if (allocated) {
     memory->destroy(setflag);
@@ -56,7 +56,7 @@ PairLJ_AB_MDF::~PairLJ_AB_MDF()
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::compute(int eflag, int vflag)
+void PairLennardMDF::compute(int eflag, int vflag)
 {
   int i,j,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
@@ -153,7 +153,7 @@ void PairLJ_AB_MDF::compute(int eflag, int vflag)
    allocate all arrays
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::allocate()
+void PairLennardMDF::allocate()
 {
   allocated = 1;
   int n = atom->ntypes;
@@ -180,7 +180,7 @@ void PairLJ_AB_MDF::allocate()
    global settings
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::settings(int narg, char **arg)
+void PairLennardMDF::settings(int narg, char **arg)
 {
   if (narg != 2) error->all(FLERR,"Illegal pair_style command");
 
@@ -207,7 +207,7 @@ void PairLJ_AB_MDF::settings(int narg, char **arg)
    set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::coeff(int narg, char **arg)
+void PairLennardMDF::coeff(int narg, char **arg)
 {
   if (narg != 4 && narg != 6)
     error->all(FLERR,"Incorrect args for pair coefficients");
@@ -248,7 +248,7 @@ void PairLJ_AB_MDF::coeff(int narg, char **arg)
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
 
-double PairLJ_AB_MDF::init_one(int i, int j)
+double PairLennardMDF::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
 
@@ -259,6 +259,7 @@ double PairLJ_AB_MDF::init_one(int i, int j)
   lj3[i][j] = aparm[i][j];
   lj4[i][j] = bparm[i][j];
 
+  cut[j][i] = cut[i][j];
   cut_inner[j][i] = cut_inner[i][j];
   cut_inner_sq[j][i] = cut_inner_sq[i][j];
   lj1[j][i] = lj1[i][j];
@@ -273,7 +274,7 @@ double PairLJ_AB_MDF::init_one(int i, int j)
   proc 0 writes to restart file
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::write_restart(FILE *fp)
+void PairLennardMDF::write_restart(FILE *fp)
 {
   write_restart_settings(fp);
 
@@ -294,7 +295,7 @@ void PairLJ_AB_MDF::write_restart(FILE *fp)
   proc 0 reads from restart file, bcasts
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::read_restart(FILE *fp)
+void PairLennardMDF::read_restart(FILE *fp)
 {
   read_restart_settings(fp);
   allocate();
@@ -324,7 +325,7 @@ void PairLJ_AB_MDF::read_restart(FILE *fp)
   proc 0 writes to restart file
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::write_restart_settings(FILE *fp)
+void PairLennardMDF::write_restart_settings(FILE *fp)
 {
   fwrite(&mix_flag,sizeof(int),1,fp);
 }
@@ -333,7 +334,7 @@ void PairLJ_AB_MDF::write_restart_settings(FILE *fp)
   proc 0 reads from restart file, bcasts
 ------------------------------------------------------------------------- */
 
-void PairLJ_AB_MDF::read_restart_settings(FILE *fp)
+void PairLennardMDF::read_restart_settings(FILE *fp)
 {
   int me = comm->me;
   if (me == 0) {
@@ -342,9 +343,32 @@ void PairLJ_AB_MDF::read_restart_settings(FILE *fp)
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }
 
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairLennardMDF::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d %g %g\n",i,aparm[i][i],bparm[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairLennardMDF::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d %g %g %g %g\n",
+              i,j,aparm[i][j],bparm[i][j],
+              cut_inner[i][j],cut[i][j]);
+}
+
 /* ---------------------------------------------------------------------- */
 
-double PairLJ_AB_MDF::single(int /*i*/, int /*j*/, int itype, int jtype,
+double PairLennardMDF::single(int /*i*/, int /*j*/, int itype, int jtype,
                              double rsq,
                              double /*factor_coul*/, double factor_lj,
                              double &fforce)
@@ -368,7 +392,7 @@ double PairLJ_AB_MDF::single(int /*i*/, int /*j*/, int itype, int jtype,
     dt = 30.* d*d * dd*dd * rr / dp;
 
     forcelj = forcelj*tt + philj*dt;
-    philj *= dt;
+    philj *= tt;
   }
 
   fforce = factor_lj*forcelj*r2inv;
@@ -378,7 +402,7 @@ double PairLJ_AB_MDF::single(int /*i*/, int /*j*/, int itype, int jtype,
 
 /* ---------------------------------------------------------------------- */
 
-void *PairLJ_AB_MDF::extract(const char *str, int &dim)
+void *PairLennardMDF::extract(const char *str, int &dim)
 {
   dim = 2;
   if (strcmp(str,"a") == 0) return (void *) aparm;
diff --git a/src/USER-MISC/pair_lennard_mdf.h b/src/USER-MISC/pair_lennard_mdf.h
index c4ffc80cd6..19325429ea 100644
--- a/src/USER-MISC/pair_lennard_mdf.h
+++ b/src/USER-MISC/pair_lennard_mdf.h
@@ -13,21 +13,21 @@
 
 #ifdef PAIR_CLASS
 
-PairStyle(lennard/mdf,PairLJ_AB_MDF)
+PairStyle(lennard/mdf,PairLennardMDF)
 
 #else
 
-#ifndef LMP_PAIR_LJ_AB_MDF_H
-#define LMP_PAIR_LJ_AB_MDF_H
+#ifndef LMP_PAIR_LENNARD_MDF_H
+#define LMP_PAIR_LENNARD_MDF_H
 
 #include "pair.h"
 
 namespace LAMMPS_NS {
 
-class PairLJ_AB_MDF : public Pair {
+class PairLennardMDF : public Pair {
  public:
-  PairLJ_AB_MDF(class LAMMPS *);
-  virtual ~PairLJ_AB_MDF();
+  PairLennardMDF(class LAMMPS *);
+  virtual ~PairLennardMDF();
 
   virtual void compute(int, int);
   void settings(int, char **);
@@ -37,6 +37,8 @@ class PairLJ_AB_MDF : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
   double single(int, int, int, int, double, double, double, double &);
   void *extract(const char *, int &);
 
diff --git a/src/USER-OMP/reaxc_init_md_omp.cpp b/src/USER-OMP/reaxc_init_md_omp.cpp
index 4ef214bf47..b0cf4a5fb8 100644
--- a/src/USER-OMP/reaxc_init_md_omp.cpp
+++ b/src/USER-OMP/reaxc_init_md_omp.cpp
@@ -36,6 +36,7 @@
 #include "reaxc_lookup.h"
 #include "reaxc_tool_box.h"
 #include "error.h"
+#include "fmt/format.h"
 
 // Functions defined in reaxc_init_md.cpp
 extern int Init_MPI_Datatypes(reax_system*, storage*, mpi_datatypes*, MPI_Comm, char*);
@@ -55,37 +56,38 @@ int Init_ListsOMP(reax_system *system, control_params *control,
   int mincap = system->mincap;
   double safezone = system->safezone;
   double saferzone = system->saferzone;
+  LAMMPS_NS::Error *error = system->error_ptr;
 
-  bond_top = (int*) calloc( system->total_cap, sizeof(int) );
-  hb_top = (int*) calloc( system->local_cap, sizeof(int) );
-  Estimate_Storages( system, control, lists,
-                     &Htop, hb_top, bond_top, &num_3body );
+  bond_top = (int*) calloc(system->total_cap, sizeof(int));
+  hb_top = (int*) calloc(system->local_cap, sizeof(int));
+  Estimate_Storages(system, control, lists,
+                     &Htop, hb_top, bond_top, &num_3body);
 
   if (control->hbond_cut > 0) {
     /* init H indexes */
     total_hbonds = 0;
-    for( i = 0; i < system->n; ++i ) {
+    for(i = 0; i < system->n; ++i) {
       system->my_atoms[i].num_hbonds = hb_top[i];
       total_hbonds += hb_top[i];
     }
     total_hbonds = (int)(MAX(total_hbonds*saferzone,mincap*system->minhbonds));
 
-    if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
-                    *lists+HBONDS ) ) {
-      system->error_ptr->one( FLERR, "Not enough space for hbonds list. Terminating!" );
+    if(!Make_List(system->Hcap, total_hbonds, TYP_HBOND,
+                    *lists+HBONDS)) {
+      error->one(FLERR, "Not enough space for hbonds list. Terminating!");
     }
   }
 
   total_bonds = 0;
-  for( i = 0; i < system->N; ++i ) {
+  for(i = 0; i < system->N; ++i) {
     system->my_atoms[i].num_bonds = bond_top[i];
     total_bonds += bond_top[i];
   }
-  bond_cap = (int)(MAX( total_bonds*safezone, mincap*MIN_BONDS ));
+  bond_cap = (int)(MAX(total_bonds*safezone, mincap*MIN_BONDS));
 
-  if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
-                  *lists+BONDS ) ) {
-    system->error_ptr->one( FLERR, "Not enough space for bonds list. Terminating!\n" );
+  if(!Make_List(system->total_cap, bond_cap, TYP_BOND,
+                  *lists+BONDS)) {
+    error->one(FLERR, "Not enough space for bonds list. Terminating!\n");
   }
 
   int nthreads = control->nthreads;
@@ -93,18 +95,18 @@ int Init_ListsOMP(reax_system *system, control_params *control,
 
   for (i = 0; i < bonds->num_intrs; ++i)
     bonds->select.bond_list[i].bo_data.CdboReduction =
-      (double*) smalloc(system->error_ptr, sizeof(double)*nthreads, "CdboReduction");
+      (double*) smalloc(error, sizeof(double)*nthreads, "CdboReduction");
 
   /* 3bodies list */
-  cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
-  if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
-                  *lists+THREE_BODIES ) ){
+  cap_3body = (int)(MAX(num_3body*safezone, MIN_3BODIES));
+  if(!Make_List(bond_cap, cap_3body, TYP_THREE_BODY,
+                  *lists+THREE_BODIES)){
 
-    system->error_ptr->one( FLERR, "Problem in initializing angles list. Terminating!" );
+    error->one(FLERR, "Problem in initializing angles list. Terminating!");
   }
 
-  free( hb_top );
-  free( bond_top );
+  free(hb_top);
+  free(bond_top);
 
   return SUCCESS;
 }
@@ -112,60 +114,43 @@ int Init_ListsOMP(reax_system *system, control_params *control,
 /* ---------------------------------------------------------------------- */
 
 // The only difference with the MPI-only function is calls to Init_ListsOMP and Init_Force_FunctionsOMP().
-void InitializeOMP( reax_system *system, control_params *control,
+void InitializeOMP(reax_system *system, control_params *control,
                  simulation_data *data, storage *workspace,
                  reax_list **lists, output_controls *out_control,
-                 mpi_datatypes *mpi_data, MPI_Comm comm )
+                 mpi_datatypes *mpi_data, MPI_Comm comm)
 {
   char msg[MAX_STR];
-  char errmsg[512];
+  LAMMPS_NS::Error *error = system->error_ptr;
 
+  if (Init_MPI_Datatypes(system,workspace,mpi_data,comm,msg) == FAILURE)
+    error->one(FLERR,"init_mpi_datatypes: could not create datatypes. "
+               "Mpi_data could not be initialized! Terminating.");
 
-  if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
-    system->error_ptr->one( FLERR, "init_mpi_datatypes: could not create datatypes. "
-    "Mpi_data couldn't be initialized! Terminating.");
-  }
+  if (Init_System(system,control,msg) == FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. System could not be "
+                                  "initialized! Terminating.",msg));
 
-  if (Init_System(system, control, msg) == FAILURE) {
-    snprintf( errmsg, 512, "Error on: %s. "
-    "System could not be initialized! Terminating.", msg );
-    system->error_ptr->one(FLERR, errmsg);
-  }
+  if (Init_Simulation_Data(system,control,data,msg) == FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. Sim_data could not be "
+                                  "initialized! Terminating.",msg));
 
-  if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
-    snprintf( errmsg, 512, "Error on: %s. "
-    "Sim_data couldn't be initialized! Terminating.", msg );
-    system->error_ptr->one(FLERR, errmsg);
-  }
+  if (Init_Workspace(system,control,workspace,msg) == FAILURE)
+    error->one(FLERR,"init_workspace: not enough memory. "
+               "Workspace could not be initialized. Terminating.");
 
-  if (Init_Workspace( system, control, workspace, msg ) ==
-      FAILURE) {
-    system->error_ptr->one(FLERR, "init_workspace: not enough memory. "
-    "Workspace couldn't be initialized! Terminating.");
-  }
+  if (Init_ListsOMP(system,control,data,workspace,lists,mpi_data,msg) == FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. System could not be "
+                                  "initialized. Terminating.",msg));
 
-  if (Init_ListsOMP( system, control, data, workspace, lists, mpi_data, msg ) ==
-      FAILURE) {
-      snprintf( errmsg, 512, "Error on: %s. "
-      "System could not be initialized! Terminating.", msg );
-      system->error_ptr->one(FLERR, errmsg);
-    }
+  if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. Could not open output files! "
+                                  "Terminating.",msg));
 
-  if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
-    snprintf( errmsg, 512, "Error on: %s"
-    "Could not open output files! Terminating.", msg );
-    system->error_ptr->one(FLERR, errmsg);
-  }
+  if (control->tabulate)
+    if (Init_Lookup_Tables(system,control,workspace,mpi_data,msg) == FAILURE)
+      error->one(FLERR,fmt::format("Error on: {}. Could not create lookup "
+                                    "table. Terminating.",msg));
 
-  if (control->tabulate) {
-    if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
-      snprintf( errmsg, 512, "Error on: %s."
-      " Couldn't create lookup table! Terminating.", msg );
-      system->error_ptr->one(FLERR, errmsg);
-    }
-  }
-
-
-  Init_Force_FunctionsOMP( control );
+  Init_Force_FunctionsOMP(control);
 }
 
diff --git a/src/USER-REAXC/reaxc_init_md.cpp b/src/USER-REAXC/reaxc_init_md.cpp
index 73210ef029..66585dad4a 100644
--- a/src/USER-REAXC/reaxc_init_md.cpp
+++ b/src/USER-REAXC/reaxc_init_md.cpp
@@ -39,8 +39,9 @@
 #include "reaxc_tool_box.h"
 
 #include "error.h"
+#include "fmt/format.h"
 
-int Init_System( reax_system *system, control_params *control, char * /*msg*/ )
+int Init_System(reax_system *system, control_params *control, char * /*msg*/)
 {
   int i;
   reax_atom *atom;
@@ -51,32 +52,32 @@ int Init_System( reax_system *system, control_params *control, char * /*msg*/ )
 
   // determine the local and total capacity
 
-  system->local_cap = MAX( (int)(system->n * safezone), mincap);
-  system->total_cap = MAX( (int)(system->N * safezone), mincap);
+  system->local_cap = MAX((int)(system->n * safezone), mincap);
+  system->total_cap = MAX((int)(system->N * safezone), mincap);
 
   /* estimate numH and Hcap */
   system->numH = 0;
   if (control->hbond_cut > 0)
-    for( i = 0; i < system->n; ++i ) {
+    for(i = 0; i < system->n; ++i) {
       atom = &(system->my_atoms[i]);
       if (system->reax_param.sbp[ atom->type ].p_hbond == 1 && atom->type >= 0)
         atom->Hindex = system->numH++;
       else atom->Hindex = -1;
     }
-  system->Hcap = (int)(MAX( system->numH * saferzone, mincap ));
+  system->Hcap = (int)(MAX(system->numH * saferzone, mincap));
 
   return SUCCESS;
 }
 
 
-int Init_Simulation_Data( reax_system *system, control_params *control,
-                          simulation_data *data, char * /*msg*/ )
+int Init_Simulation_Data(reax_system *system, control_params *control,
+                          simulation_data *data, char * /*msg*/)
 {
-  Reset_Simulation_Data( data, control->virial );
+  Reset_Simulation_Data(data, control->virial);
 
   /* initialize the timer(s) */
   if (system->my_rank == MASTER_NODE) {
-    data->timing.start = Get_Time( );
+    data->timing.start = Get_Time();
   }
 
   data->step = data->prev_steps = 0;
@@ -84,81 +85,82 @@ int Init_Simulation_Data( reax_system *system, control_params *control,
   return SUCCESS;
 }
 
-void Init_Taper( control_params *control,  storage *workspace )
+void Init_Taper(control_params *control,  storage *workspace)
 {
   double d1, d7;
   double swa, swa2, swa3;
   double swb, swb2, swb3;
+  LAMMPS_NS::Error *error = control->error_ptr;
 
   swa = control->nonb_low;
   swb = control->nonb_cut;
 
-  if (fabs( swa ) > 0.01 && control->me == 0)
-    control->error_ptr->warning( FLERR, "Non-zero lower Taper-radius cutoff" );
+  if (fabs(swa) > 0.01 && control->me == 0)
+    error->warning(FLERR, "Non-zero lower Taper-radius cutoff");
 
   if (swb < 0) {
-    control->error_ptr->all(FLERR,"Negative upper Taper-radius cutoff");
+    error->all(FLERR,"Negative upper Taper-radius cutoff");
   }
-  else if( swb < 5 && control->me == 0) {
+  else if(swb < 5 && control->me == 0) {
     char errmsg[256];
-    snprintf(errmsg, 256, "Very low Taper-radius cutoff: %f", swb );
-    control->error_ptr->warning( FLERR, errmsg );
+    snprintf(errmsg, 256, "Very low Taper-radius cutoff: %f", swb);
+    error->warning(FLERR, errmsg);
   }
 
   d1 = swb - swa;
-  d7 = pow( d1, 7.0 );
-  swa2 = SQR( swa );
-  swa3 = CUBE( swa );
-  swb2 = SQR( swb );
-  swb3 = CUBE( swb );
+  d7 = pow(d1, 7.0);
+  swa2 = SQR(swa);
+  swa3 = CUBE(swa);
+  swb2 = SQR(swb);
+  swb3 = CUBE(swb);
 
   workspace->Tap[7] =  20.0 / d7;
   workspace->Tap[6] = -70.0 * (swa + swb) / d7;
   workspace->Tap[5] =  84.0 * (swa2 + 3.0*swa*swb + swb2) / d7;
-  workspace->Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3 ) / d7;
-  workspace->Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3 ) / d7;
+  workspace->Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3) / d7;
+  workspace->Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3) / d7;
   workspace->Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7;
   workspace->Tap[1] = 140.0 * swa3 * swb3 / d7;
   workspace->Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 -
-                     7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7;
+                     7.0*swa*swb3*swb3 + swb3*swb3*swb) / d7;
 }
 
 
-int Init_Workspace( reax_system *system, control_params *control,
-                    storage *workspace, char *msg )
+int Init_Workspace(reax_system *system, control_params *control,
+                    storage *workspace, char *msg)
 {
   int ret;
 
-  ret = Allocate_Workspace( system, control, workspace,
-                            system->local_cap, system->total_cap, msg );
+  ret = Allocate_Workspace(system, control, workspace,
+                            system->local_cap, system->total_cap, msg);
   if (ret != SUCCESS)
     return ret;
 
-  memset( &(workspace->realloc), 0, sizeof(reallocate_data) );
-  Reset_Workspace( system, workspace );
+  memset(&(workspace->realloc), 0, sizeof(reallocate_data));
+  Reset_Workspace(system, workspace);
 
   /* Initialize the Taper function */
-  Init_Taper( control, workspace);
+  Init_Taper(control, workspace);
 
   return SUCCESS;
 }
 
 
 /************** setup communication data structures  **************/
-int Init_MPI_Datatypes( reax_system *system, storage * /*workspace*/,
-                        mpi_datatypes *mpi_data, MPI_Comm comm, char * /*msg*/ )
+int Init_MPI_Datatypes(reax_system *system, storage * /*workspace*/,
+                        mpi_datatypes *mpi_data, MPI_Comm comm, char * /*msg*/)
 {
 
   /* setup the world */
   mpi_data->world = comm;
-  MPI_Comm_size( comm, &(system->wsize) );
+  MPI_Comm_size(comm, &(system->wsize));
 
   return SUCCESS;
 }
 
-int  Init_Lists( reax_system *system, control_params *control,
+int  Init_Lists(reax_system *system, control_params *control,
                  simulation_data * /*data*/, storage * /*workspace*/, reax_list **lists,
-                 mpi_datatypes * /*mpi_data*/, char * /*msg*/ )
+                 mpi_datatypes * /*mpi_data*/, char * /*msg*/)
 {
   int i, total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop;
   int *hb_top, *bond_top;
@@ -166,94 +168,93 @@ int  Init_Lists( reax_system *system, control_params *control,
   int mincap = system->mincap;
   double safezone = system->safezone;
   double saferzone = system->saferzone;
+  LAMMPS_NS::Error *error = system->error_ptr;
 
-  bond_top = (int*) calloc( system->total_cap, sizeof(int) );
-  hb_top = (int*) calloc( system->local_cap, sizeof(int) );
-  Estimate_Storages( system, control, lists,
+  bond_top = (int*) calloc(system->total_cap, sizeof(int));
+  hb_top = (int*) calloc(system->local_cap, sizeof(int));
+  Estimate_Storages(system, control, lists,
                      &Htop, hb_top, bond_top, &num_3body);
 
   if (control->hbond_cut > 0) {
     /* init H indexes */
     total_hbonds = 0;
-    for( i = 0; i < system->n; ++i ) {
+    for(i = 0; i < system->n; ++i) {
       system->my_atoms[i].num_hbonds = hb_top[i];
       total_hbonds += hb_top[i];
     }
     total_hbonds = (int)(MAX(total_hbonds*saferzone,mincap*system->minhbonds));
 
-    if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
-                    *lists+HBONDS ) ) {
-      control->error_ptr->one(FLERR, "Not enough space for hbonds list.");
-    }
+    if(!Make_List(system->Hcap, total_hbonds, TYP_HBOND,
+                    *lists+HBONDS))
+      error->one(FLERR, "Not enough space for hbonds list.");
+
     (*lists+HBONDS)->error_ptr = system->error_ptr;
   }
 
   total_bonds = 0;
-  for( i = 0; i < system->N; ++i ) {
+  for(i = 0; i < system->N; ++i) {
     system->my_atoms[i].num_bonds = bond_top[i];
     total_bonds += bond_top[i];
   }
-  bond_cap = (int)(MAX( total_bonds*safezone, mincap*MIN_BONDS ));
+  bond_cap = (int)(MAX(total_bonds*safezone, mincap*MIN_BONDS));
+
+  if(!Make_List(system->total_cap, bond_cap, TYP_BOND,
+                  *lists+BONDS))
+    error->one(FLERR, "Not enough space for bonds list.");
 
-  if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
-                  *lists+BONDS ) ) {
-    control->error_ptr->one(FLERR, "Not enough space for bonds list.");
-  }
   (*lists+BONDS)->error_ptr = system->error_ptr;
 
   /* 3bodies list */
-  cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
-  if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
-                  *lists+THREE_BODIES ) ){
-    control->error_ptr->one(FLERR,"Problem in initializing angles list.");
-  }
+  cap_3body = (int)(MAX(num_3body*safezone, MIN_3BODIES));
+  if(!Make_List(bond_cap, cap_3body, TYP_THREE_BODY,
+                  *lists+THREE_BODIES))
+    error->one(FLERR,"Problem in initializing angles list.");
+
   (*lists+THREE_BODIES)->error_ptr = system->error_ptr;
 
-  free( hb_top );
-  free( bond_top );
+  free(hb_top);
+  free(bond_top);
 
   return SUCCESS;
 }
 
-void Initialize( reax_system *system, control_params *control,
+void Initialize(reax_system *system, control_params *control,
                  simulation_data *data, storage *workspace,
                  reax_list **lists, output_controls *out_control,
-                 mpi_datatypes *mpi_data, MPI_Comm comm )
+                 mpi_datatypes *mpi_data, MPI_Comm comm)
 {
   char msg[MAX_STR];
+  LAMMPS_NS::Error *error = system->error_ptr;
 
-  if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
-    control->error_ptr->one(FLERR,"Could not create datatypes");
-  }
+  if (Init_MPI_Datatypes(system,workspace,mpi_data,comm,msg) == FAILURE)
+    error->one(FLERR,"init_mpi_datatypes: could not create datatypes. "
+               "Mpi_data could not be initialized! Terminating.");
 
-  if (Init_System(system, control, msg) == FAILURE) {
-    control->error_ptr->one(FLERR,"System could not be initialized");
-  }
+  if (Init_System(system,control,msg) == FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. System could not be "
+                                  "initialized! Terminating.",msg));
 
-  if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
-    control->error_ptr->one(FLERR,"Sim_data could not be initialized");
-  }
+  if (Init_Simulation_Data( system,control,data,msg) == FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. Sim_data could not be "
+                                  "initialized! Terminating.",msg));
 
-  if (Init_Workspace( system, control, workspace, msg ) ==
-      FAILURE) {
-    control->error_ptr->one(FLERR,"Workspace could not be initialized");
-  }
+  if (Init_Workspace( system,control,workspace,msg) == FAILURE)
+    error->one(FLERR,"init_workspace: not enough memory. "
+               "Workspace could not be initialized. Terminating.");
 
-  if (Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
-      FAILURE) {
-    control->error_ptr->one(FLERR,"Lists could not be initialized");
-    }
+  if (Init_Lists(system, control, data, workspace, lists, mpi_data, msg) ==FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. System could not be "
+                                  "initialized. Terminating.",msg));
 
-  if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
-    control->error_ptr->one(FLERR,"Could not open output files");
-  }
+  if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE)
+    error->one(FLERR,fmt::format("Error on: {}. Could not open output files! "
+                                  "Terminating.",msg));
 
-  if (control->tabulate) {
-    if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
-    control->error_ptr->one(FLERR,"Lookup table could not be created");
-    }
-  }
+  if (control->tabulate)
+    if (Init_Lookup_Tables(system,control,workspace,mpi_data,msg) == FAILURE)
+      error->one(FLERR,fmt::format("Error on: {}. Could not create lookup "
+                                    "table. Terminating.",msg));
 
 
-  Init_Force_Functions( control );
+  Init_Force_Functions(control);
 }
diff --git a/src/balance.cpp b/src/balance.cpp
index 8c1e7a4ed7..49cd6ecad5 100644
--- a/src/balance.cpp
+++ b/src/balance.cpp
@@ -113,7 +113,7 @@ void Balance::command(int narg, char **arg)
   if (domain->box_exist == 0)
     error->all(FLERR,"Balance command before simulation box is defined");
 
-  if (me == 0 && screen) fprintf(screen,"Balancing ...\n");
+  if (me == 0) utils::logmesg(lmp,"Balancing ...\n");
 
   // parse required arguments
 
@@ -573,8 +573,8 @@ int *Balance::bisection(int sortflag)
   int triclinic = domain->triclinic;
 
   double *boxlo,*boxhi,*prd;
-  
-  if (triclinic == 0) { 
+
+  if (triclinic == 0) {
     boxlo = domain->boxlo;
     boxhi = domain->boxhi;
     prd = domain->prd;
@@ -587,7 +587,7 @@ int *Balance::bisection(int sortflag)
   // shrink-wrap simulation box around atoms for input to RCB
   // leads to better-shaped sub-boxes when atoms are far from box boundaries
   // if triclinic, do this in lamda coords
-  
+
   double shrink[6],shrinkall[6];
 
   shrink[0] = boxhi[0]; shrink[1] = boxhi[1]; shrink[2] = boxhi[2];
@@ -597,7 +597,7 @@ int *Balance::bisection(int sortflag)
   int nlocal = atom->nlocal;
 
   if (triclinic) domain->x2lamda(nlocal);
-  
+
   for (int i = 0; i < nlocal; i++) {
     shrink[0] = MIN(shrink[0],x[i][0]);
     shrink[1] = MIN(shrink[1],x[i][1]);
diff --git a/src/change_box.cpp b/src/change_box.cpp
index f50ba23db4..ff7658a01f 100644
--- a/src/change_box.cpp
+++ b/src/change_box.cpp
@@ -49,7 +49,7 @@ void ChangeBox::command(int narg, char **arg)
     error->all(FLERR,"Change_box command before simulation box is defined");
   if (narg < 2) error->all(FLERR,"Illegal change_box command");
 
-  if (comm->me == 0 && screen) fprintf(screen,"Changing box ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"Changing box ...\n");
 
   // group
 
diff --git a/src/comm.h b/src/comm.h
index d308a26367..1743436e22 100644
--- a/src/comm.h
+++ b/src/comm.h
@@ -253,6 +253,22 @@ be generated and LAMMPS may lose atoms or use incorrect periodic images of atoms
 interaction lists.  To avoid, either define pair style zero with a suitable cutoff
 or use comm_modify cutoff.
 
+W: Communication cutoff is shorter than a bond length based estimate. This may lead to errors.
+
+Since LAMMPS stores topology data with individual atoms, all atoms comprising
+a bond, angle, dihedral or improper must be present on any sub-domain that
+"owns" the atom with the information, either as a local or a ghost atom. The
+communication cutoff is what determines up to what distance from a sub-domain
+boundary ghost atoms are created.  The communication cutoff is by default the
+largest non-bonded cutoff plus the neighbor skin distance, but for short or
+non-bonded cutoffs and/or long bonds, this may not be sufficient. This warning
+indicates that there is an increased risk of a simulation stopping unexpectedly
+because of Bond/Angle/Dihedral/Improper atoms missing.  It can be silenced by
+manually setting the communication cutoff via comm_modify cutoff.  However,
+since the heuristic used to determine the estimate is not always accurate, it
+is not changed automatically and the warning may be ignored depending on the
+specific system being simulated.
+
 UNDOCUMENTED
 
 U: OMP_NUM_THREADS environment is not set.
diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp
index de02823ec8..b7878ea86a 100644
--- a/src/comm_tiled.cpp
+++ b/src/comm_tiled.cpp
@@ -98,7 +98,7 @@ void CommTiled::init()
   Comm::init();
 
   nswap = 2 * domain->dimension;
-  
+
   int bufextra_old = bufextra;
   init_exchange();
   if (bufextra > bufextra_old) grow_send(maxsend+bufextra,2);
@@ -192,7 +192,7 @@ void CommTiled::setup()
     cut = MIN(prd[0],prd[1]);
     if (dimension == 3) cut = MIN(cut,prd[2]);
     cut *= EPSILON*EPSILON;
-    cutghost[0] = cutghost[1] = cutghost[2] = cut; 
+    cutghost[0] = cutghost[1] = cutghost[2] = cut;
   }
 
   // setup forward/reverse communication
@@ -325,16 +325,16 @@ void CommTiled::setup()
           sbox[3] = MIN(oboxhi[0],hi1[0]);
           sbox[4] = MIN(oboxhi[1],hi1[1]);
           sbox[5] = MIN(oboxhi[2],hi1[2]);
-	  
+
         } else {
           pbc_flag[iswap][i] = 1;
           if (idir == 0) pbc[iswap][i][idim] = 1;
           else pbc[iswap][i][idim] = -1;
-	  if (triclinic) {
+          if (triclinic) {
             if (idim == 1) pbc[iswap][i][5] = pbc[iswap][i][idim];
             if (idim == 2) pbc[iswap][i][4] = pbc[iswap][i][3] = pbc[iswap][i][idim];
           }
-	  
+
           sbox[0] = MAX(oboxlo[0],lo2[0]);
           sbox[1] = MAX(oboxlo[1],lo2[1]);
           sbox[2] = MAX(oboxlo[2],lo2[2]);
@@ -346,7 +346,7 @@ void CommTiled::setup()
         if (idir == 0) {
           sbox[idim] = sublo[idim];
           if (i < noverlap1)
-	    sbox[3+idim] = MIN(sbox[3+idim]+cutghost[idim],subhi[idim]);
+            sbox[3+idim] = MIN(sbox[3+idim]+cutghost[idim],subhi[idim]);
           else sbox[3+idim] = MIN(sbox[3+idim]-prd[idim]+cutghost[idim],subhi[idim]);
         } else {
           if (i < noverlap1) sbox[idim] = MAX(sbox[idim]-cutghost[idim],sublo[idim]);
@@ -370,7 +370,7 @@ void CommTiled::setup()
     }
   }
 
-  
+
   // setup exchange communication = subset of forward/reverse comm procs
   // loop over dimensions
   // determine which procs I will exchange with in each dimension
@@ -717,10 +717,10 @@ void CommTiled::exchange()
           buf_send[nsend++] = proc;
           nsend += avec->pack_exchange(i,&buf_send[nsend]);
         } else {
-	  // DEBUG statment
-	  // error->warning(FLERR,"Losing atom in CommTiled::exchange() send, "
-	  // 		    "likely bad dynamics");
-	}
+          // DEBUG statment
+          // error->warning(FLERR,"Losing atom in CommTiled::exchange() send, "
+          //                "likely bad dynamics");
+        }
         avec->copy(nlocal-1,i,1);
         nlocal--;
       } else i++;
@@ -769,10 +769,10 @@ void CommTiled::exchange()
         if (value >= lo && value < hi) {
           m += avec->unpack_exchange(&buf_recv[m]);
           continue;
-	} else {
-	  // DEBUG statment
-	  // error->warning(FLERR,"Losing atom in CommTiled::exchange() recv");
-	}
+        } else {
+          // DEBUG statment
+          // error->warning(FLERR,"Losing atom in CommTiled::exchange() recv");
+        }
       }
       m += static_cast<int> (buf_recv[m]);
     }
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index ba8940de95..03d5d2f65d 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -282,8 +282,7 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
   int group1bit = group->bitmask[igroup1];
   int group2bit = group->bitmask[igroup2];
 
-  if (comm->me == 0 && screen)
-    fprintf(screen,"System init for delete_atoms ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"System init for delete_atoms ...\n");
 
   // request a full neighbor list for use by this command
 
diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp
index 99c0f1b951..7b9260bc52 100644
--- a/src/delete_bonds.cpp
+++ b/src/delete_bonds.cpp
@@ -51,11 +51,10 @@ void DeleteBonds::command(int narg, char **arg)
   // init entire system since comm->borders is done
   // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
 
-  if (comm->me == 0 && screen)
-    fprintf(screen,"System init for delete_bonds ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"System init for delete_bonds ...\n");
   lmp->init();
 
-  if (comm->me == 0 && screen) fprintf(screen,"Deleting bonds ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"Deleting bonds ...\n");
 
   // identify group
 
diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp
index 748adca811..82b7bf9492 100644
--- a/src/displace_atoms.cpp
+++ b/src/displace_atoms.cpp
@@ -34,6 +34,7 @@
 #include "math_extra.h"
 #include "memory.h"
 #include "error.h"
+#include "utils.h"
 #include "fmt/format.h"
 
 using namespace LAMMPS_NS;
@@ -68,7 +69,7 @@ void DisplaceAtoms::command(int narg, char **arg)
     error->all(FLERR,"Cannot displace_atoms after "
                "reading restart file with per-atom info");
 
-  if (comm->me == 0 && screen) fprintf(screen,"Displacing atoms ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"Displacing atoms ...\n");
 
   // group and style
 
diff --git a/src/read_dump.cpp b/src/read_dump.cpp
index 0426323a21..f47eda3262 100644
--- a/src/read_dump.cpp
+++ b/src/read_dump.cpp
@@ -120,7 +120,7 @@ void ReadDump::command(int narg, char **arg)
 
   // find the snapshot and read/bcast/process header info
 
-  if (me == 0 && screen) fprintf(screen,"Scanning dump file ...\n");
+  if (me == 0) utils::logmesg(lmp,"Scanning dump file ...\n");
 
   bigint ntimestep = seek(nstep,1);
   if (ntimestep < 0)
@@ -135,8 +135,7 @@ void ReadDump::command(int narg, char **arg)
 
   // read in the snapshot and reset system
 
-  if (me == 0 && screen)
-    fprintf(screen,"Reading snapshot from dump file ...\n");
+  if (me == 0) utils::logmesg(lmp,"Reading snapshot from dump file ...\n");
 
   bigint natoms_prev = atom->natoms;
   atoms();
diff --git a/src/read_restart.cpp b/src/read_restart.cpp
index 3bbae596d1..157992bce3 100644
--- a/src/read_restart.cpp
+++ b/src/read_restart.cpp
@@ -103,7 +103,7 @@ void ReadRestart::command(int narg, char **arg)
   // open single restart file or base file for multiproc case
 
   if (me == 0) {
-    if (screen) fprintf(screen,"Reading restart file ...\n");
+    utils::logmesg(lmp,"Reading restart file ...\n");
     std::string hfile = file;
     if (multiproc) {
       hfile.replace(hfile.find("%"),1,"base");
@@ -620,10 +620,9 @@ void ReadRestart::header()
 
     if (flag == VERSION) {
       char *version = read_string();
-      if (me == 0) {
-        if (screen) fprintf(screen,"  restart file = %s, LAMMPS = %s\n",
-                            version,universe->version);
-      }
+      if (me == 0)
+        utils::logmesg(lmp,fmt::format("  restart file = {}, LAMMPS = {}\n",
+                                       version,universe->version));
       delete [] version;
 
       // we have no forward compatibility, thus exit with error
diff --git a/src/replicate.cpp b/src/replicate.cpp
index 51ef75c37b..a9c0dbf7b3 100644
--- a/src/replicate.cpp
+++ b/src/replicate.cpp
@@ -48,7 +48,7 @@ void Replicate::command(int narg, char **arg)
   int me = comm->me;
   int nprocs = comm->nprocs;
 
-  if (me == 0 && screen) fprintf(screen,"Replicating atoms ...\n");
+  if (me == 0) utils::logmesg(lmp,"Replicating atoms ...\n");
 
   // nrep = total # of replications
 
@@ -349,12 +349,14 @@ void Replicate::command(int narg, char **arg)
     int size_buf_all = 0;
     MPI_Allreduce(&n, &size_buf_all, 1, MPI_INT, MPI_SUM, world);
 
-    if (me == 0 && screen) {
-      fmt::print(screen,"  bounding box image = ({} {} {}) to ({} {} {})\n",
-              _imagelo[0],_imagelo[1],_imagelo[2],
-              _imagehi[0],_imagehi[1],_imagehi[2]);
-      fmt::print(screen,"  bounding box extra memory = {:.2f} MB\n",
-              (double)size_buf_all*sizeof(double)/1024/1024);
+    if (me == 0) {
+      auto mesg = fmt::format("  bounding box image = ({:.8} {:.8} {:.8}) "
+                              "to ({:.8} {:.8} {:.8})\n",
+                              _imagelo[0],_imagelo[1],_imagelo[2],
+                              _imagehi[0],_imagehi[1],_imagehi[2]);
+      mesg += fmt::format("  bounding box extra memory = {:.2f} MB\n",
+                          (double)size_buf_all*sizeof(double)/1024/1024);
+      utils::logmesg(lmp,mesg);
     }
 
     // rnk offsets
@@ -632,11 +634,10 @@ void Replicate::command(int narg, char **arg)
     int sum = 0;
     MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world);
     double avg = (double) sum / nprocs;
-    if (me == 0 && screen)
-      fprintf(screen,"  average # of replicas added to proc = %.2f "
-              "out of %i (%.2f %%)\n",
-              avg,nx*ny*nz,avg/(nx*ny*nz)*100.0);
-
+    if (me == 0)
+      utils::logmesg(lmp,fmt::format("  average # of replicas added to proc ="
+                                     " {:.2f} out of {} ({:.2f}%)\n",
+                                     avg,nx*ny*nz,avg/(nx*ny*nz)*100.0));
   } else {
 
     for (int iproc = 0; iproc < nprocs; iproc++) {
diff --git a/src/set.cpp b/src/set.cpp
index 8527f93df7..01176f68c0 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -82,7 +82,7 @@ void Set::command(int narg, char **arg)
   // loop over keyword/value pairs
   // call appropriate routine to reset attributes
 
-  if (comm->me == 0 && screen) fprintf(screen,"Setting atom values ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"Setting atom values ...\n");
 
   int allcount,origarg;
 
@@ -1264,7 +1264,7 @@ void Set::topology(int keyword)
   // init entire system since comm->exchange is done
   // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
 
-  if (comm->me == 0 && screen) fprintf(screen,"  system init for set ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"  system init for set ...\n");
   lmp->init();
 
   if (domain->triclinic) domain->x2lamda(atom->nlocal);
diff --git a/src/utils.cpp b/src/utils.cpp
index ce6a64f22f..3369f684ba 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -82,7 +82,7 @@ using namespace LAMMPS_NS;
  *  even with char * type variables.
  *  Example: utils::strmatch(text, std::string("^") + charptr)
  */
-bool utils::strmatch(std::string text, std::string pattern)
+bool utils::strmatch(const std::string &text, const std::string &pattern)
 {
   const int pos = re_match(text.c_str(),pattern.c_str());
   return (pos >= 0);
diff --git a/src/utils.h b/src/utils.h
index 3919f8d4d5..9762087ca4 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -35,7 +35,7 @@ namespace LAMMPS_NS {
      *  \param pattern the search pattern, which may contain regexp markers
      *  \return true if the pattern matches, false if not
      */
-    bool strmatch(std::string text, std::string pattern);
+    bool strmatch(const std::string &text, const std::string &pattern);
 
     /** \brief Send message to screen and logfile, if available
      *
@@ -155,7 +155,7 @@ namespace LAMMPS_NS {
      * \param line string that should be trimmed
      * \return new string without comment (string)
      */
-    std::string trim_comment(const std::string & line);
+    std::string trim_comment(const std::string &line);
 
     /**
      * \brief Count words in string
diff --git a/src/write_data.cpp b/src/write_data.cpp
index 1ec24a99cd..511552cf69 100644
--- a/src/write_data.cpp
+++ b/src/write_data.cpp
@@ -105,8 +105,7 @@ void WriteData::command(int narg, char **arg)
   //     e.g. pair hybrid coeffs, dpd ghost-atom velocity setting
 
   if (noinit == 0) {
-    if (comm->me == 0 && screen)
-      fputs("System init for write_data ...\n",screen);
+    if (comm->me == 0) utils::logmesg(lmp,"System init for write_data ...\n");
     lmp->init();
 
     // move atoms to new processors before writing file
diff --git a/src/write_restart.cpp b/src/write_restart.cpp
index 9129126aae..062add3b3c 100644
--- a/src/write_restart.cpp
+++ b/src/write_restart.cpp
@@ -87,8 +87,7 @@ void WriteRestart::command(int narg, char **arg)
   // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
 
   if (noinit == 0) {
-    if (comm->me == 0 && screen)
-      fprintf(screen,"System init for write_restart ...\n");
+    if (comm->me == 0) utils::logmesg(lmp,"System init for write_restart ...\n");
     lmp->init();
 
     // move atoms to new processors before writing file
diff --git a/unittest/force-styles/tests/atomic-pair-born.yaml b/unittest/force-styles/tests/atomic-pair-born.yaml
new file mode 100644
index 0000000000..89c43d44e9
--- /dev/null
+++ b/unittest/force-styles/tests/atomic-pair-born.yaml
@@ -0,0 +1,89 @@
+---
+lammps_version: 21 Jul 2020
+date_generated: Thu Jul 30 03:30:06 202
+epsilon: 5e-13
+prerequisites: ! |
+  pair born
+pre_commands: ! ""
+post_commands: ! ""
+input_file: in.metal
+pair_style: born 8.0
+pair_coeff: ! |
+  * * 6.08 0.317 2.340 24.18 11.51
+extract: ! ""
+natoms: 32
+init_vdwl: 836.893711811487
+init_coul: 0
+init_stress: ! |2-
+   2.1907762159508770e+03  2.1533146615415390e+03  2.1384803888060319e+03  3.1458958046100697e+01 -5.2331768978302460e+00 -1.2410855173406260e+01
+init_forces: ! |2
+    1  6.1952278656609625e+00  1.6628887543469570e+01  2.3406241273969628e+00
+    2 -7.1498478250272974e+00 -1.9650729468423989e+01  1.3392331597463706e+01
+    3 -8.4302759991532952e+00  4.5935465074977714e+01 -1.1339628764311961e+01
+    4  6.6013347699419995e+00 -2.2770126898611998e+00 -1.4793531675390774e+01
+    5  8.8132464681297229e+00  3.2091565756514420e+01  1.0510805603350610e+01
+    6  1.0484038408463748e+01  9.4085458810652689e+00 -1.1561120732403498e+00
+    7 -1.4142880457618997e+01  2.1179702223374814e+01 -1.6277650983857441e+00
+    8 -1.6212550653790164e+01 -1.5776737518534750e+01  8.2546279217585035e+00
+    9  2.2842749772666266e+00 -1.5405998927734670e+01 -6.1169286615136649e+00
+   10 -5.6088561559285921e+00 -2.9980200955395699e+01 -2.8361114533155039e+01
+   11 -1.0796936404967482e+01 -1.4085775254255909e+01  1.0322820000470387e+01
+   12  4.9935473695250936e+00 -2.0508984017848242e+01 -3.9212800338387188e-01
+   13 -3.8342965060925644e+00  5.8453344255223003e+00 -1.1781310586322938e+01
+   14  1.2960559403161112e+01 -5.7879190752621907e+00  1.6324972340830911e+01
+   15 -7.5631698457336611e+00  7.5962132139231553e+00 -7.5008652746855180e+00
+   16  2.7227909176641056e+01 -2.2159697717922022e+01  9.0018115646343890e+00
+   17  9.6743623705665449e+00  3.1763667293227797e+01  1.7437438229370321e+00
+   18 -3.3175129754412701e+01 -5.7620879906320663e+00 -1.6758037406277673e+01
+   19 -3.4363169461589663e+01 -1.6361294228706516e+01  9.5257528857146188e+00
+   20  6.3606452908246656e+00  4.8320683138296658e+00 -7.1004418004298362e+00
+   21  3.7194448979025765e+01  1.0842684015015067e+01  3.9593216918575974e-01
+   22  4.0813085193467058e-01 -2.1849731691352865e+01 -1.5393240154689587e+00
+   23  9.0194539312342243e+00  2.7646519773287814e+01 -1.8893395925821068e+01
+   24  2.7633552340313585e+01  2.1176334996344043e+01  1.9050374400981202e+01
+   25  1.8215010599705714e+01  7.8343723407555865e-01 -1.1516192459397450e+01
+   26 -9.6728258529328315e+00 -3.3921248222992517e+01  2.9174727715069260e+01
+   27  1.0237099691004113e+01  1.1774417628369038e+00  9.0059613217903767e-01
+   28  5.7207612737003188e+00  1.2091973176094161e+01  1.6885831883966489e+01
+   29  8.0920152886696819e+00  2.3976025150348377e+01 -1.0543989422319544e+00
+   30 -1.2005631838887281e+01 -6.4092999695138637e+00 -1.9850411539058847e+01
+   31 -5.5566952298732328e+00 -2.8420402803389738e+01  1.7818115115845121e+01
+   32 -4.3603353069761802e+01 -1.4618745302080386e+01 -5.8614805227083417e+00
+run_vdwl: 808.581991307929
+run_coul: 0
+run_stress: ! |2-
+   2.1352577876164714e+03  2.1048336346087822e+03  2.0866093207888316e+03  1.6710708502154183e+01 -7.2805172515122214e+00 -5.1007546794256298e+00
+run_forces: ! |2
+    1  4.3996728759831116e+00  1.3617808952853350e+01  3.9365269387608910e+00
+    2 -5.8870109187129227e+00 -1.6653786974892565e+01  1.0458442718468866e+01
+    3 -2.5384607642764081e+00  2.7976453386304165e+01 -7.9538384346985778e+00
+    4  2.6011652570379979e+00  3.4444982822291809e-01 -1.4078044299688214e+01
+    5  6.2739400167978490e+00  2.3321353952206238e+01  9.8524679762955127e+00
+    6  9.7615972917437634e+00  9.8314066166309821e+00  9.1579317972042351e-01
+    7 -8.4994423410616129e+00  1.7133945132359301e+01 -2.1335415687539081e+00
+    8 -1.4401879593786829e+01 -1.0717410919398073e+01  6.2772265942246497e+00
+    9  2.4581092844990748e+00 -1.5139809581001368e+01 -6.0486224458444502e+00
+   10  5.7295535185278368e-01 -2.2473696014663325e+01 -1.8469640259236243e+01
+   11 -8.3438231201666824e+00 -1.1773385047512056e+01  8.0911869578298905e+00
+   12  5.6018202198927955e+00 -1.3190833667603776e+01  3.2425752936252150e-01
+   13 -5.1665403937031575e+00  3.8641569820964516e+00 -8.0278457568777704e+00
+   14  1.1059028830337486e+01 -3.5867141219142817e-01  8.3359673389637585e+00
+   15 -7.0225257033738311e+00  3.4315031664638429e+00 -5.9505539807290697e+00
+   16  1.5086691796848593e+01 -1.7105542466895855e+01  6.8923022058483863e+00
+   17  7.7552246185118392e+00  2.3712928096926632e+01 -1.1565629751691666e+00
+   18 -2.2415417000611026e+01 -5.6849910964810348e+00 -1.0167777125296530e+01
+   19 -2.2914820778853272e+01 -9.2271207471412673e+00  6.1677237945523231e+00
+   20  6.5866441066676380e+00  6.2266458578595580e+00 -6.7918108679194686e+00
+   21  2.5378769975471076e+01  6.5502175694342872e+00  2.7113928300966954e+00
+   22  1.4722121850540351e+00 -1.6174770388907241e+01 -3.6146790093186110e+00
+   23  3.9177677709068743e+00  1.6640247776368867e+01 -1.5217886244976798e+01
+   24  1.6292888686152281e+01  1.6539457853286358e+01  1.3246494291709737e+01
+   25  1.2399183978889422e+01  2.2450939828300065e+00 -8.4542270100046313e+00
+   26 -5.5487323330521612e+00 -2.0228255850187249e+01  1.8528693239502584e+01
+   27  7.5279483917695700e+00 -3.0072122179470506e+00  6.9304022265424459e-01
+   28 -3.8508457491974468e-01  9.5889969664907753e+00  1.3098728886480401e+01
+   29  5.7371217847825058e+00  1.5830079463104170e+01 -9.7700534543667272e-01
+   30 -7.4069396659644475e+00 -5.7323481111133043e+00 -1.2021778057070120e+01
+   31 -3.9192345419264893e+00 -2.1509286285601934e+01  1.3203747851309895e+01
+   32 -3.0432830692790130e+01 -7.8776248019003354e+00 -1.6701791747605144e+00
+...
diff --git a/unittest/force-styles/tests/mol-pair-lennard_mdf.yaml b/unittest/force-styles/tests/mol-pair-lennard_mdf.yaml
new file mode 100644
index 0000000000..caa4e04edc
--- /dev/null
+++ b/unittest/force-styles/tests/mol-pair-lennard_mdf.yaml
@@ -0,0 +1,99 @@
+---
+lammps_version: 21 Jul 2020
+date_generated: Thu Jul 30 11:43:01 202
+epsilon: 5e-14
+prerequisites: ! |
+  atom full
+  pair lennard/mdf
+pre_commands: ! ""
+post_commands: ! ""
+input_file: in.fourmol
+pair_style: lennard/mdf 6.5 8.0
+pair_coeff: ! |
+  1 1 4768.37158203125 19.53125
+  1 2 33.0002002739906 1.148916015625
+  1 3 22973.5604278016 42.87055906125
+  1 4 16088.5533432005 33.3863396081125
+  1 5 16088.5533432005 33.3863396081125
+  2 2 0.02 0.02
+  2 3 294.233100455466 3.43064484
+  2 4 4.8828125e-06 0.0003125
+  2 5 190.825740216415 2.57106933958137
+  3 3 92233.7203685478 85.89934592
+  3 4 66122.1953250058 67.6836765582081
+  3 5 66122.1953250058 67.6836765582081
+  4 4 47259.767027313 53.25022086
+  4 5 47259.767027313 53.25022086
+  5 5 47259.767027313 53.25022086
+extract: ! |
+  a 2
+natoms: 29
+init_vdwl: 749.238128769917
+init_coul: 0
+init_stress: ! |2-
+   2.1793844857971999e+03  2.1988910807167222e+03  4.6653956464347129e+03 -7.5956623392527081e+02  2.4752218664716132e+01  6.6652254832700407e+02
+init_forces: ! |2
+    1 -2.3333456445729308e+01  2.6994541715130174e+02  3.3272836403600138e+02
+    2  1.5828554631286457e+02  1.3025008843150334e+02 -1.8629682358090739e+02
+    3 -1.3528903744071818e+02 -3.8704313350789715e+02 -1.4568978426110166e+02
+    4 -7.8711102183152653e+00  2.1350433612475337e+00 -5.5954483189689057e+00
+    5 -2.5176757267276115e+00 -4.0521510680612840e+00  1.2152704057983788e+01
+    6 -8.3190657145323701e+02  9.6394157585967673e+02  1.1509101058588756e+03
+    7  5.8203391152120780e+01 -3.3609018510471799e+02 -1.7179623897115050e+03
+    8  1.4451392646293473e+02 -1.0927476052490431e+02  3.9990594285329564e+02
+    9  7.9156945404948118e+01  8.5273009664329706e+01  3.5032175683823232e+02
+   10  5.3118875274635047e+02 -6.1040990309251686e+02 -1.8355872279721592e+02
+   11 -2.3530151263226644e+00 -5.9077525030323867e+00 -9.6590611210259372e+00
+   12  1.7527040556419166e+01  1.0633031637151417e+01 -7.9253531066396325e+00
+   13  8.0986456357509482e+00 -3.2098047359651849e+00 -1.4896595191878603e-01
+   14 -3.3852718245345126e+00  6.8637005694452435e-01 -8.7507089611400559e+00
+   15 -2.0454759617429094e-01  8.4846189340026168e+00  3.0131599324912615e+00
+   16  4.6326331422454814e+02 -3.3087729037420274e+02 -1.1893030196750653e+03
+   17 -4.5334322511776759e+02  3.1554322933650400e+02  1.2058421140909506e+03
+   18 -1.8856205920908232e-02 -3.3657380846641591e-02  3.1345559228195895e-02
+   19  3.2119183384465589e-04 -2.4442950403900178e-04  1.7531402321548882e-03
+   20 -9.9797590977484174e-04 -1.0244014514626872e-03  3.9741350958193128e-04
+   21 -7.1566239135890754e+01 -8.1615637646874490e+01  2.2589575706544480e+02
+   22 -1.0808841110845918e+02 -2.6193795811011739e+01 -1.6957912582510812e+02
+   23  1.7964463249215564e+02  1.0782103141601063e+02 -5.6305809828670455e+01
+   24  3.6591395075232896e+01 -2.1181605371471514e+02  1.1218279667490484e+02
+   25 -1.4851496073494792e+02  2.3907129094414373e+01 -1.2485640713592730e+02
+   26  1.1191134635354324e+02  1.8789782100385878e+02  1.2650134235513654e+01
+   27  5.1810629413951034e+01 -2.2705426061698427e+02  9.0848881987846823e+01
+   28 -1.8041314536249808e+02  7.7534100857114538e+01 -1.2206963127965679e+02
+   29  1.2861063445049953e+02  1.4952718810862598e+02  3.1216037810340779e+01
+run_vdwl: 146.142938771306
+run_coul: 0
+run_stress: ! |2-
+   6.2607827422370826e+02  6.7185771646959927e+02  4.8484269825432062e+02 -3.1994971119380659e+02 -2.8453599030981216e+01  1.2292331223237372e+02
+run_forces: ! |2
+    1  1.2964000396689146e+01  8.0892267678816864e+01  6.0323237031964297e+01
+    2  1.8550219067108042e+01  1.4016234288534708e+01 -2.2757514465613156e+01
+    3 -2.4823731194481803e+02  1.0508770689878342e+02  8.6172340602494984e+01
+    4 -8.6006937785047128e+00  3.1744103796518282e+00 -6.2092116768831902e+00
+    5 -1.7308376426194816e+00 -2.5575752575563615e+00  1.0273464022561958e+01
+    6  1.1683409017502134e+02 -9.8726850559552091e+01 -1.4208744781635963e+02
+    7  3.8041358692158358e+00 -1.0705007215066624e+01 -4.2399446322832375e+01
+    8 -1.9935799766053112e+01  1.9886787977710259e+01  5.5348128888062327e+01
+    9  1.1205806049421172e+01  1.3198926774269228e+01  4.9879830296866828e+01
+   10  6.4812479379225579e+01 -7.6374990958105826e+01 -7.5660080738756506e+01
+   11 -6.2988692484756150e+00 -8.0952467316475403e+00 -1.6391401801081937e+01
+   12  1.4369329449673979e+01  1.2193215509952948e+01 -8.2527090649331516e+00
+   13  7.3377104327473912e+00 -2.7954976458396401e+00 -1.0299311972585752e-01
+   14 -2.7952420267096518e+00  4.3889411992317034e-01 -7.1100894724883545e+00
+   15 -3.9779622065865505e-01  8.1106030137207288e+00  3.0711098878236167e+00
+   16  5.2787474086118777e+01 -6.7733534176990759e+01  1.5088298239601448e+01
+   17 -1.4633177660642183e+01  1.0014392510695107e+01  4.0799972741104774e+01
+   18 -1.8058839496734813e-02 -3.1846480759107010e-02  2.9996583084227670e-02
+   19  2.5788523138616643e-04 -2.4689081230698301e-04  1.6839267797863844e-03
+   20 -8.8284959462040197e-04 -9.1548924391997993e-04  3.8796476197999974e-04
+   21 -9.9361952309492949e+00 -1.2017775020353746e+01  3.8871115643964693e+01
+   22 -1.8317763935337126e+01 -4.4399886451888442e+00 -2.9440610409402939e+01
+   23  2.8242484527873138e+01  1.6470730701095345e+01 -9.4179941582837472e+00
+   24  1.2413262420490618e+01 -3.6423612222048455e+01  2.3794298924499856e+01
+   25 -3.0031415092604043e+01  5.3796434394870465e+00 -2.5778439006629586e+01
+   26  1.7604732980064579e+01  3.1032191967689151e+01  1.9586447754505298e+00
+   27  3.4727551332751974e+00 -3.8994887375138603e+01  1.1702853254878343e+01
+   28 -2.5472489346144421e+01  1.1754711219486996e+01 -1.7174787246988220e+01
+   29  2.2007795730451438e+01  2.7247258188486931e+01  5.4673625160789685e+00
+...
diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp
index 8e10fa0253..630ff1ba7f 100644
--- a/unittest/formats/test_atom_styles.cpp
+++ b/unittest/formats/test_atom_styles.cpp
@@ -33,6 +33,8 @@
 #if !defined(_FORTIFY_SOURCE) || (_FORTIFY_SOURCE == 0)
 #if defined(__INTEL_COMPILER)
 #define _do_nothing
+#elif defined(__clang__)
+#pragma clang optimize off
 #elif defined(__GNUC__)
 #if (__GNUC__ > 4) || ((__GNUC__ == 4) && (__GNUC_MINOR__ >= 9))
 #pragma GCC optimize("no-var-tracking-assignments", "O0")