diff --git a/lib/gpu/Makefile.linux b/lib/gpu/Makefile.linux
index 1e689a3..d72c0ba 100644
--- a/lib/gpu/Makefile.linux
+++ b/lib/gpu/Makefile.linux
@@ -7,7 +7,7 @@
 
 EXTRAMAKE = Makefile.lammps.standard
 
-ifeq($(CUDA_HOME),)
+ifeq ($(CUDA_HOME),)
 CUDA_HOME = /usr/local/cuda
 endif
 
diff --git a/lib/gpu/lal_base_three.cpp b/lib/gpu/lal_base_three.cpp
index c41aad7..14f642e 100644
--- a/lib/gpu/lal_base_three.cpp
+++ b/lib/gpu/lal_base_three.cpp
@@ -12,7 +12,7 @@
     begin                : Tue April 2, 2013
     email                : brownw@ornl.gov
  ***************************************************************************/
- 
+
 #include "lal_base_three.h"
 using namespace LAMMPS_AL;
 #define BaseThreeT BaseThree<numtyp, acctyp>
@@ -45,7 +45,7 @@ int BaseThreeT::bytes_per_atom_atomic(const int max_nbors) const {
   #ifdef THREE_CONCURRENT
   b+=ans2->bytes_per_atom();
   #endif
-  return b;     
+  return b;
 }
 
 template <class numtyp, class acctyp>
