Merge branch 'collected-small-changes' into add_set_time

cmake/presets/windows.cmake

@@ -18,7 +18,6 @@ set(WIN_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
-  ELECTRODE
   EFF
   EXTRA-COMPUTE
   EXTRA-DUMP

doc/src/compute_tally.rst

@@ -198,7 +198,9 @@ potentials only include the pair potential portion of the EAM
 interaction when used by this compute, not the embedding term.  Also
 bonded or Kspace interactions do not contribute to this compute.
 
-The computes in this package are not compatible with dynamic groups.
+When used with dynamic groups, a :doc:`run 0 <run>` command needs to
+be inserted in order to initialize the dynamic groups before accessing
+the computes.
 
 Related commands
 """"""""""""""""

lib/atc/LammpsInterface.cpp

@@ -519,8 +519,8 @@ double LammpsInterface::domain_yz() const { return lammps_->domain->yz; }
 int LammpsInterface::domain_triclinic() const { return lammps_->domain->triclinic; }
 
 void LammpsInterface::box_periodicity(int & xperiodic,
-                                          int & yperiodic,
+                                      int & yperiodic,
-                                          int & zperiodic) const
+                                      int & zperiodic) const
 {
   xperiodic = lammps_->domain->xperiodic;
   yperiodic = lammps_->domain->yperiodic;
@@ -546,6 +546,7 @@ bool LammpsInterface::region_bounds(const char * regionName,
   double & xscale, double & yscale, double & zscale) const
 {
   int iRegion = region_id(regionName);
+  if (iRegion < 0) throw ATC_Error("Unknown region " + to_string(regionName));
   xscale = region_xscale(iRegion);
   yscale = region_yscale(iRegion);
   zscale = region_zscale(iRegion);

lib/gpu/lal_device.cpp

@@ -201,9 +201,9 @@ int DeviceT::init_device(MPI_Comm world, MPI_Comm replica, const int ngpu,
     unsigned best_cus = gpu->cus(0);
     bool type_match = (gpu->device_type(0) == type);
     for (int i = 1; i < gpu->num_devices(); i++) {
-      if (type_match && gpu->device_type(i)!=type)
+      if (type_match && (gpu->device_type(i) != type))
         continue;
-      if (type_match && gpu->device_type(i) == type) {
+      if (!type_match && (gpu->device_type(i) == type)) {
         type_match = true;
         best_cus = gpu->cus(i);
         best_device = i;

src/ELECTRODE/pppm_electrode.cpp

@@ -670,7 +670,7 @@ void PPPMElectrode::compute_matrix(bigint *imat, double **matrix, bool timer_fla
   compute(1, 0);
 
   // fft green's function k -> r
-  std::vector<double> greens_real(nz_pppm * ny_pppm * nx_pppm, 0.);
+  std::vector<double> greens_real((std::size_t) nz_pppm * ny_pppm * nx_pppm, 0.0);
   for (int i = 0, n = 0; i < nfft; i++) {
     work2[n++] = greensfn[i];
     work2[n++] = ZEROF;

src/EXTRA-FIX/fix_electron_stopping.cpp

@@ -106,6 +106,7 @@ FixElectronStopping::FixElectronStopping(LAMMPS *lmp, int narg, char **arg) :
 
 FixElectronStopping::~FixElectronStopping()
 {
+  delete[] idregion;
   memory->destroy(elstop_ranges);
 }
 

src/GRANULAR/pair_granular.cpp

@@ -188,7 +188,7 @@ void PairGranular::compute(int eflag, int vflag)
   double *history,*allhistory,**firsthistory;
 
   bool touchflag = false;
-  const bool historyupdate = update->setupflag != 0;
+  const bool historyupdate = update->setupflag == 0;
 
   ev_init(eflag,vflag);
 

src/INTEL/pppm_disp_intel.cpp

@@ -825,8 +825,6 @@ void PPPMDispIntel::make_rho_c(IntelBuffers<flt_t,acc_t> * /*buffers*/)
   // (dx,dy,dz) = distance to "lower left" grid pt
   // (mx,my,mz) = global coords of moving stencil pt
 
-  //double *q = atom->q;
-  //double **x = atom->x;
   int nlocal = atom->nlocal;
   int nthr = comm->nthreads;
 
@@ -847,7 +845,6 @@ void PPPMDispIntel::make_rho_c(IntelBuffers<flt_t,acc_t> * /*buffers*/)
     const flt_t xi = delxinv;
     const flt_t yi = delyinv;
     const flt_t zi = delzinv;
-    const flt_t fshift = shift;
     const flt_t fshiftone = shiftone;
     const flt_t fdelvolinv = delvolinv;
 
@@ -1011,7 +1008,6 @@ void PPPMDispIntel::make_rho_g(IntelBuffers<flt_t,acc_t> * /*buffers*/)
     const flt_t xi = delxinv_6;
     const flt_t yi = delyinv_6;
     const flt_t zi = delzinv_6;
-    const flt_t fshift = shift_6;
     const flt_t fshiftone = shiftone_6;
     const flt_t fdelvolinv = delvolinv_6;
 
@@ -1150,20 +1146,13 @@ void PPPMDispIntel::make_rho_a(IntelBuffers<flt_t,acc_t> * /*buffers*/)
 {
   // clear 3d density array
 
-  memset(&(density_brick_a0[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
+  memset(&(density_brick_a0[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-         ngrid_6*sizeof(FFT_SCALAR));
+  memset(&(density_brick_a1[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a1[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
+  memset(&(density_brick_a2[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-         ngrid_6*sizeof(FFT_SCALAR));
+  memset(&(density_brick_a3[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a2[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
+  memset(&(density_brick_a4[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-         ngrid_6*sizeof(FFT_SCALAR));
+  memset(&(density_brick_a5[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a3[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
+  memset(&(density_brick_a6[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,ngrid_6*sizeof(FFT_SCALAR));
-         ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a4[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
-         ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a5[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
-         ngrid_6*sizeof(FFT_SCALAR));
-  memset(&(density_brick_a6[nzlo_out_6][nylo_out_6][nxlo_out_6]),0,
-         ngrid_6*sizeof(FFT_SCALAR));
 
   // loop over my charges, add their contribution to nearby grid points
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@@ -1171,117 +1160,112 @@ void PPPMDispIntel::make_rho_a(IntelBuffers<flt_t,acc_t> * /*buffers*/)
   // (mx,my,mz) = global coords of moving stencil pt
 
   int nlocal = atom->nlocal;
+  double **x = atom->x;
 
-    double **x = atom->x;
+  const flt_t lo0 = boxlo[0];
+  const flt_t lo1 = boxlo[1];
+  const flt_t lo2 = boxlo[2];
+  const flt_t xi = delxinv_6;
+  const flt_t yi = delyinv_6;
+  const flt_t zi = delzinv_6;
+  const flt_t fshiftone = shiftone_6;
+  const flt_t fdelvolinv = delvolinv_6;
 
-    const int nix = nxhi_out_6 - nxlo_out_6 + 1;
+  for (int i = 0; i < nlocal; i++) {
-    const int niy = nyhi_out_6 - nylo_out_6 + 1;
 
-    const flt_t lo0 = boxlo[0];
+    int nx = part2grid_6[i][0];
-    const flt_t lo1 = boxlo[1];
+    int ny = part2grid_6[i][1];
-    const flt_t lo2 = boxlo[2];
+    int nz = part2grid_6[i][2];
-    const flt_t xi = delxinv_6;
-    const flt_t yi = delyinv_6;
-    const flt_t zi = delzinv_6;
-    const flt_t fshift = shift_6;
-    const flt_t fshiftone = shiftone_6;
-    const flt_t fdelvolinv = delvolinv_6;
 
-    for (int i = 0; i < nlocal; i++) {
+    int nxsum = nx + nlower_6;
+    int nysum = ny + nlower_6;
+    int nzsum = nz + nlower_6;
 
-      int nx = part2grid_6[i][0];
+    FFT_SCALAR dx = nx+fshiftone - (x[i][0]-lo0)*xi;
-      int ny = part2grid_6[i][1];
+    FFT_SCALAR dy = ny+fshiftone - (x[i][1]-lo1)*yi;
-      int nz = part2grid_6[i][2];
+    FFT_SCALAR dz = nz+fshiftone - (x[i][2]-lo2)*zi;
 
-      int nxsum = nx + nlower_6;
+    _alignvar(flt_t rho[3][INTEL_P3M_ALIGNED_MAXORDER], 64) = {0};
-      int nysum = ny + nlower_6;
-      int nzsum = nz + nlower_6;
 
-      FFT_SCALAR dx = nx+fshiftone - (x[i][0]-lo0)*xi;
+    if (use_table) {
-      FFT_SCALAR dy = ny+fshiftone - (x[i][1]-lo1)*yi;
+      dx = dx*half_rho_scale + half_rho_scale_plus;
-      FFT_SCALAR dz = nz+fshiftone - (x[i][2]-lo2)*zi;
+      int idx = dx;
-
+      dy = dy*half_rho_scale + half_rho_scale_plus;
-      _alignvar(flt_t rho[3][INTEL_P3M_ALIGNED_MAXORDER], 64) = {0};
+      int idy = dy;
-
+      dz = dz*half_rho_scale + half_rho_scale_plus;
-      if (use_table) {
+      int idz = dz;
-        dx = dx*half_rho_scale + half_rho_scale_plus;
+#if defined(LMP_SIMD_COMPILER)
-        int idx = dx;
-        dy = dy*half_rho_scale + half_rho_scale_plus;
-        int idy = dy;
-        dz = dz*half_rho_scale + half_rho_scale_plus;
-        int idz = dz;
-        #if defined(LMP_SIMD_COMPILER)
 #if defined(USE_OMP_SIMD)
-        #pragma omp simd
+#pragma omp simd
 #else
-        #pragma simd
+#pragma simd
 #endif
-        #endif
-        for (int k = 0; k < INTEL_P3M_ALIGNED_MAXORDER; k++) {
-          rho[0][k] = rho6_lookup[idx][k];
-          rho[1][k] = rho6_lookup[idy][k];
-          rho[2][k] = rho6_lookup[idz][k];
-        }
-      } else {
-        #if defined(LMP_SIMD_COMPILER)
-#if defined(USE_OMP_SIMD)
-        #pragma omp simd
-#else
-        #pragma simd
 #endif
-        #endif
+      for (int k = 0; k < INTEL_P3M_ALIGNED_MAXORDER; k++) {
-        for (int k = nlower_6; k <= nupper_6; k++) {
+        rho[0][k] = rho6_lookup[idx][k];
-          FFT_SCALAR r1,r2,r3;
+        rho[1][k] = rho6_lookup[idy][k];
-          r1 = r2 = r3 = ZEROF;
+        rho[2][k] = rho6_lookup[idz][k];
-
-          for (int l = order_6-1; l >= 0; l--) {
-            r1 = rho_coeff_6[l][k] + r1*dx;
-            r2 = rho_coeff_6[l][k] + r2*dy;
-            r3 = rho_coeff_6[l][k] + r3*dz;
-          }
-          rho[0][k-nlower_6] = r1;
-          rho[1][k-nlower_6] = r2;
-          rho[2][k-nlower_6] = r3;
-        }
       }
-
+    } else {
-      const int type = atom->type[i];
+#if defined(LMP_SIMD_COMPILER)
-      FFT_SCALAR z0 = fdelvolinv;
-
-      #if defined(LMP_SIMD_COMPILER)
-      #pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
-      #endif
-      for (int n = 0; n < order_6; n++) {
-        int mz = n + nzsum;
-        FFT_SCALAR y0 = z0*rho[2][n];
-        #if defined(LMP_SIMD_COMPILER)
-        #pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
-        #endif
-        for (int m = 0; m < order_6; m++) {
-          int my = m + nysum;
-          FFT_SCALAR x0 = y0*rho[1][m];
-          #if defined(LMP_SIMD_COMPILER)
 #if defined(USE_OMP_SIMD)
-          #pragma omp simd
+#pragma omp simd
 #else
-          #pragma simd
+#pragma simd
 #endif
-          #pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
+#endif
-          #endif
+      for (int k = nlower_6; k <= nupper_6; k++) {
-          for (int l = 0; l < order; l++) {
+        FFT_SCALAR r1,r2,r3;
-            int mx = l + nxsum;
+        r1 = r2 = r3 = ZEROF;
-            FFT_SCALAR w = x0*rho[0][l];
+
-            density_brick_a0[mz][my][mx] += w*B[7*type];
+        for (int l = order_6-1; l >= 0; l--) {
-            density_brick_a1[mz][my][mx] += w*B[7*type+1];
+          r1 = rho_coeff_6[l][k] + r1*dx;
-            density_brick_a2[mz][my][mx] += w*B[7*type+2];
+          r2 = rho_coeff_6[l][k] + r2*dy;
-            density_brick_a3[mz][my][mx] += w*B[7*type+3];
+          r3 = rho_coeff_6[l][k] + r3*dz;
-            density_brick_a4[mz][my][mx] += w*B[7*type+4];
+        }
-            density_brick_a5[mz][my][mx] += w*B[7*type+5];
+        rho[0][k-nlower_6] = r1;
-            density_brick_a6[mz][my][mx] += w*B[7*type+6];
+        rho[1][k-nlower_6] = r2;
-          }
+        rho[2][k-nlower_6] = r3;
+      }
+    }
+
+    const int type = atom->type[i];
+    FFT_SCALAR z0 = fdelvolinv;
+
+#if defined(LMP_SIMD_COMPILER)
+#pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
+#endif
+    for (int n = 0; n < order_6; n++) {
+      int mz = n + nzsum;
+      FFT_SCALAR y0 = z0*rho[2][n];
+#if defined(LMP_SIMD_COMPILER)
+#pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
+#endif
+      for (int m = 0; m < order_6; m++) {
+        int my = m + nysum;
+        FFT_SCALAR x0 = y0*rho[1][m];
+#if defined(LMP_SIMD_COMPILER)
+#if defined(USE_OMP_SIMD)
+#pragma omp simd
+#else
+#pragma simd
+#endif
+#pragma loop_count min(2), max(INTEL_P3M_ALIGNED_MAXORDER), avg(7)
+#endif
+        for (int l = 0; l < order; l++) {
+          int mx = l + nxsum;
+          FFT_SCALAR w = x0*rho[0][l];
+          density_brick_a0[mz][my][mx] += w*B[7*type];
+          density_brick_a1[mz][my][mx] += w*B[7*type+1];
+          density_brick_a2[mz][my][mx] += w*B[7*type+2];
+          density_brick_a3[mz][my][mx] += w*B[7*type+3];
+          density_brick_a4[mz][my][mx] += w*B[7*type+4];
+          density_brick_a5[mz][my][mx] += w*B[7*type+5];
+          density_brick_a6[mz][my][mx] += w*B[7*type+6];
         }
       }
     }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1310,7 +1294,6 @@ void PPPMDispIntel::make_rho_none(IntelBuffers<flt_t,acc_t> * /*buffers*/)
     shared(nthr, nlocal, global_density) if (!_use_lrt)
   #endif
   {
-    int type;
     double **x = atom->x;
 
     const int nix = nxhi_out_6 - nxlo_out_6 + 1;
@@ -1322,7 +1305,6 @@ void PPPMDispIntel::make_rho_none(IntelBuffers<flt_t,acc_t> * /*buffers*/)
     const flt_t xi = delxinv_6;
     const flt_t yi = delyinv_6;
     const flt_t zi = delzinv_6;
-    const flt_t fshift = shift_6;
     const flt_t fshiftone = shiftone_6;
     const flt_t fdelvolinv = delvolinv_6;
 
@@ -1391,7 +1373,6 @@ void PPPMDispIntel::make_rho_none(IntelBuffers<flt_t,acc_t> * /*buffers*/)
         }
       }
 
-      type = atom->type[i];
       FFT_SCALAR z0 = fdelvolinv;
 
       #if defined(LMP_SIMD_COMPILER)
@@ -1416,7 +1397,6 @@ void PPPMDispIntel::make_rho_none(IntelBuffers<flt_t,acc_t> * /*buffers*/)
           #endif
           for (int l = 0; l < order; l++) {
             int mzyx = l + mzy;
-            FFT_SCALAR w0 = x0*rho[0][l];
             for (int k = 0; k < nsplit; k++)
               my_density[mzyx + k*ngrid_6] += x0*rho[0][l];
           }

src/KOKKOS/verlet_kokkos.h

@@ -35,6 +35,7 @@ class VerletKokkos : public Verlet {
   void setup(int) override;
   void setup_minimal(int) override;
   void run(int) override;
+  void force_clear() override;
 
   KOKKOS_INLINE_FUNCTION
   void operator() (const int& i) const {
@@ -43,13 +44,9 @@ class VerletKokkos : public Verlet {
     f(i,2) += f_merge_copy(i,2);
   }
 
-
  protected:
   DAT::t_f_array f_merge_copy,f;
-
-  void force_clear() override;
 };
-
 }
 
 #endif

src/OPENMP/pair_lubricate_poly_omp.cpp

@@ -38,8 +38,8 @@ enum{EDGE,CONSTANT,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
 
-PairLubricatePolyOMP::PairLubricatePolyOMP(LAMMPS *lmp) :
+PairLubricatePolyOMP::PairLubricatePolyOMP(LAMMPS *_lmp) :
-  PairLubricatePoly(lmp), ThrOMP(lmp, THR_PAIR)
+  PairLubricatePoly(_lmp), ThrOMP(_lmp, THR_PAIR)
 {
   suffix_flag |= Suffix::OMP;
   respa_enable = 0;

src/OPENMP/pair_lubricate_poly_omp.h

@@ -17,7 +17,7 @@
 
 #ifdef PAIR_CLASS
 // clang-format off
-PairStyle(lubricate/poly/omp,PairLubricateOMP);
+PairStyle(lubricate/poly/omp,PairLubricatePolyOMP);
 // clang-format on
 #else
 

src/OPENMP/respa_omp.cpp

@@ -117,7 +117,7 @@ void RespaOMP::setup(int flag)
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
-    force_clear(newton[ilevel]);
+    force_clear();
     modify->setup_pre_force_respa(vflag,ilevel);
 
     if (nhybrid_styles > 0) {
@@ -211,7 +211,7 @@ void RespaOMP::setup_minimal(int flag)
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
-    force_clear(newton[ilevel]);
+    force_clear();
     modify->setup_pre_force_respa(vflag,ilevel);
 
     if (nhybrid_styles > 0) {
@@ -343,7 +343,7 @@ void RespaOMP::recurse(int ilevel)
     // so that any order dependencies are the same
     // when potentials are invoked at same level
 
-    force_clear(newton[ilevel]);
+    force_clear();
     if (modify->n_pre_force_respa) {
       timer->stamp();
       modify->pre_force_respa(vflag,ilevel,iloop);

src/OPT/pair_ilp_graphene_hbn_opt.cpp

@@ -64,7 +64,7 @@ static bool check_vdw(tagint itag, tagint jtag, double *xi, double *xj);
 
 PairILPGrapheneHBNOpt::PairILPGrapheneHBNOpt(LAMMPS *lmp) :
     PairILPGrapheneHBN(lmp), layered_neigh(nullptr), first_layered_neigh(nullptr),
-    num_intra(nullptr), num_inter(nullptr), num_vdw(nullptr), special_type(nullptr)
+    special_type(nullptr), num_intra(nullptr), num_inter(nullptr), num_vdw(nullptr)
 {
   if (lmp->citeme) lmp->citeme->add(cite_ilp_cur);
 

src/SMTBQ/pair_smatb.cpp

@@ -404,7 +404,7 @@ double PairSMATB::init_one(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATB::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+int PairSMATB::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
 {
   int i, j, m;
 
@@ -465,11 +465,10 @@ void PairSMATB::write_restart_settings(FILE *fp)
 void PairSMATB::read_restart_settings(FILE *fp)
 {
   int me = comm->me;
-  size_t result;
   if (me == 0) {
-    result = fread(&offset_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
-    result = fread(&mix_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
-    result = fread(&tail_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error);
   }
   MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
   MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
@@ -504,15 +503,13 @@ void PairSMATB::write_restart(FILE *fp)
 void PairSMATB::read_restart(FILE *fp)
 {
   read_restart_settings(fp);
-
   allocate();
-  size_t result;
 
   int i, j;
   int me = comm->me;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
+      if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
       MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
       if (setflag[i][j]) {
         if (me == 0) {

src/SMTBQ/pair_smatb_single.cpp

@@ -87,7 +87,6 @@ void PairSMATBSingle::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  int *type = atom->type;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
@@ -237,16 +236,14 @@ void PairSMATBSingle::settings(int narg, char **)
 
 void PairSMATBSingle::allocate()
 {
-  int n = atom->ntypes;
+  int np1 = atom->ntypes + 1;
-  int natoms = atom->natoms;
 
-  memory->create(setflag, n + 1, n + 1, "pair_smatb:setflag");
+  memory->create(setflag, np1, np1, "pair_smatb:setflag");
-  for (int i = 1; i <= n; i++) {
+  for (int i = 1; i < np1; i++) {
-    for (int j = i; j <= n; j++) { setflag[i][j] = 0; }
+    for (int j = i; j < np1; j++) { setflag[i][j] = 0; }
   }
 
-  memory->create(cutsq, n + 1, n + 1, "pair_smatb:cutsq");
+  memory->create(cutsq, np1, np1, "pair_smatb:cutsq");
-
   allocated = 1;
 }
 
@@ -299,58 +296,40 @@ void PairSMATBSingle::init_style()
 double PairSMATBSingle::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
+
   //calculating the polynomial linking to zero
   double es = cutOffEnd - cutOffStart;
   double es2 = es * es;
   double es3 = es2 * es;
 
   //variables for poly for p and A
-  double expp = A * exp(p * (1. - cutOffStart / r0));
+  double expp = A * exp(p * (1.0 - cutOffStart / r0));
-  double ap = -1. / es3;
+  double ap = -1.0 / es3;
   double bp = p / (r0 * es2);
   double cp = -(p * p) / (es * r0 * r0);
 
-  a5 = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
+  a5 = expp * (12.0 * ap + 6.0 * bp + cp) / (2.0 * es2);
-  a4 = expp * (15. * ap + 7. * bp + cp) / es;
+  a4 = expp * (15.0 * ap + 7.0 * bp + cp) / es;
-  a3 = expp * (20. * ap + 8. * bp + cp) / 2.;
+  a3 = expp * (20.0 * ap + 8.0 * bp + cp) / 2.0;
 
   //variables for poly for q and qsi
-  double expq = QSI * exp(q * (1. - cutOffStart / r0));
+  double expq = QSI * exp(q * (1.0 - cutOffStart / r0));
   double aq = -1 / es3;
   double bq = q / (es2 * r0);
   double cq = -(q * q) / (es * r0 * r0);
 
-  x5 = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
+  x5 = expq * (12.0 * aq + 6.0 * bq + cq) / (2.0 * es2);
-  x4 = expq * (15. * aq + 7. * bq + cq) / es;
+  x4 = expq * (15.0 * aq + 7.0 * bq + cq) / es;
-  x3 = expq * (20. * aq + 8. * bq + cq) / 2.;
+  x3 = expq * (20.0 * aq + 8.0 * bq + cq) / 2.0;
 
   cutOffEnd2 = cutOffEnd * cutOffEnd;
-  if (i != j) {
-    setflag[j][i] = 1;
-    cutOffEnd2 = cutOffEnd2;
-
-    r0 = r0;
-    p = p;
-    q = q;
-    A = A;
-    QSI = QSI;
-    cutOffStart = cutOffStart;
-    cutOffEnd = cutOffEnd;
-
-    a3 = a3;
-    a4 = a4;
-    a5 = a5;
-    x3 = x3;
-    x4 = x4;
-    x5 = x5;
-  }
-
   return cutOffEnd;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
+                                       int * /*pbc*/)
 {
   int i, j, m;
 
@@ -412,11 +391,10 @@ void PairSMATBSingle::write_restart_settings(FILE *fp)
 void PairSMATBSingle::read_restart_settings(FILE *fp)
 {
   int me = comm->me;
-  size_t result;
   if (me == 0) {
-    result = fread(&offset_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
-    result = fread(&mix_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
-    result = fread(&tail_flag, sizeof(int), 1, fp);
+    utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error);
   }
   MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
   MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
@@ -451,15 +429,13 @@ void PairSMATBSingle::write_restart(FILE *fp)
 void PairSMATBSingle::read_restart(FILE *fp)
 {
   read_restart_settings(fp);
-
   allocate();
-  size_t result;
 
   int i, j;
   int me = comm->me;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
+      if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
       MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
       if (setflag[i][j]) {
         if (me == 0) {

src/TALLY/compute_force_tally.cpp

@@ -39,7 +39,6 @@ ComputeForceTally::ComputeForceTally(LAMMPS *lmp, int narg, char **arg) : Comput
   vector_flag = 0;
   peratom_flag = 1;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   comm_reverse = size_peratom_cols = 3;
   extscalar = 1;

src/TALLY/compute_heat_flux_tally.cpp

@@ -37,7 +37,6 @@ ComputeHeatFluxTally::ComputeHeatFluxTally(LAMMPS *lmp, int narg, char **arg) :
 
   vector_flag = 1;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   comm_reverse = 7;
   extvector = 1;

src/TALLY/compute_heat_flux_virial_tally.cpp

@@ -40,7 +40,6 @@ ComputeHeatFluxVirialTally::ComputeHeatFluxVirialTally(LAMMPS *lmp, int narg, ch
   vector_flag = 0;
   peratom_flag = 1;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   comm_reverse = size_peratom_cols = 3;
   extscalar = 1;

src/TALLY/compute_pe_mol_tally.cpp

@@ -36,7 +36,6 @@ ComputePEMolTally::ComputePEMolTally(LAMMPS *lmp, int narg, char **arg) : Comput
   vector_flag = 1;
   size_vector = 4;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   extvector = 1;
   peflag = 1;    // we need Pair::ev_tally() to be run

src/TALLY/compute_pe_tally.cpp

@@ -38,7 +38,6 @@ ComputePETally::ComputePETally(LAMMPS *lmp, int narg, char **arg) : Compute(lmp,
   vector_flag = 0;
   peratom_flag = 1;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   comm_reverse = size_peratom_cols = 2;
   extscalar = 1;

src/TALLY/compute_stress_tally.cpp

@@ -39,7 +39,6 @@ ComputeStressTally::ComputeStressTally(LAMMPS *lmp, int narg, char **arg) : Comp
   vector_flag = 0;
   peratom_flag = 1;
   timeflag = 1;
-  dynamic_group_allow = 0;
 
   comm_reverse = size_peratom_cols = 6;
   extscalar = 0;

src/integrate.h

@@ -26,6 +26,7 @@ class Integrate : protected Pointers {
   virtual void setup(int flag) = 0;
   virtual void setup_minimal(int) = 0;
   virtual void run(int) = 0;
+  virtual void force_clear() = 0;
   virtual void cleanup() {}
   virtual void reset_dt() {}
   virtual double memory_usage() { return 0; }

src/min.h

@@ -34,6 +34,7 @@ class Min : protected Pointers {
   virtual void setup(int flag = 1);
   virtual void setup_minimal(int);
   virtual void run(int);
+  virtual void force_clear();
   void cleanup();
   int request(class Pair *, int, double);
   virtual double memory_usage() { return 0; }
@@ -138,7 +139,6 @@ class Min : protected Pointers {
   int neigh_every, neigh_delay, neigh_dist_check;    // neighboring params
 
   virtual double energy_force(int);
-  virtual void force_clear();
 
   void ev_setup();
   void ev_set(bigint);

src/region_intersect.cpp

@@ -117,7 +117,7 @@ void RegIntersect::init()
 {
   Region::init();
 
-  // re-build list of sub-regions in case other regions were deleted
+  // re-build list of sub-regions in case regions were changed
   // error if a sub-region was deleted
 
   for (int ilist = 0; ilist < nregion; ilist++) {

src/respa.cpp

@@ -413,7 +413,7 @@ void Respa::setup(int flag)
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
-    force_clear(newton[ilevel]);
+    force_clear();
     modify->setup_pre_force_respa(vflag, ilevel);
 
     if (nhybrid_styles > 0) {
@@ -481,7 +481,7 @@ void Respa::setup_minimal(int flag)
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
-    force_clear(newton[ilevel]);
+    force_clear();
     modify->setup_pre_force_respa(vflag, ilevel);
 
     if (nhybrid_styles > 0) {
@@ -644,7 +644,7 @@ void Respa::recurse(int ilevel)
     // so that any order dependencies are the same
     // when potentials are invoked at same level
 
-    force_clear(newton[ilevel]);
+    force_clear();
     if (modify->n_pre_force_respa) {
       timer->stamp();
       modify->pre_force_respa(vflag, ilevel, iloop);
@@ -717,7 +717,7 @@ void Respa::recurse(int ilevel)
    clear other arrays as needed
 ------------------------------------------------------------------------- */
 
-void Respa::force_clear(int /*newtonflag*/)
+void Respa::force_clear()
 {
   if (external_force_clear) return;
 

src/respa.h

@@ -51,6 +51,7 @@ class Respa : public Integrate {
   void setup(int) override;
   void setup_minimal(int) override;
   void run(int) override;
+  void force_clear() override;
   void cleanup() override;
   void reset_dt() override;
 
@@ -65,7 +66,6 @@ class Respa : public Integrate {
   class FixRespa *fix_respa;    // Fix to store the force level array
 
   virtual void recurse(int);
-  void force_clear(int);
   void sum_flevel_f();
   void set_compute_flags(int ilevel);
 };

src/verlet.h

@@ -27,18 +27,16 @@ namespace LAMMPS_NS {
 class Verlet : public Integrate {
  public:
   Verlet(class LAMMPS *, int, char **);
-
   void init() override;
   void setup(int flag) override;
   void setup_minimal(int) override;
   void run(int) override;
+  void force_clear() override;
   void cleanup() override;
 
  protected:
   int triclinic;    // 0 if domain is orthog, 1 if triclinic
   int torqueflag, extraflag;
-
-  virtual void force_clear();
 };
 
 }    // namespace LAMMPS_NS

unittest/commands/CMakeLists.txt

@@ -21,6 +21,10 @@ add_executable(test_groups test_groups.cpp)
 target_link_libraries(test_groups PRIVATE lammps GTest::GMock)
 add_test(NAME Groups COMMAND test_groups)
 
+add_executable(test_regions test_regions.cpp)
+target_link_libraries(test_regions PRIVATE lammps GTest::GMock)
+add_test(NAME Regions COMMAND test_regions)
+
 add_executable(test_delete_atoms test_delete_atoms.cpp)
 target_link_libraries(test_delete_atoms PRIVATE lammps GTest::GMock)
 add_test(NAME DeleteAtoms COMMAND test_delete_atoms)

unittest/commands/test_regions.cpp

@@ -0,0 +1,270 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "lammps.h"
+
+#include "atom.h"
+#include "domain.h"
+#include "group.h"
+#include "info.h"
+#include "input.h"
+#include "region.h"
+
+#include "../testing/core.h"
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+#include <cstring>
+#include <vector>
+
+// whether to print verbose output (i.e. not capturing LAMMPS screen output).
+bool verbose = false;
+
+using LAMMPS_NS::utils::split_words;
+
+namespace LAMMPS_NS {
+using ::testing::ExitedWithCode;
+using ::testing::StrEq;
+
+class RegionTest : public LAMMPSTest {
+protected:
+    Atom *atom;
+    Domain *domain;
+    Group *group;
+
+    void SetUp() override
+    {
+        testbinary = "RegionTest";
+        LAMMPSTest::SetUp();
+        atom   = lmp->atom;
+        domain = lmp->domain;
+        group  = lmp->group;
+    }
+
+    void atomic_system()
+    {
+        BEGIN_HIDE_OUTPUT();
+        command("units real");
+        command("lattice sc 1.0 origin 0.125 0.125 0.125");
+        command("region box block -10 10 -10 10 -10 10");
+        command("create_box 8 box");
+        command("create_atoms 1 box");
+        command("mass * 1.0");
+        END_HIDE_OUTPUT();
+    }
+};
+
+TEST_F(RegionTest, NoBox)
+{
+    auto list = domain->get_region_list();
+    ASSERT_EQ(list.size(), 0);
+    BEGIN_HIDE_OUTPUT();
+    command("region reg1 block 0 1 0 1 0 1");
+    command("region reg2 cone x 0 0 2 1 0 1 units box");
+    command("region reg3 plane 0 0 0 0 0 1 side out");
+    command("region reg4 prism 0 1 0 1 0 1 0.1 0.2 0.3");
+    command("region reg5 sphere 0 0 0 1");
+    command("region reg6 union 3 reg1 reg2 reg3");
+    command("region reg7 intersect 3 reg1 reg2 reg4");
+    command("region reg8 ellipsoid 0 0 0 2 1 2");
+    END_HIDE_OUTPUT();
+    list = domain->get_region_list();
+    EXPECT_EQ(list.size(), 8);
+
+    auto reg = domain->get_region_by_id("reg1");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg2");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 0);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg3");
+    EXPECT_EQ(reg->interior, 0);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 0);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg4");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg5");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg6");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 0);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg7");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    reg = domain->get_region_by_id("reg8");
+    EXPECT_EQ(reg->interior, 1);
+    EXPECT_EQ(reg->scaleflag, 1);
+    EXPECT_EQ(reg->bboxflag, 1);
+    EXPECT_EQ(reg->varshape, 0);
+    EXPECT_EQ(reg->dynamic, 0);
+    EXPECT_EQ(reg->moveflag, 0);
+    EXPECT_EQ(reg->rotateflag, 0);
+    EXPECT_EQ(reg->openflag, 0);
+
+    BEGIN_HIDE_OUTPUT();
+    command("region reg3 delete");
+    command("region reg5 delete");
+    command("region reg6 delete");
+    command("region reg1 delete");
+    END_HIDE_OUTPUT();
+    list = domain->get_region_list();
+    EXPECT_EQ(list.size(), 4);
+
+    reg = domain->get_region_by_id("reg7");
+    TEST_FAILURE(".*ERROR: Region intersect region reg1 does not exist.*", reg->init(););
+
+    TEST_FAILURE(".*ERROR: Unrecognized region style 'xxx'.*", command("region new1 xxx"););
+    TEST_FAILURE(".*ERROR: Illegal region command.*", command("region new1 block 0 1"););
+    TEST_FAILURE(".*ERROR: Delete region reg3 does not exist.*", command("region reg3 delete"););
+}
+
+TEST_F(RegionTest, Counts)
+{
+    atomic_system();
+
+    BEGIN_HIDE_OUTPUT();
+    command("region reg1 block 0 5 0 5 -5 5");
+    command("region reg2 cone y 0 0 10 0 -5 5 units box");
+    command("region reg3 plane 0 0 0 0 1 0 side out");
+    command("region reg4 prism -5 5 -5 5 -5 5 0.1 0.2 0.3 units box");
+    command("region reg5 sphere 1 0 0 6");
+    command("region reg6 union 3 reg1 reg2 reg3");
+    command("region reg7 intersect 3 reg1 reg2 reg4");
+    command("region reg8 ellipsoid 0 0 0 5 2 3");
+    command("region reg9 ellipsoid 1 0 0 6 6 6");           // same as sphere
+    command("region reg10 prism 0 5 0 5 -5 5 0.0 0.0 0.0"); // same as block
+    END_HIDE_OUTPUT();
+
+    auto x    = atom->x;
+    auto reg1 = domain->get_region_by_id("reg1");
+    auto reg2 = domain->get_region_by_id("reg2");
+    auto reg3 = domain->get_region_by_id("reg3");
+    auto reg4 = domain->get_region_by_id("reg4");
+    auto reg5 = domain->get_region_by_id("reg5");
+    auto reg6 = domain->get_region_by_id("reg6");
+    auto reg7 = domain->get_region_by_id("reg7");
+    auto reg8 = domain->get_region_by_id("reg8");
+    auto reg9 = domain->get_region_by_id("reg9");
+    auto reg10 = domain->get_region_by_id("reg10");
+    int count1, count2, count3, count4, count5, count6, count7, count8, count9, count10;
+    count1 = count2 = count3 = count4 = count5 = count6 = count7 = count8 = count9 = count10 = 0;
+    reg1->prematch();
+    reg2->prematch();
+    reg3->prematch();
+    reg4->prematch();
+    reg5->prematch();
+    reg6->prematch();
+    reg7->prematch();
+    reg8->prematch();
+    reg9->prematch();
+    reg10->prematch();
+
+    for (int i = 0; i < atom->nlocal; ++i) {
+        if (reg1->match(x[i][0], x[i][1], x[i][2])) ++count1;
+        if (reg2->match(x[i][0], x[i][1], x[i][2])) ++count2;
+        if (reg3->match(x[i][0], x[i][1], x[i][2])) ++count3;
+        if (reg4->match(x[i][0], x[i][1], x[i][2])) ++count4;
+        if (reg5->match(x[i][0], x[i][1], x[i][2])) ++count5;
+        if (reg6->match(x[i][0], x[i][1], x[i][2])) ++count6;
+        if (reg7->match(x[i][0], x[i][1], x[i][2])) ++count7;
+        if (reg8->match(x[i][0], x[i][1], x[i][2])) ++count8;
+        if (reg9->match(x[i][0], x[i][1], x[i][2])) ++count9;
+        if (reg10->match(x[i][0], x[i][1], x[i][2])) ++count10;
+    }
+    EXPECT_EQ(count1, 250);
+    EXPECT_EQ(count2, 1132);
+    EXPECT_EQ(count3, 4000);
+    EXPECT_EQ(count4, 1000);
+    EXPECT_EQ(count5, 907);
+    EXPECT_EQ(count6, 4314);
+    EXPECT_EQ(count7, 86);
+    EXPECT_EQ(count8, 122);
+    EXPECT_EQ(count9, count5);
+    EXPECT_EQ(count10, count1);
+}
+} // namespace LAMMPS_NS
+
+int main(int argc, char **argv)
+{
+    MPI_Init(&argc, &argv);
+    ::testing::InitGoogleMock(&argc, argv);
+
+    if (platform::mpi_vendor() == "Open MPI" && !LAMMPS_NS::Info::has_exceptions())
+        std::cout << "Warning: using OpenMPI without exceptions. Death tests will be skipped\n";
+
+    // handle arguments passed via environment variable
+    if (const char *var = getenv("TEST_ARGS")) {
+        std::vector<std::string> env = split_words(var);
+        for (auto arg : env) {
+            if (arg == "-v") {
+                verbose = true;
+            }
+        }
+    }
+
+    if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;
+
+    int rv = RUN_ALL_TESTS();
+    MPI_Finalize();
+    return rv;
+}