 	
 	__syncthreads();
 	
+	bool is_valid = (bin_pairs[threadIdx.x].bin < Mx*My*Mz);
+	
 	// now: each unique start point does it's own atomicAdd to find the starting offset
-	__shared__ unsigned int start_offset[block_size];
-	
-	if (is_unique)
-		start_offset[unique[threadIdx.x]] = atomicAdd(&d_bin_size[bin_pairs[threadIdx.x].bin], start[unique[threadIdx.x]+1] - start[unique[threadIdx.x]]);
+	// the is_valid check is to prevent writing to out of bounds memory at the tail end of the array
+	if (is_unique && is_valid)
+		bin_pairs[unique[threadIdx.x]].start_offset = atomicAdd(&d_bin_size[bin_pairs[threadIdx.x].bin], start[unique[threadIdx.x]+1] - start[unique[threadIdx.x]]);
 	
 	__syncthreads();
 	
 	// finally! we can write out all the particles
-	unsigned int offset = start_offset[unique[threadIdx.x]];
-	if (offset + threadIdx.x - start[unique[threadIdx.x]] < Nmax)
+	// the is_valid check is to prevent writing to out of bounds memory at the tail end of the array
+	unsigned int offset = bin_pairs[unique[threadIdx.x]].start_offset;
+	if (offset + threadIdx.x - start[unique[threadIdx.x]] < Nmax && is_valid)
 		d_idxlist[bin_pairs[threadIdx.x].bin*Nmax + offset + threadIdx.x - start[unique[threadIdx.x]]] = bin_pairs[threadIdx.x].id;
 	}
 	
-void rebin_particles_simple_sort(unsigned int *idxlist, unsigned int *bin_size, float4 *pos, unsigned int N, float Lx, float Ly, float Lz, unsigned int Mx, unsigned int My, unsigned int Mz, unsigned int Nmax)
+void rebin_particles_simple_sort(unsigned int *idxlist, unsigned int *bin_size, float4 *pos, unsigned int N, float Lx, float Ly, float Lz, unsigned int Mx, unsigned int My, unsigned int Mz, unsigned int Nmax, unsigned int block_size)
 	{
 	// run one particle per thread
-	const int block_size = 64;
 	int n_blocks = (int)ceil(float(N)/(float)block_size);
 
 	// make even bin dimensions

	scan_naive<unsigned int, block_size>(unique);
	
	// determine start location of each unique value in the array
	__shared__ unsigned int start[block_size+1];
	
	if (is_unique)
		start[unique[threadIdx.x]] = threadIdx.x;

	unsigned int tid = threadIdx.x;
	
	// Parallel bitonic sort.
	#pragma unroll
	for (int k = 2; k <= block_size; k *= 2)
		{
		// Bitonic merge:
		#pragma unroll
		for (int j = k / 2; j>0; j /= 2)
			{
			int ixj = tid ^ j;