@@ -62,6 +62,7 @@ int BaseThreeT::init_three(const int nlocal, const int nall,
     gpu_nbor=1;
   else if (device->gpu_mode()==Device<numtyp,acctyp>::GPU_HYB_NEIGH)
     gpu_nbor=2;
+  _gpu_nbor=gpu_nbor;
 
   int _gpu_host=0;
   int host_nlocal=hd_balancer.first_host_count(nlocal,gpu_split,gpu_nbor);
@@ -76,7 +77,7 @@ int BaseThreeT::init_three(const int nlocal, const int nall,
     _nbor_data=&(nbor->dev_nbor);
   if (_threads_per_atom*_threads_per_atom>device->warp_size())
     return -10;
-    
+
   int success=device->init(*ans,false,false,nlocal,host_nlocal,nall,nbor,
                            maxspecial,_gpu_host,max_nbors,cell_size,false,
                            _threads_per_atom);
@@ -93,7 +94,7 @@ int BaseThreeT::init_three(const int nlocal, const int nall,
     return -3;
   ans2->cq(_end_command_queue);
   #endif
-    
+
   _block_pair=device->pair_block_size();
   _block_size=device->block_ellipse();
   compile_kernels(*ucl_device,pair_program,k_two,k_three_center,k_three_end);
@@ -111,7 +112,7 @@ int BaseThreeT::init_three(const int nlocal, const int nall,
   #ifdef THREE_CONCURRENT
   _max_an_bytes+=ans2->gpu_bytes();
   #endif
-  
+
   return 0;
 }
 
@@ -158,7 +159,7 @@ void BaseThreeT::clear_atomic() {
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
 int * BaseThreeT::reset_nbors(const int nall, const int inum, const int nlist,
-                              int *ilist, int *numj, int **firstneigh, 
+                              int *ilist, int *numj, int **firstneigh,
                               bool &success) {
   success=true;
 
@@ -168,7 +169,12 @@ int * BaseThreeT::reset_nbors(const int nall, const int inum, const int nlist,
   if (!success)
     return NULL;
 
-  nbor->get_host3(nall,nlist,ilist,numj,firstneigh,block_size());
+  // originally the requirement that nall == nlist was enforced
+  // to allow direct indexing neighbors of neighbors after re-arrangement
+//  nbor->get_host3(nall,nlist,ilist,numj,firstneigh,block_size());
+
+  // now the requirement is removed, allowing to work within pair hybrid
+  nbor->get_host(nlist,ilist,numj,firstneigh,block_size());
 
   double bytes=ans->gpu_bytes()+nbor->gpu_bytes();
   #ifdef THREE_CONCURRENT
@@ -176,7 +182,7 @@ int * BaseThreeT::reset_nbors(const int nall, const int inum, const int nlist,
   #endif
   if (bytes>_max_an_bytes)
     _max_an_bytes=bytes;
-  
+
   return ilist;
 }
 
@@ -185,11 +191,11 @@ int * BaseThreeT::reset_nbors(const int nall, const int inum, const int nlist,
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
 inline int BaseThreeT::build_nbor_list(const int inum, const int host_inum,
-                                         const int nall, double **host_x,
-                                         int *host_type, double *sublo,
-                                         double *subhi, tagint *tag,
-                                         int **nspecial, tagint **special,
-                                         bool &success) {
+                                       const int nall, double **host_x,
+                                       int *host_type, double *sublo,
+                                       double *subhi, tagint *tag,
+                                       int **nspecial, tagint **special,
+                                       bool &success) {
   success=true;
   resize_atom(inum,nall,success);
   resize_local(nall,host_inum,nbor->max_nbors(),success);
@@ -214,11 +220,11 @@ inline int BaseThreeT::build_nbor_list(const int inum, const int host_inum,
 // Copy nbor list from host if necessary and then calculate forces, virials,..
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
-void BaseThreeT::compute(const int f_ago, const int nlocal, const int nall, 
+void BaseThreeT::compute(const int f_ago, const int inum_full, const int nall,
                          const int nlist, double **host_x, int *host_type,
-                         int *ilist, int *numj, int **firstneigh, 
+                         int *ilist, int *numj, int **firstneigh,
                          const bool eflag, const bool vflag, const bool eatom,
-                         const bool vatom, int &host_start, 
+                         const bool vatom, int &host_start,
                          const double cpu_time, bool &success) {
   acc_timers();
   if (nlist==0) {
@@ -228,9 +234,9 @@ void BaseThreeT::compute(const int f_ago, const int nlocal, const int nall,
     zero_timers();
     return;
   }
-  
+
   int ago=hd_balancer.ago_first(f_ago);
-  int inum=hd_balancer.balance(ago,nlocal,cpu_time);
+  int inum=hd_balancer.balance(ago,inum_full,cpu_time);
   ans->inum(inum);
   #ifdef THREE_CONCURRENT
   ans2->inum(inum);
@@ -270,7 +276,7 @@ template <class numtyp, class acctyp>
 int ** BaseThreeT::compute(const int ago, const int inum_full,
                                  const int nall, double **host_x, int *host_type,
                                  double *sublo, double *subhi, tagint *tag,
-                                 int **nspecial, tagint **special, const bool eflag, 
+                                 int **nspecial, tagint **special, const bool eflag,
                                  const bool vflag, const bool eatom,
                                  const bool vatom, int &host_start,
                                  int **ilist, int **jnum,
@@ -283,7 +289,7 @@ int ** BaseThreeT::compute(const int ago, const int inum_full,
     zero_timers();
     return NULL;
   }
-  
+
   hd_balancer.balance(cpu_time);
   int inum=hd_balancer.get_gpu_count(ago,inum_full);
   ans->inum(inum);
@@ -291,7 +297,7 @@ int ** BaseThreeT::compute(const int ago, const int inum_full,
   ans2->inum(inum);
   #endif
   host_start=inum;
- 
+
   // Build neighbor list on GPU if necessary
   if (ago==0) {
     build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
@@ -321,7 +327,7 @@ int ** BaseThreeT::compute(const int ago, const int inum_full,
   device->add_ans_object(ans2);
   #endif
   hd_balancer.stop_timer();
-  
+
   return nbor->host_jlist.begin()-host_start;
 }
 
@@ -352,7 +358,7 @@ void BaseThreeT::compile_kernels(UCL_Device &dev, const void *pair_str,
   k_three_end.cq(ucl_device->cq(_end_command_queue));
   k_three_end_vatom.cq(ucl_device->cq(_end_command_queue));
   #endif
-  
+
   _compiled=true;
 }
 
diff --git a/lib/gpu/lal_base_three.h b/lib/gpu/lal_base_three.h
index 0af2904..4f27ecd 100644
--- a/lib/gpu/lal_base_three.h
+++ b/lib/gpu/lal_base_three.h
@@ -44,7 +44,7 @@ class BaseThree {
     * \param gpu_split fraction of particles handled by device
     * \param k_two name for the kernel for 2-body force calculation
     * \param k_three name for the kernel for 3-body force calculation
-    * 
+    *
     * Returns:
     * -  0 if successfull
     * - -1 if fix gpu not found
@@ -53,8 +53,8 @@ class BaseThree {
     * - -5 Double precision is not supported on card
     * - -10 if invalid thread_per_atom setting **/
   int init_three(const int nlocal, const int nall, const int max_nbors,
-                 const int maxspecial, const double cell_size, 
-                 const double gpu_split, FILE *screen, 
+                 const int maxspecial, const double cell_size,
+                 const double gpu_split, FILE *screen,
                  const void *pair_program, const char *k_two,
                  const char *k_three_center, const char *k_three_end);
 
@@ -88,7 +88,7 @@ class BaseThree {
     * \note host_inum is 0 if the host is performing neighboring
     * \note nlocal+host_inum=total number local particles
     * \note olist_size=0 **/
-  inline void resize_local(const int inum, const int host_inum, 
+  inline void resize_local(const int inum, const int host_inum,
                            const int max_nbors, bool &success) {
     nbor->resize(inum,host_inum,max_nbors,success);
   }
@@ -133,33 +133,33 @@ class BaseThree {
   /// Build neighbor list on device
   int build_nbor_list(const int inum, const int host_inum,
                        const int nall, double **host_x, int *host_type,
-                       double *sublo, double *subhi, tagint *tag, int **nspecial, 
+                       double *sublo, double *subhi, tagint *tag, int **nspecial,
                        tagint **special, bool &success);
 
   /// Pair loop with host neighboring
-  void compute(const int f_ago, const int inum_full, const int nall, 
+  void compute(const int f_ago, const int inum_full, const int nall,
                const int nlist, double **host_x, int *host_type,
                int *ilist, int *numj, int **firstneigh, const bool eflag,
                const bool vflag, const bool eatom, const bool vatom,
                int &host_start, const double cpu_time, bool &success);
 
   /// Pair loop with device neighboring
-  int * compute(const int ago, const int inum_full, const int nall, 
+  int * compute(const int ago, const int inum_full, const int nall,
                 double **host_x, int *host_type, double *sublo,
                 double *subhi, tagint *tag, int **nspecial,
-                tagint **special, const bool eflag, const bool vflag, 
-                const bool eatom, const bool vatom, int &host_start, 
+                tagint **special, const bool eflag, const bool vflag,
+                const bool eatom, const bool vatom, int &host_start,
                 const double cpu_time, bool &success);
 
   /// Pair loop with device neighboring
   int ** compute(const int ago, const int inum_full,
                  const int nall, double **host_x, int *host_type, double *sublo,
                  double *subhi, tagint *tag, int **nspecial,
-                 tagint **special, const bool eflag, const bool vflag, 
-                 const bool eatom, const bool vatom, int &host_start, 
+                 tagint **special, const bool eflag, const bool vflag,
+                 const bool eatom, const bool vatom, int &host_start,
                  int **ilist, int **numj, const double cpu_time, bool &success);
 
-  // -------------------------- DEVICE DATA ------------------------- 
+  // -------------------------- DEVICE DATA -------------------------
 
   /// Device Properties and Atom and Neighbor storage
   Device<numtyp,acctyp> *device;
@@ -186,7 +186,7 @@ class BaseThree {
   Answer<numtyp,acctyp> *ans;
   #ifdef THREE_CONCURRENT
   Answer<numtyp,acctyp> *ans2;
-  #endif  
+  #endif
 
   // --------------------------- NBOR DATA ----------------------------
 
@@ -205,15 +205,16 @@ class BaseThree {
  protected:
   bool _compiled;
   int _block_pair, _block_size, _threads_per_atom, _end_command_queue;
+  int _gpu_nbor;
   double _max_bytes, _max_an_bytes;
   double _gpu_overhead, _driver_overhead;
   UCL_D_Vec<int> *_nbor_data;
 
-  void compile_kernels(UCL_Device &dev, const void *pair_string, 
+  void compile_kernels(UCL_Device &dev, const void *pair_string,
                        const char *k_two, const char *k_three_center,
                        const char *k_three_end);
 
-  virtual void loop(const bool _eflag, const bool _vflag, 
+  virtual void loop(const bool _eflag, const bool _vflag,
                     const int evatom) = 0;
 };
 
diff --git a/lib/gpu/lal_neighbor.cpp b/lib/gpu/lal_neighbor.cpp
index 074eaa8..0a9933a 100644
--- a/lib/gpu/lal_neighbor.cpp
+++ b/lib/gpu/lal_neighbor.cpp
@@ -10,7 +10,7 @@
     This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
  __________________________________________________________________________
 
-    begin                : 
+    begin                :
     email                : brownw@ornl.gov, penwang@nvidia.com
  ***************************************************************************/
 
@@ -32,13 +32,13 @@ int Neighbor::bytes_per_atom(const int max_nbors) const {
 }
 
 bool Neighbor::init(NeighborShared *shared, const int inum,
-                       const int host_inum, const int max_nbors, 
-                       const int maxspecial, UCL_Device &devi, 
-                       const int gpu_nbor, const int gpu_host, 
+                       const int host_inum, const int max_nbors,
+                       const int maxspecial, UCL_Device &devi,
+                       const int gpu_nbor, const int gpu_host,
                        const bool pre_cut, const int block_cell_2d,
                        const int block_cell_id, const int block_nbor_build,
                        const int threads_per_atom, const int warp_size,
-                       const bool time_device, 
+                       const bool time_device,
                        const std::string compile_flags) {
   clear();
 
@@ -56,10 +56,10 @@ bool Neighbor::init(NeighborShared *shared, const int inum,
     _gpu_host=false;
   else if (gpu_host==1)
     _gpu_host=true;
-  else 
+  else
     // Not yet implemented
     assert(0==1);
-  
+
   if (pre_cut || gpu_nbor==0)
     _alloc_packed=true;
   else
@@ -71,7 +71,7 @@ bool Neighbor::init(NeighborShared *shared, const int inum,
     _packed_permissions=UCL_READ_ONLY;
 
   bool success=true;
-    
+
   // Initialize timers for the selected GPU
   _nbor_time_avail=false;
   time_nbor.init(*dev);
@@ -88,7 +88,7 @@ bool Neighbor::init(NeighborShared *shared, const int inum,
   _max_atoms=static_cast<int>(static_cast<double>(inum)*1.10);
   if (_max_atoms==0)
     _max_atoms=1000;
-    
+
   _max_host=static_cast<int>(static_cast<double>(host_inum)*1.10);
   _max_nbors=(max_nbors/threads_per_atom+1)*threads_per_atom;
 
@@ -102,21 +102,21 @@ bool Neighbor::init(NeighborShared *shared, const int inum,
   alloc(success);
   if (!success)
     return false;
-    
+
   if (_use_packing==false)
     _shared->compile_kernels(devi,gpu_nbor,compile_flags);
 
   return success;
 }
 
-void Neighbor::alloc(bool &success) { 
+void Neighbor::alloc(bool &success) {
   dev_nbor.clear();
   host_acc.clear();
   int nt=_max_atoms+_max_host;
-  if (_use_packing==false || _gpu_nbor>0) 
-    success=success && 
+  if (_use_packing==false || _gpu_nbor>0)
+    success=success &&
             (dev_nbor.alloc((_max_nbors+2)*_max_atoms,*dev)==UCL_SUCCESS);
-  else 
+  else
     success=success && (dev_nbor.alloc(3*_max_atoms,*dev,
                                        UCL_READ_ONLY)==UCL_SUCCESS);
   success=success && (host_acc.alloc(nt*2,*dev,
@@ -127,14 +127,17 @@ void Neighbor::alloc(bool &success) {
     dev_packed.clear();
     success=success && (dev_packed.alloc((_max_nbors+2)*_max_atoms,*dev,
                                          _packed_permissions)==UCL_SUCCESS);
-    _c_bytes+=dev_packed.row_bytes();                                         
-  } 
+    dev_acc.clear();
+    success=success && (dev_acc.alloc(_max_atoms,*dev,
+                                      UCL_READ_WRITE)==UCL_SUCCESS);
+    _c_bytes+=dev_packed.row_bytes()+dev_acc.row_bytes();
+  }
   if (_max_host>0) {
     nbor_host.clear();
     dev_numj_host.clear();
     host_ilist.clear();
     host_jlist.clear();
-    
+
     success=(nbor_host.alloc(_max_nbors*_max_host,*dev,UCL_READ_WRITE,
                              UCL_READ_WRITE)==UCL_SUCCESS) && success;
     success=success && (dev_numj_host.alloc(_max_host,*dev,
@@ -152,7 +155,7 @@ void Neighbor::alloc(bool &success) {
     for (int i=0; i<_max_host; i++) {
       host_jlist[i]=ptr;
       ptr+=_max_nbors;
-    }                                                 
+    }
     _c_bytes+=nbor_host.device.row_bytes()+dev_numj_host.row_bytes();
   } else {
     // Some OpenCL implementations return errors for NULL pointers as args
@@ -176,7 +179,7 @@ void Neighbor::alloc(bool &success) {
 
   _allocated=true;
 }
-  
+
 void Neighbor::clear() {
   _gpu_bytes=0.0;
   _cell_bytes=0.0;
@@ -194,6 +197,7 @@ void Neighbor::clear() {
 
     host_packed.clear();
     host_acc.clear();
+    dev_acc.clear();
     dev_nbor.clear();
     nbor_host.clear();
     dev_packed.clear();
@@ -219,13 +223,13 @@ double Neighbor::host_memory_usage() const {
              host_ilist.row_bytes()+host_jlist.row_bytes();
     else
       return 0;
-  } else 
+  } else
     return host_packed.row_bytes()*host_packed.rows()+host_acc.row_bytes()+
            sizeof(Neighbor);
 }
 
 void Neighbor::get_host(const int inum, int *ilist, int *numj,
-                           int **firstneigh, const int block_size) {  
+                        int **firstneigh, const int block_size) {
   _nbor_time_avail=true;
   time_nbor.start();
 
@@ -242,7 +246,7 @@ void Neighbor::get_host(const int inum, int *ilist, int *numj,
   int dev_count=0;
   int *h_ptr=host_packed.begin();
   _nbor_pitch=inum;
-  
+
   for (int ii=0; ii<inum; ii++) {
     int i=ilist[ii];
     int nj=numj[i];
@@ -250,13 +254,13 @@ void Neighbor::get_host(const int inum, int *ilist, int *numj,
     host_acc[ii+inum]=acc_count;
 
     acc_count+=nj;
-    
+
     int *jlist=firstneigh[i];
     for (int jj=0; jj<nj; jj++) {
       *h_ptr=jlist[jj];
       h_ptr++;
       ij_count++;
-       
+
       if (ij_count==IJ_SIZE) {
         dev_nbor.sync();
         host_offset.view_offset(IJ_SIZE*(copy_count%2),host_packed,IJ_SIZE);
@@ -278,8 +282,17 @@ void Neighbor::get_host(const int inum, int *ilist, int *numj,
   UCL_D_Vec<int> acc_view;
   acc_view.view_offset(inum,dev_nbor,inum*2);
   ucl_copy(acc_view,host_acc,true);
+
+  UCL_H_Vec<int> host_view;
+  host_view.alloc(_max_atoms,*dev,UCL_READ_WRITE);
+  for (int ii=0; ii<inum; ii++) {
+    int i=ilist[ii];
+    host_view[i] = ii;
+  }
+  ucl_copy(dev_acc,host_view,true);
+
   time_nbor.stop();
-  
+
   if (_use_packing==false) {
     time_kernel.start();
     int GX=static_cast<int>(ceil(static_cast<double>(inum)*_threads_per_atom/
@@ -294,7 +307,7 @@ void Neighbor::get_host(const int inum, int *ilist, int *numj,
 // inum=nlocal is forced to be true to allow direct indexing of neighbors of
 // neighbors
 void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
-                         int **firstneigh, const int block_size) {  
+                         int **firstneigh, const int block_size) {
   _nbor_time_avail=true;
   time_nbor.start();
 
@@ -311,7 +324,7 @@ void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
   int dev_count=0;
   int *h_ptr=host_packed.begin();
   _nbor_pitch=inum;
-  
+
   if (nlist!=inum)
     host_acc.zero(inum);
 
@@ -322,7 +335,7 @@ void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
     host_acc[i+inum]=acc_count;
     acc_count+=nj;
   }
-  
+
   for (int i=0; i<inum; i++) {
     int nj=host_acc[i];
     int *jlist=firstneigh[i];
@@ -330,7 +343,7 @@ void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
       *h_ptr=jlist[jj];
       h_ptr++;
       ij_count++;
-       
+
       if (ij_count==IJ_SIZE) {
         dev_nbor.sync();
         host_offset.view_offset(IJ_SIZE*(copy_count%2),host_packed,IJ_SIZE);
@@ -353,7 +366,7 @@ void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
   acc_view.view_offset(inum,dev_nbor,inum*2);
   ucl_copy(acc_view,host_acc,true);
   time_nbor.stop();
-  
+
   if (_use_packing==false) {
     time_kernel.start();
     int GX=static_cast<int>(ceil(static_cast<double>(inum)*_threads_per_atom/
@@ -366,7 +379,7 @@ void Neighbor::get_host3(const int inum, const int nlist, int *ilist, int *numj,
 
 template <class numtyp, class acctyp>
 void Neighbor::resize_max_neighbors(const int maxn, bool &success) {
-  if (maxn>_max_nbors) {  
+  if (maxn>_max_nbors) {
     int mn=static_cast<int>(static_cast<double>(maxn)*1.10);
     mn=(mn/_threads_per_atom+1)*_threads_per_atom;
     success=success && (dev_nbor.resize((mn+1)*_max_atoms)==UCL_SUCCESS);
@@ -377,7 +390,7 @@ void Neighbor::resize_max_neighbors(const int maxn, bool &success) {
       for (int i=0; i<_max_host; i++) {
         host_jlist[i]=ptr;
         ptr+=mn;
-      }                                                 
+      }
       _gpu_bytes+=nbor_host.row_bytes();
     } else {
       nbor_host.device.view(dev_nbor);
@@ -393,8 +406,8 @@ void Neighbor::resize_max_neighbors(const int maxn, bool &success) {
 
 template <class numtyp, class acctyp>
 void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
-                               const int nall, Atom<numtyp,acctyp> &atom, 
-                               double *sublo, double *subhi, tagint *tag, 
+                               const int nall, Atom<numtyp,acctyp> &atom,
+                               double *sublo, double *subhi, tagint *tag,
                                int **nspecial, tagint **special, bool &success,
                                int &mn) {
   _nbor_time_avail=true;
@@ -409,7 +422,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
   ncell_3d = ncellx * ncelly * ncellz;
   if (ncell_3d+1>_ncells) {
     cell_counts.clear();
-  
+
     if (_gpu_nbor==2) {
       if (_ncells>0)
         delete [] cell_iter;
@@ -419,7 +432,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
       cell_counts.device.clear();
       cell_counts.device.alloc(ncell_3d+1,dev_nbor);
     }
-    
+
     _ncells=ncell_3d+1;
     _cell_bytes=cell_counts.device.row_bytes();
   }
@@ -445,17 +458,17 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
     const int g2x=static_cast<int>(ceil(static_cast<double>(_maxspecial)/b2x));
     const int g2y=static_cast<int>(ceil(static_cast<double>(nt)/b2y));
     _shared->k_transpose.set_size(g2x,g2y,b2x,b2y);
-    _shared->k_transpose.run(&dev_special,&dev_special_t,&_maxspecial,&nt);        
+    _shared->k_transpose.run(&dev_special,&dev_special_t,&_maxspecial,&nt);
     time_transpose.stop();
   }
-  
+
   // If binning on CPU, do this now
   if (_gpu_nbor==2) {
     double stime = MPI_Wtime();
     int *cell_id=atom.host_cell_id.begin();
     int *particle_id=atom.host_particle_id.begin();
-    
-    // Build cell list on CPU                               
+
+    // Build cell list on CPU
     cell_counts.host.zero();
     double i_cell_size=1.0/_cell_size;
 
@@ -475,12 +488,12 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
       int iz = static_cast<int>(pz*i_cell_size+1);
       iz = std::max(iz,_cells_in_cutoff);
       iz = std::min(iz,ncellz-offset_hi);
-    
+
       int id = ix+iy*ncellx+iz*ncellx*ncelly;
       cell_id[i] = id;
       cell_counts[id+1]++;
     }
-    
+
     for (int i=nt; i<nall; i++) {
       double px, py, pz;
       px=x[i][0]-sublo[0];
@@ -496,12 +509,12 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
       int iz = static_cast<int>(pz*i_cell_size+1);
       iz = std::max(iz,0);
       iz = std::min(iz,ncellz-1);
-    
+
       int id = ix+iy*ncellx+iz*ncellx*ncelly;
       cell_id[i] = id;
       cell_counts[id+1]++;
     }
-    
+
     mn=0;
     for (int i=0; i<_ncells; i++)
       mn=std::max(mn,cell_counts[i]);
@@ -531,7 +544,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
     ucl_copy(atom.dev_particle_id,atom.host_particle_id,true);
     time_hybrid2.stop();
     _bin_time+=MPI_Wtime()-stime;
-  }        
+  }
 
   time_kernel.start();
 
@@ -547,7 +560,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
     const numtyp sublo1=static_cast<numtyp>(sublo[1]);
     const numtyp sublo2=static_cast<numtyp>(sublo[2]);
     _shared->k_cell_id.set_size(GX,neigh_block);
-    _shared->k_cell_id.run(&atom.x, &atom.dev_cell_id, 
+    _shared->k_cell_id.run(&atom.x, &atom.dev_cell_id,
                            &atom.dev_particle_id, &sublo0, &sublo1,
                            &sublo2, &i_cell_size, &ncellx, &ncelly, &ncellz,
                            &nt, &nall, &_cells_in_cutoff);
@@ -556,10 +569,10 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
 
     /* calculate cell count */
     _shared->k_cell_counts.set_size(GX,neigh_block);
-    _shared->k_cell_counts.run(&atom.dev_cell_id, &cell_counts, &nall, 
+    _shared->k_cell_counts.run(&atom.dev_cell_id, &cell_counts, &nall,
                                &ncell_3d);
-  } 
-  
+  }
+
   /* build the neighbor list */
   const int cell_block=_block_nbor_build;
   _shared->k_build_nbor.set_size(ncellx-ghost_cells,(ncelly-ghost_cells)*
@@ -579,7 +592,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
     host_offset.view_offset(inum,host_acc,nt-inum);
     ucl_copy(host_offset,dev_numj_host,nt-inum,true);
   }
-  
+
   if (_gpu_nbor!=2) {
     host_acc.sync();
     mn=host_acc[0];
@@ -587,7 +600,7 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
       mn=std::max(mn,host_acc[i]);
     set_nbor_block_size(mn);
 
-    if (mn>_max_nbors) {  
+    if (mn>_max_nbors) {
       resize_max_neighbors<numtyp,acctyp>(mn,success);
       if (!success)
         return;
@@ -599,13 +612,13 @@ void Neighbor::build_nbor_list(double **x, const int inum, const int host_inum,
       return;
     }
   }
-  
+
   if (_maxspecial>0) {
     const int GX2=static_cast<int>(ceil(static_cast<double>
                                           (nt*_threads_per_atom)/cell_block));
     _shared->k_special.set_size(GX2,cell_block);
     _shared->k_special.run(&dev_nbor, &nbor_host, &dev_numj_host,
-                           &atom.dev_tag, &dev_nspecial, &dev_special, 
+                           &atom.dev_tag, &dev_nspecial, &dev_special,
                            &inum, &nt, &_max_nbors, &_threads_per_atom);
   }
   time_kernel.stop();
diff --git a/lib/gpu/lal_neighbor.h b/lib/gpu/lal_neighbor.h
index 7653291..0516883 100644
--- a/lib/gpu/lal_neighbor.h
+++ b/lib/gpu/lal_neighbor.h
@@ -10,7 +10,7 @@
     This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
  __________________________________________________________________________
 
-    begin                : 
+    begin                :
     email                : brownw@ornl.gov, penwang@nvidia.com
  ***************************************************************************/
 
@@ -28,12 +28,12 @@ class Neighbor {
  public:
   Neighbor() : _allocated(false), _use_packing(false), _ncells(0) {}
   ~Neighbor() { clear(); }
- 
+
   /// Determine whether neighbor unpacking should be used
-  /** If false, twice as much memory is reserved to allow unpacking neighbors by 
+  /** If false, twice as much memory is reserved to allow unpacking neighbors by
     * atom for coalesced access. **/
   void packing(const bool use_packing) { _use_packing=use_packing; }
-  
+
   /// Clear any old data and setup for new LAMMPS run
   /** \param inum Initial number of particles whose neighbors stored on device
     * \param host_inum Initial number of particles whose nbors copied to host
@@ -45,20 +45,20 @@ class Neighbor {
     *                 1 if gpu_nbor is true, and host needs a half nbor list,
     *                 2 if gpu_nbor is true, and host needs a full nbor list
     * \param pre_cut True if cutoff test will be performed in separate kernel
-    *                than the force kernel 
+    *                than the force kernel
     * \param threads_per_atom Number of threads used per atom for force
-    *                         calculation 
+    *                         calculation
     * \param compile_flags Flags for JIT compiling **/
   bool init(NeighborShared *shared, const int inum, const int host_inum,
             const int max_nbors, const int maxspecial, UCL_Device &dev,
             const int gpu_nbor, const int gpu_host, const bool pre_cut,
-            const int block_cell_2d, const int block_cell_id, 
+            const int block_cell_2d, const int block_cell_id,
             const int block_nbor_build, const int threads_per_atom,
-            const int warp_size, const bool time_device, 
+            const int warp_size, const bool time_device,
             const std::string compile_flags);
 
   /// Set the size of the cutoff+skin
-  inline void cell_size(const double size, const double cutoff) { 
+  inline void cell_size(const double size, const double cutoff) {
     _cell_size=size;
     _cutoff=cutoff;
     if (cutoff>size)
@@ -66,7 +66,7 @@ class Neighbor {
     else
       _cells_in_cutoff=1;
   }
-  
+
   /// Get the size of the cutoff+skin
   inline double cell_size() const { return _cell_size; }
 
@@ -88,7 +88,7 @@ class Neighbor {
     * \param host_inum Number of particles whose nbors will be copied to host
     * \param max_nbor Current max number of neighbors for a particle
     * \param success False if insufficient memory **/
-  inline void resize(const int inum, const int host_inum, const int max_nbor, 
+  inline void resize(const int inum, const int host_inum, const int max_nbor,
                      bool &success) {
     if (inum>_max_atoms || max_nbor>_max_nbors || host_inum>_max_host) {
       _max_atoms=static_cast<int>(static_cast<double>(inum)*1.10);
@@ -124,40 +124,40 @@ class Neighbor {
 
   /// Free all memory on host and device
   void clear();
- 
+
   /// Bytes per atom used on device
   int bytes_per_atom(const int max_nbors) const;
-  
+
   /// Total host memory used by class
   double host_memory_usage() const;
-  
+
   /// Returns the type of neighboring:
   /** - 0 if neighboring will be performed on host
     * - 1 if neighboring will be performed on device
     * - 2 if binning on host and neighboring on device **/
   inline int gpu_nbor() const { return _gpu_nbor; }
-  
+
   /// Make a copy of unpacked nbor lists in the packed storage area (for gb)
-  inline void copy_unpacked(const int inum, const int maxj) 
+  inline void copy_unpacked(const int inum, const int maxj)
     { ucl_copy(dev_packed,dev_nbor,inum*(maxj+2),true); }
 
-  /// Copy neighbor list from host (first time or from a rebuild)  
-  void get_host(const int inum, int *ilist, int *numj, 
+  /// Copy neighbor list from host (first time or from a rebuild)
+  void get_host(const int inum, int *ilist, int *numj,
                 int **firstneigh, const int block_size);
-  
-  /// Copy neighbor list from host for 3-body (first time or from a rebuild)  
-  void get_host3(const int inum, const int nlist, int *ilist, int *numj, 
+
+  /// Copy neighbor list from host for 3-body (first time or from a rebuild)
+  void get_host3(const int inum, const int nlist, int *ilist, int *numj,
                  int **firstneigh, const int block_size);
-  
+
   /// Return the stride in elements for each nbor row
   inline int nbor_pitch() const { return _nbor_pitch; }
-  
+
   /// Return the maximum number of atoms that can currently be stored
   inline int max_atoms() const { return _max_atoms; }
 
   /// Return the maximum number of nbors for a particle based on current alloc
   inline int max_nbors() const { return _max_nbors; }
-  
+
   /// Return the time spent binning on the CPU for hybrid neighbor builds
   inline double bin_time() const { return _bin_time; }
 
@@ -171,9 +171,9 @@ class Neighbor {
 
   /// Build nbor list on the device
   template <class numtyp, class acctyp>
-  void build_nbor_list(double **x, const int inum, const int host_inum, 
+  void build_nbor_list(double **x, const int inum, const int host_inum,
                        const int nall, Atom<numtyp,acctyp> &atom, double *sublo,
-                       double *subhi, tagint *tag, int **nspecial, tagint **special, 
+                       double *subhi, tagint *tag, int **nspecial, tagint **special,
                        bool &success, int &max_nbors);
 
   /// Return the number of bytes used on device
@@ -184,7 +184,7 @@ class Neighbor {
 
     return res;
   }
-  
+
   // ------------------------------- Data -------------------------------
 
   /// Device neighbor matrix
@@ -199,6 +199,8 @@ class Neighbor {
   UCL_H_Vec<int> host_packed;
   /// Host storage for nbor counts (row 1) & accumulated neighbor counts (row2)
   UCL_H_Vec<int> host_acc;
+  /// Device storage for accessing atom indices from the neighbor list (3-body)
+  UCL_D_Vec<int> dev_acc;
 
   // ----------------- Data for GPU Neighbor Calculation ---------------
 
@@ -219,7 +221,7 @@ class Neighbor {
 
   /// Device timers
   UCL_Timer time_nbor, time_kernel, time_hybrid1, time_hybrid2, time_transpose;
-  
+
  private:
   NeighborShared *_shared;
   UCL_Device *dev;
@@ -231,14 +233,14 @@ class Neighbor {
 
   double _gpu_bytes, _c_bytes, _cell_bytes;
   void alloc(bool &success);
-  
+
   int _block_cell_2d, _block_cell_id, _max_block_nbor_build, _block_nbor_build;
   int _ncells, _threads_per_atom, _total_atoms;
   int _cells_in_cutoff;
 
   template <class numtyp, class acctyp>
   inline void resize_max_neighbors(const int maxn, bool &success);
-  
+
   int _warp_size;
   inline void set_nbor_block_size(const int mn) {
     int desired=mn/(2*_warp_size);
diff --git a/lib/gpu/lal_sw.cpp b/lib/gpu/lal_sw.cpp
index 1f68616..3492d70 100644
--- a/lib/gpu/lal_sw.cpp
+++ b/lib/gpu/lal_sw.cpp
@@ -33,10 +33,10 @@ SWT::SW() : BaseThree<numtyp,acctyp>(), _allocated(false) {
 }
 
 template <class numtyp, class acctyp>
-SWT::~SW() { 
+SWT::~SW() {
   clear();
 }
- 
+
 template <class numtyp, class acctyp>
 int SWT::bytes_per_atom(const int max_nbors) const {
   return this->bytes_per_atom_atomic(max_nbors);
@@ -45,7 +45,7 @@ int SWT::bytes_per_atom(const int max_nbors) const {
 template <class numtyp, class acctyp>
 int SWT::init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
            const double cell_size, const double gpu_split, FILE *_screen,
-           int* host_map, const int nelements, int*** host_elem2param, const int nparams, 
+           int* host_map, const int nelements, int*** host_elem2param, const int nparams,
            const double* epsilon, const double* sigma,
            const double* lambda, const double* gamma,
            const double* costheta, const double* biga,
@@ -76,41 +76,41 @@ int SWT::init(const int ntypes, const int nlocal, const int nall, const int max_
                              UCL_WRITE_ONLY);
 
   for (int i=0; i<nparams; i++) {
-    dview[i].x=(numtyp)0; 
+    dview[i].x=(numtyp)0;
     dview[i].y=(numtyp)0;
-    dview[i].z=(numtyp)0; 
+    dview[i].z=(numtyp)0;
     dview[i].w=(numtyp)0;
   }
 
   // pack coefficients into arrays
   sw1.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
-  
+
   for (int i=0; i<nparams; i++) {
     dview[i].x=static_cast<numtyp>(epsilon[i]);
     dview[i].y=static_cast<numtyp>(sigma[i]);
     dview[i].z=static_cast<numtyp>(lambda[i]);
     dview[i].w=static_cast<numtyp>(gamma[i]);
   }
-  
+
   ucl_copy(sw1,dview,false);
   sw1_tex.get_texture(*(this->pair_program),"sw1_tex");
   sw1_tex.bind_float(sw1,4);
 
   sw2.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
-  
+
   for (int i=0; i<nparams; i++) {
     dview[i].x=static_cast<numtyp>(biga[i]);
     dview[i].y=static_cast<numtyp>(bigb[i]);
     dview[i].z=static_cast<numtyp>(powerp[i]);
     dview[i].w=static_cast<numtyp>(powerq[i]);
   }
-  
+
   ucl_copy(sw2,dview,false);
   sw2_tex.get_texture(*(this->pair_program),"sw2_tex");
   sw2_tex.bind_float(sw2,4);
 
   sw3.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
-  
+
   for (int i=0; i<nparams; i++) {
     double sw_cut = cut[i];
     double sw_cutsq = cutsq[i];
@@ -121,7 +121,7 @@ int SWT::init(const int ntypes, const int nlocal, const int nall, const int max_
     dview[i].z=static_cast<numtyp>(costheta[i]);
     dview[i].w=(numtyp)0;
   }
-  
+
   ucl_copy(sw3,dview,false);
   sw3_tex.get_texture(*(this->pair_program),"sw3_tex");
   sw3_tex.bind_float(sw3,4);
@@ -192,31 +192,32 @@ void SWT::loop(const bool _eflag, const bool _vflag, const int evatom) {
     vflag=1;
   else
     vflag=0;
-  
+
   int GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
                                (BX/this->_threads_per_atom)));
 
-  // this->_nbor_data == nbor->dev_packed for gpu_nbor == 0 and tpa == 1 
-  // this->_nbor_data == nbor->dev_nbor for gpu_nbor == 1
+  // this->_nbor_data == nbor->dev_packed for gpu_nbor == 0 and tpa > 1
+  // this->_nbor_data == nbor->dev_nbor for gpu_nbor == 1 or tpa == 1
   int ainum=this->ans->inum();
   int nbor_pitch=this->nbor->nbor_pitch();
   this->time_pair.start();
+
   this->k_pair.set_size(GX,BX);
-  this->k_pair.run(&this->atom->x, &sw1, &sw2, &sw3, 
+  this->k_pair.run(&this->atom->x, &sw1, &sw2, &sw3,
                    &map, &elem2param, &_nelements,
                    &this->nbor->dev_nbor, &this->_nbor_data->begin(),
-                   &this->ans->force, &this->ans->engv, 
-                   &eflag, &vflag, &ainum, &nbor_pitch, 
+                   &this->ans->force, &this->ans->engv,
+                   &eflag, &vflag, &ainum, &nbor_pitch,
                    &this->_threads_per_atom);
 
   BX=this->block_size();
   GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
-                           (BX/(KTHREADS*JTHREADS)))); 
+                           (BX/(KTHREADS*JTHREADS))));
   this->k_three_center.set_size(GX,BX);
-  this->k_three_center.run(&this->atom->x, &sw1, &sw2, &sw3, 
+  this->k_three_center.run(&this->atom->x, &sw1, &sw2, &sw3,
                            &map, &elem2param, &_nelements,
-                           &this->nbor->dev_nbor, &this->_nbor_data->begin(), 
-                           &this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum, 
+                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                           &this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum,
                            &nbor_pitch, &this->_threads_per_atom, &evatom);
 
   Answer<numtyp,acctyp> *end_ans;
@@ -227,21 +228,24 @@ void SWT::loop(const bool _eflag, const bool _vflag, const int evatom) {
   #endif
   if (evatom!=0) {
     this->k_three_end_vatom.set_size(GX,BX);
-    this->k_three_end_vatom.run(&this->atom->x, &sw1, &sw2, &sw3, 
-                          &map, &elem2param, &_nelements, 
-                          &this->nbor->dev_nbor, &this->_nbor_data->begin(), 
+    this->k_three_end_vatom.run(&this->atom->x, &sw1, &sw2, &sw3,
+                          &map, &elem2param, &_nelements,
+                          &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
 
   } else {
     this->k_three_end.set_size(GX,BX);
-    this->k_three_end.run(&this->atom->x, &sw1, &sw2, &sw3, 
-                          &map, &elem2param, &_nelements, 
-                          &this->nbor->dev_nbor, &this->_nbor_data->begin(), 
-                          &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum, 
-                          &nbor_pitch, &this->_threads_per_atom);
+    this->k_three_end.run(&this->atom->x, &sw1, &sw2, &sw3,
+                          &map, &elem2param, &_nelements,
+                          &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
+                          &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
 
   }
+
   this->time_pair.stop();
 }
 
diff --git a/lib/gpu/lal_sw.cu b/lib/gpu/lal_sw.cu
index 1e358fb..63019d1 100644
--- a/lib/gpu/lal_sw.cu
+++ b/lib/gpu/lal_sw.cu
@@ -138,16 +138,16 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
                    const __global int *restrict map,
                    const __global int *restrict elem2param,
                    const int nelements,
-                   const __global int * dev_nbor, 
-                   const __global int * dev_packed, 
-                   __global acctyp4 *restrict ans, 
-                   __global acctyp *restrict engv, 
-                   const int eflag, const int vflag, const int inum, 
+                   const __global int * dev_nbor,
+                   const __global int * dev_packed,
+                   __global acctyp4 *restrict ans,
+                   __global acctyp *restrict engv,
+                   const int eflag, const int vflag, const int inum,
                    const int nbor_pitch, const int t_per_atom) {
   __local int n_stride;
   int tid, ii, offset;
   atom_info(t_per_atom,ii,tid,offset);
-  
+
   acctyp energy=(acctyp)0;
   acctyp4 f;
   f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
@@ -166,9 +166,9 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
     numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
     int itype=ix.w;
     itype=map[itype];
-    
+
     for ( ; nbor<nbor_end; nbor+=n_stride) {
-  
+
       int j=dev_packed[nbor];
       j &= NEIGHMASK;
 
@@ -183,7 +183,7 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
       numtyp dely = ix.y-jx.y;
       numtyp delz = ix.z-jx.z;
       numtyp rsq = delx*delx+dely*dely+delz*delz;
-        
+
       if (rsq<sw3[ijparam].y) { // sw_cutsq = sw3[ijparam].y
         numtyp4 sw1_ijparam; fetch4(sw1_ijparam,ijparam,sw1_tex);
         numtyp sw_epsilon=sw1_ijparam.x;
@@ -218,13 +218,13 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
         numtyp force = (pre_sw_c1*rp-pre_sw_c2*rq +
                        (pre_sw_c3*rp-pre_sw_c4*rq) * rainv)*
                        expsrainv*ucl_recip(rsq);
-      
+
         f.x+=delx*force;
         f.y+=dely*force;
         f.z+=delz*force;
 
-        if (eflag>0) 
-          energy+=(pre_sw_c5*rp - pre_sw_c6*rq) * expsrainv; 
+        if (eflag>0)
+          energy+=(pre_sw_c5*rp - pre_sw_c6*rq) * expsrainv;
 
         if (vflag>0) {
           virial[0] += delx*delx*force;
@@ -329,19 +329,19 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
   fjz = delr1z*(frad1+csfac1)-delr2z*facang12;                               \
 }
 
-__kernel void k_sw_three_center(const __global numtyp4 *restrict x_, 
+__kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
                                 const __global numtyp4 *restrict sw1,
                                 const __global numtyp4 *restrict sw2,
                                 const __global numtyp4 *restrict sw3,
                                 const __global int *restrict map,
                                 const __global int *restrict elem2param,
                                 const int nelements,
-                                const __global int * dev_nbor, 
-                                const __global int * dev_packed, 
-                                __global acctyp4 *restrict ans, 
-                                __global acctyp *restrict engv, 
-                                const int eflag, const int vflag, 
-                                const int inum,  const int nbor_pitch, 
+                                const __global int * dev_nbor,
+                                const __global int * dev_packed,
+                                __global acctyp4 *restrict ans,
+                                __global acctyp *restrict engv,
+                                const int eflag, const int vflag,
+                                const int inum,  const int nbor_pitch,
                                 const int t_per_atom, const int evatom) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
@@ -351,7 +351,7 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
 
   int tid, ii, offset;
   atom_info(tpa_sq,ii,tid,offset);
-  
+
   acctyp energy=(acctyp)0;
   acctyp4 f;
   f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
@@ -360,7 +360,7 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
     virial[i]=(acctyp)0;
 
   __syncthreads();
-  
+
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end;
 
@@ -370,11 +370,11 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
     int offset_k=tid & (t_per_atom-1);
 
     numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
-    int itype=ix.w; 
+    int itype=ix.w;
     itype=map[itype];
 
     for ( ; nbor_j<nbor_end; nbor_j+=n_stride) {
-  
+
       int j=dev_packed[nbor_j];
       j &= NEIGHMASK;
 
@@ -439,7 +439,7 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
           f.x -= fjx + fkx;
           f.y -= fjy + fky;
           f.z -= fjz + fkz;
-        }  
+        }
       }
     } // for nbor
 
@@ -458,20 +458,21 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
   } // if ii
 }
 
-__kernel void k_sw_three_end(const __global numtyp4 *restrict x_, 
+__kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
                              const __global numtyp4 *restrict sw1,
                              const __global numtyp4 *restrict sw2,
                              const __global numtyp4 *restrict sw3,
                              const __global int *restrict map,
                              const __global int *restrict elem2param,
                              const int nelements,
-                             const __global int * dev_nbor, 
-                             const __global int * dev_packed, 
-                             __global acctyp4 *restrict ans, 
-                             __global acctyp *restrict engv, 
-                             const int eflag, const int vflag, 
-                             const int inum,  const int nbor_pitch, 
-                             const int t_per_atom) {
+                             const __global int * dev_nbor,
+                             const __global int * dev_packed,
+                             const __global int * dev_acc,
+                             __global acctyp4 *restrict ans,
+                             __global acctyp *restrict engv,
+                             const int eflag, const int vflag,
+                             const int inum,  const int nbor_pitch,
+                             const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp sw_epsilon, sw_sigma, sw_lambda, sw_gamma;
@@ -480,7 +481,7 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
 
   int tid, ii, offset;
   atom_info(tpa_sq,ii,tid,offset);
-  
+
   acctyp energy=(acctyp)0;
   acctyp4 f;
   f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
@@ -489,7 +490,7 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
     virial[i]=(acctyp)0;
 
   __syncthreads();
-  
+
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
 
@@ -527,13 +528,17 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
       sw_sigma_gamma_ij=sw1_ijparam.y*sw1_ijparam.w; //sw_sigma*sw_gamma;
       sw_cut_ij=sw3_ijparam.x;
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k]; 
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
@@ -605,13 +610,14 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
                              const __global int *restrict map,
                              const __global int *restrict elem2param,
                              const int nelements,
-                             const __global int * dev_nbor, 
-                             const __global int * dev_packed, 
-                             __global acctyp4 *restrict ans, 
-                             __global acctyp *restrict engv, 
-                             const int eflag, const int vflag, 
-                             const int inum,  const int nbor_pitch, 
-                             const int t_per_atom) {
+                             const __global int * dev_nbor,
+                             const __global int * dev_packed,
+                             const __global int * dev_acc,
+                             __global acctyp4 *restrict ans,
+                             __global acctyp *restrict engv,
+                             const int eflag, const int vflag,
+                             const int inum,  const int nbor_pitch,
+                             const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp sw_epsilon, sw_sigma, sw_lambda, sw_gamma;
@@ -620,7 +626,7 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
 
   int tid, ii, offset;
   atom_info(tpa_sq,ii,tid,offset);
-  
+
   acctyp energy=(acctyp)0;
   acctyp4 f;
   f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
@@ -629,7 +635,7 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
     virial[i]=(acctyp)0;
 
   __syncthreads();
-  
+
   if (ii<inum) {
     int i, numj, nbor_j, nbor_end, k_end;
 
@@ -666,14 +672,18 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
       sw_gamma=sw1_ijparam.w;
       sw_sigma_gamma_ij=sw1_ijparam.y*sw1_ijparam.w; //sw_sigma*sw_gamma;
       sw_cut_ij=sw3_ijparam.x;
-        
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k]; 
+
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
diff --git a/lib/gpu/lal_sw.h b/lib/gpu/lal_sw.h
index 66b36a9..3546f02 100644
--- a/lib/gpu/lal_sw.h
+++ b/lib/gpu/lal_sw.h
@@ -24,28 +24,28 @@ template <class numtyp, class acctyp>
 class SW : public BaseThree<numtyp, acctyp> {
  public:
   SW();
-  ~SW(); 
+  ~SW();
 
   /// Clear any previous data and set up for a new LAMMPS run
   /** \param max_nbors initial number of rows in the neighbor matrix
     * \param cell_size cutoff + skin
     * \param gpu_split fraction of particles handled by device
-    * 
+    *
     * Returns:
     * -  0 if successfull
     * - -1 if fix gpu not found
     * - -3 if there is an out of memory error
     * - -4 if the GPU library was not compiled for GPU
     * - -5 Double precision is not supported on card **/
-  int init(const int ntypes, const int nlocal, const int nall, const int max_nbors, 
+  int init(const int ntypes, const int nlocal, const int nall, const int max_nbors,
            const double cell_size, const double gpu_split, FILE *screen,
-           int* host_map, const int nelements, int*** host_elem2param, const int nparams, 
+           int* host_map, const int nelements, int*** host_elem2param, const int nparams,
            const double* epsilon, const double* sigma,
            const double* lambda, const double* gamma,
            const double* costheta, const double* biga,
            const double* bigb, const double* powerp,
            const double* powerq, const double* cut, const double* cutsq);
-           
+
   /// Clear all host and device data
   /** \note This is called at the beginning of the init() routine **/
   void clear();
@@ -61,7 +61,7 @@ class SW : public BaseThree<numtyp, acctyp> {
   /// If atom type constants fit in shared memory, use fast kernels
   bool shared_types;
 
-  /// Number of atom types   
+  /// Number of atom types
   int _lj_types;
 
   /// sw1.x = epsilon, sw1.y = sigma, sw1.z = lambda, sw1.w = gamma
diff --git a/lib/gpu/lal_tersoff.cpp b/lib/gpu/lal_tersoff.cpp
index bc89c53..8805a6f 100644
--- a/lib/gpu/lal_tersoff.cpp
+++ b/lib/gpu/lal_tersoff.cpp
@@ -437,16 +437,18 @@ void TersoffT::loop(const bool _eflag, const bool _vflag, const int evatom) {
     this->k_three_end_vatom.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
 
   } else {
     this->k_three_end.set_size(GX,BX);
     this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
   }
 
   this->time_pair.stop();
diff --git a/lib/gpu/lal_tersoff.cu b/lib/gpu/lal_tersoff.cu
index e98a454..4a71518 100644
--- a/lib/gpu/lal_tersoff.cu
+++ b/lib/gpu/lal_tersoff.cu
@@ -597,11 +597,12 @@ __kernel void k_tersoff_three_end(const __global numtyp4 *restrict x_,
                                   const __global acctyp4 *restrict zetaij,
                                   const __global int * dev_nbor,
                                   const __global int * dev_packed,
+                                  const __global int * dev_acc,
                                   __global acctyp4 *restrict ans,
                                   __global acctyp *restrict engv,
                                   const int eflag, const int vflag,
                                   const int inum,  const int nbor_pitch,
-                                  const int t_per_atom) {
+                                  const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;
@@ -666,13 +667,17 @@ __kernel void k_tersoff_three_end(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
@@ -818,11 +823,12 @@ __kernel void k_tersoff_three_end_vatom(const __global numtyp4 *restrict x_,
                                         const __global acctyp4 *restrict zetaij,
                                         const __global int * dev_nbor,
                                         const __global int * dev_packed,
+                                        const __global int * dev_acc,
                                         __global acctyp4 *restrict ans,
                                         __global acctyp *restrict engv,
                                         const int eflag, const int vflag,
                                         const int inum,  const int nbor_pitch,
-                                        const int t_per_atom) {
+                                        const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;
@@ -887,13 +893,17 @@ __kernel void k_tersoff_three_end_vatom(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
diff --git a/lib/gpu/lal_tersoff_mod.cpp b/lib/gpu/lal_tersoff_mod.cpp
index bfcc9c3..04baa4e 100644
--- a/lib/gpu/lal_tersoff_mod.cpp
+++ b/lib/gpu/lal_tersoff_mod.cpp
@@ -437,16 +437,18 @@ void TersoffMT::loop(const bool _eflag, const bool _vflag, const int evatom) {
     this->k_three_end_vatom.run(&this->atom->x, &ts1, &ts2, &ts4, &ts5, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
 
   } else {
     this->k_three_end.set_size(GX,BX);
     this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts4, &ts5, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
   }
 
   this->time_pair.stop();
diff --git a/lib/gpu/lal_tersoff_mod.cu b/lib/gpu/lal_tersoff_mod.cu
index ba4ad32..13c8147 100644
--- a/lib/gpu/lal_tersoff_mod.cu
+++ b/lib/gpu/lal_tersoff_mod.cu
@@ -605,11 +605,12 @@ __kernel void k_tersoff_mod_three_end(const __global numtyp4 *restrict x_,
                                   const __global acctyp4 *restrict zetaij,
                                   const __global int * dev_nbor,
                                   const __global int * dev_packed,
+                                  const __global int * dev_acc,
                                   __global acctyp4 *restrict ans,
                                   __global acctyp *restrict engv,
                                   const int eflag, const int vflag,
                                   const int inum,  const int nbor_pitch,
-                                  const int t_per_atom) {
+                                  const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c1, c2, c3, c4, c5, h;
@@ -676,13 +677,17 @@ __kernel void k_tersoff_mod_three_end(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
@@ -835,11 +840,12 @@ __kernel void k_tersoff_mod_three_end_vatom(const __global numtyp4 *restrict x_,
                                         const __global acctyp4 *restrict zetaij,
                                         const __global int * dev_nbor,
                                         const __global int * dev_packed,
+                                        const __global int * dev_acc,
                                         __global acctyp4 *restrict ans,
                                         __global acctyp *restrict engv,
                                         const int eflag, const int vflag,
                                         const int inum,  const int nbor_pitch,
-                                        const int t_per_atom) {
+                                        const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c1, c2, c3, c4, c5, h;
@@ -906,13 +912,17 @@ __kernel void k_tersoff_mod_three_end_vatom(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
diff --git a/lib/gpu/lal_tersoff_zbl.cpp b/lib/gpu/lal_tersoff_zbl.cpp
index 57688f5..9193d9a 100644
--- a/lib/gpu/lal_tersoff_zbl.cpp
+++ b/lib/gpu/lal_tersoff_zbl.cpp
@@ -463,16 +463,18 @@ void TersoffZT::loop(const bool _eflag, const bool _vflag, const int evatom) {
     this->k_three_end_vatom.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
 
   } else {
     this->k_three_end.set_size(GX,BX);
     this->k_three_end.run(&this->atom->x, &ts1, &ts2, &ts4, &cutsq,
                           &map, &elem2param, &_nelements, &_nparams, &_zetaij,
                           &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+                          &this->nbor->dev_acc,
                           &end_ans->force, &end_ans->engv, &eflag, &vflag, &ainum,
-                          &nbor_pitch, &this->_threads_per_atom);
+                          &nbor_pitch, &this->_threads_per_atom, &this->_gpu_nbor);
   }
 
   this->time_pair.stop();
diff --git a/lib/gpu/lal_tersoff_zbl.cu b/lib/gpu/lal_tersoff_zbl.cu
index 0d6c5a3..97eb72b 100644
--- a/lib/gpu/lal_tersoff_zbl.cu
+++ b/lib/gpu/lal_tersoff_zbl.cu
@@ -617,11 +617,12 @@ __kernel void k_tersoff_zbl_three_end(const __global numtyp4 *restrict x_,
                                   const __global acctyp4 *restrict zetaij,
                                   const __global int * dev_nbor,
                                   const __global int * dev_packed,
+                                  const __global int * dev_acc,
                                   __global acctyp4 *restrict ans,
                                   __global acctyp *restrict engv,
                                   const int eflag, const int vflag,
                                   const int inum,  const int nbor_pitch,
-                                  const int t_per_atom) {
+                                  const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;
@@ -686,13 +687,17 @@ __kernel void k_tersoff_zbl_three_end(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
@@ -838,11 +843,12 @@ __kernel void k_tersoff_zbl_three_end_vatom(const __global numtyp4 *restrict x_,
                                         const __global acctyp4 *restrict zetaij,
                                         const __global int * dev_nbor,
                                         const __global int * dev_packed,
+                                        const __global int * dev_acc,
                                         __global acctyp4 *restrict ans,
                                         __global acctyp *restrict engv,
                                         const int eflag, const int vflag,
                                         const int inum,  const int nbor_pitch,
-                                        const int t_per_atom) {
+                                        const int t_per_atom, const int gpu_nbor) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
   numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;
@@ -907,13 +913,17 @@ __kernel void k_tersoff_zbl_three_end_vatom(const __global numtyp4 *restrict x_,
       mdelr1[1] = -delr1[1];
       mdelr1[2] = -delr1[2];
 
-      int nbor_k=j+nbor_pitch;
-      int numk=dev_nbor[nbor_k];
+      int nbor_k,numk;
       if (dev_nbor==dev_packed) {
+        if (gpu_nbor) nbor_k=j+nbor_pitch;
+        else nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
         nbor_k+=offset_k;
       } else {
+        nbor_k=dev_acc[j]+nbor_pitch;
+        numk=dev_nbor[nbor_k];
         nbor_k+=nbor_pitch;
         nbor_k=dev_nbor[nbor_k];
         k_end=nbor_k+numk;
diff --git a/src/GPU/pair_sw_gpu.cpp b/src/GPU/pair_sw_gpu.cpp
index 67faf98..d2d21a1 100644
--- a/src/GPU/pair_sw_gpu.cpp
+++ b/src/GPU/pair_sw_gpu.cpp
@@ -115,7 +115,7 @@ void PairSWGPU::compute(int eflag, int vflag)
     numneigh = list->numneigh;
     firstneigh = list->firstneigh;
 
-    sw_gpu_compute(neighbor->ago, atom->nlocal, nall, inum+list->gnum,
+    sw_gpu_compute(neighbor->ago, inum, nall, inum+list->gnum,
                    atom->x, atom->type, ilist, numneigh, firstneigh, eflag,
                    vflag, eflag_atom, vflag_atom, host_start, cpu_time,
                    success);
diff --git a/src/GPU/pair_tersoff_gpu.cpp b/src/GPU/pair_tersoff_gpu.cpp
index 7b225de..41d8e6b 100644
--- a/src/GPU/pair_tersoff_gpu.cpp
+++ b/src/GPU/pair_tersoff_gpu.cpp
@@ -119,7 +119,7 @@ void PairTersoffGPU::compute(int eflag, int vflag)
     numneigh = list->numneigh;
     firstneigh = list->firstneigh;
 
-    tersoff_gpu_compute(neighbor->ago, atom->nlocal, nall, inum+list->gnum,
+    tersoff_gpu_compute(neighbor->ago, inum, nall, inum+list->gnum,
                    atom->x, atom->type, ilist, numneigh, firstneigh, eflag,
                    vflag, eflag_atom, vflag_atom, host_start, cpu_time,
                    success);
diff --git a/src/GPU/pair_tersoff_mod_gpu.cpp b/src/GPU/pair_tersoff_mod_gpu.cpp
index 3313904..25ea2cb 100644
--- a/src/GPU/pair_tersoff_mod_gpu.cpp
+++ b/src/GPU/pair_tersoff_mod_gpu.cpp
@@ -112,7 +112,7 @@ void PairTersoffMODGPU::compute(int eflag, int vflag)
     numneigh = list->numneigh;
     firstneigh = list->firstneigh;
 
-    tersoff_mod_gpu_compute(neighbor->ago, atom->nlocal, nall, inum+list->gnum,
+    tersoff_mod_gpu_compute(neighbor->ago, inum, nall, inum+list->gnum,
                    atom->x, atom->type, ilist, numneigh, firstneigh, eflag,
                    vflag, eflag_atom, vflag_atom, host_start, cpu_time,
                    success);
diff --git a/src/GPU/pair_tersoff_zbl_gpu.cpp b/src/GPU/pair_tersoff_zbl_gpu.cpp
index 86f02e4..b4b679b 100644
--- a/src/GPU/pair_tersoff_zbl_gpu.cpp
+++ b/src/GPU/pair_tersoff_zbl_gpu.cpp
@@ -120,7 +120,7 @@ void PairTersoffZBLGPU::compute(int eflag, int vflag)
     numneigh = list->numneigh;
     firstneigh = list->firstneigh;
 
-    tersoff_zbl_gpu_compute(neighbor->ago, atom->nlocal, nall, inum+list->gnum,
+    tersoff_zbl_gpu_compute(neighbor->ago, inum, nall, inum+list->gnum,
                    atom->x, atom->type, ilist, numneigh, firstneigh, eflag,
                    vflag, eflag_atom, vflag_atom, host_start, cpu_time,
                    success);
