//	calculates the normalized force for N particles
//	from the shielded interaction with all N-1 other particles
//	the constant factor Q1 Q2 / 4 pi epsilon_0 must be applied externally

	__device__ __constant__ float kdebye;	// inverse shielding length
	__device__ __constant__ float eps2;	// square of softening distance
//------------------------------------------------------------------------
	__device__ float3 pairInteraction(float4 bi, float4 bj,float3 ai)
	{
		float3 r;
		r.x = bj.x - bi.x;
		r.y = bj.y - bi.y;
		r.z = bj.z - bi.z;
		float dist2 = r.x*r.x + r.y*r.y + r.z*r.z + eps2;
		float rdist = rsqrtf(dist2);
		float kappa = kdebye/rdist;
		float s = rdist*rdist*rdist*(1.f+kappa)*__expf(-kappa);
		ai.x += r.x * s;
		ai.y += r.y * s;
		ai.z += r.z * s;
		return ai;
	}
//---------------------------------------------------------------------------
	__device__ float3 evalTile(float4 iPos, float3 acc)
	{
		extern __shared__ float4 sharedPos[];
		for(int j=0; j<blockDim.x; j++)
		{
			float4 jPos=sharedPos[j];
			acc = pairInteraction(iPos,jPos,acc);
		}
		return acc;
	}
//---------------------------------------------------------------------------
	__global__ void calcForce(void *d_pos, void *d_acc, int N)
	{
		extern __shared__ float4 sharedPos[];
		float4 *globalPos = (float4 *)d_pos;
		float4 *globalAcc = (float4 *)d_acc;
		float3 accel = {0.f,0.f,0.f};
		float4 iPos;
		int i,tile;
		int gtid = blockIdx.x*blockDim.x+threadIdx.x;
		iPos = globalPos[gtid];
		for(i=0, tile=0; i<N; i+=blockDim.x, tile++)
		{
			int gidx = tile*blockDim.x + threadIdx.x;
			sharedPos[threadIdx.x] = globalPos[gidx];
			__syncthreads();
			accel = evalTile(iPos,accel);
			__syncthreads();
		}
		float4 acc4 = {accel.x,accel.y,accel.z,0.0f};
		globalAcc[gtid] = acc4;
	}