computPartGridForce wraps about RB center, not about grid origin; may result in performance loss

src/ComputeGridGrid.cu

@@ -76,7 +76,7 @@ __global__
 void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForce,
 				const int num, const int* __restrict__ particleIds, 
 				const RigidBodyGrid* __restrict__ u,
-				const Matrix3 basis_u_inv, const Vector3 origin_u,
+				const Matrix3 basis_u_inv, const Vector3 center_u, const Vector3 origin_u,
 				ForceEnergy* __restrict__ retForceTorque, float* __restrict__ energy, bool get_energy, int scheme, BaseGrid* sys_d) {
 
 	extern __shared__ ForceEnergy s[];
@@ -90,9 +90,7 @@ void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForc
 	torque[tid] = ForceEnergy(0.f,0.f);
 	if (i < num) {
 		const int id = particleIds[i];
-		//Vector3 p = pos[id] - origin_u;
-		Vector3 p = sys_d->wrapDiff(pos[id]-origin_u); /* TODO: wrap about RB center, not origin */
-		// TODO: wrap to center of u
+		Vector3 p = sys_d->wrapDiff(pos[id]-center_u) + center_u - origin_u
 		const Vector3 u_ijk_float = basis_u_inv.transform( p );
 
                 ForceEnergy fe;

src/ComputeGridGrid.cuh

@@ -17,7 +17,7 @@ extern __global__
 void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForce,
 				const int num, const int* __restrict__ particleIds,
 				const RigidBodyGrid* __restrict__ u,
-				const Matrix3 basis_u_inv, const Vector3 origin_u,
+				const Matrix3 basis_u_inv, const Vector3 center_u, const Vector3 origin_u,
 				ForceEnergy* __restrict__ retForceTorque, float* energy, bool get_energy, int scheme, BaseGrid* sys_d);
 
 extern __global__

src/RigidBody.cu

@@ -184,7 +184,7 @@ void RigidBody::callGridParticleForceKernel(Vector3* pos_d, Vector3* force_d, in
 		computePartGridForce<<< nb, NUMTHREADS, NUMTHREADS*2*sizeof(ForceEnergy), stream >>>(
 			pos_d, force_d, numParticles[i], particles_d[i],
 			t->RBC->grids_d+idx,
-			B, c, forcestorques_d+forcestorques_offset[fto_idx++], energy, get_energy, scheme, sys_d);
+			B, getPosition(), c, forcestorques_d+forcestorques_offset[fto_idx++], energy, get_energy, scheme, sys_d);
 	}
 }