diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index f6837260..1d766ef6 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -80,6 +80,9 @@ extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad,
 extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int Np);
 
 // Density functional hydrodynamics LBM
+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 Fx, double Fy, double Fz, int start, int finish, int Np);
@@ -92,6 +95,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *NeighborList, double *Aq, double *Bq,
 
 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);
+
 // BOUNDARY CONDITION ROUTINES
 
 //extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, double din,
diff --git a/cpu/dfh.cpp b/cpu/dfh.cpp
index 6543de3b..acc139da 100644
--- a/cpu/dfh.cpp
+++ b/cpu/dfh.cpp
@@ -1,6 +1,40 @@
 #include <math.h>
 #include <stdio.h>
 
+extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){
+	int idx,n;
+	double phi,nA,nB;
+
+	for (idx=start; idx<finish; idx++){
+		phi = Phi[idx];
+		if (phi > 0.f){
+			nA = 1.0; nB = 0.f;
+		}
+		else{
+			nB = 1.0; nA = 0.f;
+		}
+		Den[idx] = nA;
+		Den[Np+idx] = nB;
+
+		Aq[idx]=0.3333333333333333*nA;
+		Aq[Np+idx]=0.1111111111111111*nA;
+		Aq[2*Np+idx]=0.1111111111111111*nA;
+		Aq[3*Np+idx]=0.1111111111111111*nA;
+		Aq[4*Np+idx]=0.1111111111111111*nA;
+		Aq[5*Np+idx]=0.1111111111111111*nA;
+		Aq[6*Np+idx]=0.1111111111111111*nA;
+
+		Bq[idx]=0.3333333333333333*nB;
+		Bq[Np+idx]=0.1111111111111111*nB;
+		Bq[2*Np+idx]=0.1111111111111111*nB;
+		Bq[3*Np+idx]=0.1111111111111111*nB;
+		Bq[4*Np+idx]=0.1111111111111111*nB;
+		Bq[5*Np+idx]=0.1111111111111111*nB;
+		Bq[6*Np+idx]=0.1111111111111111*nB;
+	}
+}
+
+
 // 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,
@@ -1378,3 +1412,71 @@ extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, doub
 		Phi[n] = (nA-nB)/(nA+nB); 	
 	}	
 }
+
+extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, int start, int finish, int Np){
+	
+	int n,nn;
+	// distributions
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double m3,m5,m7;
+	double nx,ny,nz;
+
+	// non-conserved moments
+	// additional variables needed for computations
+
+	for (n=start; n<finish; n++){
+		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));
+		
+		//...........Normalize the Color Gradient.................................
+		//	C = sqrt(nx*nx+ny*ny+nz*nz);
+		//	nx = nx/C;
+		//	ny = ny/C;
+		//	nz = nz/C;
+		//...Store the Color Gradient....................
+		ColorGrad[n] = nx;
+		ColorGrad[Np+n] = ny;
+		ColorGrad[2*Np+n] = nz;
+		//...............................................
+	}
+}
diff --git a/gpu/Extras.cu b/gpu/Extras.cu
index f84aac51..064f89ab 100644
--- a/gpu/Extras.cu
+++ b/gpu/Extras.cu
@@ -13,10 +13,10 @@ extern "C" int ScaLBL_SetDevice(int rank){
 }
 
 extern "C" void ScaLBL_AllocateDeviceMemory(void** address, size_t size){
-       	cudaMalloc(address,size);
+	cudaMalloc(address,size);
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
-	   printf("Error in cudaMalloc: %s \n",cudaGetErrorString(err));
+		printf("Error in cudaMalloc: %s \n",cudaGetErrorString(err));
 	}	
 }
 
