                     for (uint j = startIndex; j < endIndex; j++) {
                         Real3 posRadB = mR3(sortedPosRad[j]);
                         Real3 rij = Distance(posRadA, posRadB);
-                        Real d = length(rij);
-                        if (d > SuppRadii || sortedRhoPreMu[j].w < -1.5)
+                        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, h_i);
                         Real3 grw_vj = grad_i_wij * paramsD.volume0;

                        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

    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);
