diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 38878067..f48419be 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -86,18 +86,18 @@ extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, doubl
 extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int Np);
 
 extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
-		double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Phi, double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *NeighborList, double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np);
 
-extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *NeighborList, double *Phi, double *ColorGrad, int start, int finish, int Np);
+extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *NeighborList, double *Phi, double *ColorGrad, double *Potential, int start, int finish, int Np);
 
 // BOUNDARY CONDITION ROUTINES
 
diff --git a/cpu/dfh.cpp b/cpu/dfh.cpp
index 990b0292..0a85936a 100644
--- a/cpu/dfh.cpp
+++ b/cpu/dfh.cpp
@@ -59,7 +59,7 @@ extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq
 
 // LBM based on density functional hydrodynamics
 extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
-		double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 
 	int ijk,nn,n;
@@ -105,60 +105,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double
 		rlx_setA = 1.f/tau;
 		rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 
-		//					COMPUTE THE COLOR GRADIENT
-		//........................................................................
-		//.................Read Phase Indicator Values............................
-		//........................................................................
-		nn = neighborList[n+Np]%Np;
-		m1 = Phi[nn];
-		nn = neighborList[n]%Np;
-		m2 = Phi[nn];
-		nn = neighborList[n+3*Np]%Np;
-		m3 = Phi[nn];
-		nn = neighborList[n+2*Np]%Np;
-		m4 = Phi[nn];		
-		nn = neighborList[n+5*Np]%Np;
-		m5 = Phi[nn];
-		nn = neighborList[n+4*Np]%Np;
-		m6 = Phi[nn];		
-		nn = neighborList[n+7*Np]%Np;
-		m7 = Phi[nn];
-		nn = neighborList[n+6*Np]%Np;
-		m8 = Phi[nn];		
-		nn = neighborList[n+9*Np]%Np;
-		m9 = Phi[nn];
-		nn = neighborList[n+8*Np]%Np;
-		m10 = Phi[nn];		
-		nn = neighborList[n+11*Np]%Np;
-		m11 = Phi[nn];
-		nn = neighborList[n+10*Np]%Np;
-		m12 = Phi[nn];		
-		nn = neighborList[n+13*Np]%Np;
-		m13 = Phi[nn];
-		nn = neighborList[n+12*Np]%Np;
-		m14 = Phi[nn];		
-		nn = neighborList[n+15*Np]%Np;
-		m15 = Phi[nn];
-		nn = neighborList[n+14*Np]%Np;
-		m16 = Phi[nn];		
-		nn = neighborList[n+17*Np]%Np;
-		m17 = Phi[nn];
-		nn = neighborList[n+16*Np]%Np;
-		m18 = Phi[nn];					
-		//............Compute the Color Gradient...................................
-		//............Compute the wn fluid Gradient...................................
-		nx = (m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-		ny = (m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-		nz = (m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));		
-		// .... read the solid potential gradient.....................
-		m1 = SolidPotential[n];
-		m2 = SolidPotential[n+Np];
-		m3 = SolidPotential[n+2*Np];
-		nx += m1;
-		ny += m2;
-		nz += m3;
-
-		//...........Normalize the Color Gradient.................................
+		//...........Read the Color Gradient.................................
+		nx = Gradient[n];
+		ny = Gradient[n+Np];
+		bz = Gradient[n+2*Np];
 		C = sqrt(nx*nx+ny*ny+nz*nz);
 		if (C==0.0) C=1.0;
 		nx = nx/C;
@@ -626,11 +576,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double
 	
 }
 
-//extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Velocity,
-//		double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
-//		double Fx, double Fy, double Fz, int start, int finish, int Np){
 extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Phi, double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 	
 	int n,nn,ijk,nread;
@@ -679,64 +626,15 @@ extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *
 		rlx_setA = 1.f/tau;
 		rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 		