@@ -33,11 +33,12 @@ extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size)
 }
 
 extern "C" void ScaLBL_AllocateZeroCopy(void** address, size_t size){
-    cudaMallocHost(address,size);
-        cudaError_t err = cudaGetLastError();
-        if (cudaSuccess != err){
-           printf("Error in cudaMallocHost: %s \n",cudaGetErrorString(err));
-        }
+	//cudaMallocHost(address,size);
+	cudaMalloc(address,size);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("Error in cudaMallocHost: %s \n",cudaGetErrorString(err));
+	}
 }
 
 extern "C" void ScaLBL_CopyToZeroCopy(void* dest, const void* source, size_t size){
diff --git a/gpu/dfh.cu b/gpu/dfh.cu
index b9d66e4b..b932db7f 100644
--- a/gpu/dfh.cu
+++ b/gpu/dfh.cu
@@ -5,6 +5,46 @@
 #define NBLOCKS 1024
 #define NTHREADS 256
 
+__global__ void dvc_ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){
+	int idx,n;
+	double phi,nA,nB;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		idx =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x+start;
+		if (idx<finish) {
+
+			phi = Phi[idx];
+			if (phi > 0.f){
+				nA = 1.0; nB = 0.f;
+			}
+			else{
+				nB = 1.0; nA = 0.f;
+			}
+			Den[idx] = nA;
+			Den[Np+idx] = nB;
+
+			Aq[idx]=0.3333333333333333*nA;
+			Aq[Np+idx]=0.1111111111111111*nA;
+			Aq[2*Np+idx]=0.1111111111111111*nA;
+			Aq[3*Np+idx]=0.1111111111111111*nA;
+			Aq[4*Np+idx]=0.1111111111111111*nA;
+			Aq[5*Np+idx]=0.1111111111111111*nA;
+			Aq[6*Np+idx]=0.1111111111111111*nA;
+
+			Bq[idx]=0.3333333333333333*nB;
+			Bq[Np+idx]=0.1111111111111111*nB;
+			Bq[2*Np+idx]=0.1111111111111111*nB;
+			Bq[3*Np+idx]=0.1111111111111111*nB;
+			Bq[4*Np+idx]=0.1111111111111111*nB;
+			Bq[5*Np+idx]=0.1111111111111111*nB;
+			Bq[6*Np+idx]=0.1111111111111111*nB;
+		}
+	}
+}
+
+
 // 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,
@@ -1359,6 +1399,75 @@ __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){
+	int n,nn;
+	// distributions
+	double m1,m2,m3,m4,m5,m6,m7,m8,m9;
+	double m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double nx,ny,nz;
+
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+		if (n<finish) {
+			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];					
+			
+			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));
+			//...............................................
+			//...Store the Color Gradient....................
+			ColorGrad[n] = nx;
+			ColorGrad[Np+n] = ny;
+			ColorGrad[2*Np+n] = nz;
+			//...............................................
+		}
+	}
+}
+extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){
+	dvc_ScaLBL_DFH_Init<<<NBLOCKS,NTHREADS >>>(Phi, Den, Aq, Bq, start, finish, Np); 
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_DFH_Init: %s \n",cudaGetErrorString(err));
+	}
+}
+
 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 Fx, double Fy, double Fz, int start, int finish, int Np){
@@ -1419,3 +1528,14 @@ 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){
+
+	int strideY=Nx;
+	int strideZ=Nx*Ny;
+	dvc_ScaLBL_D3Q19_Gradient_DFH<<<NBLOCKS,NTHREADS >>>(NeighborList, Phi, ColorGrad, 0, Np, 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 d1d7a682..4f875a52 100644
--- a/tests/lbpm_dfh_simulator.cpp
+++ b/tests/lbpm_dfh_simulator.cpp
@@ -81,6 +81,7 @@ int main(int argc, char **argv)
 		double tauA, tauB, rhoA,rhoB;
 		double Fx,Fy,Fz,tol,err;
 		double alpha, beta;
+		double bns,bws,cns,cws;
 		int BoundaryCondition;
 		int InitialCondition;
 		//	bool pBC,Restart;
@@ -112,6 +113,10 @@ int main(int argc, char **argv)
 				input >> rhoB;			// density wetting
 				input >> alpha;			// Surface Tension parameter
 				input >> beta;			// Width of the interface
+				input >> cws;			// solid interaction coefficients
+				input >> bws;			// solid interaction coefficients
+				input >> cns;			// solid interaction coefficients
+				input >> bns;			// solid interaction coefficients
 				// Line 2:  External force components (Fx,Fy, Fz)
 				input >> Fx;
 				input >> Fy;
@@ -179,6 +184,10 @@ int main(int argc, char **argv)
 		MPI_Bcast(&rhoB,1,MPI_DOUBLE,0,comm);
 		MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
 		MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&cns,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&cws,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&bns,1,MPI_DOUBLE,0,comm);
+		MPI_Bcast(&bws,1,MPI_DOUBLE,0,comm);
 		MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
 		MPI_Bcast(&InitialCondition,1,MPI_INT,0,comm);
 		MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
@@ -496,8 +505,8 @@ int main(int argc, char **argv)
 		
 		// Compute the solid interaction potential and copy result to device
 		if (rank==0) printf("Computing solid interaction potential \n");
-		int *Tmp;
-		Tmp=new int[3*Np];
+		double *Tmp;
+		Tmp=new double[3*Np];
 		Averages->UpdateMeshValues(); // this computes the gradient of distance field (among other things)
 		for (k=1; k<Nz-1; k++){
 			for (j=1; j<Ny-1; j++){
@@ -509,6 +518,13 @@ int main(int argc, char **argv)
 						double dx = Averages->SDs_x(n);
 						double dy = Averages->SDs_y(n);
 						double dz = Averages->SDs_z(n);
+						double value=cns*exp(-bns*fabs(d))-cws*exp(-bns*fabs(d));
+						Tmp[idx] = value*dx;
+						Tmp[idx+Np] = value*dy;
+						Tmp[idx+2*Np] = value*dz;
+						// initialize fluid phases
+						if (Dm.id[n] == 1)	PhaseLabel[idx] = 1.0;
+						else 				PhaseLabel[idx] = -1.0;
 					}
 				}
 			}
