 //gridindex is which grid they belong in
 //radii are atom radii
 //grids are the output and are assumed to be zeroed
-__global__ void binary_set(float3 origin, int dim, float resolution, int n, float3 *coords, short *gridindex, float *radii, float *grids) 
+template<bool Binary> __global__ void gpu_grid_set(float3 origin, int dim, float resolution, float rmult, int n, float3 *coords, short *gridindex, float *radii, float *grids) 
 {
 	//figure out what grid point we are 
 	unsigned xi = threadIdx.x + blockIdx.x*blockDim.x;

		short which = gridindex[i];
		if(which >= 0) { //because of hydrogens on ligands
			float r = radii[i];
			r *= r; //square radius
			float d = sqDistance(coord, x,y,z);
			
			if(d < r)
			{
				//set gridpoint to 1
				unsigned goffset = which*gsize;
				unsigned off = (xi*dim+yi)*dim+zi;
				//printf("%f,%f,%f %d,%d,%d  %d  %d %d\n",x,y,z, xi,yi,zi, which, goffset,off);
				grids[goffset+off] = 1.0;
			}
		}
	}

	CUDA_CHECK(cudaMemset(grids, 0, ngrids*dim*dim*dim*sizeof(float))); //TODO: see if faster to do in kernel
	if(binary)
	{
		gpu_grid_set<true><<<blocks,threads>>>(origin, dim, resolution, radiusmultiple, natoms, coords, gridindex, radii, grids);
		CUDA_CHECK(cudaPeekAtLastError() );
	}
	else 
	{
		gpu_grid_set<false><<<blocks,threads>>>(origin, dim, resolution, radiusmultiple, natoms, coords, gridindex, radii, grids);
		CUDA_CHECK(cudaPeekAtLastError() );
	}
}
