diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H
index 9bcaaeeda3a..51bac0c0820 100644
--- a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H
+++ b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H
@@ -43,20 +43,22 @@ void UpdateMomentumPerezElastic (
     T_R const dt, amrex::RandomEngine const& engine )
 {
 
-    T_PR const diffx = amrex::Math::abs(u1x-u2x);
-    T_PR const diffy = amrex::Math::abs(u1y-u2y);
-    T_PR const diffz = amrex::Math::abs(u1z-u2z);
-    T_PR const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz);
-    T_PR const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z);
-    // If g = u1 - u2 = 0, do not collide.
-    // Or if the relative difference is less than 1.0e-10.
-    if ( diffm < std::numeric_limits<T_PR>::min() || diffm/summm < 1.0e-10 ) { return; }
-
     T_PR constexpr inv_c2 = T_PR(1.0) / ( PhysConst::c * PhysConst::c );
 
     // Compute Lorentz factor gamma
-    T_PR const g1 = std::sqrt( T_PR(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 );
-    T_PR const g2 = std::sqrt( T_PR(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 );
+    T_PR const gb1sq = (u1x*u1x + u1y*u1y + u1z*u1z)*inv_c2;
+    T_PR const gb2sq = (u2x*u2x + u2y*u2y + u2z*u2z)*inv_c2;
+    T_PR const g1 = std::sqrt( T_PR(1.0) + gb1sq );
+    T_PR const g2 = std::sqrt( T_PR(1.0) + gb2sq );
+
+    T_PR const diffx = u1x-u2x;
+    T_PR const diffy = u1y-u2y;
+    T_PR const diffz = u1z-u2z;
+    T_PR const diffm = std::sqrt((diffx*diffx+diffy*diffy+diffz*diffz)*inv_c2);
+    T_PR const summm = std::sqrt(gb1sq) + std::sqrt(gb2sq);
+    // If g = u1 - u2 = 0, do not collide.
+    // Or if the relative difference is less than 1.0e-10.
+    if ( diffm < std::numeric_limits<T_PR>::min() || diffm/summm < 1.0e-10 ) { return; }
 
     // Compute momenta
     T_PR const p1x = u1x * m1;