@@ -520,26 +536,13 @@ int main(int argc, char **argv)
 		// copy the neighbor list 
 		ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
 		// initialize phi based on PhaseLabel (include solid component labels)
-		ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double));
-		//ScaLBL_CopyToDevice(Distance, PhaseLabel, N*sizeof(double));
+		ScaLBL_CopyToDevice(Phi, PhaseLabel, Np*sizeof(double));
 		//...........................................................................
 
 		if (rank==0)	printf ("Initializing distributions \n");
 		ScaLBL_D3Q19_Init(fq, Np);
 		if (rank==0)	printf ("Initializing phase field \n");
-		ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, Np);
-		if (BoundaryCondition >0 ){
-			if (Dm.kproc==0){
-				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
-				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
-			}
-			if (Dm.kproc == nprocz-1){
-				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
-				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
-			}
-		}
+		ScaLBL_DFH_Init(Phi, Den, Aq, Bq, Np);
 
 		//.......................................................................
 		// Once phase has been initialized, map solid to account for 'smeared' interface
@@ -691,6 +694,15 @@ int main(int argc, char **argv)
 
 		// ************************************************************************
     	
+		// Copy back final phase indicator field and convert to regular layout
+		DoubleArray PhaseField(Nx,Ny,Nz);
+        ScaLBL_Comm.RegularLayout(Map,Phi,PhaseField);
+    	FILE *OUTFILE;
+		sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
+		OUTFILE = fopen(LocalRankFilename,"wb");
+    	fwrite(PhaseField.data(),8,N,OUTFILE);
+    	fclose(OUTFILE);
+    	
 		PROFILE_STOP("Main");
 		PROFILE_SAVE("lbpm_color_simulator",1);
 		// ****************************************************
diff --git a/tests/lbpm_porenetwork_pp.cpp b/tests/lbpm_porenetwork_pp.cpp
index c0f705fe..6e9de20f 100644
--- a/tests/lbpm_porenetwork_pp.cpp
+++ b/tests/lbpm_porenetwork_pp.cpp
@@ -208,12 +208,12 @@ int main(int argc, char **argv)
 					double zk = double(k) - z;
 					// value of s along center line {x=alpha*s, y = beta*s, z=gamma*s} that is closest to i,j,k
 					double s = (alpha*xi + beta*yj + gamma*zk)/(alpha*alpha + beta*beta + gamma*gamma);
-					double distance=0.f;
+					double distance=Averages.SDs(i,j,k);
 					if (s > length){
-						distance = radius - sqrt((xi-X)*(xi-X) + (yj-Y)*(yj-Y) + (zk-Z)*(zk-Z));
+						//distance = radius - sqrt((xi-X)*(xi-X) + (yj-Y)*(yj-Y) + (zk-Z)*(zk-Z));
 					}
 					else if (s<0){
-						distance = radius - sqrt((xi-x)*(xi-x) + (yj-y)*(yj-y) + (zk-z)*(zk-z));
+						//distance = radius - sqrt((xi-x)*(xi-x) + (yj-y)*(yj-y) + (zk-z)*(zk-z));
 					}
 					else{
 						double xs = alpha*s;