#include <math.h>
#include <stdio.h>

#define STOKES

extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np)
{
	int n;
    double p = 1.0;//NOTE: take initial pressure p=1.0
    double chem;
    double cg_x,cg_y,cg_z;

	for (n=0; n<Np; n++){
        chem = mu_phi[n];//chemical potential
        cg_x = ColorGrad[0*Np+n];
        cg_y = ColorGrad[1*Np+n];
        cg_z = ColorGrad[2*Np+n];
        
		gqbar[0*Np+n]  = 0.3333333333333333*p;
		gqbar[1*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_x+Fx));		//double(100*n)+1.f;
		gqbar[2*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_x-Fx));	//double(100*n)+2.f;
		gqbar[3*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_y+Fy));	//double(100*n)+3.f;
		gqbar[4*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_y-Fy));	//double(100*n)+4.f;
		gqbar[5*Np+n]  = 0.055555555555555555*(p - 0.5*(chem*cg_z+Fz));	//double(100*n)+5.f;
		gqbar[6*Np+n]  = 0.055555555555555555*(p - 0.5*(-chem*cg_z-Fz));	//double(100*n)+6.f;

		gqbar[7*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_y)+Fx+Fy));   //double(100*n)+7.f;
		gqbar[8*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_y)-Fx-Fy));   //double(100*n)+8.f;
		gqbar[9*Np+n]  = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_y)+Fx-Fy));   //double(100*n)+9.f;
		gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_y)-Fx+Fy));  //double(100*n)+10.f;

		gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_z)+Fx+Fz));  //double(100*n)+11.f;
		gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_z)-Fx-Fz));  //double(100*n)+12.f;
		gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_z)+Fx-Fz));  //double(100*n)+13.f;
		gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_z)-Fx+Fz));  //double(100*n)+14.f;
		
        gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y+cg_z)+Fy+Fz)); ;  //double(100*n)+15.f;
		gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y-cg_z)-Fy-Fz));;  //double(100*n)+16.f;
		gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y-cg_z)+Fy-Fz)); ;  //double(100*n)+17.f;
		gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y+cg_z)-Fy+Fz));;  //double(100*n)+18.f;
	}
}

extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np)
{
	int n;
    double p = 1.0;//NOTE: take initial pressure p=1.0

	for (n=0; n<Np; n++){
        
		gqbar[0*Np+n]  = 0.3333333333333333*p;
		gqbar[1*Np+n]  = 0.055555555555555555*(p - 0.5*(Fx));		//double(100*n)+1.f;
		gqbar[2*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fx));	//double(100*n)+2.f;
		gqbar[3*Np+n]  = 0.055555555555555555*(p - 0.5*(Fy));	//double(100*n)+3.f;
		gqbar[4*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fy));	//double(100*n)+4.f;
		gqbar[5*Np+n]  = 0.055555555555555555*(p - 0.5*(Fz));	//double(100*n)+5.f;
		gqbar[6*Np+n]  = 0.055555555555555555*(p - 0.5*(-Fz));	//double(100*n)+6.f;

		gqbar[7*Np+n]  = 0.0277777777777778*(p-0.5*(Fx+Fy));   //double(100*n)+7.f;
		gqbar[8*Np+n]  = 0.0277777777777778*(p-0.5*(-Fx-Fy));   //double(100*n)+8.f;
		gqbar[9*Np+n]  = 0.0277777777777778*(p-0.5*(Fx-Fy));   //double(100*n)+9.f;
		gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fy));  //double(100*n)+10.f;

		gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(Fx+Fz));  //double(100*n)+11.f;
		gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(-Fx-Fz));  //double(100*n)+12.f;
		gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(Fx-Fz));  //double(100*n)+13.f;
		gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fz));  //double(100*n)+14.f;
		
        gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(Fy+Fz)); ;  //double(100*n)+15.f;
		gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(-Fy-Fz));;  //double(100*n)+16.f;
		gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(Fy-Fz)); ;  //double(100*n)+17.f;
		gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(-Fy+Fz));;  //double(100*n)+18.f;
	}
	
	
}

extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, 
                                                    double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
	int idx,n;
	double phi;
    double nx,ny,nz,cg_mag;
    double theta;//anti-diffusion term
    double cs2_inv = 4.5;//inverse of speed of sound for D3Q7
    double M = 1.0/cs2_inv*(tauM-0.5);//diffusivity (or mobility)
    double factor = 1.0;

	for (idx=start; idx<finish; idx++){

		n = Map[idx];
		phi = Phi[n];
		if (phi > 1.f)   phi = 1.0;
		if (phi < -1.f)  phi = -1.0;
		Den[idx] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);

        //compute unit normal of color gradient
        nx = ColorGrad[idx+0*Np];
        ny = ColorGrad[idx+1*Np];
        nz = ColorGrad[idx+2*Np];
        cg_mag = sqrt(nx*nx+ny*ny+nz*nz);
		double ColorMag_temp = cg_mag;
		if (cg_mag==0.0) ColorMag_temp=1.0;
		nx = nx/ColorMag_temp;
		ny = ny/ColorMag_temp;
		nz = nz/ColorMag_temp;		

        //theta = M*cs2_inv*(1-factor*phi*phi)/W;
        theta = 4.5*M*2.0*(1-phi*phi)/W;
        //theta = 0; // try more diffusive initial condition
        
		hq[0*Np+idx]=0.3333333333333333*(phi);
		hq[1*Np+idx]=0.1111111111111111*(phi+theta*nx);
		hq[2*Np+idx]=0.1111111111111111*(phi-theta*nx);
		hq[3*Np+idx]=0.1111111111111111*(phi+theta*ny);
		hq[4*Np+idx]=0.1111111111111111*(phi-theta*ny);
		hq[5*Np+idx]=0.1111111111111111*(phi+theta*nz);
		hq[6*Np+idx]=0.1111111111111111*(phi-theta*nz);
		
	}
}

extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, 
                                                          double rhoA, double rhoB, int start, int finish, int Np){

	int idx,nread;
	double fq,phi;

	for (int n=start; n<finish; n++){

		// q=0
		fq = hq[n];
		phi = fq;

		// q=1
		nread = neighborList[n]; 
		fq = hq[nread];
		phi += fq;
		
		// q=2
		nread = neighborList[n+Np]; 
		fq = hq[nread];  
		phi += fq;

		// q=3
		nread = neighborList[n+2*Np]; 
		fq = hq[nread];
		phi += fq;

		// q = 4
		nread = neighborList[n+3*Np]; 
		fq = hq[nread];
		phi += fq;

		// q=5
		nread = neighborList[n+4*Np];
		fq = hq[nread];
		phi += fq;

		// q = 6
		nread = neighborList[n+5*Np];
		fq = hq[nread];
		phi += fq;
		
		// save the number densities
		Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
		
		// save the phase indicator field
		idx = Map[n];
		Phi[idx] = phi; 
		
	}
}



extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, 
			                                               double rhoA, double rhoB, int start, int finish, int Np){
	int idx;
	double fq,phi;
	for (int n=start; n<finish; n++){
		
		// q=0
		fq = hq[n];
		phi = fq;
		
		// q=1
		fq = hq[2*Np+n];
		phi += fq;

		// f2 = hq[10*Np+n];
		fq = hq[1*Np+n];
		phi += fq;

		// q=3
		fq = hq[4*Np+n];
		phi += fq;

		// q = 4
		fq = hq[3*Np+n];
		phi += fq;

		// q=5
		fq = hq[6*Np+n];
		phi += fq;

		// q = 6
		fq = hq[5*Np+n];
		phi += fq;

		// save the number densities
		Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
		
		// save the phase indicator field
		idx = Map[n];
		Phi[idx] = phi; 	
	}	
}

extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
                                                          double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){

	int idx,nr1,nr2,nr3,nr4,nr5,nr6;
	double h0,h1,h2,h3,h4,h5,h6;
	double nx,ny,nz,C;
	double ux,uy,uz;
	double phi;
    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
    double factor = 1.0;
    double theta;

	for (int n=start; n<finish; n++){

		/* load phase indicator field */
		idx = Map[n];
		phi = Phi[idx];
        theta = 4.5*M*2.0*(1-phi*phi)/W;

		/* velocity */
		ux = Vel[0*Np+n];
		uy = Vel[1*Np+n];
		uz = Vel[2*Np+n];

        /*color gradient */
		nx = ColorGrad[0*Np+n];
		ny = ColorGrad[1*Np+n];
		nz = ColorGrad[2*Np+n];
		
		//Normalize the Color Gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		double ColorMag = C;
		if (C < 1.0e-12) ColorMag=1.0;
		nx = nx/ColorMag;
		ny = ny/ColorMag;
		nz = nz/ColorMag;		

		// q=1
		nr1 = neighborList[n]; 
		nr2 = neighborList[n+Np]; 
		nr3 = neighborList[n+2*Np]; 
		nr4 = neighborList[n+3*Np];
		nr5 = neighborList[n+4*Np];
		nr6 = neighborList[n+5*Np];
		
		//q=0
		h0 = hq[n];
		//q=1
		h1 = hq[nr1]; 
		//q=2
		h2 = hq[nr2];  
		//q=3
		h3 = hq[nr3];
		//q=4
		h4 = hq[nr4];
		//q=5
		h5 = hq[nr5];
		//q=6
		h6 = hq[nr6];

        //-------------------------------- BGK collison for phase field ---------------------------------//
		h0 -= (h0 - 0.3333333333333333*phi)/tauM;
//		h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
//		h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
//		h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
//		h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
//		h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//		h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
		h1 -= (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
		h2 -= (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
		h3 -= (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
		h4 -= (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
		h5 -= (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
		h6 -= (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM; 
		//........................................................................
		
		/*Update the distributions */
		// q = 0
		hq[n] = h0;
		hq[nr2] = h1;
		hq[nr1] = h2;
		hq[nr4] = h3;
		hq[nr3] = h4;
		hq[nr6] = h5;
		hq[nr5] = h6;
		//........................................................................

	}
}



extern "C" void ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
		double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){

	int idx,n;
	double h0,h1,h2,h3,h4,h5,h6;
	double nx,ny,nz,C;
	double ux,uy,uz;
	double phi;
    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
    double factor = 1.0;
    double theta;

	for (int n=start; n<finish; n++){

		/* load phase indicator field */
		idx = Map[n];
		phi = Phi[idx];
        theta = 4.5*M*2.0*(1-phi*phi)/W;

		/* velocity */
		ux = Vel[0*Np+n];
		uy = Vel[1*Np+n];
		uz = Vel[2*Np+n];

		/*color gradient */
		nx = ColorGrad[0*Np+n];
		ny = ColorGrad[1*Np+n];
		nz = ColorGrad[2*Np+n];
		//Normalize the Color Gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		double ColorMag = C;
		if (C < 1.0e-12) ColorMag=1.0;
		nx = nx/ColorMag;
		ny = ny/ColorMag;
		nz = nz/ColorMag;

		h0 = hq[n];
		h1 = hq[2*Np+n]; 
		h2 = hq[Np+n];  
		h3 = hq[4*Np+n];
		h4 = hq[3*Np+n];
		h5 = hq[6*Np+n];
		h6 = hq[5*Np+n];

		//-------------------------------- BGK collison for phase field ---------------------------------//
		h0 -= (h0 - 0.3333333333333333*phi)/tauM;
//		h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
//		h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
//		h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
//		h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
//		h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//		h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
		h1 -= (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
		h2 -= (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
		h3 -= (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
		h4 -= (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
		h5 -= (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
		h6 -= (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM; 
		//........................................................................

		/*Update the distributions */
		// q = 0
		hq[n] = h0;
		hq[Np+n] = h1;
		hq[2*Np+n] = h2;
		hq[3*Np+n] = h3;
		hq[4*Np+n] = h4;
		hq[5*Np+n] = h5;
		hq[6*Np+n] = h6;
		//........................................................................
		
	}
}

extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map,  double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
	int idx,n;
	double h0,h1,h2,h3,h4,h5,h6;
	double phi;

	for (int n=start; n<finish; n++){

		h0 = hq[n];
		h1 = hq[1*Np+n]; 
		h2 = hq[2*Np+n];  
		h3 = hq[3*Np+n];
		h4 = hq[4*Np+n];
		h5 = hq[5*Np+n];
		h6 = hq[6*Np+n];

		phi = h0+h1+h2+h3+h4+h5+h6;

		// save the number densities
		Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);

		// save the phase indicator field
		idx = Map[n];
		Phi[idx] = phi; 

	}
}


extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, 
                                                double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
                                                int strideY, int strideZ, int start, int finish, int Np){
	
	int nn,nn2x,ijk;
	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
    double chem;//chemical potential
    double phi; //phase field
    double rho0;//fluid density
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;
	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
	double mm3,mm5,mm7;
    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
    double nx,ny,nz;//normal color gradient
    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor

	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
    //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
	double tau;//position dependent LB relaxation time for fluid
    //double C,theta;
    // double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7

	for (int n=start; n<finish; n++){

		rho0 = Den[n];//load density

		// Get the 1D index based on regular data layout
		ijk = Map[n];
        phi = Phi[ijk];// load phase field

		// local relaxation time
		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);

		//					COMPUTE THE COLOR GRADIENT
		//........................................................................
		//.................Read Phase Indicator Values............................
		//........................................................................
		nn = ijk-1;							// neighbor index (get convention)
		m1 = Phi[nn];						// get neighbor for phi - 1
		//........................................................................
		nn = ijk+1;							// neighbor index (get convention)
		m2 = Phi[nn];						// get neighbor for phi - 2
		//........................................................................
		nn = ijk-strideY;							// neighbor index (get convention)
		m3 = Phi[nn];					// get neighbor for phi - 3
		//........................................................................
		nn = ijk+strideY;							// neighbor index (get convention)
		m4 = Phi[nn];					// get neighbor for phi - 4
		//........................................................................
		nn = ijk-strideZ;						// neighbor index (get convention)
		m5 = Phi[nn];					// get neighbor for phi - 5
		//........................................................................
		nn = ijk+strideZ;						// neighbor index (get convention)
		m6 = Phi[nn];					// get neighbor for phi - 6
		//........................................................................
		nn = ijk-strideY-1;						// neighbor index (get convention)
		m7 = Phi[nn];					// get neighbor for phi - 7
		//........................................................................
		nn = ijk+strideY+1;						// neighbor index (get convention)
		m8 = Phi[nn];					// get neighbor for phi - 8
		//........................................................................
		nn = ijk+strideY-1;						// neighbor index (get convention)
		m9 = Phi[nn];					// get neighbor for phi - 9
		//........................................................................
		nn = ijk-strideY+1;						// neighbor index (get convention)
		m10 = Phi[nn];					// get neighbor for phi - 10
		//........................................................................
		nn = ijk-strideZ-1;						// neighbor index (get convention)
		m11 = Phi[nn];					// get neighbor for phi - 11
		//........................................................................
		nn = ijk+strideZ+1;						// neighbor index (get convention)
		m12 = Phi[nn];					// get neighbor for phi - 12
		//........................................................................
		nn = ijk+strideZ-1;						// neighbor index (get convention)
		m13 = Phi[nn];					// get neighbor for phi - 13
		//........................................................................
		nn = ijk-strideZ+1;						// neighbor index (get convention)
		m14 = Phi[nn];					// get neighbor for phi - 14
		//........................................................................
		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
		m15 = Phi[nn];					// get neighbor for phi - 15
		//........................................................................
		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
		m16 = Phi[nn];					// get neighbor for phi - 16
		//........................................................................
		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
		m17 = Phi[nn];					// get neighbor for phi - 17
		//........................................................................
		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
		m18 = Phi[nn];					// get neighbor for phi - 18

        // compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
		//........................................................................
		nn2x = ijk-2;							// neighbor index (get convention)
		mm1 = Phi[nn2x];						// get neighbor for phi - 1
        mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
		//........................................................................
		nn2x = ijk+2;							// neighbor index (get convention)
		mm2 = Phi[nn2x];						// get neighbor for phi - 2
        mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
		//........................................................................
		nn2x = ijk-strideY*2;							// neighbor index (get convention)
		mm3 = Phi[nn2x];					// get neighbor for phi - 3
        mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
		//........................................................................
		nn2x = ijk+strideY*2;							// neighbor index (get convention)
		mm4 = Phi[nn2x];					// get neighbor for phi - 4
        mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
		//........................................................................
		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
		mm5 = Phi[nn2x];					// get neighbor for phi - 5
        mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
		//........................................................................
		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
		mm6 = Phi[nn2x];					// get neighbor for phi - 6
        mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
		//........................................................................
		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
		mm7 = Phi[nn2x];					// get neighbor for phi - 7
        mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
		//........................................................................
		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
		mm8 = Phi[nn2x];					// get neighbor for phi - 8
        mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
		//........................................................................
		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
		mm9 = Phi[nn2x];					// get neighbor for phi - 9
        mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
		//........................................................................
		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
		mm10 = Phi[nn2x];					// get neighbor for phi - 10
        mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
		//........................................................................
		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
		mm11 = Phi[nn2x];					// get neighbor for phi - 11
        mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
		//........................................................................
		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
		mm12 = Phi[nn2x];					// get neighbor for phi - 12
        mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
		//........................................................................
		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
		mm13 = Phi[nn2x];					// get neighbor for phi - 13
        mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
		//........................................................................
		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
		mm14 = Phi[nn2x];					// get neighbor for phi - 14
        mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
		//........................................................................
		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
		mm15 = Phi[nn2x];					// get neighbor for phi - 15
        mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
		//........................................................................
		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
		mm16 = Phi[nn2x];					// get neighbor for phi - 16
        mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
		//........................................................................
		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
		mm17 = Phi[nn2x];					// get neighbor for phi - 17
        mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
		//........................................................................
		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
		mm18 = Phi[nn2x];					// get neighbor for phi - 18
        mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);


		//............Compute the Color Gradient...................................
		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
		//............Compute the Chemical Potential...............................
        chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
        chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
		//............Compute the Mixed Gradient...................................
		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
		
		// q=0
		m0 = dist[n];
		// q=1
		nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
		m1 = dist[nr1]; // reading the f1 data into register fq

		nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
		m2 = dist[nr2];  // reading the f2 data into register fq

		// q=3
		nr3 = neighborList[n+2*Np]; // neighbor 4
		m3 = dist[nr3];

		// q = 4
		nr4 = neighborList[n+3*Np]; // neighbor 3
		m4 = dist[nr4];

		// q=5
		nr5 = neighborList[n+4*Np];
		m5 = dist[nr5];

		// q = 6
		nr6 = neighborList[n+5*Np];
		m6 = dist[nr6];
		
		// q=7
		nr7 = neighborList[n+6*Np];
		m7 = dist[nr7];

		// q = 8
		nr8 = neighborList[n+7*Np];
		m8 = dist[nr8];

		// q=9
		nr9 = neighborList[n+8*Np];
		m9 = dist[nr9];

		// q = 10
		nr10 = neighborList[n+9*Np];
		m10 = dist[nr10];

		// q=11
		nr11 = neighborList[n+10*Np];
		m11 = dist[nr11];

		// q=12
		nr12 = neighborList[n+11*Np];
		m12 = dist[nr12];

		// q=13
		nr13 = neighborList[n+12*Np];
		m13 = dist[nr13];

		// q=14
		nr14 = neighborList[n+13*Np];
		m14 = dist[nr14];

		// q=15
		nr15 = neighborList[n+14*Np];
		m15 = dist[nr15];

		// q=16
		nr16 = neighborList[n+15*Np];
		m16 = dist[nr16];

		// q=17
		nr17 = neighborList[n+16*Np];
		m17 = dist[nr17];

		// q=18
		nr18 = neighborList[n+17*Np];
		m18 = dist[nr18];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
                  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);

        //compute equilibrium distributions
            feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
     0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
      (-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
            feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
      (2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
            feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
         chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
            feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
            feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
         chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
            feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
            feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
     0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
            feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
         2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
      (2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
            feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
      (2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
            feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
            feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
     (mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
        0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz)))); 

		// q = 1
		dist[nr2] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));

		// q=2
		dist[nr1] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));

		// q = 3
		dist[nr4] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 4
		dist[nr3] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 5
		dist[nr6] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 6
		dist[nr5] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 7
		dist[nr8] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (0.2222222222222222 + (ux + uy)*(ux + uy) - 
        0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 8
		dist[nr7] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uy)*(ux + uy) + 
        0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 9
		dist[nr10] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 10
		dist[nr9] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 11
		dist[nr12] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (0.2222222222222222 + (ux + uz)*(ux + uz) - 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 12
		dist[nr11] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
       + (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 13
		dist[nr14] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q= 14
		dist[nr13] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 15
		dist[nr16] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
      (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
       + (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));

		// q = 16
		dist[nr15] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
       + (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 17
		dist[nr18] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
       + (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 18
		dist[nr17] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
        //----------------------------------------------------------------------------------------------------------------------------------------//

        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
        //Update chemical potential on device
        mu_phi[n] = chem;
        //Update color gradient on device
		ColorGrad[0*Np+n] = nx;
		ColorGrad[1*Np+n] = ny;
		ColorGrad[2*Np+n] = nz;

	}
}

extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
                                                double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
                                                int strideY, int strideZ, int start, int finish, int Np){
	
	int nn,nn2x,ijk;
	//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
    double chem;//chemical potential
    double phi; //phase field
    double rho0;//fluid density
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;
	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
	double mm3,mm5,mm7;
    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
    double nx,ny,nz;//normal color gradient
    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor

	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
    //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
	double tau;//position dependent LB relaxation time for fluid
    //double C,theta;
    //double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7

	for (int n=start; n<finish; n++){

		rho0 = Den[n];//load density

		// Get the 1D index based on regular data layout
		ijk = Map[n];
        phi = Phi[ijk];// load phase field

		// local relaxation time
		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);

		//					COMPUTE THE COLOR GRADIENT
		//........................................................................
		//.................Read Phase Indicator Values............................
		//........................................................................
		nn = ijk-1;							// neighbor index (get convention)
		m1 = Phi[nn];						// get neighbor for phi - 1
		//........................................................................
		nn = ijk+1;							// neighbor index (get convention)
		m2 = Phi[nn];						// get neighbor for phi - 2
		//........................................................................
		nn = ijk-strideY;							// neighbor index (get convention)
		m3 = Phi[nn];					// get neighbor for phi - 3
		//........................................................................
		nn = ijk+strideY;							// neighbor index (get convention)
		m4 = Phi[nn];					// get neighbor for phi - 4
		//........................................................................
		nn = ijk-strideZ;						// neighbor index (get convention)
		m5 = Phi[nn];					// get neighbor for phi - 5
		//........................................................................
		nn = ijk+strideZ;						// neighbor index (get convention)
		m6 = Phi[nn];					// get neighbor for phi - 6
		//........................................................................
		nn = ijk-strideY-1;						// neighbor index (get convention)
		m7 = Phi[nn];					// get neighbor for phi - 7
		//........................................................................
		nn = ijk+strideY+1;						// neighbor index (get convention)
		m8 = Phi[nn];					// get neighbor for phi - 8
		//........................................................................
		nn = ijk+strideY-1;						// neighbor index (get convention)
		m9 = Phi[nn];					// get neighbor for phi - 9
		//........................................................................
		nn = ijk-strideY+1;						// neighbor index (get convention)
		m10 = Phi[nn];					// get neighbor for phi - 10
		//........................................................................
		nn = ijk-strideZ-1;						// neighbor index (get convention)
		m11 = Phi[nn];					// get neighbor for phi - 11
		//........................................................................
		nn = ijk+strideZ+1;						// neighbor index (get convention)
		m12 = Phi[nn];					// get neighbor for phi - 12
		//........................................................................
		nn = ijk+strideZ-1;						// neighbor index (get convention)
		m13 = Phi[nn];					// get neighbor for phi - 13
		//........................................................................
		nn = ijk-strideZ+1;						// neighbor index (get convention)
		m14 = Phi[nn];					// get neighbor for phi - 14
		//........................................................................
		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
		m15 = Phi[nn];					// get neighbor for phi - 15
		//........................................................................
		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
		m16 = Phi[nn];					// get neighbor for phi - 16
		//........................................................................
		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
		m17 = Phi[nn];					// get neighbor for phi - 17
		//........................................................................
		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
		m18 = Phi[nn];					// get neighbor for phi - 18

        // compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
		//........................................................................
		nn2x = ijk-2;							// neighbor index (get convention)
		mm1 = Phi[nn2x];						// get neighbor for phi - 1
        mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
		//........................................................................
		nn2x = ijk+2;							// neighbor index (get convention)
		mm2 = Phi[nn2x];						// get neighbor for phi - 2
        mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
		//........................................................................
		nn2x = ijk-strideY*2;							// neighbor index (get convention)
		mm3 = Phi[nn2x];					// get neighbor for phi - 3
        mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
		//........................................................................
		nn2x = ijk+strideY*2;							// neighbor index (get convention)
		mm4 = Phi[nn2x];					// get neighbor for phi - 4
        mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
		//........................................................................
		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
		mm5 = Phi[nn2x];					// get neighbor for phi - 5
        mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
		//........................................................................
		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
		mm6 = Phi[nn2x];					// get neighbor for phi - 6
        mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
		//........................................................................
		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
		mm7 = Phi[nn2x];					// get neighbor for phi - 7
        mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
		//........................................................................
		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
		mm8 = Phi[nn2x];					// get neighbor for phi - 8
        mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
		//........................................................................
		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
		mm9 = Phi[nn2x];					// get neighbor for phi - 9
        mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
		//........................................................................
		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
		mm10 = Phi[nn2x];					// get neighbor for phi - 10
        mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
		//........................................................................
		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
		mm11 = Phi[nn2x];					// get neighbor for phi - 11
        mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
		//........................................................................
		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
		mm12 = Phi[nn2x];					// get neighbor for phi - 12
        mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
		//........................................................................
		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
		mm13 = Phi[nn2x];					// get neighbor for phi - 13
        mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
		//........................................................................
		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
		mm14 = Phi[nn2x];					// get neighbor for phi - 14
        mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
		//........................................................................
		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
		mm15 = Phi[nn2x];					// get neighbor for phi - 15
        mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
		//........................................................................
		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
		mm16 = Phi[nn2x];					// get neighbor for phi - 16
        mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
		//........................................................................
		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
		mm17 = Phi[nn2x];					// get neighbor for phi - 17
        mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
		//........................................................................
		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
		mm18 = Phi[nn2x];					// get neighbor for phi - 18
        mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);


		//............Compute the Color Gradient...................................
		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
		//............Compute the Chemical Potential...............................
        chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
        chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
		//............Compute the Mixed Gradient...................................
		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
		
		// q=0
		m0 = dist[n];
		// q=1
		m1 = dist[2*Np+n]; 

        // q=2
		m2 = dist[1*Np+n];  

		// q=3
		m3 = dist[4*Np+n];

		// q = 4
		m4 = dist[3*Np+n];

		// q=5
		m5 = dist[6*Np+n];

		// q = 6
		m6 = dist[5*Np+n];
		
		// q=7
		m7 = dist[8*Np+n];

		// q = 8
		m8 = dist[7*Np+n];

		// q=9
		m9 = dist[10*Np+n];

		// q = 10
		m10 = dist[9*Np+n];

		// q=11
		m11 = dist[12*Np+n];

		// q=12
		m12 = dist[11*Np+n];

		// q=13
		m13 = dist[14*Np+n];

		// q=14
		m14 = dist[13*Np+n];

		// q=15
		m15 = dist[16*Np+n];

		// q=16
		m16 = dist[15*Np+n];

		// q=17
		m17 = dist[18*Np+n];

		// q=18
		m18 = dist[17*Np+n];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);

        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
                  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);

        //compute equilibrium distributions
            feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
     0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
      (-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
            feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
      (2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
            feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
         chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
            feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
            feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
         chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
            feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
            feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
     0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
            feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
         2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
      (2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
            feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
      (2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
            feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
            feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
     (mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
        0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));

		// q = 1
		dist[1*Np+n] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));

		// q=2
		dist[2*Np+n] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));

		// q = 3
		dist[3*Np+n] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 4
		dist[4*Np+n] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 5
		dist[5*Np+n] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 6
		dist[6*Np+n] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 7
		dist[7*Np+n] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (0.2222222222222222 + (ux + uy)*(ux + uy) - 
        0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 8
		dist[8*Np+n] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uy)*(ux + uy) + 
        0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 9
		dist[9*Np+n] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 10
		dist[10*Np+n] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 11
		dist[11*Np+n] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (0.2222222222222222 + (ux + uz)*(ux + uz) - 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 12
		dist[12*Np+n] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
       + (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 13
		dist[13*Np+n] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q= 14
		dist[14*Np+n] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 15
		dist[15*Np+n] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
      (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
       + (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));

		// q = 16
		dist[16*Np+n] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
       + (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 17
		dist[17*Np+n] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
       + (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 18
		dist[18*Np+n] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
        //----------------------------------------------------------------------------------------------------------------------------------------//

        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
        //Update chemical potential on device
        mu_phi[n] = chem;
        //Update color gradient on device
		ColorGrad[0*Np+n] = nx;
		ColorGrad[1*Np+n] = ny;
		ColorGrad[2*Np+n] = nz;

	}
}

extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(int *neighborList, int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, 
                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
                                                int strideY, int strideZ, int start, int finish, int Np){
	
	int n,nn,nn2x,ijk;
	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
    double chem;//chemical potential
    double phi; //phase field
    double rho0;//fluid density
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;
	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
	double mm3,mm5,mm7;
    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
    double nx,ny,nz;//normal color gradient
    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor

	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
    double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
	double tau;//position dependent LB relaxation time for fluid
    double C,theta;
    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
    double phi_temp;

	for (int n=start; n<finish; n++){

		rho0 = Den[n];//load density

		// Get the 1D index based on regular data layout
		ijk = Map[n];
        phi = Phi[ijk];// load phase field
        phi_temp = phi;
        if (phi>1.f) phi_temp=1.0;
        if (phi<-1.f) phi_temp=-1.0;

		// local relaxation time
		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);

		//					COMPUTE THE COLOR GRADIENT
		//........................................................................
		//.................Read Phase Indicator Values............................
		//........................................................................
		nn = ijk-1;							// neighbor index (get convention)
		m1 = Phi[nn];						// get neighbor for phi - 1
		//........................................................................
		nn = ijk+1;							// neighbor index (get convention)
		m2 = Phi[nn];						// get neighbor for phi - 2
		//........................................................................
		nn = ijk-strideY;							// neighbor index (get convention)
		m3 = Phi[nn];					// get neighbor for phi - 3
		//........................................................................
		nn = ijk+strideY;							// neighbor index (get convention)
		m4 = Phi[nn];					// get neighbor for phi - 4
		//........................................................................
		nn = ijk-strideZ;						// neighbor index (get convention)
		m5 = Phi[nn];					// get neighbor for phi - 5
		//........................................................................
		nn = ijk+strideZ;						// neighbor index (get convention)
		m6 = Phi[nn];					// get neighbor for phi - 6
		//........................................................................
		nn = ijk-strideY-1;						// neighbor index (get convention)
		m7 = Phi[nn];					// get neighbor for phi - 7
		//........................................................................
		nn = ijk+strideY+1;						// neighbor index (get convention)
		m8 = Phi[nn];					// get neighbor for phi - 8
		//........................................................................
		nn = ijk+strideY-1;						// neighbor index (get convention)
		m9 = Phi[nn];					// get neighbor for phi - 9
		//........................................................................
		nn = ijk-strideY+1;						// neighbor index (get convention)
		m10 = Phi[nn];					// get neighbor for phi - 10
		//........................................................................
		nn = ijk-strideZ-1;						// neighbor index (get convention)
		m11 = Phi[nn];					// get neighbor for phi - 11
		//........................................................................
		nn = ijk+strideZ+1;						// neighbor index (get convention)
		m12 = Phi[nn];					// get neighbor for phi - 12
		//........................................................................
		nn = ijk+strideZ-1;						// neighbor index (get convention)
		m13 = Phi[nn];					// get neighbor for phi - 13
		//........................................................................
		nn = ijk-strideZ+1;						// neighbor index (get convention)
		m14 = Phi[nn];					// get neighbor for phi - 14
		//........................................................................
		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
		m15 = Phi[nn];					// get neighbor for phi - 15
		//........................................................................
		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
		m16 = Phi[nn];					// get neighbor for phi - 16
		//........................................................................
		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
		m17 = Phi[nn];					// get neighbor for phi - 17
		//........................................................................
		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
		m18 = Phi[nn];					// get neighbor for phi - 18

        // compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
		//........................................................................
		nn2x = ijk-2;							// neighbor index (get convention)
		mm1 = Phi[nn2x];						// get neighbor for phi - 1
        mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
		//........................................................................
		nn2x = ijk+2;							// neighbor index (get convention)
		mm2 = Phi[nn2x];						// get neighbor for phi - 2
        mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
		//........................................................................
		nn2x = ijk-strideY*2;							// neighbor index (get convention)
		mm3 = Phi[nn2x];					// get neighbor for phi - 3
        mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
		//........................................................................
		nn2x = ijk+strideY*2;							// neighbor index (get convention)
		mm4 = Phi[nn2x];					// get neighbor for phi - 4
        mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
		//........................................................................
		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
		mm5 = Phi[nn2x];					// get neighbor for phi - 5
        mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
		//........................................................................
		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
		mm6 = Phi[nn2x];					// get neighbor for phi - 6
        mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
		//........................................................................
		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
		mm7 = Phi[nn2x];					// get neighbor for phi - 7
        mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
		//........................................................................
		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
		mm8 = Phi[nn2x];					// get neighbor for phi - 8
        mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
		//........................................................................
		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
		mm9 = Phi[nn2x];					// get neighbor for phi - 9
        mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
		//........................................................................
		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
		mm10 = Phi[nn2x];					// get neighbor for phi - 10
        mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
		//........................................................................
		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
		mm11 = Phi[nn2x];					// get neighbor for phi - 11
        mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
		//........................................................................
		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
		mm12 = Phi[nn2x];					// get neighbor for phi - 12
        mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
		//........................................................................
		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
		mm13 = Phi[nn2x];					// get neighbor for phi - 13
        mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
		//........................................................................
		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
		mm14 = Phi[nn2x];					// get neighbor for phi - 14
        mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
		//........................................................................
		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
		mm15 = Phi[nn2x];					// get neighbor for phi - 15
        mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
		//........................................................................
		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
		mm16 = Phi[nn2x];					// get neighbor for phi - 16
        mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
		//........................................................................
		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
		mm17 = Phi[nn2x];					// get neighbor for phi - 17
        mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
		//........................................................................
		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
		mm18 = Phi[nn2x];					// get neighbor for phi - 18
        mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);


		//............Compute the Color Gradient...................................
		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
		//............Compute the Chemical Potential...............................
        //chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
        //chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
        chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi_temp+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi_temp));//intermediate var, i.e. the laplacian
        chem = 4.0*beta*phi_temp*(phi_temp+1.0)*(phi_temp-1.0)-kappa*chem;
		//............Compute the Mixed Gradient...................................
		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
		
        //de-noise color gradient and mixed gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		if (C<1.0e-12) nx=ny=nz=0.0;
		double mg_mag = sqrt(mgx*mgx+mgy*mgy+mgz*mgz);
		if (mg_mag<1.0e-12) mgx=mgy=mgz=0.0;
        //maybe you can also de-noise chemical potential ? within the bulk phase chem should be ZERO
        if (fabs(chem)<1.0e-12) chem=0.0;

		// q=0
		m0 = dist[n];
		// q=1
		nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
		m1 = dist[nr1]; // reading the f1 data into register fq

		nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
		m2 = dist[nr2];  // reading the f2 data into register fq

		// q=3
		nr3 = neighborList[n+2*Np]; // neighbor 4
		m3 = dist[nr3];

		// q = 4
		nr4 = neighborList[n+3*Np]; // neighbor 3
		m4 = dist[nr4];

		// q=5
		nr5 = neighborList[n+4*Np];
		m5 = dist[nr5];

		// q = 6
		nr6 = neighborList[n+5*Np];
		m6 = dist[nr6];
		
		// q=7
		nr7 = neighborList[n+6*Np];
		m7 = dist[nr7];

		// q = 8
		nr8 = neighborList[n+7*Np];
		m8 = dist[nr8];

		// q=9
		nr9 = neighborList[n+8*Np];
		m9 = dist[nr9];

		// q = 10
		nr10 = neighborList[n+9*Np];
		m10 = dist[nr10];

		// q=11
		nr11 = neighborList[n+10*Np];
		m11 = dist[nr11];

		// q=12
		nr12 = neighborList[n+11*Np];
		m12 = dist[nr12];

		// q=13
		nr13 = neighborList[n+12*Np];
		m13 = dist[nr13];

		// q=14
		nr14 = neighborList[n+13*Np];
		m14 = dist[nr14];

		// q=15
		nr15 = neighborList[n+14*Np];
		m15 = dist[nr15];

		// q=16
		nr16 = neighborList[n+15*Np];
		m16 = dist[nr16];

		// q=17
		nr17 = neighborList[n+16*Np];
		m17 = dist[nr17];

		// q=18
		nr18 = neighborList[n+17*Np];
		m18 = dist[nr18];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
                  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);

        //compute equilibrium distributions
            feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
     0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
      (-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
            feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
      (2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
            feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
         chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
            feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
            feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
         chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
            feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
            feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
     0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
            feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
         2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
      (2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
            feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
      (2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
            feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
            feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
     (mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
        0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz)))); 

		// q = 1
		dist[nr2] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));

		// q=2
		dist[nr1] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));

		// q = 3
		dist[nr4] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 4
		dist[nr3] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 5
		dist[nr6] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 6
		dist[nr5] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 7
		dist[nr8] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (0.2222222222222222 + (ux + uy)*(ux + uy) - 
        0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 8
		dist[nr7] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uy)*(ux + uy) + 
        0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 9
		dist[nr10] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 10
		dist[nr9] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 11
		dist[nr12] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (0.2222222222222222 + (ux + uz)*(ux + uz) - 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 12
		dist[nr11] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
       + (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 13
		dist[nr14] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q= 14
		dist[nr13] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 15
		dist[nr16] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
      (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
       + (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));

		// q = 16
		dist[nr15] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
       + (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 17
		dist[nr18] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
       + (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 18
		dist[nr17] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
        //----------------------------------------------------------------------------------------------------------------------------------------//


        // ----------------------------- compute phase field evolution ----------------------------------------
		//Normalize the Color Gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		double ColorMag = C;
		if (C==0.0) ColorMag=1.0;
		nx = nx/ColorMag;
		ny = ny/ColorMag;
		nz = nz/ColorMag;		
        //compute surface tension-related parameter
        //theta = 4.5*M*2.0*(1-phi*phi)/W;
        theta = 4.5*M*2.0*(1-phi_temp*phi_temp)/W;

        //load distributions of phase field
		//q=0
		h0 = hq[n];
		//q=1
		h1 = hq[nr1]; 

		//q=2
		h2 = hq[nr2];  

		//q=3
		h3 = hq[nr3];

		//q=4
		h4 = hq[nr4];

		//q=5
		h5 = hq[nr5];

		//q=6
		h6 = hq[nr6];

        //-------------------------------- BGK collison for phase field ---------------------------------//
		// q = 0
		hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;

		// q = 1
		hq[nr2] = h1 - (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;

		// q = 2
		hq[nr1] = h2 - (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;

		// q = 3
		hq[nr4] = h3 - (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;

		// q = 4
		hq[nr3] = h4 - (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;

		// q = 5
		hq[nr6] = h5 - (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;

		// q = 6
		hq[nr5] = h6 - (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM;
		//........................................................................

        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
        //Update chemical potential on device
        mu_phi[n] = chem;
        //Update color gradient on device
		ColorGrad[0*Np+n] = nx;
		ColorGrad[1*Np+n] = ny;
		ColorGrad[2*Np+n] = nz;

	}
}

extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(int *Map, double *dist, double *hq, double *Den,	double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
                                                double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 
                                                int strideY, int strideZ, int start, int finish, int Np){
	
	int n,nn,nn2x,ijk;
	//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
    double chem;//chemical potential
    double phi; //phase field
    double rho0;//fluid density
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;
	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
	double mm3,mm5,mm7;
    double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
    double nx,ny,nz;//normal color gradient
    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor

	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
    double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
	double tau;//position dependent LB relaxation time for fluid
    double C,theta;
    double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
    double phi_temp;

	for (int n=start; n<finish; n++){

		rho0 = Den[n];//load density

		// Get the 1D index based on regular data layout
		ijk = Map[n];
        phi = Phi[ijk];// load phase field
        phi_temp = phi;
        if (phi>1.f) phi_temp=1.0;
        if (phi<-1.f) phi_temp=-1.0;

		// local relaxation time
		tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);

		//					COMPUTE THE COLOR GRADIENT
		//........................................................................
		//.................Read Phase Indicator Values............................
		//........................................................................
		nn = ijk-1;							// neighbor index (get convention)
		m1 = Phi[nn];						// get neighbor for phi - 1
		//........................................................................
		nn = ijk+1;							// neighbor index (get convention)
		m2 = Phi[nn];						// get neighbor for phi - 2
		//........................................................................
		nn = ijk-strideY;							// neighbor index (get convention)
		m3 = Phi[nn];					// get neighbor for phi - 3
		//........................................................................
		nn = ijk+strideY;							// neighbor index (get convention)
		m4 = Phi[nn];					// get neighbor for phi - 4
		//........................................................................
		nn = ijk-strideZ;						// neighbor index (get convention)
		m5 = Phi[nn];					// get neighbor for phi - 5
		//........................................................................
		nn = ijk+strideZ;						// neighbor index (get convention)
		m6 = Phi[nn];					// get neighbor for phi - 6
		//........................................................................
		nn = ijk-strideY-1;						// neighbor index (get convention)
		m7 = Phi[nn];					// get neighbor for phi - 7
		//........................................................................
		nn = ijk+strideY+1;						// neighbor index (get convention)
		m8 = Phi[nn];					// get neighbor for phi - 8
		//........................................................................
		nn = ijk+strideY-1;						// neighbor index (get convention)
		m9 = Phi[nn];					// get neighbor for phi - 9
		//........................................................................
		nn = ijk-strideY+1;						// neighbor index (get convention)
		m10 = Phi[nn];					// get neighbor for phi - 10
		//........................................................................
		nn = ijk-strideZ-1;						// neighbor index (get convention)
		m11 = Phi[nn];					// get neighbor for phi - 11
		//........................................................................
		nn = ijk+strideZ+1;						// neighbor index (get convention)
		m12 = Phi[nn];					// get neighbor for phi - 12
		//........................................................................
		nn = ijk+strideZ-1;						// neighbor index (get convention)
		m13 = Phi[nn];					// get neighbor for phi - 13
		//........................................................................
		nn = ijk-strideZ+1;						// neighbor index (get convention)
		m14 = Phi[nn];					// get neighbor for phi - 14
		//........................................................................
		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
		m15 = Phi[nn];					// get neighbor for phi - 15
		//........................................................................
		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
		m16 = Phi[nn];					// get neighbor for phi - 16
		//........................................................................
		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
		m17 = Phi[nn];					// get neighbor for phi - 17
		//........................................................................
		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
		m18 = Phi[nn];					// get neighbor for phi - 18

        // compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
		//........................................................................
		nn2x = ijk-2;							// neighbor index (get convention)
		mm1 = Phi[nn2x];						// get neighbor for phi - 1
        mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
		//........................................................................
		nn2x = ijk+2;							// neighbor index (get convention)
		mm2 = Phi[nn2x];						// get neighbor for phi - 2
        mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
		//........................................................................
		nn2x = ijk-strideY*2;							// neighbor index (get convention)
		mm3 = Phi[nn2x];					// get neighbor for phi - 3
        mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
		//........................................................................
		nn2x = ijk+strideY*2;							// neighbor index (get convention)
		mm4 = Phi[nn2x];					// get neighbor for phi - 4
        mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
		//........................................................................
		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
		mm5 = Phi[nn2x];					// get neighbor for phi - 5
        mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
		//........................................................................
		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
		mm6 = Phi[nn2x];					// get neighbor for phi - 6
        mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
		//........................................................................
		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
		mm7 = Phi[nn2x];					// get neighbor for phi - 7
        mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
		//........................................................................
		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
		mm8 = Phi[nn2x];					// get neighbor for phi - 8
        mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
		//........................................................................
		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
		mm9 = Phi[nn2x];					// get neighbor for phi - 9
        mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
		//........................................................................
		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
		mm10 = Phi[nn2x];					// get neighbor for phi - 10
        mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
		//........................................................................
		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
		mm11 = Phi[nn2x];					// get neighbor for phi - 11
        mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
		//........................................................................
		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
		mm12 = Phi[nn2x];					// get neighbor for phi - 12
        mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
		//........................................................................
		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
		mm13 = Phi[nn2x];					// get neighbor for phi - 13
        mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
		//........................................................................
		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
		mm14 = Phi[nn2x];					// get neighbor for phi - 14
        mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
		//........................................................................
		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
		mm15 = Phi[nn2x];					// get neighbor for phi - 15
        mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
		//........................................................................
		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
		mm16 = Phi[nn2x];					// get neighbor for phi - 16
        mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
		//........................................................................
		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
		mm17 = Phi[nn2x];					// get neighbor for phi - 17
        mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
		//........................................................................
		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
		mm18 = Phi[nn2x];					// get neighbor for phi - 18
        mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);


		//............Compute the Color Gradient...................................
		nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
		ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
		nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
		//............Compute the Chemical Potential...............................
        //chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
        //chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
        chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi_temp+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi_temp));//intermediate var, i.e. the laplacian
        chem = 4.0*beta*phi_temp*(phi_temp+1.0)*(phi_temp-1.0)-kappa*chem;
		//............Compute the Mixed Gradient...................................
		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
		
        //de-noise color gradient and mixed gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		if (C<1.0e-12) nx=ny=nz=0.0;
		double mg_mag = sqrt(mgx*mgx+mgy*mgy+mgz*mgz);
		if (mg_mag<1.0e-12) mgx=mgy=mgz=0.0;
        //maybe you can also de-noise chemical potential ? within the bulk phase chem should be ZERO
        if (fabs(chem)<1.0e-12) chem=0.0;

		// q=0
		m0 = dist[n];
		// q=1
		m1 = dist[2*Np+n]; 

        // q=2
		m2 = dist[1*Np+n];  

		// q=3
		m3 = dist[4*Np+n];

		// q = 4
		m4 = dist[3*Np+n];

		// q=5
		m5 = dist[6*Np+n];

		// q = 6
		m6 = dist[5*Np+n];
		
		// q=7
		m7 = dist[8*Np+n];

		// q = 8
		m8 = dist[7*Np+n];

		// q=9
		m9 = dist[10*Np+n];

		// q = 10
		m10 = dist[9*Np+n];

		// q=11
		m11 = dist[12*Np+n];

		// q=12
		m12 = dist[11*Np+n];

		// q=13
		m13 = dist[14*Np+n];

		// q=14
		m14 = dist[13*Np+n];

		// q=15
		m15 = dist[16*Np+n];

		// q=16
		m16 = dist[15*Np+n];

		// q=17
		m17 = dist[18*Np+n];

		// q=18
		m18 = dist[17*Np+n];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
                  +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);

        //compute equilibrium distributions
            feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
     0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
      (-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
            feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
      (2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
            feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
     0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux + 
         chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
            feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
            feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) + 
         rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy + 
         chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz))); 
            feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz)))); 
            feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
     0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
      (-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) + 
         rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz + 
         chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz)))); 
            feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) + 
         2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))); 
            feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy + 
       0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) + 
         2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))); 
            feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
      (2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))); 
            feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz + 
       0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
      (2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) + 
         2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))); 
            feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
      (2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
       0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))); 
            feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
      (2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))); 
            feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) - 
     0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))); 
            feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
      (-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) - 
     0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz + 
       0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
      (2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) + 
         2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
        (4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))); 

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 
     (mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) + 
        0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));

		// q = 1
		dist[1*Np+n] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));

		// q=2
		dist[2*Np+n] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
     (mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) + 
        0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));

		// q = 3
		dist[3*Np+n] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 4
		dist[4*Np+n] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) + 
        0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 5
		dist[5*Np+n] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 6
		dist[6*Np+n] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
        0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) + 
        0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 7
		dist[7*Np+n] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (0.2222222222222222 + (ux + uy)*(ux + uy) - 
        0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 8
		dist[8*Np+n] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uy)*(ux + uy) + 
        0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 9
		dist[9*Np+n] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
     (mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));

		// q = 10
		dist[10*Np+n] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (ux - uy)*(ux - uy) + 
        0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
     (mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));

		// q = 11
		dist[11*Np+n] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (0.2222222222222222 + (ux + uz)*(ux + uz) - 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 12
		dist[12*Np+n] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
       + (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q = 13
		dist[13*Np+n] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
     (mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));

		// q= 14
		dist[14*Np+n] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
      (-0.2222222222222222 - (ux - uz)*(ux - uz) + 
        0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
     (mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
      (-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));

		// q = 15
		dist[15*Np+n] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
      (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
       + (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));

		// q = 16
		dist[16*Np+n] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
      (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
       + (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
      (-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 17
		dist[17*Np+n] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
       + (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));

		// q = 18
		dist[18*Np+n] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
      (-0.2222222222222222 - (uy - uz)*(uy - uz) + 
        0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
     (mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
      (-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
         (-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
        0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
        //----------------------------------------------------------------------------------------------------------------------------------------//


        // ----------------------------- compute phase field evolution ----------------------------------------
		//Normalize the Color Gradient
		C = sqrt(nx*nx+ny*ny+nz*nz);
		double ColorMag = C;
		if (C==0.0) ColorMag=1.0;
		nx = nx/ColorMag;
		ny = ny/ColorMag;
		nz = nz/ColorMag;		
        //compute surface tension-related parameter
        //theta = 4.5*M*2.0*(1-phi*phi)/W;
        theta = 4.5*M*2.0*(1-phi_temp*phi_temp)/W;

        //load distributions of phase field
		//q=0
		h0 = hq[n];
		//q=1
		h1 = hq[2*Np+n]; 

		//q=2
		h2 = hq[1*Np+n];  

		//q=3
		h3 = hq[4*Np+n];

		//q=4
		h4 = hq[3*Np+n];

		//q=5
		h5 = hq[6*Np+n];

		//q=6
		h6 = hq[5*Np+n];

        //-------------------------------- BGK collison for phase field ---------------------------------//
		// q = 0
		hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;

		// q = 1
		hq[1*Np+n] = h1 - (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;

		// q = 2
		hq[2*Np+n] = h2 - (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;

		// q = 3
		hq[3*Np+n] = h3 - (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;

		// q = 4
		hq[4*Np+n] = h4 - (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;

		// q = 5
		hq[5*Np+n] = h5 - (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;

		// q = 6
		hq[6*Np+n] = h6 - (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM;
		//........................................................................

        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
        //Update chemical potential on device
        mu_phi[n] = chem;
        //Update color gradient on device
		ColorGrad[0*Np+n] = nx;
		ColorGrad[1*Np+n] = ny;
		ColorGrad[2*Np+n] = nz;

	}
}


extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,  
                                                                double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
	
	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;

	for (int n=start; n<finish; n++){

		// q=0
		m0 = dist[n];
		// q=1
		nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
		m1 = dist[nr1]; // reading the f1 data into register fq

		nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
		m2 = dist[nr2];  // reading the f2 data into register fq

		// q=3
		nr3 = neighborList[n+2*Np]; // neighbor 4
		m3 = dist[nr3];

		// q = 4
		nr4 = neighborList[n+3*Np]; // neighbor 3
		m4 = dist[nr4];

		// q=5
		nr5 = neighborList[n+4*Np];
		m5 = dist[nr5];

		// q = 6
		nr6 = neighborList[n+5*Np];
		m6 = dist[nr6];
		
		// q=7
		nr7 = neighborList[n+6*Np];
		m7 = dist[nr7];

		// q = 8
		nr8 = neighborList[n+7*Np];
		m8 = dist[nr8];

		// q=9
		nr9 = neighborList[n+8*Np];
		m9 = dist[nr9];

		// q = 10
		nr10 = neighborList[n+9*Np];
		m10 = dist[nr10];

		// q=11
		nr11 = neighborList[n+10*Np];
		m11 = dist[nr11];

		// q=12
		nr12 = neighborList[n+11*Np];
		m12 = dist[nr12];

		// q=13
		nr13 = neighborList[n+12*Np];
		m13 = dist[nr13];

		// q=14
		nr14 = neighborList[n+13*Np];
		m14 = dist[nr14];

		// q=15
		nr15 = neighborList[n+14*Np];
		m15 = dist[nr15];

		// q=16
		nr16 = neighborList[n+15*Np];
		m16 = dist[nr16];

		// q=17
		nr17 = neighborList[n+16*Np];
		m17 = dist[nr17];

		// q=18
		nr18 = neighborList[n+17*Np];
		m18 = dist[nr18];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
   (m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
       (-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
    tau;

		// q = 1
		dist[nr2] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
      0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
   (m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
      0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;

		// q=2
		dist[nr1] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
      0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
   (m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
      0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;

		// q = 3
		dist[nr4] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
      0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
   (m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 4
		dist[nr3] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
      0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
   (m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 5
		dist[nr6] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
      0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
   (m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 6
		dist[nr5] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
      0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
   (m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
      0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q = 7
		dist[nr8] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
    (0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
    - (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 8
		dist[nr7] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
    - (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 9
		dist[nr10] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
     - (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 10
		dist[nr9] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
    (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
    - (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 11
		dist[nr12] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
    (0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
    - (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 12
		dist[nr11] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
    - (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q = 13
		dist[nr14] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
    - (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q= 14
		dist[nr13] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
    (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
    - (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 15
		dist[nr16] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
    (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
    - (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 16
		dist[nr15] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
    - (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;

		// q = 17
		dist[nr18] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
    (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
    - (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;

		// q = 18
		dist[nr17] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
    (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
    - (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
        //----------------------------------------------------------------------------------------------------------------------------------------//


        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
	}
}

extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, 
                                                                 double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
	
    double ux,uy,uz;//fluid velocity 
    double p;//pressure
	// distribution functions
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m0,m3,m5,m7;

	for (int n=start; n<finish; n++){
		
		// q=0
		m0 = dist[n];
		// q=1
		m1 = dist[2*Np+n]; 

        // q=2
		m2 = dist[1*Np+n];  

		// q=3
		m3 = dist[4*Np+n];

		// q = 4
		m4 = dist[3*Np+n];

		// q=5
		m5 = dist[6*Np+n];

		// q = 6
		m6 = dist[5*Np+n];
		
		// q=7
		m7 = dist[8*Np+n];

		// q = 8
		m8 = dist[7*Np+n];

		// q=9
		m9 = dist[10*Np+n];

		// q = 10
		m10 = dist[9*Np+n];

		// q=11
		m11 = dist[12*Np+n];

		// q=12
		m12 = dist[11*Np+n];

		// q=13
		m13 = dist[14*Np+n];

		// q=14
		m14 = dist[13*Np+n];

		// q=15
		m15 = dist[16*Np+n];

		// q=16
		m16 = dist[15*Np+n];

		// q=17
		m17 = dist[18*Np+n];

		// q=18
		m18 = dist[17*Np+n];

        //compute fluid velocity
        ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
        uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
        uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
        //compute pressure
        p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);

        //------------------------------------------------- BCK collison ------------------------------------------------------------//
		// q=0
		dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - 
   (m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
       (-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
    tau; 

		// q = 1
		dist[1*Np+n] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
      0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - 
   (m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + 
      0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;

		// q=2
		dist[2*Np+n] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + 
      0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - 
   (m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + 
      0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;

		// q = 3
		dist[3*Np+n] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + 
      0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - 
   (m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 4
		dist[4*Np+n] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + 
      0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - 
   (m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 5
		dist[5*Np+n] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + 
      0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - 
   (m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 6
		dist[6*Np+n] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + 
      0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - 
   (m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
       (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + 
      0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q = 7
		dist[7*Np+n] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
    (0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
    - (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 8
		dist[8*Np+n] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
    - (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 9
		dist[9*Np+n] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
     - (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 10
		dist[10*Np+n] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
    (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
    - (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 
      0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + 
         0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;

		// q = 11
		dist[11*Np+n] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
    (0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
    - (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 12
		dist[12*Np+n] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
    - (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q = 13
		dist[13*Np+n] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
    - (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;

		// q= 14
		dist[14*Np+n] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
    (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
    - (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + 
         0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 15
		dist[15*Np+n] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
    (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
    - (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;

		// q = 16
		dist[16*Np+n] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
    (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
    - (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;

		// q = 17
		dist[17*Np+n] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
    (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
    - (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 
      0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
       (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;

		// q = 18
		dist[18*Np+n] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
    (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
    - (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
       (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 
      0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
       (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + 
         0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
        //----------------------------------------------------------------------------------------------------------------------------------------//

        //Update velocity on device
		Vel[0*Np+n] = ux;
		Vel[1*Np+n] = uy;
		Vel[2*Np+n] = uz;
        //Update pressure on device
        Pressure[n] = p;
	}
}

extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np){

	int nn,nn2x,ijk;
	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
	double m3,m5,m7;
	double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
	double mm3,mm5,mm7;
    //double nx,ny,nz;//normal color gradient
    double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
    double phi;

	for (int n=start; n<finish; n++){

		// Get the 1D index based on regular data layout
		ijk = Map[n];
        phi = Phi[ijk];// load phase field
		//					COMPUTE THE COLOR GRADIENT
		//........................................................................
		//.................Read Phase Indicator Values............................
		//........................................................................
		nn = ijk-1;							// neighbor index (get convention)
		m1 = Phi[nn];						// get neighbor for phi - 1
		//........................................................................
		nn = ijk+1;							// neighbor index (get convention)
		m2 = Phi[nn];						// get neighbor for phi - 2
		//........................................................................
		nn = ijk-strideY;							// neighbor index (get convention)
		m3 = Phi[nn];					// get neighbor for phi - 3
		//........................................................................
		nn = ijk+strideY;							// neighbor index (get convention)
		m4 = Phi[nn];					// get neighbor for phi - 4
		//........................................................................
		nn = ijk-strideZ;						// neighbor index (get convention)
		m5 = Phi[nn];					// get neighbor for phi - 5
		//........................................................................
		nn = ijk+strideZ;						// neighbor index (get convention)
		m6 = Phi[nn];					// get neighbor for phi - 6
		//........................................................................
		nn = ijk-strideY-1;						// neighbor index (get convention)
		m7 = Phi[nn];					// get neighbor for phi - 7
		//........................................................................
		nn = ijk+strideY+1;						// neighbor index (get convention)
		m8 = Phi[nn];					// get neighbor for phi - 8
		//........................................................................
		nn = ijk+strideY-1;						// neighbor index (get convention)
		m9 = Phi[nn];					// get neighbor for phi - 9
		//........................................................................
		nn = ijk-strideY+1;						// neighbor index (get convention)
		m10 = Phi[nn];					// get neighbor for phi - 10
		//........................................................................
		nn = ijk-strideZ-1;						// neighbor index (get convention)
		m11 = Phi[nn];					// get neighbor for phi - 11
		//........................................................................
		nn = ijk+strideZ+1;						// neighbor index (get convention)
		m12 = Phi[nn];					// get neighbor for phi - 12
		//........................................................................
		nn = ijk+strideZ-1;						// neighbor index (get convention)
		m13 = Phi[nn];					// get neighbor for phi - 13
		//........................................................................
		nn = ijk-strideZ+1;						// neighbor index (get convention)
		m14 = Phi[nn];					// get neighbor for phi - 14
		//........................................................................
		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
		m15 = Phi[nn];					// get neighbor for phi - 15
		//........................................................................
		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
		m16 = Phi[nn];					// get neighbor for phi - 16
		//........................................................................
		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
		m17 = Phi[nn];					// get neighbor for phi - 17
		//........................................................................
		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
		m18 = Phi[nn];					// get neighbor for phi - 18

        // compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
		//........................................................................
		nn2x = ijk-2;							// neighbor index (get convention)
		mm1 = Phi[nn2x];						// get neighbor for phi - 1
        mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
		//........................................................................
		nn2x = ijk+2;							// neighbor index (get convention)
		mm2 = Phi[nn2x];						// get neighbor for phi - 2
        mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
		//........................................................................
		nn2x = ijk-strideY*2;							// neighbor index (get convention)
		mm3 = Phi[nn2x];					// get neighbor for phi - 3
        mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
		//........................................................................
		nn2x = ijk+strideY*2;							// neighbor index (get convention)
		mm4 = Phi[nn2x];					// get neighbor for phi - 4
        mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
		//........................................................................
		nn2x = ijk-strideZ*2;						// neighbor index (get convention)
		mm5 = Phi[nn2x];					// get neighbor for phi - 5
        mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
		//........................................................................
		nn2x = ijk+strideZ*2;						// neighbor index (get convention)
		mm6 = Phi[nn2x];					// get neighbor for phi - 6
        mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
		//........................................................................
		nn2x = ijk-strideY*2-2;						// neighbor index (get convention)
		mm7 = Phi[nn2x];					// get neighbor for phi - 7
        mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
		//........................................................................
		nn2x = ijk+strideY*2+2;						// neighbor index (get convention)
		mm8 = Phi[nn2x];					// get neighbor for phi - 8
        mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
		//........................................................................
		nn2x = ijk+strideY*2-2;						// neighbor index (get convention)
		mm9 = Phi[nn2x];					// get neighbor for phi - 9
        mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
		//........................................................................
		nn2x = ijk-strideY*2+2;						// neighbor index (get convention)
		mm10 = Phi[nn2x];					// get neighbor for phi - 10
        mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
		//........................................................................
		nn2x = ijk-strideZ*2-2;						// neighbor index (get convention)
		mm11 = Phi[nn2x];					// get neighbor for phi - 11
        mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
		//........................................................................
		nn2x = ijk+strideZ*2+2;						// neighbor index (get convention)
		mm12 = Phi[nn2x];					// get neighbor for phi - 12
        mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
		//........................................................................
		nn2x = ijk+strideZ*2-2;						// neighbor index (get convention)
		mm13 = Phi[nn2x];					// get neighbor for phi - 13
        mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
		//........................................................................
		nn2x = ijk-strideZ*2+2;						// neighbor index (get convention)
		mm14 = Phi[nn2x];					// get neighbor for phi - 14
        mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
		//........................................................................
		nn2x = ijk-strideZ*2-strideY*2;					// neighbor index (get convention)
		mm15 = Phi[nn2x];					// get neighbor for phi - 15
        mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
		//........................................................................
		nn2x = ijk+strideZ*2+strideY*2;					// neighbor index (get convention)
		mm16 = Phi[nn2x];					// get neighbor for phi - 16
        mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
		//........................................................................
		nn2x = ijk+strideZ*2-strideY*2;					// neighbor index (get convention)
		mm17 = Phi[nn2x];					// get neighbor for phi - 17
        mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
		//........................................................................
		nn2x = ijk-strideZ*2+strideY*2;					// neighbor index (get convention)
		mm18 = Phi[nn2x];					// get neighbor for phi - 18
        mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);


		//............Compute the Color Gradient...................................
		//nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
		//ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
		//nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
		//............Compute the Mixed Gradient...................................
		mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
		mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
		mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));

		ColorGrad[0*Np+n] = mgx;
		ColorGrad[1*Np+n] = mgy;
		ColorGrad[2*Np+n] = mgz;
    }
}