diff --git a/cpu/Color.cpp b/cpu/Color.cpp
index 4fd64f6a..52f7a9b3 100644
--- a/cpu/Color.cpp
+++ b/cpu/Color.cpp
@@ -1109,7 +1109,7 @@ extern "C" void MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd, do
 	//	int sendLoc;
 
 	double f0,f1,f2,f3,f4,f5,f6;
-	double na,nb;		// density values
+	double na,nb,nab;		// density values
 	double ux,uy,uz;	// flow velocity
 	double nx,ny,nz,C;	// color gradient components
 	double a1,a2,b1,b2;
@@ -1156,6 +1156,7 @@ extern "C" void MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd, do
 			f3 = B_even[2*N+n];
 			f5 = B_even[3*N+n];
 			nb = f0+f1+f2+f3+f4+f5+f6;
+			nab = 1.0/(na+nb);
 			//........................................................................
 			//....Instantiate the density distributions
 			// Generate Equilibrium Distributions and stream
@@ -1173,17 +1174,13 @@ extern "C" void MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd, do
 			//...............................................
 			// q = 0,2,4
 			// Cq = {1,0,0}, {0,1,0}, {0,0,1}
-			a1 = na*(0.1111111111111111*(1+4.5*ux));
-			b1 = nb*(0.1111111111111111*(1+4.5*ux));
-			a2 = na*(0.1111111111111111*(1-4.5*ux));
-			b2 = nb*(0.1111111111111111*(1-4.5*ux));
-			if (na*nb > 0.0){
-				delta = na*nb/(na+nb)*0.1111111111111111*nx;
-				a1 += beta*delta;
-				a2 -= beta*delta;
-				b1 -= beta*delta;
-				b2 += beta*delta;
-			}
+			delta = beta*na*nb*nab*0.1111111111111111*nx;
+			if (!(na*nb*nab>0)) delta=0;
+			a1 = na*(0.1111111111111111*(1+4.5*ux))+delta;
+			b1 = nb*(0.1111111111111111*(1+4.5*ux))-delta;
+			a2 = na*(0.1111111111111111*(1-4.5*ux))-delta;
+			b2 = nb*(0.1111111111111111*(1-4.5*ux))+delta;
+
 			A_odd[n] 	= a1;
 			A_even[N+n] = a2;
 			B_odd[n] 	= b1;
@@ -1191,17 +1188,13 @@ extern "C" void MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd, do
 			//...............................................
 			// q = 2
 			// Cq = {0,1,0}
-			a1 = na*(0.1111111111111111*(1+4.5*uy));
-			b1 = nb*(0.1111111111111111*(1+4.5*uy));
-			a2 = na*(0.1111111111111111*(1-4.5*uy));
-			b2 = nb*(0.1111111111111111*(1-4.5*uy));
-			if (na*nb > 0.0){
-				delta = na*nb/(na+nb)*0.1111111111111111*ny;
-				a1 += beta*delta;
-				a2 -= beta*delta;
-				b1 -= beta*delta;
-				b2 += beta*delta;
-			}
+			delta = beta*na*nb*nab*0.1111111111111111*ny;
+			if (!(na*nb*nab>0)) delta=0;
+			a1 = na*(0.1111111111111111*(1+4.5*uy))+delta;
+			b1 = nb*(0.1111111111111111*(1+4.5*uy))-delta;
+			a2 = na*(0.1111111111111111*(1-4.5*uy))-delta;
+			b2 = nb*(0.1111111111111111*(1-4.5*uy))+delta;
+
 			A_odd[N+n] 	= a1;
 			A_even[2*N+n] = a2;
 			B_odd[N+n] 	= b1;
@@ -1209,21 +1202,12 @@ extern "C" void MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd, do
 			//...............................................
 			// q = 4
 			// Cq = {0,0,1}
-			a1 = na*(0.1111111111111111*(1+4.5*uz));
-			b1 = nb*(0.1111111111111111*(1+4.5*uz));
-			a2 = na*(0.1111111111111111*(1-4.5*uz));
-			b2 = nb*(0.1111111111111111*(1-4.5*uz));
-			if (na*nb > 0.0){
-				delta = na*nb/(na+nb)*0.1111111111111111*nz;
-				a1 += beta*delta;
-				a2 -= beta*delta;
-				b1 -= beta*delta;
-				b2 += beta*delta;
-			}
-			A_odd[2*N+n] = a1;
-			A_even[3*N+n] = a2;
-			B_odd[2*N+n] = b1;
-			B_even[3*N+n] = b2;
+			delta = beta*na*nb*nab*0.1111111111111111*nz;
+			if (!(na*nb*nab>0)) delta=0;
+			a1 = na*(0.1111111111111111*(1+4.5*uz))+delta;
+			b1 = nb*(0.1111111111111111*(1+4.5*uz))-delta;
+			a2 = na*(0.1111111111111111*(1-4.5*uz))-delta;
+			b2 = nb*(0.1111111111111111*(1-4.5*uz)+delta);
 			//...............................................
 
 	/*		// Construction and streaming for the components