-		//					COMPUTE THE COLOR GRADIENT
-		//........................................................................
-		//.................Read Phase Indicator Values............................
-		//........................................................................
-		nn = neighborList[n+Np]%Np;
-		m1 = Phi[nn];
-		nn = neighborList[n]%Np;
-		m2 = Phi[nn];
-		nn = neighborList[n+3*Np]%Np;
-		m3 = Phi[nn];
-		nn = neighborList[n+2*Np]%Np;
-		m4 = Phi[nn];		
-		nn = neighborList[n+5*Np]%Np;
-		m5 = Phi[nn];
-		nn = neighborList[n+4*Np]%Np;
-		m6 = Phi[nn];		
-		nn = neighborList[n+7*Np]%Np;
-		m7 = Phi[nn];
-		nn = neighborList[n+6*Np]%Np;
-		m8 = Phi[nn];		
-		nn = neighborList[n+9*Np]%Np;
-		m9 = Phi[nn];
-		nn = neighborList[n+8*Np]%Np;
-		m10 = Phi[nn];		
-		nn = neighborList[n+11*Np]%Np;
-		m11 = Phi[nn];
-		nn = neighborList[n+10*Np]%Np;
-		m12 = Phi[nn];		
-		nn = neighborList[n+13*Np]%Np;
-		m13 = Phi[nn];
-		nn = neighborList[n+12*Np]%Np;
-		m14 = Phi[nn];		
-		nn = neighborList[n+15*Np]%Np;
-		m15 = Phi[nn];
-		nn = neighborList[n+14*Np]%Np;
-		m16 = Phi[nn];		
-		nn = neighborList[n+17*Np]%Np;
-		m17 = Phi[nn];
-		nn = neighborList[n+16*Np]%Np;
-		m18 = Phi[nn];					
-		//............Compute the Color Gradient...................................
-		//............Compute the wn fluid Gradient...................................
-		nx = (m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-		ny = (m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-		nz = (m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-		// .... read the solid potential gradient.....................
-		m1 = SolidPotential[n];
-		m2 = SolidPotential[n+Np];
-		m3 = SolidPotential[n+2*Np];
-		nx += m1;
-		ny += m2;
-		nz += m3;
-		//...........Normalize the Color Gradient.................................
+		//...........Read the Color Gradient.................................
+		nx = Gradient[n];
+		ny = Gradient[n+Np];
+		bz = Gradient[n+2*Np];
 		C = sqrt(nx*nx+ny*ny+nz*nz);
 		if (C==0.0) C=1.0;
 		nx = nx/C;
 		ny = ny/C;
-		nz = nz/C;		
+		nz = nz/C;			
 
 		// q=0
 		fq = dist[n];
@@ -1433,7 +1331,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, doub
 	}	
 }
 
-extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, double *SolidPotential, int start, int finish, int Np){
 	
 	int n,nn;
 	// distributions
@@ -1488,6 +1386,14 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double
 		ny = (m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
 		nz = (m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
 		
+		// .... read the solid potential gradient.....................
+		m1 = SolidPotential[n];
+		m2 = SolidPotential[n+Np];
+		m3 = SolidPotential[n+2*Np];
+		nx += m1;
+		ny += m2;
+		nz += m3;
+		
 		//...........Normalize the Color Gradient.................................
 		//	C = sqrt(nx*nx+ny*ny+nz*nz);
 		//	nx = nx/C;
diff --git a/gpu/dfh.cu b/gpu/dfh.cu
index 50835a18..86b281e7 100644
--- a/gpu/dfh.cu
+++ b/gpu/dfh.cu
@@ -75,7 +75,7 @@ __global__ void dvc_ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double
 
 // LBM based on density functional hydrodynamics
 __global__ void dvc_ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
-			double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+			double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 			double Fx, double Fy, double Fz, int start, int finish, int Np){
 	int nn,n;
 	double fq;
@@ -123,61 +123,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, dou
 				rlx_setA = 1.f/tau;
 				rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 
-				//					COMPUTE THE COLOR GRADIENT
-				//........................................................................
-				//.................Read Phase Indicator Values............................
-				//........................................................................
-				nn = neighborList[n+Np]%Np;
-				m1 = Phi[nn];
-				nn = neighborList[n]%Np;
-				m2 = Phi[nn];
-				nn = neighborList[n+3*Np]%Np;
-				m3 = Phi[nn];
-				nn = neighborList[n+2*Np]%Np;
-				m4 = Phi[nn];		
-				nn = neighborList[n+5*Np]%Np;
-				m5 = Phi[nn];
-				nn = neighborList[n+4*Np]%Np;
-				m6 = Phi[nn];		
-				nn = neighborList[n+7*Np]%Np;
-				m7 = Phi[nn];
-				nn = neighborList[n+6*Np]%Np;
-				m8 = Phi[nn];		
-				nn = neighborList[n+9*Np]%Np;
-				m9 = Phi[nn];
-				nn = neighborList[n+8*Np]%Np;
-				m10 = Phi[nn];		
-				nn = neighborList[n+11*Np]%Np;
-				m11 = Phi[nn];
-				nn = neighborList[n+10*Np]%Np;
-				m12 = Phi[nn];		
-				nn = neighborList[n+13*Np]%Np;
-				m13 = Phi[nn];
-				nn = neighborList[n+12*Np]%Np;
-				m14 = Phi[nn];		
-				nn = neighborList[n+15*Np]%Np;
-				m15 = Phi[nn];
-				nn = neighborList[n+14*Np]%Np;
-				m16 = Phi[nn];		
-				nn = neighborList[n+17*Np]%Np;
-				m17 = Phi[nn];
-				nn = neighborList[n+16*Np]%Np;
-				m18 = Phi[nn];					
-				
-				//............Compute the wn fluid Gradient...................................
-				nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-				ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-				nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-				
-				// .... read the solid potential gradient.....................
-				m1 = SolidPotential[n];
-				m2 = SolidPotential[n+Np];
-				m3 = SolidPotential[n+2*Np];
-				nx += m1;
-				ny += m2;
-				nz += m3;
-
-				//...........Normalize the Color Gradient.................................
+				//...........Read the Color Gradient.................................
+				nx = Gradient[n];
+				ny = Gradient[n+Np];
+				bz = Gradient[n+2*Np];
 				C = sqrt(nx*nx+ny*ny+nz*nz);
 				if (C==0.0) C=1.0;
 				nx = nx/C;
@@ -647,7 +596,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, dou
 
 
 __global__ void dvc_ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Phi, double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 
 	int n,nn,nread;
@@ -699,60 +648,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, doub
 			rlx_setA = 1.f/tau;
 			rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 			
-			//					COMPUTE THE COLOR GRADIENT
-			//........................................................................
-			//.................Read Phase Indicator Values............................
-			//........................................................................
-			nn = neighborList[n+Np]%Np;
-			m1 = Phi[nn];
-			nn = neighborList[n]%Np;
-			m2 = Phi[nn];
-			nn = neighborList[n+3*Np]%Np;
-			m3 = Phi[nn];
-			nn = neighborList[n+2*Np]%Np;
-			m4 = Phi[nn];		
-			nn = neighborList[n+5*Np]%Np;
-			m5 = Phi[nn];
-			nn = neighborList[n+4*Np]%Np;
-			m6 = Phi[nn];		
-			nn = neighborList[n+7*Np]%Np;
-			m7 = Phi[nn];
-			nn = neighborList[n+6*Np]%Np;
-			m8 = Phi[nn];		
-			nn = neighborList[n+9*Np]%Np;
-			m9 = Phi[nn];
-			nn = neighborList[n+8*Np]%Np;
-			m10 = Phi[nn];		
-			nn = neighborList[n+11*Np]%Np;
-			m11 = Phi[nn];
-			nn = neighborList[n+10*Np]%Np;
-			m12 = Phi[nn];		
-			nn = neighborList[n+13*Np]%Np;
-			m13 = Phi[nn];
-			nn = neighborList[n+12*Np]%Np;
-			m14 = Phi[nn];		
-			nn = neighborList[n+15*Np]%Np;
-			m15 = Phi[nn];
-			nn = neighborList[n+14*Np]%Np;
-			m16 = Phi[nn];		
-			nn = neighborList[n+17*Np]%Np;
-			m17 = Phi[nn];
-			nn = neighborList[n+16*Np]%Np;
-			m18 = Phi[nn];					
-			
-			//............Compute the Color Gradient...................................
-			nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
-			ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
-			nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
-			// .... read the solid potential gradient.....................
-			m1 = SolidPotential[n];
-			m2 = SolidPotential[n+Np];
-			m3 = SolidPotential[n+2*Np];
-			nx += m1;
-			ny += m2;
-			nz += m3;
-			
-			//...........Normalize the Color Gradient.................................
+			//...........Read the Color Gradient.................................
+			nx = Gradient[n];
+			ny = Gradient[n+Np];
+			bz = Gradient[n+2*Np];
 			C = sqrt(nx*nx+ny*ny+nz*nz);
 			if (C==0.0) C=1.0;
 			nx = nx/C;
@@ -1427,7 +1326,8 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den,
 	}
 }
 
-__global__ void dvc_ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, int start, int finish, int Np){
+__global__ void dvc_ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, double *SolidPotential, int start, int finish, int Np){
+
 	int n,nn;
 	// distributions
 	double m1,m2,m3,m4,m5,m6,m7,m8,m9;
@@ -1508,30 +1408,29 @@ extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq
 }
 
 extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
