diff --git a/hip/Poisson.cu b/hip/Poisson.cu
index 34975f58..e869a17b 100644
--- a/hip/Poisson.cu
+++ b/hip/Poisson.cu
@@ -104,7 +104,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -122,8 +122,11 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //rho_e = Den_charge[n];
+            //rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
             
@@ -184,7 +187,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -201,8 +204,11 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //rho_e = Den_charge[n];
+            //rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
 
@@ -293,10 +299,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
@@ -305,10 +311,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
diff --git a/hip/Stokes.cu b/hip/Stokes.cu
index a6a05fba..79ea9c4a 100644
--- a/hip/Stokes.cu
+++ b/hip/Stokes.cu
@@ -6,7 +6,7 @@
 #define NTHREADS 256
 
 
-__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
 
     int n;
     double fq;
@@ -47,9 +47,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -511,7 +514,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC,int start, int finish, int Np){
 
     int n;
     double fq;
@@ -551,9 +554,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -970,10 +976,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
@@ -982,10 +988,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){