diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
index b10a9256b82a446830bb6af0644a4820df245bee..d0ff5d290dfd0cfcbfc08d2f8dcd6a0bb165c75b 100644
--- a/src/CORESHELL/compute_temp_cs.cpp
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -309,7 +309,7 @@ void ComputeTempCS::vcm_pairs()
 
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxatom) {
+  if (atom->nmax > maxatom) {
     memory->destroy(vint);
     maxatom = atom->nmax;
     memory->create(vint,maxatom,3,"temp/cs:vint");
diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp
index 0e2e6038c1a873dac7e9cb0ba648002211188597..4d343367f1029fef74b6023004d2e4067a1f4b24 100644
--- a/src/KOKKOS/atom_kokkos.cpp
+++ b/src/KOKKOS/atom_kokkos.cpp
@@ -122,7 +122,7 @@ void AtomKokkos::sort()
 
   // reallocate per-atom vectors if needed
 
-  if (nlocal > maxnext) {
+  if (atom->nmax > maxnext) {
     memory->destroy(next);
     memory->destroy(permute);
     maxnext = atom->nmax;
@@ -265,4 +265,4 @@ AtomVec *AtomKokkos::new_avec(const char *style, int trysuffix, int &sflag)
   if (!avec->kokkosable)
     error->all(FLERR,"KOKKOS package requires a kokkos enabled atom_style");
   return avec;
-}
\ No newline at end of file
+}
diff --git a/src/KOKKOS/fix_langevin_kokkos.cpp b/src/KOKKOS/fix_langevin_kokkos.cpp
index de8920b55690d5a437a92982dd270819e0548218..43af4168f3996e06dea2625e9952f9cafa481d94 100644
--- a/src/KOKKOS/fix_langevin_kokkos.cpp
+++ b/src/KOKKOS/fix_langevin_kokkos.cpp
@@ -667,7 +667,7 @@ void FixLangevinKokkos<DeviceType>::compute_target()
         error->one(FLERR,"Fix langevin variable returned negative temperature");
       tsqrt = sqrt(t_target);
     } else {
-      if (nlocal > maxatom2) {
+      if (atom->nmax > maxatom2) {
         maxatom2 = atom->nmax;
         memory->destroy_kokkos(k_tforce,tforce);
         memory->create_kokkos(k_tforce,tforce,maxatom2,"langevin:tforce");
diff --git a/src/KOKKOS/fix_setforce_kokkos.cpp b/src/KOKKOS/fix_setforce_kokkos.cpp
index 862aab1e21761233b1f2ba004af2f5037cc9d446..93e5e19377dfb798089a115a7c83436d03022ef2 100755
--- a/src/KOKKOS/fix_setforce_kokkos.cpp
+++ b/src/KOKKOS/fix_setforce_kokkos.cpp
@@ -94,7 +94,7 @@ void FixSetForceKokkos<DeviceType>::post_force(int vflag)
 
   // reallocate sforce array if necessary
 
-  if (varflag == ATOM && nlocal > maxatom) {
+  if (varflag == ATOM && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     memory->destroy_kokkos(k_sforce,sforce);
     memory->create_kokkos(k_sforce,sforce,maxatom,3,"setforce:sforce");
diff --git a/src/KOKKOS/neighbor_kokkos.cpp b/src/KOKKOS/neighbor_kokkos.cpp
index c018eebb1bc013843d136d9fe36f539d02d197d3..55fc0f1e9b55c2a6d3b34ca614c58cfa009c429e 100644
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@@ -430,7 +430,7 @@ void NeighborKokkos::build_kokkos(int topoflag)
     int nlocal = atom->nlocal;
     if (includegroup) nlocal = atom->nfirst;
     int maxhold_kokkos = xhold.view<DeviceType>().dimension_0();
-    if (nlocal > maxhold || maxhold_kokkos < maxhold) {
+    if (atom->nmax > maxhold || maxhold_kokkos < maxhold) {
       maxhold = atom->nmax;
       xhold = DAT::tdual_x_array("neigh:xhold",maxhold);
     }
@@ -464,12 +464,12 @@ void NeighborKokkos::build_kokkos(int topoflag)
   // else only invoke grow() if nlocal exceeds previous list size
   // only for lists with growflag set and which are perpetual (glist)
 
-  if (anyghostlist && atom->nlocal+atom->nghost > maxatom) {
+  if (anyghostlist && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     for (i = 0; i < nglist; i++)
       if (lists[glist[i]]) lists[glist[i]]->grow(maxatom);
       else init_list_grow_kokkos(glist[i]);
-  } else if (atom->nlocal > maxatom) {
+  } else if (atom->nmax > maxatom) {
     maxatom = atom->nmax;
     for (i = 0; i < nglist; i++)
       if (lists[glist[i]]) lists[glist[i]]->grow(maxatom);
diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index f3a4188de9ae041309a454f621026c514f2d6151..39a98901cb266d5388cffef8ea08fd93c6ed3976 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -368,7 +368,7 @@ void Ewald::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(ek);
     memory->destroy3d_offset(cs,-kmax_created);
     memory->destroy3d_offset(sn,-kmax_created);
diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp
index c1a8d81e00ee53a1bf09d39860ec8ef93e6b20a9..c567ce01682337c4f1fea0bb76a79e270d79aed4 100644
--- a/src/KSPACE/ewald_disp.cpp
+++ b/src/KSPACE/ewald_disp.cpp
@@ -379,7 +379,7 @@ void EwaldDisp::reallocate()
 void EwaldDisp::reallocate_atoms()
 {
   if (eflag_atom || vflag_atom)
-    if (atom->nlocal > nmax) {
+    if (atom->nmax > nmax) {
       deallocate_peratom();
       allocate_peratom();
       nmax = atom->nmax;
diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp
index 7341688bb2c5a6a726d8d8f48151f84f7ea14115..b9d894a00adc863516ca25bf811d4ce59af8d2f4 100644
--- a/src/KSPACE/msm.cpp
+++ b/src/KSPACE/msm.cpp
@@ -488,7 +488,7 @@ void MSM::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"msm:part2grid");
diff --git a/src/KSPACE/msm_cg.cpp b/src/KSPACE/msm_cg.cpp
index c7af52ef89384cfa92302522f468bc5974a33a16..31f70b1934189ee4b2850e469039ac054a8c793a 100644
--- a/src/KSPACE/msm_cg.cpp
+++ b/src/KSPACE/msm_cg.cpp
@@ -101,7 +101,7 @@ void MSMCG::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     memory->destroy(is_charged);
     nmax = atom->nmax;
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index a4cc6f89dfccec0873411f920746fac539d51401..9b18ad88bb64fbd3abe8e41356256db4268232b5 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -638,7 +638,7 @@ void PPPM::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"pppm:part2grid");
diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp
index 6fcdb438a8d7eb94e08cc198419f5d0d9b62b3bd..3e5e4431a603d216d7cbe12191fd4b5976d435c7 100644
--- a/src/KSPACE/pppm_cg.cpp
+++ b/src/KSPACE/pppm_cg.cpp
@@ -101,7 +101,7 @@ void PPPMCG::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     memory->destroy(is_charged);
     nmax = atom->nmax;
diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp
index 2c065e926f0a64086652367ba9af53b3503c2d51..b7ec188150e51752eb962c9b90f66a6d51339751 100755
--- a/src/KSPACE/pppm_disp.cpp
+++ b/src/KSPACE/pppm_disp.cpp
@@ -900,7 +900,7 @@ void PPPMDisp::compute(int eflag, int vflag)
   }
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
 
     if (function[0]) memory->destroy(part2grid);
     if (function[1] + function[2] + function[3]) memory->destroy(part2grid_6);
diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp
index 7cd719ec73c1b0d71704b07716bb57522b573b10..a16c0d474b80f53f5578c3e54f47738300a76649 100755
--- a/src/KSPACE/pppm_stagger.cpp
+++ b/src/KSPACE/pppm_stagger.cpp
@@ -144,7 +144,7 @@ void PPPMStagger::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"pppm:part2grid");
diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index 49728d6aa61fc782051306f5c87f037c7a766d16..d4e57c2c6cb62eda855b0cc0e9556bae83f3e3bb 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -56,6 +56,9 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
   ljflag = torflag = 1;
   morseflag = 0;
 
+  nextra = 3;
+  pvector = new double[nextra];
+
   maxlocal = 0;
   REBO_numneigh = NULL;
   REBO_firstneigh = NULL;
@@ -78,6 +81,7 @@ PairAIREBO::~PairAIREBO()
   delete [] ipage;
   memory->destroy(nC);
   memory->destroy(nH);
+  delete [] pvector;
 
   if (allocated) {
     memory->destroy(setflag);
@@ -99,6 +103,7 @@ void PairAIREBO::compute(int eflag, int vflag)
 {
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = vflag_atom = 0;
+  pvector[0] = pvector[1] = pvector[2] = 0.0;
 
   REBO_neigh();
   FREBO(eflag,vflag);
@@ -490,7 +495,7 @@ void PairAIREBO::FREBO(int eflag, int vflag)
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
 
-      if (eflag) evdwl = VR + bij*VA;
+      if (eflag) pvector[0] += evdwl = VR + bij*VA;
       if (evflag) ev_tally(i,j,nlocal,newton_pair,
                            evdwl,0.0,fpair,delx,dely,delz);
     }
@@ -789,7 +794,7 @@ void PairAIREBO::FLJ(int eflag, int vflag)
       f[j][1] -= delij[1]*fpair;
       f[j][2] -= delij[2]*fpair;
 
-      if (eflag) evdwl = VA*Stb + (1.0-Str)*cij*VLJ;
+      if (eflag) pvector[1] += evdwl = VA*Stb + (1.0-Str)*cij*VLJ;
       if (evflag) ev_tally(i,j,nlocal,newton_pair,
                            evdwl,0.0,fpair,delij[0],delij[1],delij[2]);
 
@@ -1032,7 +1037,7 @@ void PairAIREBO::TORSION(int eflag, int vflag)
           Ec = 256.0*ekijl/405.0;
           Vtors = (Ec*(powint(cw2,5)))-(ekijl/10.0);
 
-          if (eflag) evdwl = Vtors*w21*w23*w34*(1.0-tspjik)*(1.0-tspijl);
+          if (eflag) pvector[2] += evdwl = Vtors*w21*w23*w34*(1.0-tspjik)*(1.0-tspijl);
 
           dndij[0] = (cross234[1]*del21[2])-(cross234[2]*del21[1]);
           dndij[1] = (cross234[2]*del21[0])-(cross234[0]*del21[2]);
@@ -1739,9 +1744,9 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
         cos321 = MAX(cos321,-1.0);
         Sp2(cos321,thmin,thmax,dcut321);
         sin321 = sqrt(1.0 - cos321*cos321);
-        sink2i = 1.0/(sin321*sin321);
-        rik2i = 1.0/(r21mag*r21mag);
-        if (sin321 != 0.0) {
+        if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0
+          sink2i = 1.0/(sin321*sin321);
+          rik2i = 1.0/(r21mag*r21mag);
           rr = (r23mag*r23mag)-(r21mag*r21mag);
           rjk[0] = r21[0]-r23[0];
           rjk[1] = r21[1]-r23[1];
@@ -1776,10 +1781,10 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
               cos234 = MIN(cos234,1.0);
               cos234 = MAX(cos234,-1.0);
               sin234 = sqrt(1.0 - cos234*cos234);
-              sinl2i = 1.0/(sin234*sin234);
-              rjl2i = 1.0/(r34mag*r34mag);
 
-              if (sin234 != 0.0) {
+              if ((sin234 > TOL) && (r34mag > TOL)) {  // XXX was sin234 != 0.0
+                sinl2i = 1.0/(sin234*sin234);
+                rjl2i = 1.0/(r34mag*r34mag);
                 w34 = Sp(r34mag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dw34);
                 rr = (r23mag*r23mag)-(r34mag*r34mag);
                 ril[0] = r23[0]+r34[0];
@@ -2667,10 +2672,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag,
           cos321 = MIN(cos321,1.0);
           cos321 = MAX(cos321,-1.0);
           sin321 = sqrt(1.0 - cos321*cos321);
-          sink2i = 1.0/(sin321*sin321);
-          rik2i = 1.0/(r21mag*r21mag);
 
-          if (sin321 != 0.0) {
+          if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0
+            sink2i = 1.0/(sin321*sin321);
+            rik2i = 1.0/(r21mag*r21mag);
             rr = (rijmag*rijmag)-(r21mag*r21mag);
             rjk[0] = r21[0]-rij[0];
             rjk[1] = r21[1]-rij[1];
@@ -2704,10 +2709,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag,
                 cos234 = MIN(cos234,1.0);
                 cos234 = MAX(cos234,-1.0);
                 sin234 = sqrt(1.0 - cos234*cos234);
-                sinl2i = 1.0/(sin234*sin234);
-                rjl2i = 1.0/(r34mag*r34mag);
 
-                if (sin234 != 0.0) {
+                if ((sin234 > TOL) && (r34mag > TOL)) { // XXX was sin234 != 0.0
+                  sinl2i = 1.0/(sin234*sin234);
+                  rjl2i = 1.0/(r34mag*r34mag);
                   w34 = Sp(r34mag,rcmin[jtype][ltype],
                            rcmaxp[jtype][ltype],dw34);
                   rr = (r23mag*r23mag)-(r34mag*r34mag);
diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp
index 689aa0a168d3a6c43ce19861a76915fb1cc73237..41a54a64fff2232bc2ff47e33155e14eb9b3d966 100644
--- a/src/MC/fix_atom_swap.cpp
+++ b/src/MC/fix_atom_swap.cpp
@@ -584,7 +584,7 @@ void FixAtomSwap::update_semi_grand_atoms_list()
   int nlocal = atom->nlocal;
   double **x = atom->x;
 
-  if (nlocal > atom_swap_nmax) {
+  if (atom->nmax > atom_swap_nmax) {
     memory->sfree(local_swap_atom_list);
     atom_swap_nmax = atom->nmax;
     local_swap_atom_list = (int *) memory->smalloc(atom_swap_nmax*sizeof(int),
@@ -639,7 +639,7 @@ void FixAtomSwap::update_swap_atoms_list()
   int *type = atom->type;
   double **x = atom->x;
 
-  if (nlocal > atom_swap_nmax) {
+  if (atom->nmax > atom_swap_nmax) {
     memory->sfree(local_swap_iatom_list);
     memory->sfree(local_swap_jatom_list);
     atom_swap_nmax = atom->nmax;
diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp
index 0c36cecb3099002d0e020a908267897222470cfd..31194ce12ec0286b7eb955d01ffb29eaf54887c5 100644
--- a/src/MC/fix_bond_swap.cpp
+++ b/src/MC/fix_bond_swap.cpp
@@ -240,7 +240,7 @@ void FixBondSwap::post_integrate()
   // randomize list of my owned atoms that are in fix group
   // grow atom list if necessary
 
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(alist);
     nmax = atom->nmax;
     memory->create(alist,nmax,"bondswap:alist");
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index c2a69399cef550feb8a180d19cd19fd6de11c019..aea029e5ba2afb00dd0a05ef4abd00d1b75c40d5 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -2197,7 +2197,7 @@ void FixGCMC::update_gas_atoms_list()
   tagint *molecule = atom->molecule;
   double **x = atom->x;
 
-  if (nlocal > gcmc_nmax) {
+  if (atom->nmax > gcmc_nmax) {
     memory->sfree(local_gas_list);
     gcmc_nmax = atom->nmax;
     local_gas_list = (int *) memory->smalloc(gcmc_nmax*sizeof(int),
diff --git a/src/MISC/fix_efield.cpp b/src/MISC/fix_efield.cpp
index d729b9b74a73bc58550d958aa248177ba069f7b9..2151a05092e8e8432f06670dc21e8e502fe4ba3c 100644
--- a/src/MISC/fix_efield.cpp
+++ b/src/MISC/fix_efield.cpp
@@ -254,7 +254,7 @@ void FixEfield::post_force(int vflag)
 
   // reallocate efield array if necessary
 
-  if (varflag == ATOM && nlocal > maxatom) {
+  if (varflag == ATOM && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     memory->destroy(efield);
     memory->create(efield,maxatom,4,"efield:efield");
diff --git a/src/MISC/fix_evaporate.cpp b/src/MISC/fix_evaporate.cpp
index 873227449d7226979744d7d64be0a0de7297c46e..42f3b844db47caed7f8f434c0f791c16dcda0282 100644
--- a/src/MISC/fix_evaporate.cpp
+++ b/src/MISC/fix_evaporate.cpp
@@ -172,7 +172,7 @@ void FixEvaporate::pre_exchange()
 
   // grow list and mark arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(list);
     memory->destroy(mark);
     nmax = atom->nmax;
diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp
index cf8cf35cd393e706bd510734faba45394e7bb029..305ab4e7166e41794a76cab78527fd3ecc99bb98 100644
--- a/src/PERI/compute_damage_atom.cpp
+++ b/src/PERI/compute_damage_atom.cpp
@@ -77,7 +77,7 @@ void ComputeDamageAtom::compute_peratom()
 
   // grow damage array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(damage);
     nmax = atom->nmax;
     memory->create(damage,nmax,"damage/atom:damage");
diff --git a/src/PERI/compute_dilatation_atom.cpp b/src/PERI/compute_dilatation_atom.cpp
index 20993fab47b2d0f8f3cf76542e151813f519b363..d79b80f4d459ad3d60ee43d9795f898860024633 100644
--- a/src/PERI/compute_dilatation_atom.cpp
+++ b/src/PERI/compute_dilatation_atom.cpp
@@ -94,7 +94,7 @@ void ComputeDilatationAtom::compute_peratom()
 
   // grow dilatation array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(dilatation);
     nmax = atom->nmax;
     memory->create(dilatation,nmax,"dilatation/atom:dilatation");
diff --git a/src/PERI/compute_plasticity_atom.cpp b/src/PERI/compute_plasticity_atom.cpp
index 2c47f18b1fc57219d9d7ef27ebb8ff74d754b169..34fd23fe6a663bccb36257724ff09a4f46b146ab 100644
--- a/src/PERI/compute_plasticity_atom.cpp
+++ b/src/PERI/compute_plasticity_atom.cpp
@@ -82,7 +82,7 @@ void ComputePlasticityAtom::compute_peratom()
 
   // grow damage array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(plasticity);
     nmax = atom->nmax;
     memory->create(plasticity,nmax,"plasticity/atom:plasticity");
diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp
index 5985c756ad1fff8d70e7e0b8e567a2e0e8002751..80d8f58f0d4ce99bce3c0f811abae7a7cf232ff6 100644
--- a/src/REPLICA/verlet_split.cpp
+++ b/src/REPLICA/verlet_split.cpp
@@ -444,7 +444,7 @@ void VerletSplit::rk_setup()
   // grow f_kspace array on master procs if necessary
 
   if (master) {
-    if (atom->nlocal > maxatom) {
+    if (atom->nmax > maxatom) {
       memory->destroy(f_kspace);
       maxatom = atom->nmax;
       memory->create(f_kspace,maxatom,3,"verlet/split:f_kspace");
diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp
index 1f7c943719d9c4da6ce97fbf9d24b60b0573579f..da948e56cca5fa526198cbd1e64b676609e2e730 100644
--- a/src/SNAP/compute_sna_atom.cpp
+++ b/src/SNAP/compute_sna_atom.cpp
@@ -183,7 +183,7 @@ void ComputeSNAAtom::compute_peratom()
 
   // grow sna array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(sna);
     nmax = atom->nmax;
     memory->create(sna,nmax,size_peratom_cols,"sna/atom:sna");
diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp
index 8767ca174574f9ae1af892938515daa69dba15d9..2164b3102917542d6fde788d7ab518aa532cbad3 100644
--- a/src/SNAP/compute_snad_atom.cpp
+++ b/src/SNAP/compute_snad_atom.cpp
@@ -179,7 +179,7 @@ void ComputeSNADAtom::compute_peratom()
 
   // grow snad array if necessary
 
-  if (ntotal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(snad);
     nmax = atom->nmax;
     memory->create(snad,nmax,size_peratom_cols,
diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp
index cbe53821b80554ba0fb791910ac41a7f979819c4..29431f29a94d1e4c21aaa62ce42f4344396b9105 100644
--- a/src/SNAP/compute_snav_atom.cpp
+++ b/src/SNAP/compute_snav_atom.cpp
@@ -182,7 +182,7 @@ void ComputeSNAVAtom::compute_peratom()
 
   // grow snav array if necessary
 
-  if (ntotal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(snav);
     nmax = atom->nmax;
     memory->create(snav,nmax,size_peratom_cols,
diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp
index 35eb21254c7787375e39528364a678e1216ed42d..ae1114732732bee6adf3cbdf303f8da77d9debf2 100644
--- a/src/SRD/fix_srd.cpp
+++ b/src/SRD/fix_srd.cpp
@@ -471,7 +471,7 @@ void FixSRD::pre_neighbor()
 
   // grow SRD per-atom bin arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     nmax = atom->nmax;
     memory->destroy(binsrd);
     memory->destroy(binnext);
diff --git a/src/USER-DPD/compute_dpd_atom.cpp b/src/USER-DPD/compute_dpd_atom.cpp
index 9629aa54358c594a378ee978643312cf5f2b6dd4..8d944a292feb641f5fea9d6193fdcee20f4b8be4 100644
--- a/src/USER-DPD/compute_dpd_atom.cpp
+++ b/src/USER-DPD/compute_dpd_atom.cpp
@@ -78,15 +78,15 @@ void ComputeDpdAtom::compute_peratom()
   double *uCond = atom->uCond;
   double *uMech = atom->uMech;
   double *dpdTheta = atom->dpdTheta;
-  int nlocal = atom->nlocal;
   int *mask = atom->mask;
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(dpdAtom);
     nmax = atom->nmax;
     memory->create(dpdAtom,nmax,size_peratom_cols,"dpd/atom:dpdAtom");
     array_atom = dpdAtom;
   }
 
+  const int nlocal = atom->nlocal;
   for (int i = 0; i < nlocal; i++){
     if (mask[i] & groupbit){
       dpdAtom[i][0] =  uCond[i];
diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp
index 93ad934d6d5a05c34081b3b7a400edb263789d8c..8d03e97b0af4d653d3d1d4809009d60b868c664a 100644
--- a/src/USER-EFF/compute_ke_atom_eff.cpp
+++ b/src/USER-EFF/compute_ke_atom_eff.cpp
@@ -75,7 +75,7 @@ void ComputeKEAtomEff::compute_peratom()
 
   // grow ke array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(ke);
     nmax = atom->nmax;
     memory->create(ke,nmax,"compute/ke/atom/eff:ke");
diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp
index 12c72479cbe6b357cd3f5b286fa8052ecb100fe3..40e8e529c0ea2baf1a0ada81edd5b0cbe06356bb 100644
--- a/src/USER-EFF/compute_temp_deform_eff.cpp
+++ b/src/USER-EFF/compute_temp_deform_eff.cpp
@@ -264,7 +264,7 @@ void ComputeTempDeformEff::remove_bias_all()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/deform/eff:vbiasall");
diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp
index 1035a3e31bbbed7c1db89c508dd7d8739776a45a..c0cd56f3a8d65c5b55d30108446cb3c58566e702 100644
--- a/src/USER-EFF/compute_temp_region_eff.cpp
+++ b/src/USER-EFF/compute_temp_region_eff.cpp
@@ -237,7 +237,7 @@ void ComputeTempRegionEff::remove_bias_all()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/region:vbiasall");
diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp
index f2ca444de204b8b933b7645018f29386fb185f99..9584623fc3970e055ae9113d76b12e62cfccb1b9 100644
--- a/src/USER-EFF/fix_langevin_eff.cpp
+++ b/src/USER-EFF/fix_langevin_eff.cpp
@@ -104,7 +104,7 @@ void FixLangevinEff::post_force_no_tally()
         error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
       tsqrt = sqrt(t_target);
     } else {
-      if (nlocal > maxatom2) {
+      if (atom->nmax > maxatom2) {
         maxatom2 = atom->nmax;
         memory->destroy(tforce);
         memory->create(tforce,maxatom2,"langevin/eff:tforce");
@@ -241,7 +241,7 @@ void FixLangevinEff::post_force_tally()
 
   // reallocate flangevin and erforcelangevin if necessary
 
-  if (atom->nlocal > maxatom1) {
+  if (atom->nmax > maxatom1) {
     memory->destroy(flangevin);
     memory->destroy(erforcelangevin);
     maxatom1 = atom->nmax;
@@ -279,7 +279,7 @@ void FixLangevinEff::post_force_tally()
         error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
       tsqrt = sqrt(t_target);
     } else {
-      if (nlocal > maxatom2) {
+      if (atom->nmax > maxatom2) {
         maxatom2 = atom->nmax;
         memory->destroy(tforce);
         memory->create(tforce,maxatom2,"langevin/eff:tforce");
diff --git a/src/USER-INTEL/pppm_intel.cpp b/src/USER-INTEL/pppm_intel.cpp
index 466b924a51e6bd9da088c9bf80c620a447ce3c8a..bd422aac828950702388cd02cded72189073bc99 100644
--- a/src/USER-INTEL/pppm_intel.cpp
+++ b/src/USER-INTEL/pppm_intel.cpp
@@ -144,7 +144,7 @@ void PPPMIntel::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"pppm:part2grid");
diff --git a/src/USER-MISC/compute_ackland_atom.cpp b/src/USER-MISC/compute_ackland_atom.cpp
index 7a60117bb84f653130043c7a2c863dd27707c4e7..ac48d75bc6dc57b46f00d22a2583ecbf9835d1e9 100644
--- a/src/USER-MISC/compute_ackland_atom.cpp
+++ b/src/USER-MISC/compute_ackland_atom.cpp
@@ -105,7 +105,7 @@ void ComputeAcklandAtom::compute_peratom()
 
   // grow structure array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(structure);
     nmax = atom->nmax;
     memory->create(structure,nmax,"compute/ackland/atom:ackland");
diff --git a/src/USER-MISC/compute_basal_atom.cpp b/src/USER-MISC/compute_basal_atom.cpp
index d323e20a1081df13f42e3d6f067c34645aad5bc3..881b25313115c0c02ebcdf336fc0a6e4fb15f522 100644
--- a/src/USER-MISC/compute_basal_atom.cpp
+++ b/src/USER-MISC/compute_basal_atom.cpp
@@ -109,7 +109,7 @@ void ComputeBasalAtom::compute_peratom()
 
   // grow structure array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(BPV);
     nmax = atom->nmax;
     memory->create(BPV,nmax,3,"basal/atom:basal");
diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp
index 574db8757c808b767ac3de50cf3e4f1cb990c5af..2210555a7c0664a4aec31e2e4f796efb3267bbf1 100644
--- a/src/USER-MISC/compute_temp_rotate.cpp
+++ b/src/USER-MISC/compute_temp_rotate.cpp
@@ -112,7 +112,7 @@ double ComputeTempRotate::compute_scalar()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/rotate:vbiasall");
@@ -175,7 +175,7 @@ void ComputeTempRotate::compute_vector()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/rotate:vbiasall");
diff --git a/src/USER-OMP/ewald_omp.cpp b/src/USER-OMP/ewald_omp.cpp
index 7ab903aab15fcebc018792db0d74828c18f246df..ad6b0bf80573f990616b53f0195ae092fa4f7e70 100644
--- a/src/USER-OMP/ewald_omp.cpp
+++ b/src/USER-OMP/ewald_omp.cpp
@@ -68,7 +68,7 @@ void EwaldOMP::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(ek);
     memory->destroy3d_offset(cs,-kmax_created);
     memory->destroy3d_offset(sn,-kmax_created);
diff --git a/src/USER-OMP/msm_cg_omp.cpp b/src/USER-OMP/msm_cg_omp.cpp
index f743abae6fd0b7ae5b25dc58842cf5db4cd1662c..9f1ebba2a8e5b78dde0a678a076c87c004cf0596 100644
--- a/src/USER-OMP/msm_cg_omp.cpp
+++ b/src/USER-OMP/msm_cg_omp.cpp
@@ -103,7 +103,7 @@ void MSMCGOMP::compute(int eflag, int vflag)
 
   // extend size of per-atom arrays if necessary
 
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(part2grid);
     memory->destroy(is_charged);
     nmax = atom->nmax;
diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp
index eeb3e134f6ac016fb795e93b640a8c252c9f4dcb..35a644007ac9060a93061d9438baf4ba1ffe434f 100644
--- a/src/USER-OMP/pair_airebo_omp.cpp
+++ b/src/USER-OMP/pair_airebo_omp.cpp
@@ -47,6 +47,8 @@ PairAIREBOOMP::PairAIREBOOMP(LAMMPS *lmp) :
 
 void PairAIREBOOMP::compute(int eflag, int vflag)
 {
+  double pv0=0.0,pv1=0.0,pv2=0.0;
+
   if (eflag || vflag) {
     ev_setup(eflag,vflag);
   } else evflag = vflag_fdotr = vflag_atom = 0;
@@ -58,7 +60,7 @@ void PairAIREBOOMP::compute(int eflag, int vflag)
   const int inum = list->inum;
 
 #if defined(_OPENMP)
-#pragma omp parallel default(none) shared(eflag,vflag)
+#pragma omp parallel default(none) shared(eflag,vflag) reduction(+:pv0,pv1,pv2)
 #endif
   {
     int ifrom, ito, tid;
@@ -68,13 +70,17 @@ void PairAIREBOOMP::compute(int eflag, int vflag)
     thr->timer(Timer::START);
     ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
 
-    FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,thr);
-    if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,thr);
-    if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,thr);
+    FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv0,thr);
+    if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv1,thr);
+    if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,&pv2,thr);
 
     thr->timer(Timer::PAIR);
     reduce_thr(this, eflag, vflag, thr);
   } // end of omp parallel region
+
+  pvector[0] = pv0;
+  pvector[1] = pv1;
+  pvector[2] = pv2;
 }
 
 /* ----------------------------------------------------------------------
@@ -586,9 +592,9 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
         cos321 = MAX(cos321,-1.0);
         Sp2(cos321,thmin,thmax,dcut321);
         sin321 = sqrt(1.0 - cos321*cos321);
-        sink2i = 1.0/(sin321*sin321);
-        rik2i = 1.0/(r21mag*r21mag);
-        if (sin321 != 0.0) {
+        if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0
+          sink2i = 1.0/(sin321*sin321);
+          rik2i = 1.0/(r21mag*r21mag);
           rr = (r23mag*r23mag)-(r21mag*r21mag);
           rjk[0] = r21[0]-r23[0];
           rjk[1] = r21[1]-r23[1];
@@ -623,10 +629,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
               cos234 = MIN(cos234,1.0);
               cos234 = MAX(cos234,-1.0);
               sin234 = sqrt(1.0 - cos234*cos234);
-              sinl2i = 1.0/(sin234*sin234);
-              rjl2i = 1.0/(r34mag*r34mag);
 
-              if (sin234 != 0.0) {
+              if ((sin234 > TOL) && (r34mag > TOL)) {  // XXX was sin234 != 0.0
+                sinl2i = 1.0/(sin234*sin234);
+                rjl2i = 1.0/(r34mag*r34mag);
                 w34 = Sp(r34mag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dw34);
                 rr = (r23mag*r23mag)-(r34mag*r34mag);
                 ril[0] = r23[0]+r34[0];
@@ -657,7 +663,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
                 cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234;
                 om1234 = cwnum/cwnom;
                 cw = om1234;
-                Etmp += ((1.0-(om1234*om1234))*w21*w34) *
+                Etmp += ((1.0-square(om1234))*w21*w34) *
                   (1.0-tspjik)*(1.0-tspijl);
 
                 dt1dik = (rik2i)-(dctik*sink2i*cos321);
@@ -684,7 +690,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
 
                 aa = (prefactor*2.0*cw/cwnom)*w21*w34 *
                   (1.0-tspjik)*(1.0-tspijl);
-                aaa1 = -prefactor*(1.0-(om1234*om1234)) *
+                aaa1 = -prefactor*(1.0-square(om1234)) *
                   (1.0-tspjik)*(1.0-tspijl);
                 aaa2 = aaa1*w21*w34;
                 at2 = aa*cwnum;
@@ -1027,8 +1033,8 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
       // evaluate splines g and derivatives dg
 
       g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);
-      Etmp = Etmp+(wjl*g*exp(lamdaijl));
-      tmp3 = tmp3+(wjl*dgdN*exp(lamdaijl));
+      Etmp += (wjl*g*exp(lamdaijl));
+      tmp3 += (wjl*dgdN*exp(lamdaijl));
       NconjtmpJ = NconjtmpJ+(kronecker(ltype,0)*wjl*Sp(Nlj,Nmin,Nmax,dS));
     }
   }
@@ -1119,7 +1125,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
                                 (crosskij[1]*crossijl[1]) +
                                 (crosskij[2]*crossijl[2])) /
                                (crosskijmag*crossijlmag));
-                Etmp += ((1.0-(omkijl*omkijl))*wik*wjl) *
+                Etmp += ((1.0-square(omkijl))*wik*wjl) *
                   (1.0-tspjik)*(1.0-tspijl);
               }
             }
@@ -1514,10 +1520,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
           cos321 = MIN(cos321,1.0);
           cos321 = MAX(cos321,-1.0);
           sin321 = sqrt(1.0 - cos321*cos321);
-          sink2i = 1.0/(sin321*sin321);
-          rik2i = 1.0/(r21mag*r21mag);
 
-          if (sin321 != 0.0) {
+          if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0
+            sink2i = 1.0/(sin321*sin321);
+            rik2i = 1.0/(r21mag*r21mag);
             rr = (rijmag*rijmag)-(r21mag*r21mag);
             rjk[0] = r21[0]-rij[0];
             rjk[1] = r21[1]-rij[1];
@@ -1551,10 +1557,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
                 cos234 = MIN(cos234,1.0);
                 cos234 = MAX(cos234,-1.0);
                 sin234 = sqrt(1.0 - cos234*cos234);
-                sinl2i = 1.0/(sin234*sin234);
-                rjl2i = 1.0/(r34mag*r34mag);
 
-                if (sin234 != 0.0) {
+                if ((sin234 > TOL) && (r34mag > TOL)) { // XXX was sin234 != 0.0
+                  sinl2i = 1.0/(sin234*sin234);
+                  rjl2i = 1.0/(r34mag*r34mag);
                   w34 = Sp(r34mag,rcmin[jtype][ltype],
                            rcmaxp[jtype][ltype],dw34);
                   rr = (r23mag*r23mag)-(r34mag*r34mag);
@@ -1586,7 +1592,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
                   cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234;
                   om1234 = cwnum/cwnom;
                   cw = om1234;
-                  Etmp += ((1.0-(om1234*om1234))*w21*w34) *
+                  Etmp += ((1.0-square(om1234))*w21*w34) *
                     (1.0-tspjik)*(1.0-tspijl);
 
                   dt1dik = (rik2i)-(dctik*sink2i*cos321);
@@ -1616,7 +1622,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
 
                   aa = (prefactor*2.0*cw/cwnom)*w21*w34 *
                     (1.0-tspjik)*(1.0-tspijl);
-                  aaa1 = -prefactor*(1.0-(om1234*om1234)) *
+                  aaa1 = -prefactor*(1.0-square(om1234)) *
                     (1.0-tspjik)*(1.0-tspijl);
                   aaa2 = aaa1*w21*w34;
                   at2 = aa*cwnum;
@@ -1836,7 +1842,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag
 ------------------------------------------------------------------------- */
 
 void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
-                              int vflag_atom, ThrData * const thr)
+                              int vflag_atom, double *pv0, ThrData * const thr)
 {
   int i,j,k,m,ii,itype,jtype;
   tagint itag,jtag;
@@ -1921,7 +1927,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
       f[j][1] -= dely*fpair;
       f[j][2] -= delz*fpair;
 
-      if (eflag) evdwl = VR + bij*VA;
+      if (eflag) *pv0 += evdwl = VR + bij*VA;
       if (evflag) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
@@ -1934,7 +1940,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
 ------------------------------------------------------------------------- */
 
 void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
-                            int vflag_atom, ThrData * const thr)
+                            int vflag_atom, double *pv1, ThrData * const thr)
 {
   int i,j,k,m,ii,jj,kk,mm,jnum,itype,jtype,ktype,mtype;
   tagint itag,jtag;
@@ -2218,7 +2224,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
       f[j][1] -= delij[1]*fpair;
       f[j][2] -= delij[2]*fpair;
 
-      if (eflag) evdwl = VA*Stb + (1.0-Str)*cij*VLJ;
+      if (eflag) *pv1 += evdwl = VA*Stb + (1.0-Str)*cij*VLJ;
       if (evflag) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
                                evdwl,0.0,fpair,delij[0],delij[1],delij[2],thr);
 
@@ -2307,8 +2313,8 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
    torsional forces and energy
 ------------------------------------------------------------------------- */
 
-void PairAIREBOOMP::TORSION_thr(int ifrom, int ito,
-                                int evflag, int eflag, ThrData * const thr)
+void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag,
+                                double *pv2, ThrData * const thr)
 {
   int i,j,k,l,ii;
   tagint itag,jtag;
@@ -2461,7 +2467,7 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito,
           Ec = 256.0*ekijl/405.0;
           Vtors = (Ec*(powint(cw2,5)))-(ekijl/10.0);
 
-          if (eflag) evdwl = Vtors*w21*w23*w34*(1.0-tspjik)*(1.0-tspijl);
+          if (eflag) *pv2 += evdwl = Vtors*w21*w23*w34*(1.0-tspjik)*(1.0-tspijl);
 
           dndij[0] = (cross234[1]*del21[2])-(cross234[2]*del21[1]);
           dndij[1] = (cross234[2]*del21[0])-(cross234[0]*del21[2]);
diff --git a/src/USER-OMP/pair_airebo_omp.h b/src/USER-OMP/pair_airebo_omp.h
index 03bc5e5f8ce6ef85b70aca3004ba57cc29630a15..5e7f4d4f584f08d4f36df8870dcd77fcf5dbcca2 100644
--- a/src/USER-OMP/pair_airebo_omp.h
+++ b/src/USER-OMP/pair_airebo_omp.h
@@ -40,10 +40,11 @@ class PairAIREBOOMP : public PairAIREBO, public ThrOMP {
                          int vflag_atom, ThrData * const thr);
 
   void FREBO_thr(int ifrom, int ito, int evflag, int eflag,
-                 int vflag_atom, ThrData * const thr);
+                 int vflag_atom, double *pv0, ThrData * const thr);
   void FLJ_thr(int ifrom, int ito, int evflag, int eflag,
-               int vflag_atom, ThrData * const thr);
-  void TORSION_thr(int ifrom, int ito, int evflag, int eflag, ThrData * const thr);
+               int vflag_atom, double *pv1, ThrData * const thr);
+  void TORSION_thr(int ifrom, int ito, int evflag, int eflag,
+                   double *pv2, ThrData * const thr);
   void REBO_neigh_thr();
 };
 
diff --git a/src/USER-REAXC/compute_spec_atom.cpp b/src/USER-REAXC/compute_spec_atom.cpp
index c719011746e1d60630f7175e22a321f889ff6607..38cf9ad6e6524f495b177f60749b66c0110b2906 100644
--- a/src/USER-REAXC/compute_spec_atom.cpp
+++ b/src/USER-REAXC/compute_spec_atom.cpp
@@ -145,7 +145,7 @@ void ComputeSpecAtom::compute_peratom()
 
   // grow vector or array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     nmax = atom->nmax;
     if (nvalues == 1) {
       memory->destroy(vector);
diff --git a/src/USER-SMD/compute_smd_contact_radius.cpp b/src/USER-SMD/compute_smd_contact_radius.cpp
index e641aea5a253721723eda1dc761b71a52b78589b..991b1a2e674316d4b832aafaa60d620e40e8012d 100644
--- a/src/USER-SMD/compute_smd_contact_radius.cpp
+++ b/src/USER-SMD/compute_smd_contact_radius.cpp
@@ -77,7 +77,7 @@ void ComputeSMDContactRadius::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(contact_radius_vector);
     nmax = atom->nmax;
     contact_radius_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:contact_radius_vector");
diff --git a/src/USER-SMD/compute_smd_damage.cpp b/src/USER-SMD/compute_smd_damage.cpp
index 6630a937352bd96a4cf4691a7268c613ddcd7ac3..7cb4d2b4b8dac1608091c72789d4456b404641d7 100644
--- a/src/USER-SMD/compute_smd_damage.cpp
+++ b/src/USER-SMD/compute_smd_damage.cpp
@@ -77,7 +77,7 @@ void ComputeSMDDamage::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(damage_vector);
     nmax = atom->nmax;
     damage_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:damage_vector");
diff --git a/src/USER-SMD/compute_smd_hourglass_error.cpp b/src/USER-SMD/compute_smd_hourglass_error.cpp
index 99998ed7ee6487731181ba0fac9a6f57a28ee2bb..cf0f67bfbb2f1e0f970aa2e31bd78bfa0c4e592e 100644
--- a/src/USER-SMD/compute_smd_hourglass_error.cpp
+++ b/src/USER-SMD/compute_smd_hourglass_error.cpp
@@ -77,7 +77,7 @@ void ComputeSMDHourglassError::compute_peratom() {
 
 	// grow output Vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->sfree(hourglass_error_vector);
 		nmax = atom->nmax;
 		hourglass_error_vector = (double *) memory->smalloc(nmax * sizeof(double), "atom:hourglass_error_vector");
diff --git a/src/USER-SMD/compute_smd_internal_energy.cpp b/src/USER-SMD/compute_smd_internal_energy.cpp
index fc869f9f93fc58b68904bd5f8bcc613af451c8c6..aabaa8906cb1bf8c7ca5a2be75f9963d98230af2 100644
--- a/src/USER-SMD/compute_smd_internal_energy.cpp
+++ b/src/USER-SMD/compute_smd_internal_energy.cpp
@@ -77,7 +77,7 @@ void ComputeSMDInternalEnergy::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(internal_energy_vector);
     nmax = atom->nmax;
     internal_energy_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:internal_energy_vector");
diff --git a/src/USER-SMD/compute_smd_plastic_strain.cpp b/src/USER-SMD/compute_smd_plastic_strain.cpp
index 7c0977beb09585774232d70e699f50e18898d857..86a5c6329f3a351e51f880473b7443dfb6446416 100644
--- a/src/USER-SMD/compute_smd_plastic_strain.cpp
+++ b/src/USER-SMD/compute_smd_plastic_strain.cpp
@@ -77,7 +77,7 @@ void ComputeSMDPlasticStrain::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(plastic_strain_vector);
     nmax = atom->nmax;
     plastic_strain_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:plastic_strain_vector");
diff --git a/src/USER-SMD/compute_smd_plastic_strain_rate.cpp b/src/USER-SMD/compute_smd_plastic_strain_rate.cpp
index 7df08bd4aca20011d191d43b700fa3731bafdc6f..899e456c8883623f2812de8a207dadb60e69d2e9 100644
--- a/src/USER-SMD/compute_smd_plastic_strain_rate.cpp
+++ b/src/USER-SMD/compute_smd_plastic_strain_rate.cpp
@@ -77,7 +77,7 @@ void ComputeSMDPlasticStrainRate::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(plastic_strain_rate_vector);
     nmax = atom->nmax;
     plastic_strain_rate_vector = (double *) memory->smalloc(nmax*sizeof(double),"atom:plastic_strain_rate_vector");
diff --git a/src/USER-SMD/compute_smd_rho.cpp b/src/USER-SMD/compute_smd_rho.cpp
index 231165dc0cec74e94d5f17e6de5ed50f6eda8e02..91b75d900341b3b4e8def9c12d480cc8730a6d11 100644
--- a/src/USER-SMD/compute_smd_rho.cpp
+++ b/src/USER-SMD/compute_smd_rho.cpp
@@ -75,7 +75,7 @@ void ComputeSMDRho::compute_peratom() {
 
 	// grow rhoVector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->sfree(rhoVector);
 		nmax = atom->nmax;
 		rhoVector = (double *) memory->smalloc(nmax * sizeof(double), "atom:rhoVector");
diff --git a/src/USER-SMD/compute_smd_tlsph_defgrad.cpp b/src/USER-SMD/compute_smd_tlsph_defgrad.cpp
index 50f849b7f1f2376f4058dee5dca4524bf8caff2a..c2ed667ef632b2abbe83aef9524f9dda803f59c3 100644
--- a/src/USER-SMD/compute_smd_tlsph_defgrad.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_defgrad.cpp
@@ -81,7 +81,7 @@ void ComputeSMDTLSPHDefgrad::compute_peratom() {
 	invoked_peratom = update->ntimestep;
 
 	// grow vector array if necessary
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(defgradVector);
 		nmax = atom->nmax;
 		memory->create(defgradVector, nmax, size_peratom_cols, "defgradVector");
diff --git a/src/USER-SMD/compute_smd_tlsph_dt.cpp b/src/USER-SMD/compute_smd_tlsph_dt.cpp
index d2bcbf246f74255503a3dcf4d82b71f16aa6d5b5..54b4f497fb2d5711cde43b40e9a378335bf72959 100644
--- a/src/USER-SMD/compute_smd_tlsph_dt.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_dt.cpp
@@ -77,7 +77,7 @@ void ComputeSMDTlsphDt::compute_peratom() {
 
 	// grow rhoVector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->sfree(dt_vector);
 		nmax = atom->nmax;
 		dt_vector = (double *) memory->smalloc(nmax * sizeof(double),
diff --git a/src/USER-SMD/compute_smd_tlsph_num_neighs.cpp b/src/USER-SMD/compute_smd_tlsph_num_neighs.cpp
index 53e7a5c116091b44349d98b52160800148d2eafa..d0df4ff3d7051d74874ef6c2b50901126321ef39 100644
--- a/src/USER-SMD/compute_smd_tlsph_num_neighs.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_num_neighs.cpp
@@ -72,7 +72,7 @@ void ComputeSMDTLSPHNumNeighs::init() {
 void ComputeSMDTLSPHNumNeighs::compute_peratom() {
     invoked_peratom = update->ntimestep;
 
-    if (atom->nlocal > nmax) {
+    if (atom->nmax > nmax) {
         memory->destroy(numNeighsRefConfigOutput);
         nmax = atom->nmax;
         memory->create(numNeighsRefConfigOutput, nmax, "tlsph/num_neighs:numNeighsRefConfigOutput");
diff --git a/src/USER-SMD/compute_smd_tlsph_shape.cpp b/src/USER-SMD/compute_smd_tlsph_shape.cpp
index 005fe6e39dd1f57ed778f4b679da68f9de516f0a..cec773b5c3151432cbcd35f42ccdcdd6e565112c 100644
--- a/src/USER-SMD/compute_smd_tlsph_shape.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_shape.cpp
@@ -82,7 +82,7 @@ void ComputeSmdTlsphShape::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(strainVector);
 		nmax = atom->nmax;
 		memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
diff --git a/src/USER-SMD/compute_smd_tlsph_strain.cpp b/src/USER-SMD/compute_smd_tlsph_strain.cpp
index 6b1b9c723154a7675c09869224f00af35e5a9c71..2037d6c3543869331f889a4ba698dc670b8a8c06 100644
--- a/src/USER-SMD/compute_smd_tlsph_strain.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_strain.cpp
@@ -83,7 +83,7 @@ void ComputeSMDTLSPHstrain::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(strainVector);
 		nmax = atom->nmax;
 		memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
diff --git a/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp b/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp
index 8ae776da39aa1750cc865a8a17af10c2c933fbdf..db91b080ba8a4460b3b65e3ab7c890459fcadf7e 100644
--- a/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp
@@ -77,7 +77,7 @@ void ComputeSMDTLSPHStrainRate::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(strain_rate_array);
 		nmax = atom->nmax;
 		memory->create(strain_rate_array, nmax, size_peratom_cols, "stresstensorVector");
diff --git a/src/USER-SMD/compute_smd_tlsph_stress.cpp b/src/USER-SMD/compute_smd_tlsph_stress.cpp
index 541cc39eff002b70dd81b93a1e0dfc0efea024d1..d1fce57f9290fdff43ab4406ecaef1539cb800d7 100644
--- a/src/USER-SMD/compute_smd_tlsph_stress.cpp
+++ b/src/USER-SMD/compute_smd_tlsph_stress.cpp
@@ -88,7 +88,7 @@ void ComputeSMDTLSPHStress::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(stress_array);
 		nmax = atom->nmax;
 		memory->create(stress_array, nmax, size_peratom_cols, "stresstensorVector");
diff --git a/src/USER-SMD/compute_smd_triangle_mesh_vertices.cpp b/src/USER-SMD/compute_smd_triangle_mesh_vertices.cpp
index 3e66e346114dcfcaa15337460dc268f968dfdf57..822d62d6db33a1f6c58609a6cb06907bb4d9b166 100644
--- a/src/USER-SMD/compute_smd_triangle_mesh_vertices.cpp
+++ b/src/USER-SMD/compute_smd_triangle_mesh_vertices.cpp
@@ -85,7 +85,7 @@ void ComputeSMDTriangleVertices::compute_peratom() {
 
     // grow vector array if necessary
 
-    if (atom->nlocal > nmax) {
+    if (atom->nmax > nmax) {
         memory->destroy(outputVector);
         nmax = atom->nmax;
         memory->create(outputVector, nmax, size_peratom_cols, "defgradVector");
diff --git a/src/USER-SMD/compute_smd_ulsph_effm.cpp b/src/USER-SMD/compute_smd_ulsph_effm.cpp
index 529ce5bcd6a7df21eeb696ed5e69c57304b01821..87de3df2562ad8c697709789049ef7666726f091 100644
--- a/src/USER-SMD/compute_smd_ulsph_effm.cpp
+++ b/src/USER-SMD/compute_smd_ulsph_effm.cpp
@@ -77,7 +77,7 @@ void ComputeSMD_Ulsph_Effm::compute_peratom() {
 
 	// grow rhoVector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->sfree(dt_vector);
 		nmax = atom->nmax;
 		dt_vector = (double *) memory->smalloc(nmax * sizeof(double),
diff --git a/src/USER-SMD/compute_smd_ulsph_num_neighs.cpp b/src/USER-SMD/compute_smd_ulsph_num_neighs.cpp
index 169eaa7dfb93dbc16fb6b91e9398682a99a2072a..efe6bba47ea16232af896ddd45092e75f682f69e 100644
--- a/src/USER-SMD/compute_smd_ulsph_num_neighs.cpp
+++ b/src/USER-SMD/compute_smd_ulsph_num_neighs.cpp
@@ -72,7 +72,7 @@ void ComputeSMDULSPHNumNeighs::init() {
 void ComputeSMDULSPHNumNeighs::compute_peratom() {
     invoked_peratom = update->ntimestep;
 
-    if (atom->nlocal > nmax) {
+    if (atom->nmax > nmax) {
         memory->destroy(numNeighsOutput);
         nmax = atom->nmax;
         memory->create(numNeighsOutput, nmax, "ulsph/num_neighs:numNeighsRefConfigOutput");
diff --git a/src/USER-SMD/compute_smd_ulsph_strain.cpp b/src/USER-SMD/compute_smd_ulsph_strain.cpp
index 395221d683775a0aab34bc513fa448bbe1bd08dc..2ecb79e6700854b22c49ef25d9399fe66d17fd2a 100644
--- a/src/USER-SMD/compute_smd_ulsph_strain.cpp
+++ b/src/USER-SMD/compute_smd_ulsph_strain.cpp
@@ -83,7 +83,7 @@ void ComputeSMDULSPHstrain::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(strainVector);
 		nmax = atom->nmax;
 		memory->create(strainVector, nmax, size_peratom_cols, "strainVector");
diff --git a/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp b/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp
index e35331c64e91c115cbd7c9303ee2b4450839391f..4ba4c6c52003cd5bfb635d6859edc28d836f16dd 100644
--- a/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp
+++ b/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp
@@ -77,7 +77,7 @@ void ComputeSMDULSPHStrainRate::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(strain_rate_array);
 		nmax = atom->nmax;
 		memory->create(strain_rate_array, nmax, size_peratom_cols, "stresstensorVector");
diff --git a/src/USER-SMD/compute_smd_ulsph_stress.cpp b/src/USER-SMD/compute_smd_ulsph_stress.cpp
index 9be8646f132f31653bb6b9163c1cc414db3ddef5..e26e0c12f23f29800fbb3f6701e28437f04ba0c2 100644
--- a/src/USER-SMD/compute_smd_ulsph_stress.cpp
+++ b/src/USER-SMD/compute_smd_ulsph_stress.cpp
@@ -88,7 +88,7 @@ void ComputeSMDULSPHStress::compute_peratom() {
 
 	// grow vector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->destroy(stress_array);
 		nmax = atom->nmax;
 		memory->create(stress_array, nmax, size_peratom_cols, "stresstensorVector");
diff --git a/src/USER-SMD/compute_smd_vol.cpp b/src/USER-SMD/compute_smd_vol.cpp
index 658ede643982345021e6aeeeaad6b82e9ab818ba..a12164d9e253eadffaf68cfcd850644488f9c0f4 100644
--- a/src/USER-SMD/compute_smd_vol.cpp
+++ b/src/USER-SMD/compute_smd_vol.cpp
@@ -76,7 +76,7 @@ void ComputeSMDVol::compute_peratom() {
 
 	// grow volVector array if necessary
 
-	if (atom->nlocal > nmax) {
+	if (atom->nmax > nmax) {
 		memory->sfree(volVector);
 		nmax = atom->nmax;
 		volVector = (double *) memory->smalloc(nmax * sizeof(double), "atom:volVector");
diff --git a/src/USER-SMD/fix_smd_setvel.cpp b/src/USER-SMD/fix_smd_setvel.cpp
index 30ebbeb332123d2e7ada65b3093802be190b16c1..568ffb79775fffaf66f17f809edf27b05dd1bd97 100644
--- a/src/USER-SMD/fix_smd_setvel.cpp
+++ b/src/USER-SMD/fix_smd_setvel.cpp
@@ -257,7 +257,7 @@ void FixSMDSetVel::post_force(int vflag) {
 
 	// reallocate sforce array if necessary
 
-	if (varflag == ATOM && nlocal > maxatom) {
+	if (varflag == ATOM && atom->nmax > maxatom) {
 		maxatom = atom->nmax;
 		memory->destroy(sforce);
 		memory->create(sforce, maxatom, 3, "setvelocity:sforce");
diff --git a/src/USER-SPH/compute_meso_e_atom.cpp b/src/USER-SPH/compute_meso_e_atom.cpp
index 51035decf80a0a0611cc2a68f0d0a6bb7d0bf3d6..dc8d0b55ad32f32d27f1bb14da002cfe57b359ec 100644
--- a/src/USER-SPH/compute_meso_e_atom.cpp
+++ b/src/USER-SPH/compute_meso_e_atom.cpp
@@ -65,7 +65,7 @@ void ComputeMesoEAtom::compute_peratom()
 
   // grow evector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(evector);
     nmax = atom->nmax;
     evector = (double *) memory->smalloc(nmax*sizeof(double),"evector/atom:evector");
diff --git a/src/USER-SPH/compute_meso_rho_atom.cpp b/src/USER-SPH/compute_meso_rho_atom.cpp
index d30d62c82e50e7b51c1d530f9c85bf45a1b6f675..9b4da8b339f55b53b07bcc0b4de11fa1b1098287 100644
--- a/src/USER-SPH/compute_meso_rho_atom.cpp
+++ b/src/USER-SPH/compute_meso_rho_atom.cpp
@@ -65,7 +65,7 @@ void ComputeMesoRhoAtom::compute_peratom()
 
   // grow rhoVector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(rhoVector);
     nmax = atom->nmax;
     rhoVector = (double *) memory->smalloc(nmax*sizeof(double),"atom:rhoVector");
diff --git a/src/USER-SPH/compute_meso_t_atom.cpp b/src/USER-SPH/compute_meso_t_atom.cpp
index f14e0c83024f8115bd0e6366a9f15aea3d136ab5..d901f73eff5feaba1490d349e6c3757e26113b69 100644
--- a/src/USER-SPH/compute_meso_t_atom.cpp
+++ b/src/USER-SPH/compute_meso_t_atom.cpp
@@ -66,7 +66,7 @@ void ComputeMesoTAtom::compute_peratom()
 
   // grow tvector array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->sfree(tvector);
     nmax = atom->nmax;
     tvector = (double *) memory->smalloc(nmax*sizeof(double),"tvector/atom:tvector");
diff --git a/src/USER-VTK/dump_custom_vtk.cpp b/src/USER-VTK/dump_custom_vtk.cpp
index 8db751e291aeaf863fff62b8a61a69ebf6bb9cfb..0e4bc459766a69b8864724454b88a4f5997547ab 100644
--- a/src/USER-VTK/dump_custom_vtk.cpp
+++ b/src/USER-VTK/dump_custom_vtk.cpp
@@ -266,7 +266,7 @@ int DumpCustomVTK::count()
   // grow choose and variable vbuf arrays if needed
 
   int nlocal = atom->nlocal;
-  if (nlocal > maxlocal) {
+  if (atom->nmax > maxlocal) {
     maxlocal = atom->nmax;
 
     memory->destroy(choose);
diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp
index ffdf5b7fadb460b80124f6aace81435e97672191..f4f48fd5e8f40f35f9a80c8ff9ee73c8c8a323d0 100644
--- a/src/VORONOI/compute_voronoi_atom.cpp
+++ b/src/VORONOI/compute_voronoi_atom.cpp
@@ -188,8 +188,7 @@ void ComputeVoronoi::compute_peratom()
   invoked_peratom = update->ntimestep;
 
   // grow per atom array if necessary
-  int nlocal = atom->nlocal;
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(voro);
     nmax = atom->nmax;
     memory->create(voro,nmax,size_peratom_cols,"voronoi/atom:voro");
@@ -199,7 +198,7 @@ void ComputeVoronoi::compute_peratom()
   // decide between occupation or per-frame tesselation modes
   if (occupation) {
     // build cells only once
-    int i, nall = nlocal + atom->nghost;
+    int i, nall = atom->nlocal + atom->nghost;
     if (con_mono==NULL && con_poly==NULL) {
       // generate the voronoi cell network for the initial structure
       buildCells();
@@ -320,7 +319,7 @@ void ComputeVoronoi::buildCells()
     if (!input->variable->atomstyle(radvar))
       error->all(FLERR,"Variable for voronoi radius is not atom style");
     // prepare destination buffer for variable evaluation
-    if (nlocal > rmax) {
+    if (atom->nmax > rmax) {
       memory->destroy(rfield);
       rmax = atom->nmax;
       memory->create(rfield,rmax,"voronoi/atom:rfield");
@@ -379,7 +378,7 @@ void ComputeVoronoi::checkOccupation()
          **x = atom->x;
 
   // prepare destination buffer for variable evaluation
-  if (nall > lmax) {
+  if (atom->nmax > lmax) {
     memory->destroy(lnext);
     lmax = atom->nmax;
     memory->create(lnext,lmax,"voronoi/atom:lnext");
diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp
index a6561f31368c5002a5f5eb6fd6a39839a9916af1..6c35ba1e8244014cd49516be139948eca11d49a6 100644
--- a/src/compute_centro_atom.cpp
+++ b/src/compute_centro_atom.cpp
@@ -107,7 +107,7 @@ void ComputeCentroAtom::compute_peratom()
 
   // grow centro array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(centro);
     nmax = atom->nmax;
     memory->create(centro,nmax,"centro/atom:centro");
diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp
index c968073f07f1809a01167e6a51915f4ef629840a..4d9a5f2275cdfff75881c71e9928357b41acd647 100644
--- a/src/compute_chunk_atom.cpp
+++ b/src/compute_chunk_atom.cpp
@@ -627,7 +627,7 @@ void ComputeChunkAtom::compute_peratom()
 
   // grow floating point chunk vector if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(chunk);
     nmax = atom->nmax;
     memory->create(chunk,nmax,"chunk/atom:chunk");
@@ -890,7 +890,7 @@ void ComputeChunkAtom::assign_chunk_ids()
 
   // grow integer chunk index vector if necessary
 
-  if (atom->nlocal > nmaxint) {
+  if (atom->nmax > nmaxint) {
     memory->destroy(ichunk);
     memory->destroy(exclude);
     nmaxint = atom->nmax;
@@ -988,7 +988,7 @@ void ComputeChunkAtom::assign_chunk_ids()
     }
 
   } else if (which == VARIABLE) {
-    if (nlocal > maxvar) {
+    if (atom->nmax > maxvar) {
       maxvar = atom->nmax;
       memory->destroy(varatom);
       memory->create(varatom,maxvar,"chunk/atom:varatom");
diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp
index ade4ce5d87ef5e84fecf86cef11ae5744ea06e2d..c81785f2152464bf293144efdea23ed1004ba0e7 100644
--- a/src/compute_cluster_atom.cpp
+++ b/src/compute_cluster_atom.cpp
@@ -104,7 +104,7 @@ void ComputeClusterAtom::compute_peratom()
 
   // grow clusterID array if necessary
 
-  if (atom->nlocal+atom->nghost > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(clusterID);
     nmax = atom->nmax;
     memory->create(clusterID,nmax,"cluster/atom:clusterID");
diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp
index 4e0d7e907bdb9cbc951d0f6451be6f368fba6be2..46c6726b0007edcf508b37522ace3c1fd166d69c 100644
--- a/src/compute_cna_atom.cpp
+++ b/src/compute_cna_atom.cpp
@@ -123,7 +123,7 @@ void ComputeCNAAtom::compute_peratom()
 
   // grow arrays if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(nearest);
     memory->destroy(nnearest);
     memory->destroy(pattern);
diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp
index 2769edc12dc5b6ab76ef086be94bf60ba2932d21..c6aa561eb6be808fe8e736bb635a85c3e16bdff3 100644
--- a/src/compute_coord_atom.cpp
+++ b/src/compute_coord_atom.cpp
@@ -125,7 +125,7 @@ void ComputeCoordAtom::compute_peratom()
 
   // grow coordination array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     if (ncol == 1) {
       memory->destroy(cvec);
       nmax = atom->nmax;
diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp
index 51d8ab22dee38452663386472d3b3cb44a72ccdb..f274c32f5db8908a45667fcd61449f7063a56949 100644
--- a/src/compute_displace_atom.cpp
+++ b/src/compute_displace_atom.cpp
@@ -110,7 +110,7 @@ void ComputeDisplaceAtom::compute_peratom()
 
   // grow local displacement array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(displace);
     nmax = atom->nmax;
     memory->create(displace,nmax,4,"displace/atom:displace");
diff --git a/src/compute_erotate_sphere_atom.cpp b/src/compute_erotate_sphere_atom.cpp
index 9c2d4d9caf1fed7c5fb4af4ed0ce9ac7fb4db3c4..ed6946e8c87abf9521654759f3564627624c4ee7 100644
--- a/src/compute_erotate_sphere_atom.cpp
+++ b/src/compute_erotate_sphere_atom.cpp
@@ -74,7 +74,7 @@ void ComputeErotateSphereAtom::compute_peratom()
 
   // grow erot array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(erot);
     nmax = atom->nmax;
     memory->create(erot,nmax,"erotate/sphere/atom:erot");
diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp
index d1257e9d8b8fee8301191a7653013a17673e0108..cf9bd6d09f3cf0dcc4b0951eed35977c16558273 100644
--- a/src/compute_hexorder_atom.cpp
+++ b/src/compute_hexorder_atom.cpp
@@ -148,7 +148,7 @@ void ComputeHexOrderAtom::compute_peratom()
 
   // grow order parameter array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(qnarray);
     nmax = atom->nmax;
     memory->create(qnarray,nmax,ncol,"hexorder/atom:qnarray");
diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp
index 5cc6f7a9395d5e6339ffd2b17146fea00edacf9f..aa893981a40542860ec4beab09f012395f24a764 100644
--- a/src/compute_ke_atom.cpp
+++ b/src/compute_ke_atom.cpp
@@ -63,7 +63,7 @@ void ComputeKEAtom::compute_peratom()
 
   // grow ke array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(ke);
     nmax = atom->nmax;
     memory->create(ke,nmax,"ke/atom:ke");
diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp
index 7a87cc2915d6695e9a8146536e4ac5f9593891c0..fb61726c025f57a312fe7f605a3bb96c7481de74 100644
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@@ -180,7 +180,7 @@ void ComputeOrientOrderAtom::compute_peratom()
 
   // grow order parameter array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(qnarray);
     nmax = atom->nmax;
     memory->create(qnarray,nmax,ncol,"orientorder/atom:qnarray");
diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp
index 54fcb9d722dc367e9485b76693d2b369e88bb2e0..c9575dba5f78a87047322f0e6ad7bc5d18dfd63a 100644
--- a/src/compute_property_atom.cpp
+++ b/src/compute_property_atom.cpp
@@ -384,7 +384,7 @@ void ComputePropertyAtom::compute_peratom()
 
   // grow vector or array if necessary
 
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     nmax = atom->nmax;
     if (nvalues == 1) {
       memory->destroy(vector);
diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index 47d85332beaf3bfb9337620df673ea6b7adb698c..a56770569f9d3264948eb570f676196df1460c0f 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -574,7 +574,7 @@ double ComputeReduce::compute_one(int m, int flag)
   // evaluate atom-style variable
 
   } else if (which[m] == VARIABLE) {
-    if (nlocal > maxatom) {
+    if (atom->nmax > maxatom) {
       maxatom = atom->nmax;
       memory->destroy(varatom);
       memory->create(varatom,maxatom,"reduce:varatom");
diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp
index c098942c75b913b14a431fa8fc556fcd2443ceb9..c867a0190ccde7331cb0060304124570891fd7f7 100644
--- a/src/compute_reduce_region.cpp
+++ b/src/compute_reduce_region.cpp
@@ -200,7 +200,7 @@ double ComputeReduceRegion::compute_one(int m, int flag)
   // evaluate atom-style variable
 
   } else if (which[m] == VARIABLE) {
-    if (nlocal > maxatom) {
+    if (atom->nmax > maxatom) {
       maxatom = atom->nmax;
       memory->destroy(varatom);
       memory->create(varatom,maxatom,"reduce/region:varatom");
diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp
index 8c1d2d85fcb81b06bdb4d948882169c2172fef63..5af995252c6f5b2df516b4fc680a17ed93bb9d47 100644
--- a/src/compute_temp_deform.cpp
+++ b/src/compute_temp_deform.cpp
@@ -231,7 +231,7 @@ void ComputeTempDeform::remove_bias_all()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/deform:vbiasall");
diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index b330495c5461b076b82803c602fce79991a5ec71..7678403d881dd85eec8d4eff6964c7088780e6d2 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -195,7 +195,7 @@ void ComputeTempPartial::remove_bias_all()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/partial:vbiasall");
diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp
index ba5ed8ca68df9caa3d1cfc47f83c7e7e36cb850d..810d6dd08bafe09d06da4c25ff235fdf9091e566 100644
--- a/src/compute_temp_ramp.cpp
+++ b/src/compute_temp_ramp.cpp
@@ -244,7 +244,7 @@ void ComputeTempRamp::remove_bias_all()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  if (nlocal > maxbias) {
+  if (atom->nmax > maxbias) {
     memory->destroy(vbiasall);
     maxbias = atom->nmax;
     memory->create(vbiasall,maxbias,3,"temp/ramp:vbiasall");
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 68097e543d1918e806ba5d5be4671e6dcce5b10b..1e7356d0a164f1bf4814396bf4c95a6db50d6dd5 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -386,7 +386,7 @@ int DumpCustom::count()
   // grow choose and variable vbuf arrays if needed
 
   int nlocal = atom->nlocal;
-  if (nlocal > maxlocal) {
+  if (atom->nmax > maxlocal) {
     maxlocal = atom->nmax;
 
     memory->destroy(choose);
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index 2ab9e84467abb4a5964d50bd30873d38c8daec8b..b9ff534c6616023dbdb8780f8486b10fcfc241c0 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -901,7 +901,7 @@ void DumpImage::create_image()
     // communicate choose flag for ghost atoms to know if they are selected
     // if bcolor/bdiam = ATOM, setup bufcopy to comm atom color/diam attributes
 
-    if (nall > maxbufcopy) {
+    if (atom->nmax > maxbufcopy) {
       maxbufcopy = atom->nmax;
       memory->destroy(chooseghost);
       memory->create(chooseghost,maxbufcopy,"dump:chooseghost");
diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp
index 0ffa7751cd4f14d3fe8dbe4bcef496b7f71a3b9c..2b5c682266aba970424e5341027bfd9ab7c0b5f5 100644
--- a/src/fix_addforce.cpp
+++ b/src/fix_addforce.cpp
@@ -251,7 +251,7 @@ void FixAddForce::post_force(int vflag)
 
   // reallocate sforce array if necessary
 
-  if ((varflag == ATOM || estyle == ATOM) && nlocal > maxatom) {
+  if ((varflag == ATOM || estyle == ATOM) && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     memory->destroy(sforce);
     memory->create(sforce,maxatom,4,"addforce:sforce");
diff --git a/src/fix_ave_chunk.cpp b/src/fix_ave_chunk.cpp
index ff6210d215384c76e4df54492f985e44aeb258e4..8390e8bbad781d0a6d101ca77288e6d328b2e431 100644
--- a/src/fix_ave_chunk.cpp
+++ b/src/fix_ave_chunk.cpp
@@ -721,7 +721,7 @@ void FixAveChunk::end_of_step()
     // evaluate atom-style variable
 
     } else if (which[m] == VARIABLE) {
-      if (nlocal > maxvar) {
+      if (atom->nmax > maxvar) {
         maxvar = atom->nmax;
         memory->destroy(varatom);
         memory->create(varatom,maxvar,"ave/chunk:varatom");
diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp
index 902934be0f1c369f0f690d4370e907b4abccaa35..9123a8561f33489c11865cd9aa732b3b111d94fc 100644
--- a/src/fix_ave_histo.cpp
+++ b/src/fix_ave_histo.cpp
@@ -754,7 +754,7 @@ void FixAveHisto::end_of_step()
         bin_vector(nvec,varvec,1);
 
       } else if (which[i] == VARIABLE && kind == PERATOM) {
-        if (atom->nlocal > maxatom) {
+        if (atom->nmax > maxatom) {
           memory->destroy(vector);
           maxatom = atom->nmax;
           memory->create(vector,maxatom,"ave/histo:vector");
diff --git a/src/fix_ave_histo_weight.cpp b/src/fix_ave_histo_weight.cpp
index 47901a516404dca2afb541d50928fb9a2f4ba385..8671a38a932a0ed9d9d18384b5c86de19c7a4ff8 100644
--- a/src/fix_ave_histo_weight.cpp
+++ b/src/fix_ave_histo_weight.cpp
@@ -256,7 +256,7 @@ void FixAveHistoWeight::end_of_step()
     weight = input->variable->compute_equal(m);
 
   } else if (which[i] == VARIABLE && kind == PERATOM) {
-    if (atom->nlocal > maxatom) {
+    if (atom->nmax > maxatom) {
       memory->destroy(vector);
       maxatom = atom->nmax;
       memory->create(vector,maxatom,"ave/histo/weight:vector");
@@ -385,7 +385,7 @@ void FixAveHistoWeight::end_of_step()
     bin_one_weights(input->variable->compute_equal(m),weight);
 
   } else if (which[i] == VARIABLE && kind == PERATOM) {
-    if (atom->nlocal > maxatom) {
+    if (atom->nmax > maxatom) {
       memory->destroy(vector);
       maxatom = atom->nmax;
       memory->create(vector,maxatom,"ave/histo/weight:vector");
diff --git a/src/fix_heat.cpp b/src/fix_heat.cpp
index 65dc960f424d7faf9f4d41683e80391ea45c753f..0609d45efb806491d3cab729ec6d66b9380f7410 100644
--- a/src/fix_heat.cpp
+++ b/src/fix_heat.cpp
@@ -153,7 +153,7 @@ void FixHeat::end_of_step()
 
   // reallocate per-atom arrays if necessary
 
-  if (hstyle == ATOM && atom->nlocal > maxatom) {
+  if (hstyle == ATOM && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     memory->destroy(vheat);
     memory->destroy(vscale);
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index 9e857339b4bbfbc86d271fdee90041bbd0b68f73..ff8d4f18408a703abe08e1245972735c8c0db490 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -511,7 +511,7 @@ void FixLangevin::post_force_untemplated
   // reallocate flangevin if necessary
 
   if (Tp_TALLY) {
-    if (atom->nlocal > maxatom1) {
+    if (atom->nmax > maxatom1) {
       memory->destroy(flangevin);
       maxatom1 = atom->nmax;
       memory->create(flangevin,maxatom1,3,"langevin:flangevin");
@@ -641,7 +641,7 @@ void FixLangevin::compute_target()
         error->one(FLERR,"Fix langevin variable returned negative temperature");
       tsqrt = sqrt(t_target);
     } else {
-      if (nlocal > maxatom2) {
+      if (atom->nmax > maxatom2) {
         maxatom2 = atom->nmax;
         memory->destroy(tforce);
         memory->create(tforce,maxatom2,"langevin:tforce");
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
index ca3727520caa90598f69c8cfcb9622247e3a7e92..592419cc1e90d0674147d3027c139b15e2ad7041 100644
--- a/src/fix_move.cpp
+++ b/src/fix_move.cpp
@@ -726,7 +726,7 @@ void FixMove::initial_integrate(int vflag)
 
     // reallocate displace and velocity arrays as necessary
 
-    if ((displaceflag || velocityflag) && nlocal > maxatom) {
+    if ((displaceflag || velocityflag) && atom->nmax > maxatom) {
       maxatom = atom->nmax;
       if (displaceflag) {
         memory->destroy(displace);
diff --git a/src/fix_setforce.cpp b/src/fix_setforce.cpp
index aff8f117c5b425494cd7e35250c6edb38988f481..d73ab0db63faf5b72169918dcceec6c136eaf02e 100644
--- a/src/fix_setforce.cpp
+++ b/src/fix_setforce.cpp
@@ -231,7 +231,7 @@ void FixSetForce::post_force(int vflag)
 
   // reallocate sforce array if necessary
 
-  if (varflag == ATOM && nlocal > maxatom) {
+  if (varflag == ATOM && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     memory->destroy(sforce);
     memory->create(sforce,maxatom,3,"setforce:sforce");
diff --git a/src/fix_store_force.cpp b/src/fix_store_force.cpp
index 521e4568912ae6db3cfef5e950959422c3e6429d..f0852092871c716a547711b96969c1e744904413 100644
--- a/src/fix_store_force.cpp
+++ b/src/fix_store_force.cpp
@@ -96,7 +96,7 @@ void FixStoreForce::min_setup(int vflag)
 
 void FixStoreForce::post_force(int vflag)
 {
-  if (atom->nlocal > nmax) {
+  if (atom->nmax > nmax) {
     nmax = atom->nmax;
     memory->destroy(foriginal);
     memory->create(foriginal,nmax,3,"store/force:foriginal");
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index c30c8a74bd42265ecd18e264991e208c890cc9d3..d0fdb1407a681ae7d08c0d344587e952c76fa76d 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1542,7 +1542,7 @@ void Neighbor::build(int topoflag)
     double **x = atom->x;
     int nlocal = atom->nlocal;
     if (includegroup) nlocal = atom->nfirst;
-    if (nlocal > maxhold) {
+    if (atom->nmax > maxhold) {
       maxhold = atom->nmax;
       memory->destroy(xhold);
       memory->create(xhold,maxhold,3,"neigh:xhold");
@@ -1577,10 +1577,10 @@ void Neighbor::build(int topoflag)
   // else only invoke grow() if nlocal exceeds previous list size
   // only for lists with growflag set and which are perpetual (glist)
 
-  if (anyghostlist && atom->nlocal+atom->nghost > maxatom) {
+  if (anyghostlist && atom->nmax > maxatom) {
     maxatom = atom->nmax;
     for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxatom);
-  } else if (atom->nlocal > maxatom) {
+  } else if (atom->nmax > maxatom) {
     maxatom = atom->nmax;
     for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxatom);
   }