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

    // Get the interaction from neighbor particles
    for (uint n = 0; n < j_num; n++) {
        uint j = j_list[n];
        Real3 posRadB = mR3(sortedPosRad[j]);
        Real3 dist3 = Distance(posRadA, posRadB);
        Real d = length(dist3);
        Real4 rhoPresMuB = sortedRhoPreMu[j];
        if (rhoPresMuA.w > -0.5 && rhoPresMuB.w > -0.5)
            continue;  // No BCE-BCE interaction

        Real invd = 1.0 / d;
        Real3 velMasB = sortedVelMas[j];
        Real3 TauXxYyZzB = sortedTauXxYyZz[j];

            uint j = j_list[n];
            Real3 posRadB = mR3(sortedPosRad[j]);
            Real3 rij = Distance(posRadA, posRadB);
            Real3 grad_i_wij = GradWh(rij, hA);
            Real3 grw_vj = grad_i_wij * paramsD.volume0;
            mGi[0] -= rij.x * grw_vj.x;
