         make_nonbonded_plan(box, r_ghost));
     }
 
+//! Kernel to construct list of ghosts to send
+__global__ void gpu_make_exchange_ghost_list_kernel(const unsigned int n_total,
+                                         const unsigned int N,
+                                         const unsigned int dir,
+                                         const unsigned char *plan,
+                                         const unsigned int *tag,
+                                         unsigned int *copy_ghosts,
+                                         unsigned int *n_copy_ghosts)
+    {
+    unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x;
+    
+    if (idx >= n_total) return;
+
+    unsigned char p = plan[idx];
+    bool do_send = (p & (1 << dir));
+
+    // do not send back ghost particles just received from the opposite direction
+    if ((dir%2) && (p & (1 << (dir-1))) && idx >= N)
+        do_send = false;
+
+    if (do_send)
+        {
+        unsigned int n = atomicInc(n_copy_ghosts, 0xffffffff);
+        copy_ghosts[n] = tag[idx];
+        }
+    }
+
+    
 //! Construct a list of particle tags to send as ghost particles
 /*! \param n_total Total number of particles to check
  * \param N number of local particles

                                  unsigned int dir,
                                  unsigned char *d_plan,
                                  unsigned int *d_global_tag,
                                  unsigned int* d_copy_ghosts,
                                  unsigned int &n_copy_ghosts)
    {
    thrust::device_ptr<unsigned char> plan_ptr(d_plan);
    thrust::device_ptr<unsigned int> global_tag_ptr(d_global_tag);
    thrust::device_ptr<unsigned int> copy_ghosts_ptr(d_copy_ghosts);

    thrust::device_ptr<unsigned int> copy_ghosts_end_ptr;

    thrust::counting_iterator<unsigned int> idx_it(0);

    copy_ghosts_end_ptr = thrust::copy_if(global_tag_ptr,
                                          global_tag_ptr+n_total,
                                          thrust::make_zip_iterator(thrust::make_tuple(
                                            idx_it,
                                           plan_ptr)),
                                          copy_ghosts_ptr,
                                          select_particle_ghost(dir,N));

    n_copy_ghosts =  copy_ghosts_end_ptr - copy_ghosts_ptr;
    }

//! Fill send buffers of particles we are sending as ghost particles with partial particle data


thrust::device_vector<unsigned int> *keys;       //!< Temporary vector of sort keys

unsigned int *d_n_copy_ghosts;     //! Counter for ghost list construction

void gpu_allocate_tmp_storage()
    {
    keys = new thrust::device_vector<unsigned int>;

    cudaMalloc(&d_n_copy_ghosts, sizeof(unsigned int));
    }

void gpu_deallocate_tmp_storage()
    {
    delete keys;

    cudaFree(d_n_copy_ghosts);
    }

//! GPU Kernel to find incomplete bonds
