 	torquingTerm.z = (-omega3.y * j1.x + omega3.x * j1.y) * omega3.x + (-omega3.y * j1.y + omega3.x * j2.x) * omega3.y
 			+ (-omega3.y * j1.z + omega3.x * j2.y) * omega3.z;
 
-	torquingTerm = solid_SPH_massD * LF_totalTorque3[rigidSphereA] - torquingTerm;
+	torquingTerm = solid_SPH_massD * LF_totalTorqueOfAcc3[rigidSphereA] - torquingTerm;
 	//*** from this point j1 and j2 will represent the j_Inverse
 	j1 = jInvD1[rigidSphereA];
 	j2 = jInvD2[rigidSphereA];

		ForceSPH(posRadD, velMasD, vel_XSPH_D, rhoPresMuD, bodyIndexD, derivVelRhoD, referenceArray, numAllMarkers, paramsH, 0.5 * delT); //?$ right now, it does not consider paramsH.gravity or other stuff on rigid bodies. they should be applied at rigid body solver
		UpdateFluid(posRadD2, velMasD2, vel_XSPH_D, rhoPresMuD2, derivVelRhoD, referenceArray, 0.5 * delT); //assumes ...D2 is a copy of ...D
		//UpdateBoundary(posRadD2, velMasD2, rhoPresMuD2, derivVelRhoD, referenceArray, 0.5 * delT);		//assumes ...D2 is a copy of ...D
		UpdateRigidBody(posRadD2, velMasD2, posRigidD2, posRadRigidCumulativeD2, velMassRigidD2, qD2, AD1_2, AD2_2, AD3_2, omegaLRF_D2, derivVelRhoD, rigidIdentifierD,
				rigidSPH_MeshPos_LRF_D, referenceArray, jD1, jD2, jInvD1, jInvD2, paramsH, numRigidBodies, startRigidMarkers, numRigid_SphMarkers, float(tStep)/stepEnd, 0.5 * delT);

		UpdateFlexibleBody(posRadD2, velMasD2,

			derivRigid *= ( .2 * paramsD.HSML / (dTD * dTD) ) / length(derivRigid);
	}
	//** end tweak
	totalAccRigid3[rigidSphereA] = derivRigid;
}
//--------------------------------------------------------------------------------------------------------------------------------
__global__ void CalcTorqueOf_SPH_Marker_Acceleration(real3* torqueMarkersD, real4* derivVelRhoD, real3* posRadD, int* rigidIdentifierD, real3* posRigidD) {
	uint index = blockIdx.x * blockDim.x + threadIdx.x;
	uint rigidMarkerIndex = index + startRigidMarkersD;
	if (index >= numRigid_SphMarkersD) {
