diff --git a/gpu/Color.cu b/gpu/Color.cu
index f53fe679..24a881f0 100644
--- a/gpu/Color.cu
+++ b/gpu/Color.cu
@@ -1110,87 +1110,98 @@ __global__  void dvc_MassColorCollideD3Q7(char *ID, double *A_even, double *A_od
 			id = ID[n];
 			if ( id != 0){
 
-			//.....Load the Color gradient.........
-			nx = ColorGrad[n];
-			ny = ColorGrad[N+n];
-			nz = ColorGrad[2*N+n];
-			C = sqrt(nx*nx+ny*ny+nz*nz);
-			if (C == 0.0) C=1.0;
-			nx = nx/C;
-			ny = ny/C;
-			nz = nz/C;
-			//....Load the flow velocity...........
-			ux = Velocity[n];
-			uy = Velocity[N+n];
-			uz = Velocity[2*N+n];
-			//........................................................................
-			//					READ THE DISTRIBUTIONS
-			//		(read from opposite array due to previous swap operation)
-			//........................................................................
-			f2 = A_odd[n];
-			f4 = A_odd[N+n];
-			f6 = A_odd[2*N+n];
-			f0 = A_even[n];
-			f1 = A_even[N+n];
-			f3 = A_even[2*N+n];
-			f5 = A_even[3*N+n];
-			na = f0+f1+f2+f3+f4+f5+f6;
-			//........................................................................
-			f2 = B_odd[n];
-			f4 = B_odd[N+n];
-			f6 = B_odd[2*N+n];
-			f0 = B_even[n];
-			f1 = B_even[N+n];
-			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);
-			//........................................................................
-			delta = beta*na*nb*nab*0.1111111111111111*nx;
-		      	if (na*nb*nab<0.0)) delta=0.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;
+            //.....Load the Color gradient.........
+            nx = ColorGrad[n];
+            ny = ColorGrad[N+n];
+            nz = ColorGrad[2*N+n];
+            C = sqrt(nx*nx+ny*ny+nz*nz);
+            if (C==0.0) C=1.0;
+            nx = nx/C;
+            ny = ny/C;
+            nz = nz/C;
+            //....Load the flow velocity...........
+            ux = Velocity[n];
+            uy = Velocity[N+n];
+            uz = Velocity[2*N+n];
+            //........................................................................
+            //					READ THE DISTRIBUTIONS
+            //		(read from opposite array due to previous swap operation)
+            //........................................................................
+            f2 = A_odd[n];
+            f4 = A_odd[N+n];
+            f6 = A_odd[2*N+n];
+            f0 = A_even[n];
+            f1 = A_even[N+n];
+            f3 = A_even[2*N+n];
+            f5 = A_even[3*N+n];
+            na = f0+f1+f2+f3+f4+f5+f6;
+            //........................................................................
+            f2 = B_odd[n];
+            f4 = B_odd[N+n];
+            f6 = B_odd[2*N+n];
+            f0 = B_even[n];
+            f1 = B_even[N+n];
+            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
+            // Stationary value - distribution 0
+            A_even[n] = 0.3333333333333333*na;
+            B_even[n] = 0.3333333333333333*nb;
+            // Non-Stationary equilibrium distributions
+            //feq[0] = 0.1111111111111111*(1+4.5*ux);
+            //feq[1] = 0.1111111111111111*(1-4.5*ux);
+            //feq[2] = 0.1111111111111111*(1+4.5*uy);
+            //feq[3] = 0.1111111111111111*(1-4.5*uy);
+            //feq[4] = 0.1111111111111111*(1+4.5*uz);
+            //feq[5] = 0.1111111111111111*(1-4.5*uz);
 
-			A_odd[n]	= a1;
-			A_even[N+n] = a2;
-			B_odd[n]  = b1;
-			B_even[N+n] = b2;
-			//........................................................................
-			
-			// q = 2
-			// Cq = {0,1,0}
-			delta = beta*na*nb*nab*0.1111111111111111*ny;
-		      	if (na*nb*nab<0.0)) delta=0.0;
+            //...............................................
+            // q = 0,2,4
+            // Cq = {1,0,0}, {0,1,0}, {0,0,1}
+            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;
 
-			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] 	= a1;
+            A_even[N+n] = a2;
+            B_odd[n] 	= b1;
+            B_even[N+n] = b2;
+            //...............................................
+            // q = 2
+            // Cq = {0,1,0}
+            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;
-			B_even[2*N+n] = b2;
-			//........................................................................
-			// q = 4
-			// Cq = {0,0,1}
-			delta = beta*na*nb*nab*0.1111111111111111*nz;
-		      	if (na*nb*nab<0.0)) delta=0.0;
+            A_odd[N+n] 	= a1;
+            A_even[2*N+n] = a2;
+            B_odd[N+n] 	= b1;
+            B_even[2*N+n] = b2;
+            //...............................................
+            // q = 4
+            // Cq = {0,0,1}
+            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;
+
+            A_odd[2*N+n] = a1;
+            A_even[3*N+n] = a2;
+            B_odd[2*N+n] = b1;
+            B_even[3*N+n] = b2;
 
-			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;
-			
-			A_odd[2*N+n] = a1;
-			A_even[3*N+n] = a2;
-			B_odd[2*N+n] = b1;
-			B_even[3*N+n] = b2;
-			//........................................................................
-			
 		}
 	}
 }
