From 876b033ea86239715bbf36b9862b7b878ba1041c Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 31 May 2016 16:33:07 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15083
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/CORESHELL/compute_temp_cs.cpp             |  2 +-
 src/KOKKOS/atom_kokkos.cpp                    |  4 +-
 src/KOKKOS/fix_langevin_kokkos.cpp            |  2 +-
 src/KOKKOS/fix_setforce_kokkos.cpp            |  2 +-
 src/KOKKOS/neighbor_kokkos.cpp                |  6 +-
 src/KSPACE/ewald.cpp                          |  2 +-
 src/KSPACE/ewald_disp.cpp                     |  2 +-
 src/KSPACE/msm.cpp                            |  2 +-
 src/KSPACE/msm_cg.cpp                         |  2 +-
 src/KSPACE/pppm.cpp                           |  2 +-
 src/KSPACE/pppm_cg.cpp                        |  2 +-
 src/KSPACE/pppm_disp.cpp                      |  2 +-
 src/KSPACE/pppm_stagger.cpp                   |  2 +-
 src/MANYBODY/pair_airebo.cpp                  | 35 +++++-----
 src/MC/fix_atom_swap.cpp                      |  4 +-
 src/MC/fix_bond_swap.cpp                      |  2 +-
 src/MC/fix_gcmc.cpp                           |  2 +-
 src/MISC/fix_efield.cpp                       |  2 +-
 src/MISC/fix_evaporate.cpp                    |  2 +-
 src/PERI/compute_damage_atom.cpp              |  2 +-
 src/PERI/compute_dilatation_atom.cpp          |  2 +-
 src/PERI/compute_plasticity_atom.cpp          |  2 +-
 src/REPLICA/verlet_split.cpp                  |  2 +-
 src/SNAP/compute_sna_atom.cpp                 |  2 +-
 src/SNAP/compute_snad_atom.cpp                |  2 +-
 src/SNAP/compute_snav_atom.cpp                |  2 +-
 src/SRD/fix_srd.cpp                           |  2 +-
 src/USER-DPD/compute_dpd_atom.cpp             |  4 +-
 src/USER-EFF/compute_ke_atom_eff.cpp          |  2 +-
 src/USER-EFF/compute_temp_deform_eff.cpp      |  2 +-
 src/USER-EFF/compute_temp_region_eff.cpp      |  2 +-
 src/USER-EFF/fix_langevin_eff.cpp             |  6 +-
 src/USER-INTEL/pppm_intel.cpp                 |  2 +-
 src/USER-MISC/compute_ackland_atom.cpp        |  2 +-
 src/USER-MISC/compute_basal_atom.cpp          |  2 +-
 src/USER-MISC/compute_temp_rotate.cpp         |  4 +-
 src/USER-OMP/ewald_omp.cpp                    |  2 +-
 src/USER-OMP/msm_cg_omp.cpp                   |  2 +-
 src/USER-OMP/pair_airebo_omp.cpp              | 66 ++++++++++---------
 src/USER-OMP/pair_airebo_omp.h                |  7 +-
 src/USER-REAXC/compute_spec_atom.cpp          |  2 +-
 src/USER-SMD/compute_smd_contact_radius.cpp   |  2 +-
 src/USER-SMD/compute_smd_damage.cpp           |  2 +-
 src/USER-SMD/compute_smd_hourglass_error.cpp  |  2 +-
 src/USER-SMD/compute_smd_internal_energy.cpp  |  2 +-
 src/USER-SMD/compute_smd_plastic_strain.cpp   |  2 +-
 .../compute_smd_plastic_strain_rate.cpp       |  2 +-
 src/USER-SMD/compute_smd_rho.cpp              |  2 +-
 src/USER-SMD/compute_smd_tlsph_defgrad.cpp    |  2 +-
 src/USER-SMD/compute_smd_tlsph_dt.cpp         |  2 +-
 src/USER-SMD/compute_smd_tlsph_num_neighs.cpp |  2 +-
 src/USER-SMD/compute_smd_tlsph_shape.cpp      |  2 +-
 src/USER-SMD/compute_smd_tlsph_strain.cpp     |  2 +-
 .../compute_smd_tlsph_strain_rate.cpp         |  2 +-
 src/USER-SMD/compute_smd_tlsph_stress.cpp     |  2 +-
 .../compute_smd_triangle_mesh_vertices.cpp    |  2 +-
 src/USER-SMD/compute_smd_ulsph_effm.cpp       |  2 +-
 src/USER-SMD/compute_smd_ulsph_num_neighs.cpp |  2 +-
 src/USER-SMD/compute_smd_ulsph_strain.cpp     |  2 +-
 .../compute_smd_ulsph_strain_rate.cpp         |  2 +-
 src/USER-SMD/compute_smd_ulsph_stress.cpp     |  2 +-
 src/USER-SMD/compute_smd_vol.cpp              |  2 +-
 src/USER-SMD/fix_smd_setvel.cpp               |  2 +-
 src/USER-SPH/compute_meso_e_atom.cpp          |  2 +-
 src/USER-SPH/compute_meso_rho_atom.cpp        |  2 +-
 src/USER-SPH/compute_meso_t_atom.cpp          |  2 +-
 src/USER-VTK/dump_custom_vtk.cpp              |  2 +-
 src/VORONOI/compute_voronoi_atom.cpp          |  9 ++-
 src/compute_centro_atom.cpp                   |  2 +-
 src/compute_chunk_atom.cpp                    |  6 +-
 src/compute_cluster_atom.cpp                  |  2 +-
 src/compute_cna_atom.cpp                      |  2 +-
 src/compute_coord_atom.cpp                    |  2 +-
 src/compute_displace_atom.cpp                 |  2 +-
 src/compute_erotate_sphere_atom.cpp           |  2 +-
 src/compute_hexorder_atom.cpp                 |  2 +-
 src/compute_ke_atom.cpp                       |  2 +-
 src/compute_orientorder_atom.cpp              |  2 +-
 src/compute_property_atom.cpp                 |  2 +-
 src/compute_reduce.cpp                        |  2 +-
 src/compute_reduce_region.cpp                 |  2 +-
 src/compute_temp_deform.cpp                   |  2 +-
 src/compute_temp_partial.cpp                  |  2 +-
 src/compute_temp_ramp.cpp                     |  2 +-
 src/dump_custom.cpp                           |  2 +-
 src/dump_image.cpp                            |  2 +-
 src/fix_addforce.cpp                          |  2 +-
 src/fix_ave_chunk.cpp                         |  2 +-
 src/fix_ave_histo.cpp                         |  2 +-
 src/fix_ave_histo_weight.cpp                  |  4 +-
 src/fix_heat.cpp                              |  2 +-
 src/fix_langevin.cpp                          |  4 +-
 src/fix_move.cpp                              |  2 +-
 src/fix_setforce.cpp                          |  2 +-
 src/fix_store_force.cpp                       |  2 +-
 src/neighbor.cpp                              |  6 +-
 96 files changed, 170 insertions(+), 159 deletions(-)

diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
index b10a9256b8..d0ff5d290d 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 0e2e6038c1..4d343367f1 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 de8920b556..43af4168f3 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 862aab1e21..93e5e19377 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 c018eebb1b..55fc0f1e9b 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 f3a4188de9..39a98901cb 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 c1a8d81e00..c567ce0168 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 7341688bb2..b9d894a00a 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 c7af52ef89..31f70b1934 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 a4cc6f89df..9b18ad88bb 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 6fcdb438a8..3e5e4431a6 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 2c065e926f..b7ec188150 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 7cd719ec73..a16c0d474b 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 49728d6aa6..d4e57c2c6c 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 689aa0a168..41a54a64ff 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 0c36cecb30..31194ce12e 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 c2a69399ce..aea029e5ba 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 d729b9b74a..2151a05092 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 873227449d..42f3b844db 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 cf8cf35cd3..305ab4e716 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 20993fab47..d79b80f4d4 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 2c47f18b1f..34fd23fe6a 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 5985c756ad..80d8f58f0d 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 1f7c943719..da948e56cc 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 8767ca1745..2164b31029 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 cbe53821b8..29431f29a9 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 35eb21254c..ae11147327 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 9629aa5435..8d944a292f 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 93ad934d6d..8d03e97b0a 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 12c72479cb..40e8e529c0 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 1035a3e31b..c0cd56f3a8 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 f2ca444de2..9584623fc3 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 466b924a51..bd422aac82 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 7a60117bb8..ac48d75bc6 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 d323e20a10..881b253131 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 574db8757c..2210555a7c 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 7ab903aab1..ad6b0bf805 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 f743abae6f..9f1ebba2a8 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 eeb3e134f6..35a644007a 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 03bc5e5f8c..5e7f4d4f58 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 c719011746..38cf9ad6e6 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 e641aea5a2..991b1a2e67 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 6630a93735..7cb4d2b4b8 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 99998ed7ee..cf0f67bfbb 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 fc869f9f93..aabaa8906c 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 7c0977beb0..86a5c6329f 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 7df08bd4ac..899e456c88 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 231165dc0c..91b75d9003 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 50f849b7f1..c2ed667ef6 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 d2bcbf246f..54b4f497fb 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 53e7a5c116..d0df4ff3d7 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 005fe6e39d..cec773b5c3 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 6b1b9c7231..2037d6c354 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 8ae776da39..db91b080ba 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 541cc39eff..d1fce57f92 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 3e66e34611..822d62d6db 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 529ce5bcd6..87de3df256 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 169eaa7dfb..efe6bba47e 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 395221d683..2ecb79e670 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 e35331c64e..4ba4c6c520 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 9be8646f13..e26e0c12f2 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 658ede6439..a12164d9e2 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 30ebbeb332..568ffb7977 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 51035decf8..dc8d0b55ad 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 d30d62c82e..9b4da8b339 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 f14e0c8302..d901f73eff 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 8db751e291..0e4bc45976 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 ffdf5b7fad..f4f48fd5e8 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 a6561f3136..6c35ba1e82 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 c968073f07..4d9a5f2275 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 ade4ce5d87..c81785f215 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 4e0d7e907b..46c6726b00 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 2769edc12d..c6aa561eb6 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 51d8ab22de..f274c32f5d 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 9c2d4d9caf..ed6946e8c8 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 d1257e9d8b..cf9bd6d09f 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 5cc6f7a939..aa893981a4 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 7a87cc2915..fb61726c02 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 54fcb9d722..c9575dba5f 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 47d85332be..a56770569f 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 c098942c75..c867a0190c 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 8c1d2d85fc..5af995252c 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 b330495c54..7678403d88 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 ba5ed8ca68..810d6dd08b 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 68097e543d..1e7356d0a1 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 2ab9e84467..b9ff534c66 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 0ffa7751cd..2b5c682266 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 ff6210d215..8390e8bbad 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 902934be0f..9123a8561f 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 47901a5164..8671a38a93 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 65dc960f42..0609d45efb 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 9e857339b4..ff8d4f1840 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 ca3727520c..592419cc1e 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 aff8f117c5..d73ab0db63 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 521e456891..f085209287 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 c30c8a74bd..d0fdb1407a 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);
   }
-- 
GitLab