     Real4 rhoPreMuA = sortedRhoPreMu[index];
     Real3 velMasA = sortedVelMas[index];
     Real SuppRadii = RESOLUTION_LENGTH_MULT * paramsD.HSML;
+    Real SqRadii = SuppRadii * SuppRadii;
 
     Real3 posRadA = mR3(sortedPosRad[index]);
     Real3 deltaV = mR3(0);

            uint j = j_list[n];
            Real3 posRadB = mR3(sortedPosRad[j]);
            Real3 rij = Distance(posRadA, posRadB);
            // Real dd = rij.x*rij.x + rij.y*rij.y + rij.z*rij.z;
            // if (dd > SqRadii || sortedRhoPreMu[j].w < -1.5)
            //     continue;
            Real3 grad_i_wij = GradWh(rij, hA);
            Real3 grw_vj = grad_i_wij * paramsD.volume0;
            mGi[0] -= rij.x * grw_vj.x;

                        B[6 * 1 + 5] += com_part * ZZ;  // 25

                        // mn=13
                        com_part = (A_i[2] * eij.x + A_i[11] * eij.y + 
                            A_i[20] * eij.z + rij.x * eij.z) * V_j;
                        B[6 * 2 + 0] += com_part * XX;  // 31
                        B[6 * 2 + 1] += com_part * XY;  // 32
                        B[6 * 2 + 2] += com_part * XZ;  // 33