diff --git a/gpu/D3Q7.cu b/gpu/D3Q7.cu
index d8de261c..79c13f32 100644
--- a/gpu/D3Q7.cu
+++ b/gpu/D3Q7.cu
@@ -109,57 +109,55 @@ __global__  void dvc_SwapD3Q7(char *ID, double *disteven, double *distodd, int N
 		if (n<N ){
 		   	id = ID[n];
 			if (id > 0){
-			//.......Back out the 3-D indices for node n..............
-			k = n/(Nx*Ny);
-			j = (n-Nx*Ny*k)/Nx;
-			i = n-Nx*Ny*k-Nz*j;
-			//........................................................................
-			// Retrieve even distributions from the local node (swap convention)
-			//		f0 = disteven[n];  // Does not particupate in streaming
-			f1 = distodd[n];
-			f3 = distodd[N+n];
-			f5 = distodd[2*N+n];
-			//........................................................................
+                //.......Back out the 3-D indices for node n..............
+                k = n/(Nx*Ny);
+                j = (n-Nx*Ny*k)/Nx;
+                i = n-Nx*Ny*k-Nz*j;
+                //........................................................................
+                // Retrieve even distributions from the local node (swap convention)
+                //		f0 = disteven[n];  // Does not particupate in streaming
+                f1 = distodd[n];
+                f3 = distodd[N+n];
+                f5 = distodd[2*N+n];
+                //........................................................................
 
-			//........................................................................
-			// Retrieve odd distributions from neighboring nodes (swap convention)
-			//........................................................................
-			nn = n+1;							// neighbor index (pull convention)
-			if (!(i+1<Nx))	nn -= Nx;			// periodic BC along the x-boundary
-			//if (i+1<Nx){
-			f2 = disteven[N+nn];					// pull neighbor for distribution 2
-			if (!(f2 < 0.0)){
-				distodd[n] = f2;
-				disteven[N+nn] = f1;
-			}
-			//}
-			//........................................................................
-			nn = n+Nx;							// neighbor index (pull convention)
-			if (!(j+1<Ny))	nn -= Nx*Ny;		// Perioidic BC along the y-boundary
-			//if (j+1<Ny){
-			f4 = disteven[2*N+nn];				// pull neighbor for distribution 4
-			if (!(f4 < 0.0)){
-				distodd[N+n] = f4;
-				disteven[2*N+nn] = f3;
-				//	}
-			}
-			//........................................................................
-			nn = n+Nx*Ny;						// neighbor index (pull convention)
-			if (!(k+1<Nz))	nn -= Nx*Ny*Nz;		// Perioidic BC along the z-boundary
-			//if (k+1<Nz){
-			f6 = disteven[3*N+nn];				// pull neighbor for distribution 6
-			if (!(f6 < 0.0)){
-				distodd[2*N+n] = f6;
-				disteven[3*N+nn] = f5;
-				//	}
-			}
+                //........................................................................
+                // Retrieve odd distributions from neighboring nodes (swap convention)
+                //........................................................................
+                nn = n+1;							// neighbor index (pull convention)
+                if (!(i+1<Nx))	nn -= Nx;			// periodic BC along the x-boundary
+                //if (i+1<Nx){
+                f2 = disteven[N+nn];					// pull neighbor for distribution 2
+                if (!(f2 < 0.0)){
+                    distodd[n] = f2;
+                    disteven[N+nn] = f1;
+                }
+                //}
+                //........................................................................
+                nn = n+Nx;							// neighbor index (pull convention)
+                if (!(j+1<Ny))	nn -= Nx*Ny;		// Perioidic BC along the y-boundary
+                //if (j+1<Ny){
+                f4 = disteven[2*N+nn];				// pull neighbor for distribution 4
+                if (!(f4 < 0.0)){
+                    distodd[N+n] = f4;
+                    disteven[2*N+nn] = f3;
+                }
+                //........................................................................
+                nn = n+Nx*Ny;						// neighbor index (pull convention)
+                if (!(k+1<Nz))	nn -= Nx*Ny*Nz;		// Perioidic BC along the z-boundary
+                //if (k+1<Nz){
+                f6 = disteven[3*N+nn];				// pull neighbor for distribution 6
+                if (!(f6 < 0.0)){
+                    distodd[2*N+n] = f6;
+                    disteven[3*N+nn] = f5;
+                }
 			}
 		}
 	}
 }
 
 //*************************************************************************
-__global__  void dvc_ComputeDensityD3Q7(char *ID, double *disteven, double *distodd, double *Den, 
+__global__  void dvc_ComputeDensityD3Q7(char *ID, double *disteven, double *distodd, double *Den,
 		int Nx, int Ny, int Nz)
 {
 	char id;
@@ -172,19 +170,19 @@ __global__  void dvc_ComputeDensityD3Q7(char *ID, double *disteven, double *dist
 		//........Get 1-D index for this thread....................
 		n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
 		if (n<N){
-		id = ID[n];
-		if ( id > 0 ){
-			// Read the distributions
-			f0 = disteven[n];
-			f2 = disteven[N+n];
-			f4 = disteven[2*N+n];
-			f6 = disteven[3*N+n];
-			f1 = distodd[n];
-			f3 = distodd[N+n];
-			f5 = distodd[2*N+n];
-			// Compute the density
-			Den[n] = f0+f1+f2+f3+f4+f5+f6;
-			}
+            id = ID[n];
+            if (id > 0 ){
+                // Read the distributions
+                f0 = disteven[n];
+                f2 = disteven[N+n];
+                f4 = disteven[2*N+n];
+                f6 = disteven[3*N+n];
+                f1 = distodd[n];
+                f3 = distodd[N+n];
+                f5 = distodd[2*N+n];
+                // Compute the density
+                Den[n] = f0+f1+f2+f3+f4+f5+f6;
+            }
 		}
 	}
 }