-		double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 
 	cudaProfilerStart();
 	cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_DFH, cudaFuncCachePreferL1);
 
-	dvc_ScaLBL_D3Q19_AAeven_DFH<<<NBLOCKS,NTHREADS >>>(neighborList, dist, Aq, Bq, Den, Phi, SolidPotential, rhoA, rhoB, tauA, tauB, 
+	dvc_ScaLBL_D3Q19_AAeven_DFH<<<NBLOCKS,NTHREADS >>>(neighborList, dist, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB, 
 			alpha, beta, Fx, Fy, Fz, start, finish, Np);
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
 		printf("CUDA error in ScaLBL_D3Q19_AAeven_DFH: %s \n",cudaGetErrorString(err));
 	}
 	cudaProfilerStop();
-
 }
 
 extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, 
-		double *Phi, double *SolidPotential, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Phi, double *Gradient, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 
 	cudaProfilerStart();
 	cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_DFH, cudaFuncCachePreferL1);
 	
-	dvc_ScaLBL_D3Q19_AAodd_DFH<<<NBLOCKS,NTHREADS >>>(neighborList,dist, Aq, Bq, Den, Phi, SolidPotential, 
+	dvc_ScaLBL_D3Q19_AAodd_DFH<<<NBLOCKS,NTHREADS >>>(neighborList,dist, Aq, Bq, Den, Phi, Gradient, 
 			rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz,  start, finish, Np);
 
 	cudaError_t err = cudaGetLastError();
@@ -1567,9 +1466,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, doub
 
 }
 
-extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *NeighborList, double *Phi, double *ColorGrad, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, double *SolidPotential, int start, int finish, int Np){
 
-	dvc_ScaLBL_D3Q19_Gradient_DFH<<<NBLOCKS,NTHREADS >>>(NeighborList, Phi, ColorGrad, 0, Np, Np);
+	dvc_ScaLBL_D3Q19_Gradient_DFH<<<NBLOCKS,NTHREADS >>>(neighborList, Phi, ColorGrad, SolidPotential, start, finish, Np);
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
 		printf("CUDA error in ScaLBL_D3Q19_Gradient_DFH: %s \n",cudaGetErrorString(err));
diff --git a/tests/lbpm_dfh_simulator.cpp b/tests/lbpm_dfh_simulator.cpp
index 7236aa3c..907836be 100644
--- a/tests/lbpm_dfh_simulator.cpp
+++ b/tests/lbpm_dfh_simulator.cpp
@@ -470,6 +470,7 @@ int main(int argc, char **argv)
 		double *Den, *Phi;
 		double *SolidPotential;
 		double *Velocity;
+		double *Gradient;
 		double *Pressure;
 		
 		//...........................................................................
@@ -482,6 +483,7 @@ int main(int argc, char **argv)
 		ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Np);		
 		ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
 		ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
+		ScaLBL_AllocateDeviceMemory((void **) &Gradient, 3*sizeof(double)*Np);
 		ScaLBL_AllocateDeviceMemory((void **) &SolidPotential, 3*sizeof(double)*Np);
 		
 		//...........................................................................
@@ -605,14 +607,16 @@ int main(int argc, char **argv)
 			ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
 			ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
 			
+			// compute the gradient 
+			ScaLBL_D3Q19_Gradient_DFH(neighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
+			ScaLBL_Comm_Regular.SendHalo(Phi);
+			ScaLBL_D3Q19_Gradient_DFH(neighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.first_interior, Np);
+			ScaLBL_Comm_Regular.RecvGrad(Gradient);
+			
 			// Perform the collision operation
 			ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
-			// Halo exchange for phase field
-			ScaLBL_Comm_Regular.SendHalo(Phi);
-
-			ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, SolidPotential, rhoA, rhoB, tauA, tauB,
+			ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
 					alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
-			ScaLBL_Comm_Regular.RecvHalo(Phi);
 			ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
 			// Set BCs
 			if (BoundaryCondition > 0){
@@ -627,7 +631,7 @@ int main(int argc, char **argv)
 				din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
 				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 			}
-			ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, SolidPotential, rhoA, rhoB, tauA, tauB,
+			ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
 					alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.next, Np);
 			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
@@ -638,14 +642,17 @@ int main(int argc, char **argv)
 			ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
 			ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
 			ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
+			
+			// compute the gradient 
+			ScaLBL_D3Q19_Gradient_DFH(neighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
+			ScaLBL_Comm_Regular.SendHalo(Phi);
+			ScaLBL_D3Q19_Gradient_DFH(neighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.first_interior, Np);
+			ScaLBL_Comm_Regular.RecvGrad(Gradient);
 
 			// Perform the collision operation
 			ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
-			// Halo exchange for phase field
-			ScaLBL_Comm_Regular.SendHalo(Phi);
-			ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, SolidPotential, rhoA, rhoB, tauA, tauB,
+			ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
 					alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
-			ScaLBL_Comm_Regular.RecvHalo(Phi);
 			ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
 			// Set boundary conditions
 			if (BoundaryCondition > 0){
@@ -660,7 +667,7 @@ int main(int argc, char **argv)
 				din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
 				ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 			}
-			ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, SolidPotential, rhoA, rhoB, tauA, tauB,
+			ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
 					alpha, beta, Fx, Fy, Fz,  0, ScaLBL_Comm.next, Np);
 			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 			//************************************************************************