 	This version dynamically avoids divergent warps
 	It is also assumed that the neighborlist is padded out to have a number of columsn a 
 	multiple of the block size so that reads can be simply done in a big block
-	\bug Performance is extremely slow due to an implemented workaround for the annoying timeout bug
-	\todo Remove limitation on excessively padded neighbor list
-	\todo Remove the warp divergence avoiding code (it is here as part of the bug workaround)
 */
 extern "C" __global__ void calcLJForces_kernel(float4 *d_forces, gpu_pdata_arrays pdata, gpu_nlist_array nlist, float r_cutsq, gpu_boxsize box)
 	{
 	// start by identifying which particle we are to handle
 	int pidx = blockIdx.x * blockDim.x + threadIdx.x;
+	
+	if (pidx >= pdata.N)
+		return;
 	
 	// load in the length of the list (each thread loads it individually, 
 	// we are going to dynamically avoid divergent warps by finding the maximum

	return retval;
	}

//! Hardcode the warp size for the divergence avoiding code
#define WARP_SIZE 32

//! Kernel for calculating lj forces
/*! This kerenel is called to calculate the lennard-jones forces on all N particles


		period = 1;
	
	// loop over neighbors
	// to prevent triggering a bug in the GPU hardware, each thread loops over the maximum
	// possible number of neighbors. Then an if in the thread body prevents computations 
	// from being performed on invalid neighbors
	for (int neigh_idx = 0; neigh_idx < nlist.height; neigh_idx++)
		{
		if (neigh_idx < n_neigh)
			{
			int cur_neigh = nlist.list[nlist.pitch*(neigh_idx+1) + pidx];
		
			// get the neighbor's position
			float4 neigh_pos = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
			neigh_pos = tex1Dfetch(pdata_pos_tex, cur_neigh);
	
			float nx = neigh_pos.x;
			float ny = neigh_pos.y;
			float nz = neigh_pos.z;
			int neigh_typ = __float_as_int(neigh_pos.w);
		
			// calculate dr (with periodic boundary conditions)
			float dx = pos.x - nx;
			float dy = pos.y - ny;
			float dz = pos.z - nz;
			
			if (period)
				{
				dx -= box.Lx * rintf(dx * box.Lxinv);
				dy -= box.Ly * rintf(dy * box.Lyinv);
				dz -= box.Lz * rintf(dz * box.Lzinv);
				}
				
			float rsq = dx*dx + dy*dy + dz*dz;
			
			float r2inv;
			if (rsq >= r_cutsq)
				r2inv = 0.0f;
			else
				r2inv = 1.0f / rsq;
	
			int ptype = __float_as_int(pos.w);
			int typ_pair = neigh_typ * MAX_NTYPES + ptype;
	
			// note that the code for calculating the force was borrowed from lammps
			float r6inv = r2inv*r2inv*r2inv;
			// the lj1 and lj2 params are encoded into the x and y components of the params array
			float lj1 = c_ljparams.lj1[typ_pair];
			float lj2 = c_ljparams.lj2[typ_pair];
			float fforce = r2inv * r6inv * (lj1  * r6inv - lj2);
			
			// add up the forces
			fx += dx * fforce;
			fy += dy * fforce;
			fz += dz * fforce;
			}
		}
	
	// now that the force calculation is complete, write out the result if we are a valid particle
	float4 force;
	force.x = fx;
	force.y = fy;
	force.z = fz;
	force.w = 0.0f;
	d_forces[pidx] = force;
	}


