diff --git a/common/Membrane.cpp b/common/Membrane.cpp
index 20f17c41..4b5b37dd 100644
--- a/common/Membrane.cpp
+++ b/common/Membrane.cpp
@@ -420,6 +420,8 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
 		}
 	}
 	
+	int Q = 7; // for D3Q7 model
+	
 	/* go through the neighborlist structure */
 	/* count & cut the links */
 	if (rank == 0) printf("   Cut membrane links... \n");
@@ -470,82 +472,84 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
 						mlink++;
 					}
 
-					neighbor=Map(i-1,j-1,k);
-					dist=Distance(i-1,j-1,k);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[6*Np+idx]=idx + 8*Np;
-					}
+					if (Q > 7){
+						neighbor=Map(i-1,j-1,k);
+						dist=Distance(i-1,j-1,k);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[6*Np+idx]=idx + 8*Np;
+						}
 
-					neighbor=Map(i+1,j+1,k);
-					dist=Distance(i+1,j+1,k);
-					if (dist*locdist < 0.0){
-						neighborList[7*Np+idx]=idx + 7*Np;
-						mlink++;
-					}
+						neighbor=Map(i+1,j+1,k);
+						dist=Distance(i+1,j+1,k);
+						if (dist*locdist < 0.0){
+							neighborList[7*Np+idx]=idx + 7*Np;
+							mlink++;
+						}
 
-					neighbor=Map(i-1,j+1,k);
-					dist=Distance(i-1,j+1,k);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[8*Np+idx]=idx + 10*Np;
-					}
+						neighbor=Map(i-1,j+1,k);
+						dist=Distance(i-1,j+1,k);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[8*Np+idx]=idx + 10*Np;
+						}
 
-					neighbor=Map(i+1,j-1,k);
-					dist=Distance(i+1,j-1,k);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[9*Np+idx]=idx + 9*Np;
-						mlink++;
-					}
+						neighbor=Map(i+1,j-1,k);
+						dist=Distance(i+1,j-1,k);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[9*Np+idx]=idx + 9*Np;
+							mlink++;
+						}
 
-					neighbor=Map(i-1,j,k-1);
-					dist=Distance(i-1,j,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[10*Np+idx]=idx + 12*Np;
-					}
+						neighbor=Map(i-1,j,k-1);
+						dist=Distance(i-1,j,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[10*Np+idx]=idx + 12*Np;
+						}
 
-					neighbor=Map(i+1,j,k+1);
-					dist=Distance(i+1,j,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[11*Np+idx]=idx + 11*Np;
-						mlink++;
-					}
+						neighbor=Map(i+1,j,k+1);
+						dist=Distance(i+1,j,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[11*Np+idx]=idx + 11*Np;
+							mlink++;
+						}
 
-					neighbor=Map(i-1,j,k+1);
-					dist=Distance(i-1,j,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[12*Np+idx]=idx + 14*Np;
-					}
+						neighbor=Map(i-1,j,k+1);
+						dist=Distance(i-1,j,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[12*Np+idx]=idx + 14*Np;
+						}
 
-					neighbor=Map(i+1,j,k-1);
-					dist=Distance(i+1,j,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[13*Np+idx]=idx + 13*Np;
-						mlink++;
-					}
+						neighbor=Map(i+1,j,k-1);
+						dist=Distance(i+1,j,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[13*Np+idx]=idx + 13*Np;
+							mlink++;
+						}
 
-					neighbor=Map(i,j-1,k-1);
-					dist=Distance(i,j-1,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[14*Np+idx]=idx + 16*Np;
-					}
+						neighbor=Map(i,j-1,k-1);
+						dist=Distance(i,j-1,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[14*Np+idx]=idx + 16*Np;
+						}
 
-					neighbor=Map(i,j+1,k+1);
-					dist=Distance(i,j+1,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[15*Np+idx]=idx + 15*Np;
-						mlink++;
-					}
+						neighbor=Map(i,j+1,k+1);
+						dist=Distance(i,j+1,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[15*Np+idx]=idx + 15*Np;
+							mlink++;
+						}
 
-					neighbor=Map(i,j-1,k+1);
-					dist=Distance(i,j-1,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[16*Np+idx]=idx + 18*Np;
-					}
+						neighbor=Map(i,j-1,k+1);
+						dist=Distance(i,j-1,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[16*Np+idx]=idx + 18*Np;
+						}
 
-					neighbor=Map(i,j+1,k-1);
-					dist=Distance(i,j+1,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						neighborList[17*Np+idx]=idx + 17*Np;
-						mlink++;
+						neighbor=Map(i,j+1,k-1);
+						dist=Distance(i,j+1,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							neighborList[17*Np+idx]=idx + 17*Np;
+							mlink++;
+						}
 					}
 				}
 			}
@@ -578,8 +582,8 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
 
 					neighbor=Map(i+1,j,k);
 					dist=Distance(i+1,j,k);
-					if (dist*locdist < 0.0){
-						if (locdist < 0.0 && !(neighbor<0)){
+					if (dist*locdist < 0.0 && !(neighbor<0)){
+						if (locdist < 0.0 ){
 							localSite = 2*mlink;
 							neighborSite = 2*mlink+1;
 						}
@@ -630,113 +634,116 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
 						mlink++;
 					}
 
-					neighbor=Map(i+1,j+1,k);
-					dist=Distance(i+1,j+1,k);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
-						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 7*Np;
-						membraneLinks[neighborSite] = neighbor+8*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					}
+					if (Q > 7){
 
-					neighbor=Map(i+1,j-1,k);
-					dist=Distance(i+1,j-1,k);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
+						neighbor=Map(i+1,j+1,k);
+						dist=Distance(i+1,j+1,k);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 7*Np;
+							membraneLinks[neighborSite] = neighbor+8*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
 						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 9*Np;
-						membraneLinks[neighborSite] = neighbor + 10*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					}
 
-					neighbor=Map(i+1,j,k+1);
-					dist=Distance(i+1,j,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
+						neighbor=Map(i+1,j-1,k);
+						dist=Distance(i+1,j-1,k);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 9*Np;
+							membraneLinks[neighborSite] = neighbor + 10*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
 						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 11*Np;
-						membraneLinks[neighborSite] = neighbor + 12*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					}
 
-					neighbor=Map(i+1,j,k-1);
-					dist=Distance(i+1,j,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
+						neighbor=Map(i+1,j,k+1);
+						dist=Distance(i+1,j,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 11*Np;
+							membraneLinks[neighborSite] = neighbor + 12*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
 						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 13*Np;
-						membraneLinks[neighborSite] = neighbor + 14*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					}
 
-					neighbor=Map(i,j+1,k+1);
-					dist=Distance(i,j+1,k+1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
+						neighbor=Map(i+1,j,k-1);
+						dist=Distance(i+1,j,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 13*Np;
+							membraneLinks[neighborSite] = neighbor + 14*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
 						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 15*Np;
-						membraneLinks[neighborSite] = neighbor + 16*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					}
 
-					neighbor=Map(i,j+1,k-1);
-					dist=Distance(i,j+1,k-1);
-					if (dist*locdist < 0.0 && !(neighbor<0)){
-						if (locdist < 0.0){
-							localSite = 2*mlink;
-							neighborSite = 2*mlink+1;
+						neighbor=Map(i,j+1,k+1);
+						dist=Distance(i,j+1,k+1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 15*Np;
+							membraneLinks[neighborSite] = neighbor + 16*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
 						}
-						else{
-							localSite = 2*mlink+1;
-							neighborSite = 2*mlink;
-						}
-						membraneLinks[localSite] = idx + 17*Np;
-						membraneLinks[neighborSite] = neighbor + 18*Np;
-						membraneDist[localSite] = locdist;
-						membraneDist[neighborSite] = dist;
-						mlink++;
-					} 
+
+						neighbor=Map(i,j+1,k-1);
+						dist=Distance(i,j+1,k-1);
+						if (dist*locdist < 0.0 && !(neighbor<0)){
+							if (locdist < 0.0){
+								localSite = 2*mlink;
+								neighborSite = 2*mlink+1;
+							}
+							else{
+								localSite = 2*mlink+1;
+								neighborSite = 2*mlink;
+							}
+							membraneLinks[localSite] = idx + 17*Np;
+							membraneLinks[neighborSite] = neighbor + 18*Np;
+							membraneDist[localSite] = locdist;
+							membraneDist[neighborSite] = dist;
+							mlink++;
+						} 
+					}
 				}
 			}
 		}
@@ -1384,13 +1391,26 @@ void Membrane::AssignCoefficients(int *Map, double *Psi, string method){
 		ThresholdMassFractionOut = 1.0;				
 		ThresholdMassFractionIn = 1.0;
 	}
-	
 	if (method == "Voltage Gated Potassium"){
 		MassFractionIn = 0.0;
 		MassFractionOut = 0.0;
 		ThresholdMassFractionOut = 0.0;				
 		ThresholdMassFractionIn = 1.0;
 	}
+	if (method == "impermeable"){
+		printf(".... impermeable membrane \n");
+		MassFractionIn = 0.0;
+		MassFractionOut = 0.0;
+		ThresholdMassFractionOut = 0.0;				
+		ThresholdMassFractionIn = 0.0;
+	}
+	if (method == "Na+"){
+		printf(".... Na+ permeable membrane \n");
+		MassFractionIn = 0.5;
+		MassFractionOut = 0.5;
+		ThresholdMassFractionOut = 0.5;				
+		ThresholdMassFractionIn = 0.5;
+	}
 	
 	ScaLBL_D3Q7_Membrane_AssignLinkCoef(MembraneLinks, Map, MembraneDistance, Psi, MembraneCoef,
 			Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp
index 915edda3..6ae2e849 100644
--- a/cpu/Ion.cpp
+++ b/cpu/Ion.cpp
@@ -144,7 +144,7 @@ extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
 		iq = membrane[2*link]; 	ip = membrane[2*link+1];
 		nq = iq%Np;				np = ip%Np;
 		fq  = dist[iq];			fp = dist[ip];
-		fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+ap*fq;
+		fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+aq*fq;
 		Cq = Den[nq];			Cp = Den[np];
 		Cq += fqq - fq;			Cp += fpp - fp;
 		Den[nq] = Cq;			Den[np] = Cp;
@@ -257,7 +257,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
     for (n = start; n < finish; n++) {
 
         //Load data
-        Ci = Den[n];
         Ex = ElectricField[n + 0 * Np];
         Ey = ElectricField[n + 1 * Np];
         Ez = ElectricField[n + 2 * Np];
@@ -290,6 +289,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
         f6 = dist[nr6];
 
         // compute diffusive flux
+        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
         flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
         flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
         flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@@ -303,6 +303,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         FluxElectrical[n + 2 * Np] = uEPz * Ci;
         
+        Den[n] = Ci;
+
         /* use logistic function to prevent negative distributions*/
         X = 4.0 * (ux + uEPx);
         Y = 4.0 * (uy + uEPy);
@@ -363,7 +365,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
     for (n = start; n < finish; n++) {
 
         //Load data
-        Ci = Den[n];
+        //Ci = Den[n];
         Ex = ElectricField[n + 0 * Np];
         Ey = ElectricField[n + 1 * Np];
         Ez = ElectricField[n + 2 * Np];
@@ -383,6 +385,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
         f6 = dist[5 * Np + n];
 
         // compute diffusive flux
+        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
         flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
         flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
         flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@@ -396,6 +399,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         FluxElectrical[n + 2 * Np] = uEPz * Ci;
         
+        Den[n] = Ci;
+        
         /* use logistic function to prevent negative distributions*/
         X = 4.0 * (ux + uEPx);
         Y = 4.0 * (uy + uEPy);
diff --git a/cuda/Ion.cu b/cuda/Ion.cu
index bacb1cfc..cdae76e8 100644
--- a/cuda/Ion.cu
+++ b/cuda/Ion.cu
@@ -213,7 +213,7 @@ __global__  void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *co
 			iq = membrane[2*link]; 	ip = membrane[2*link+1];
 			nq = iq%Np;				np = ip%Np;
 			fq  = dist[iq];			fp = dist[ip];
-			fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+ap*fq;
+			fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+aq*fq;
 			Cq = Den[nq];			Cp = Den[np];
 			Cq += fqq - fq;			Cp += fpp - fp;
 			Den[nq] = Cq;			Den[np] = Cp;
@@ -333,93 +333,96 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
 
-            //Load data
-            Ci=Den[n];
-            Ex=ElectricField[n+0*Np];
-            Ey=ElectricField[n+1*Np];
-            Ez=ElectricField[n+2*Np];
-            ux=Velocity[n+0*Np];
-            uy=Velocity[n+1*Np];
-            uz=Velocity[n+2*Np];
-            uEPx=zi*Di/Vt*Ex;
-            uEPy=zi*Di/Vt*Ey;
-            uEPz=zi*Di/Vt*Ez;
+	        //Load data
+	        Ex = ElectricField[n + 0 * Np];
+	        Ey = ElectricField[n + 1 * Np];
+	        Ez = ElectricField[n + 2 * Np];
+	        ux = Velocity[n + 0 * Np];
+	        uy = Velocity[n + 1 * Np];
+	        uz = Velocity[n + 2 * Np];
+	        uEPx = zi * Di / Vt * Ex;
+	        uEPy = zi * Di / Vt * Ey;
+	        uEPz = zi * Di / Vt * Ez;
 
-            // q=0
-            f0 = dist[n];
-            // q=1
-            nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
-            f1 = dist[nr1]; // reading the f1 data into register fq
-            // q=2
-            nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
-            f2 = dist[nr2];  // reading the f2 data into register fq
-            // q=3
-            nr3 = neighborList[n+2*Np]; // neighbor 4
-            f3 = dist[nr3];
-            // q=4
-            nr4 = neighborList[n+3*Np]; // neighbor 3
-            f4 = dist[nr4];
-            // q=5
-            nr5 = neighborList[n+4*Np];
-            f5 = dist[nr5];
-            // q=6
-            nr6 = neighborList[n+5*Np];
-            f6 = dist[nr6];
-            
-            // compute diffusive flux
-            flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
-            flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
-            flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
-            FluxDiffusive[n+0*Np] = flux_diffusive_x;
-            FluxDiffusive[n+1*Np] = flux_diffusive_y;
-            FluxDiffusive[n+2*Np] = flux_diffusive_z;
-            FluxAdvective[n+0*Np] = ux*Ci;
-            FluxAdvective[n+1*Np] = uy*Ci;
-            FluxAdvective[n+2*Np] = uz*Ci;
-            FluxElectrical[n+0*Np] = uEPx*Ci;
-            FluxElectrical[n+1*Np] = uEPy*Ci;
-            FluxElectrical[n+2*Np] = uEPz*Ci;
-            
-            /* use logistic function to prevent negative distributions*/
-            X = 4.0 * (ux + uEPx);
-            Y = 4.0 * (uy + uEPy);
-            Z = 4.0 * (uz + uEPz);
-            factor_x = X / sqrt(1 + X*X);
-            factor_y = Y / sqrt(1 + Y*Y);
-            factor_z = Z / sqrt(1 + Z*Z);
+	        // q=0
+	        f0 = dist[n];
+	        // q=1
+	        nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
+	        f1 = dist[nr1];        // reading the f1 data into register fq
+	        // q=2
+	        nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist)
+	        f2 = dist[nr2];             // reading the f2 data into register fq
+	        // q=3
+	        nr3 = neighborList[n + 2 * Np]; // neighbor 4
+	        f3 = dist[nr3];
+	        // q=4
+	        nr4 = neighborList[n + 3 * Np]; // neighbor 3
+	        f4 = dist[nr4];
+	        // q=5
+	        nr5 = neighborList[n + 4 * Np];
+	        f5 = dist[nr5];
+	        // q=6
+	        nr6 = neighborList[n + 5 * Np];
+	        f6 = dist[nr6];
 
-            // q=0
-            dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+	        // compute diffusive flux
+	        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
+	        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
+	        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
+	        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
+	        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
+	        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
+	        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
+	        FluxAdvective[n + 0 * Np] = ux * Ci;
+	        FluxAdvective[n + 1 * Np] = uy * Ci;
+	        FluxAdvective[n + 2 * Np] = uz * Ci;
+	        FluxElectrical[n + 0 * Np] = uEPx * Ci;
+	        FluxElectrical[n + 1 * Np] = uEPy * Ci;
+	        FluxElectrical[n + 2 * Np] = uEPz * Ci;
+	        
+	        Den[n] = Ci;
 
-            // q = 1
-            dist[nr2] =
-                f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
-            //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
+	        /* use logistic function to prevent negative distributions*/
+	        X = 4.0 * (ux + uEPx);
+	        Y = 4.0 * (uy + uEPy);
+	        Z = 4.0 * (uz + uEPz);
+	        factor_x = X / sqrt(1 + X*X);
+	        factor_y = Y / sqrt(1 + Y*Y);
+	        factor_z = Z / sqrt(1 + Z*Z);
+
+	        // q=0
+	        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+
+	        // q = 1
+	        dist[nr2] =
+	            f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
+	        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
 
-            // q=2
-            dist[nr1] =
-                    f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-            //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+	        // q=2
+	        dist[nr1] =
+	                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
+	        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
-            // q = 3
-            dist[nr4] =
-                    f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
-            //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+	        // q = 3
+	        dist[nr4] =
+	                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
+	        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
-            // q = 4
-            dist[nr3] =
-                    f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-            //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+	        // q = 4
+	        dist[nr3] =
+	                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
+	        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
-            // q = 5
-            dist[nr6] =
-                    f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 +  factor_z);
-            //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+	        // q = 5
+	        dist[nr6] =
+	                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 +  factor_z);
+	        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+
+	        // q = 6
+	        dist[nr5] =
+	            f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
 
-            // q = 6
-            dist[nr5] =
-                f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
 		}
 	}
 }
@@ -441,80 +444,83 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
 
-            //Load data
-            Ci=Den[n];
-            Ex=ElectricField[n+0*Np];
-            Ey=ElectricField[n+1*Np];
-            Ez=ElectricField[n+2*Np];
-            ux=Velocity[n+0*Np];
-            uy=Velocity[n+1*Np];
-            uz=Velocity[n+2*Np];
-            uEPx=zi*Di/Vt*Ex;
-            uEPy=zi*Di/Vt*Ey;
-            uEPz=zi*Di/Vt*Ez;
+		       //Load data
+		        //Ci = Den[n];
+		        Ex = ElectricField[n + 0 * Np];
+		        Ey = ElectricField[n + 1 * Np];
+		        Ez = ElectricField[n + 2 * Np];
+		        ux = Velocity[n + 0 * Np];
+		        uy = Velocity[n + 1 * Np];
+		        uz = Velocity[n + 2 * Np];
+		        uEPx = zi * Di / Vt * Ex;
+		        uEPy = zi * Di / Vt * Ey;
+		        uEPz = zi * Di / Vt * Ez;
 
-            f0 = dist[n];
-            f1 = dist[2*Np+n];
-            f2 = dist[1*Np+n];
-            f3 = dist[4*Np+n];
-            f4 = dist[3*Np+n];
-            f5 = dist[6*Np+n];
-            f6 = dist[5*Np+n];
-            
-            // compute diffusive flux
-            flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
-            flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
-            flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
-            FluxDiffusive[n+0*Np] = flux_diffusive_x;
-            FluxDiffusive[n+1*Np] = flux_diffusive_y;
-            FluxDiffusive[n+2*Np] = flux_diffusive_z;
-            FluxAdvective[n+0*Np] = ux*Ci;
-            FluxAdvective[n+1*Np] = uy*Ci;
-            FluxAdvective[n+2*Np] = uz*Ci;
-            FluxElectrical[n+0*Np] = uEPx*Ci;
-            FluxElectrical[n+1*Np] = uEPy*Ci;
-            FluxElectrical[n+2*Np] = uEPz*Ci;
+		        f0 = dist[n];
+		        f1 = dist[2 * Np + n];
+		        f2 = dist[1 * Np + n];
+		        f3 = dist[4 * Np + n];
+		        f4 = dist[3 * Np + n];
+		        f5 = dist[6 * Np + n];
+		        f6 = dist[5 * Np + n];
 
-            /* use logistic function to prevent negative distributions*/
-            X = 4.0 * (ux + uEPx);
-            Y = 4.0 * (uy + uEPy);
-            Z = 4.0 * (uz + uEPz);
-            factor_x = X / sqrt(1 + X*X);
-            factor_y = Y / sqrt(1 + Y*Y);
-            factor_z = Z / sqrt(1 + Z*Z);
+		        // compute diffusive flux
+		        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
+		        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
+		        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
+		        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
+		        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
+		        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
+		        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
+		        FluxAdvective[n + 0 * Np] = ux * Ci;
+		        FluxAdvective[n + 1 * Np] = uy * Ci;
+		        FluxAdvective[n + 2 * Np] = uz * Ci;
+		        FluxElectrical[n + 0 * Np] = uEPx * Ci;
+		        FluxElectrical[n + 1 * Np] = uEPy * Ci;
+		        FluxElectrical[n + 2 * Np] = uEPz * Ci;
+		        
+		        Den[n] = Ci;
+		        
+		        /* use logistic function to prevent negative distributions*/
+		        X = 4.0 * (ux + uEPx);
+		        Y = 4.0 * (uy + uEPy);
+		        Z = 4.0 * (uz + uEPz);
+		        factor_x = X / sqrt(1 + X*X);
+		        factor_y = Y / sqrt(1 + Y*Y);
+		        factor_z = Z / sqrt(1 + Z*Z);
 
-            // q=0
-            dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+		        // q=0
+		        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
 
-            // q = 1
-            dist[1 * Np + n] =
-                    f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
-            //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
+		        // q = 1
+		        dist[1 * Np + n] =
+		                f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
+		        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
-            // q=2
-            dist[2 * Np + n] =
-                    f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-            //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+		        // q=2
+		        dist[2 * Np + n] =
+		                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
+		        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
-            // q = 3
-            dist[3 * Np + n] =
-                    f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
-            //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+		        // q = 3
+		        dist[3 * Np + n] =
+		                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
+		        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
-            // q = 4
-            dist[4 * Np + n] =
-                    f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-            //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+		        // q = 4
+		        dist[4 * Np + n] =
+		                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
+		        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
-            // q = 5
-            dist[5 * Np + n] =
-                    f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
-            //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+		        // q = 5
+		        dist[5 * Np + n] =
+		                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
+		        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
-            // q = 6
-            dist[6 * Np + n] =
-                    f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
-            //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
+		        // q = 6
+		        dist[6 * Np + n] =
+		                f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
+		        //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
 		}
 	}
 }
diff --git a/example/SingleCell/NaCl-cell.py b/example/SingleCell/NaCl-cell.py
new file mode 100644
index 00000000..78ee6076
--- /dev/null
+++ b/example/SingleCell/NaCl-cell.py
@@ -0,0 +1,37 @@
+import numpy as np
+import matplotlib.pylab as plt
+
+D=np.ones((40,40,40),dtype="uint8")
+
+cx = 20
+cy = 20
+cz = 20
+
+for i in range(0, 40):
+    for j in range (0, 40):
+        for k in range (0,40):
+            dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
+            if (dist < 15.5 ) :
+                D[i,j,k] = 2
+
+D.tofile("cell_40x40x40.raw")
+
+                
+C1=np.zeros((40,40,40),dtype="double")
+C2=np.zeros((40,40,40),dtype="double")
+
+for i in range(0, 40):
+    for j in range (0, 40):
+        for k in range (0,40):
+            #outside the cell
+            C1[i,j,k] = 125.0e-6  # Na
+            C2[i,j,k] = 125.0e-6  # Cl
+            dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
+            # inside the cell
+            if (dist < 15.5 ) :
+                C1[i,j,k] = 5.0e-6
+                C2[i,j,k] = 5.0e-6
+
+C1.tofile("cell_concentration_Na_40x40x40.raw")
+C2.tofile("cell_concentration_Cl_40x40x40.raw")
+
diff --git a/example/SingleCell/NaCl.db b/example/SingleCell/NaCl.db
new file mode 100644
index 00000000..f9bc4999
--- /dev/null
+++ b/example/SingleCell/NaCl.db
@@ -0,0 +1,75 @@
+MultiphysController {
+    timestepMax = 200
+    num_iter_Ion_List = 2
+    analysis_interval  = 50
+    tolerance = 1.0e-9
+    visualization_interval = 100        // Frequency to write visualization data
+    analysis_interval = 50    // Frequency to perform analysis
+}
+Stokes {
+    tau = 1.0
+    F = 0, 0, 0
+    ElectricField = 0, 0, 0 //body electric field; user-input unit: [V/m]
+    nu_phys = 0.889e-6      //fluid kinematic viscosity; user-input unit: [m^2/sec]
+}
+Ions {
+    IonConcentrationFile = "cell_concentration_Na_40x40x40.raw", "double", "cell_concentration_Cl_40x40x40.raw", "double"
+    temperature = 293.15 //unit [K]
+    number_ion_species = 2  //number of ions
+    tauList = 1.0, 1.0
+    IonDiffusivityList = 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
+    IonValenceList = 1, -1 //valence charge of ions; dimensionless; positive/negative integer
+    IonConcentrationList = 1.0e-6, 1.0e-6 //user-input unit: [mol/m^3]
+    BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration
+    //SolidLabels = 0 //solid labels for assigning solid boundary condition; ONLY for BC_Solid=1
+    //SolidValues = 1.0e-5 // user-input surface ion concentration unit: [mol/m^2]; ONLY for BC_Solid=1
+    FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
+}
+Poisson {
+    epsilonR = 78.5 //fluid dielectric constant [dimensionless]
+    BC_Inlet  = 0  // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
+    BC_Outlet = 0  // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
+    //--------------------------------------------------------------------------
+    //--------------------------------------------------------------------------
+    BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density
+    SolidLabels = 0 //solid labels for assigning solid boundary condition
+    SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2]
+    WriteLog = true //write convergence log for LB-Poisson solver
+    // ------------------------------- Testing Utilities ----------------------------------------
+    // ONLY for code debugging; the followings test sine/cosine voltage BCs; disabled by default
+    TestPeriodic = false
+    TestPeriodicTime = 1.0 //unit:[sec]
+    TestPeriodicTimeConv = 0.01 //unit:[sec]
+    TestPeriodicSaveInterval = 0.2 //unit:[sec]
+    //------------------------------ advanced setting ------------------------------------
+    timestepMax = 100000 //max timestep for obtaining steady-state electrical potential
+    analysis_interval  = 200 //timestep checking steady-state convergence
+    tolerance = 1.0e-6  //stopping criterion for steady-state solution
+}
+Domain {
+    Filename = "cell_40x40x40.raw"
+    nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
+    n = 40, 40, 40      // Size of local domain (Nx,Ny,Nz)
+    N = 40, 40, 40         // size of the input image
+    voxel_length = 1.0   //resolution; user-input unit: [um]
+    BC = 0              // Boundary condition type
+    ReadType = "8bit"
+    ReadValues  = 0, 1, 2
+    WriteValues = 0, 1, 2
+}
+Analysis {
+    analysis_interval = 100
+    subphase_analysis_interval = 50    // Frequency to perform analysis
+    restart_interval = 5000    // Frequency to write restart data
+    restart_file = "Restart"    // Filename to use for restart file (will append rank)
+    N_threads    = 4            // Number of threads to use
+    load_balance = "independent" // Load balance method to use: "none", "default", "independent"
+}
+Visualization {
+    save_electric_potential = true
+    save_concentration = true
+    save_velocity = true
+}
+Membrane {
+    MembraneLabels = 2
+}
diff --git a/hip/Ion.hip b/hip/Ion.hip
index a9b55f8c..d7e863bb 100644
--- a/hip/Ion.hip
+++ b/hip/Ion.hip
@@ -212,7 +212,7 @@ __global__  void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *co
 			iq = membrane[2*link]; 	ip = membrane[2*link+1];
 			nq = iq%Np;				np = ip%Np;
 			fq  = dist[iq];			fp = dist[ip];
-			fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+ap*fq;
+			fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+aq*fq;
 			Cq = Den[nq];			Cp = Den[np];
 			Cq += fqq - fq;			Cp += fpp - fp;
 			Den[nq] = Cq;			Den[np] = Cp;
@@ -331,93 +331,96 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
 
-            //Load data
-            Ci=Den[n];
-            Ex=ElectricField[n+0*Np];
-            Ey=ElectricField[n+1*Np];
-            Ez=ElectricField[n+2*Np];
-            ux=Velocity[n+0*Np];
-            uy=Velocity[n+1*Np];
-            uz=Velocity[n+2*Np];
-            uEPx=zi*Di/Vt*Ex;
-            uEPy=zi*Di/Vt*Ey;
-            uEPz=zi*Di/Vt*Ez;
+           //Load data
+        Ex = ElectricField[n + 0 * Np];
+        Ey = ElectricField[n + 1 * Np];
+        Ez = ElectricField[n + 2 * Np];
+        ux = Velocity[n + 0 * Np];
+        uy = Velocity[n + 1 * Np];
+        uz = Velocity[n + 2 * Np];
+        uEPx = zi * Di / Vt * Ex;
+        uEPy = zi * Di / Vt * Ey;
+        uEPz = zi * Di / Vt * Ez;
 
-            // q=0
-            f0 = dist[n];
-            // q=1
-            nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
-            f1 = dist[nr1]; // reading the f1 data into register fq
-            // q=2
-            nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
-            f2 = dist[nr2];  // reading the f2 data into register fq
-            // q=3
-            nr3 = neighborList[n+2*Np]; // neighbor 4
-            f3 = dist[nr3];
-            // q=4
-            nr4 = neighborList[n+3*Np]; // neighbor 3
-            f4 = dist[nr4];
-            // q=5
-            nr5 = neighborList[n+4*Np];
-            f5 = dist[nr5];
-            // q=6
-            nr6 = neighborList[n+5*Np];
-            f6 = dist[nr6];
-            
-            // compute diffusive flux
-            flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
-            flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
-            flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
-            FluxDiffusive[n+0*Np] = flux_diffusive_x;
-            FluxDiffusive[n+1*Np] = flux_diffusive_y;
-            FluxDiffusive[n+2*Np] = flux_diffusive_z;
-            FluxAdvective[n+0*Np] = ux*Ci;
-            FluxAdvective[n+1*Np] = uy*Ci;
-            FluxAdvective[n+2*Np] = uz*Ci;
-            FluxElectrical[n+0*Np] = uEPx*Ci;
-            FluxElectrical[n+1*Np] = uEPy*Ci;
-            FluxElectrical[n+2*Np] = uEPz*Ci;
-            
-            /* use logistic function to prevent negative distributions*/
-            X = 4.0 * (ux + uEPx);
-            Y = 4.0 * (uy + uEPy);
-            Z = 4.0 * (uz + uEPz);
-            factor_x = X / sqrt(1 + X*X);
-            factor_y = Y / sqrt(1 + Y*Y);
-            factor_z = Z / sqrt(1 + Z*Z);
+        // q=0
+        f0 = dist[n];
+        // q=1
+        nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
+        f1 = dist[nr1];        // reading the f1 data into register fq
+        // q=2
+        nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist)
+        f2 = dist[nr2];             // reading the f2 data into register fq
+        // q=3
+        nr3 = neighborList[n + 2 * Np]; // neighbor 4
+        f3 = dist[nr3];
+        // q=4
+        nr4 = neighborList[n + 3 * Np]; // neighbor 3
+        f4 = dist[nr4];
+        // q=5
+        nr5 = neighborList[n + 4 * Np];
+        f5 = dist[nr5];
+        // q=6
+        nr6 = neighborList[n + 5 * Np];
+        f6 = dist[nr6];
 
-            // q=0
-            dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+        // compute diffusive flux
+        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
+        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
+        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
+        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
+        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
+        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
+        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
+        FluxAdvective[n + 0 * Np] = ux * Ci;
+        FluxAdvective[n + 1 * Np] = uy * Ci;
+        FluxAdvective[n + 2 * Np] = uz * Ci;
+        FluxElectrical[n + 0 * Np] = uEPx * Ci;
+        FluxElectrical[n + 1 * Np] = uEPy * Ci;
+        FluxElectrical[n + 2 * Np] = uEPz * Ci;
+        
+        Den[n] = Ci;
 
-            // q = 1
-            dist[nr2] =
-                f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
-            //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
+        /* use logistic function to prevent negative distributions*/
+        X = 4.0 * (ux + uEPx);
+        Y = 4.0 * (uy + uEPy);
+        Z = 4.0 * (uz + uEPz);
+        factor_x = X / sqrt(1 + X*X);
+        factor_y = Y / sqrt(1 + Y*Y);
+        factor_z = Z / sqrt(1 + Z*Z);
+
+        // q=0
+        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+
+        // q = 1
+        dist[nr2] =
+            f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
+        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
 
-            // q=2
-            dist[nr1] =
-                    f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-            //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+        // q=2
+        dist[nr1] =
+                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
+        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
-            // q = 3
-            dist[nr4] =
-                    f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
-            //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+        // q = 3
+        dist[nr4] =
+                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
+        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
-            // q = 4
-            dist[nr3] =
-                    f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-            //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+        // q = 4
+        dist[nr3] =
+                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
+        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
-            // q = 5
-            dist[nr6] =
-                    f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 +  factor_z);
-            //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+        // q = 5
+        dist[nr6] =
+                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 +  factor_z);
+        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
-            // q = 6
-            dist[nr5] =
-                f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
+        // q = 6
+        dist[nr5] =
+            f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
+   
 		}
 	}
 }
@@ -440,79 +443,82 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
 		if (n<finish) {
 
             //Load data
-            Ci=Den[n];
-            Ex=ElectricField[n+0*Np];
-            Ey=ElectricField[n+1*Np];
-            Ez=ElectricField[n+2*Np];
-            ux=Velocity[n+0*Np];
-            uy=Velocity[n+1*Np];
-            uz=Velocity[n+2*Np];
-            uEPx=zi*Di/Vt*Ex;
-            uEPy=zi*Di/Vt*Ey;
-            uEPz=zi*Di/Vt*Ez;
+        //Ci = Den[n];
+        Ex = ElectricField[n + 0 * Np];
+        Ey = ElectricField[n + 1 * Np];
+        Ez = ElectricField[n + 2 * Np];
+        ux = Velocity[n + 0 * Np];
+        uy = Velocity[n + 1 * Np];
+        uz = Velocity[n + 2 * Np];
+        uEPx = zi * Di / Vt * Ex;
+        uEPy = zi * Di / Vt * Ey;
+        uEPz = zi * Di / Vt * Ez;
 
-            f0 = dist[n];
-            f1 = dist[2*Np+n];
-            f2 = dist[1*Np+n];
-            f3 = dist[4*Np+n];
-            f4 = dist[3*Np+n];
-            f5 = dist[6*Np+n];
-            f6 = dist[5*Np+n];
-            
-            // compute diffusive flux
-            flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
-            flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
-            flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
-            FluxDiffusive[n+0*Np] = flux_diffusive_x;
-            FluxDiffusive[n+1*Np] = flux_diffusive_y;
-            FluxDiffusive[n+2*Np] = flux_diffusive_z;
-            FluxAdvective[n+0*Np] = ux*Ci;
-            FluxAdvective[n+1*Np] = uy*Ci;
-            FluxAdvective[n+2*Np] = uz*Ci;
-            FluxElectrical[n+0*Np] = uEPx*Ci;
-            FluxElectrical[n+1*Np] = uEPy*Ci;
-            FluxElectrical[n+2*Np] = uEPz*Ci;
+        f0 = dist[n];
+        f1 = dist[2 * Np + n];
+        f2 = dist[1 * Np + n];
+        f3 = dist[4 * Np + n];
+        f4 = dist[3 * Np + n];
+        f5 = dist[6 * Np + n];
+        f6 = dist[5 * Np + n];
 
-            /* use logistic function to prevent negative distributions*/
-            X = 4.0 * (ux + uEPx);
-            Y = 4.0 * (uy + uEPy);
-            Z = 4.0 * (uz + uEPz);
-            factor_x = X / sqrt(1 + X*X);
-            factor_y = Y / sqrt(1 + Y*Y);
-            factor_z = Z / sqrt(1 + Z*Z);
+        // compute diffusive flux
+        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
+        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
+        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
+        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
+        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
+        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
+        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
+        FluxAdvective[n + 0 * Np] = ux * Ci;
+        FluxAdvective[n + 1 * Np] = uy * Ci;
+        FluxAdvective[n + 2 * Np] = uz * Ci;
+        FluxElectrical[n + 0 * Np] = uEPx * Ci;
+        FluxElectrical[n + 1 * Np] = uEPy * Ci;
+        FluxElectrical[n + 2 * Np] = uEPz * Ci;
+        
+        Den[n] = Ci;
+        
+        /* use logistic function to prevent negative distributions*/
+        X = 4.0 * (ux + uEPx);
+        Y = 4.0 * (uy + uEPy);
+        Z = 4.0 * (uz + uEPz);
+        factor_x = X / sqrt(1 + X*X);
+        factor_y = Y / sqrt(1 + Y*Y);
+        factor_z = Z / sqrt(1 + Z*Z);
 
-            // q=0
-            dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+        // q=0
+        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
 
-            // q = 1
-            dist[1 * Np + n] =
-                    f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
-            //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
+        // q = 1
+        dist[1 * Np + n] =
+                f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
+        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
-            // q=2
-            dist[2 * Np + n] =
-                    f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-            //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+        // q=2
+        dist[2 * Np + n] =
+                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
+        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
-            // q = 3
-            dist[3 * Np + n] =
-                    f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
-            //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+        // q = 3
+        dist[3 * Np + n] =
+                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
+        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
-            // q = 4
-            dist[4 * Np + n] =
-                    f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-            //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+        // q = 4
+        dist[4 * Np + n] =
+                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
+        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
-            // q = 5
-            dist[5 * Np + n] =
-                    f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
-            //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+        // q = 5
+        dist[5 * Np + n] =
+                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
+        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
-            // q = 6
-            dist[6 * Np + n] =
-                    f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
-            //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
+        // q = 6
+        dist[6 * Np + n] =
+                f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
+        //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
 		}
 	}
 }
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 58ddbf0d..ab2ee2ed 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -1296,90 +1296,39 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
     //ScaLBL_Comm->Barrier(); comm.barrier();
     //auto t1 = std::chrono::system_clock::now();
 
-    
     for (size_t ic = 0; ic < number_ion_species; ic++) {
         /* set the mass transfer coefficients for the membrane */
-        IonMembrane->AssignCoefficients(dvcMap, Psi, "default");
+    	if (ic == 0)
+    		IonMembrane->AssignCoefficients(dvcMap, Psi, "Na+");
+    	else {
+    		IonMembrane->AssignCoefficients(dvcMap, Psi, "impermeable");
+    	}
         timestep = 0;
         while (timestep < timestepMax[ic]) {
             //************************************************************************/
             // *************ODD TIMESTEP*************//
             timestep++;
-            //Update ion concentration and charge density
-            IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
-            ScaLBL_D3Q7_AAodd_IonConcentration(
-                IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
-                ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-            IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
-            /* ScaLBL_Comm->Barrier();
-            //--------------------------------------- Set boundary conditions -------------------------------------//
-            if (BoundaryConditionInlet[ic] > 0) {
-                switch (BoundaryConditionInlet[ic]) {
-                case 1:
-                    ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep);
-                    break;
-                case 21:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 22:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 23:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
-                        &Velocity[2 * Np], &ElectricField[2 * Np],
-                        IonDiffusivity[ic], IonValence[ic], Vt, timestep);
-                    break;
-                }
-            }
-            if (BoundaryConditionOutlet[ic] > 0) {
-                switch (BoundaryConditionOutlet[ic]) {
-                case 1:
-                    ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep);
-                    break;
-                case 21:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 22:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 23:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], &ElectricField[2 * Np],
-                        IonDiffusivity[ic], IonValence[ic], Vt, timestep);
-                    break;
-                }
-            }
-            */
-            //----------------------------------------------------------------------------------------------------//
-            ScaLBL_D3Q7_AAodd_IonConcentration(IonMembrane->NeighborList, &fq[ic * Np * 7],
-                                               &Ci[ic * Np], 0,
-                                               ScaLBL_Comm->LastExterior(), Np);
 
             //LB-Ion collison
+            IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
             ScaLBL_D3Q7_AAodd_Ion(
                 IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
                 &FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
                 &FluxElectrical[3 * ic * Np], Velocity, ElectricField,
                 IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt,
                 ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+            
+            IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
+
             ScaLBL_D3Q7_AAodd_Ion(
                 IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
                 &FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
                 &FluxElectrical[3 * ic * Np], Velocity, ElectricField,
                 IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt, 0,
                 ScaLBL_Comm->LastExterior(), Np);
+            
+           IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]);
+
 
             if (BoundaryConditionSolid == 1) {
                 //TODO IonSolid may also be species-dependent
@@ -1390,80 +1339,26 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
 
             // *************EVEN TIMESTEP*************//            
             timestep++;
-            //Update ion concentration and charge density
-            IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
-            ScaLBL_D3Q7_AAeven_IonConcentration(
-                &fq[ic * Np * 7], &Ci[ic * Np], ScaLBL_Comm->FirstInterior(),
-                ScaLBL_Comm->LastInterior(), Np);
-            IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
-            ScaLBL_Comm->Barrier();
-            //--------------------------------------- Set boundary conditions -------------------------------------//
-            /*if (BoundaryConditionInlet[ic] > 0) {
-                switch (BoundaryConditionInlet[ic]) {
-                case 1:
-                    ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep);
-                    break;
-                case 21:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],q
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 22:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 23:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
-                        &Velocity[2 * Np], &ElectricField[2 * Np],
-                        IonDiffusivity[ic], IonValence[ic], Vt, timestep);
-                    break;
-                }
-            }
-            if (BoundaryConditionOutlet[ic] > 0) {
-                switch (BoundaryConditionOutlet[ic]) {
-                case 1:
-                    ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep);
-                    break;
-                case 21:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 22:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], timestep);
-                    break;
-                case 23:
-                    ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z(
-                        IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
-                        &Velocity[2 * Np], &ElectricField[2 * Np],
-                        IonDiffusivity[ic], IonValence[ic], Vt, timestep);
-                    break;
-                }
-            } 
-            */
-            //----------------------------------------------------------------------------------------------------//
-            ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic * Np * 7], &Ci[ic * Np],
-                                                0, ScaLBL_Comm->LastExterior(),
-                                                Np);
 
             //LB-Ion collison
+            IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
             ScaLBL_D3Q7_AAeven_Ion(
                 &fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
                 &FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
                 Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],
                 rlx[ic], Vt, ScaLBL_Comm->FirstInterior(),
                 ScaLBL_Comm->LastInterior(), Np);
+
+            IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
+
             ScaLBL_D3Q7_AAeven_Ion(
                 &fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
                 &FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
                 Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],
                 rlx[ic], Vt, 0, ScaLBL_Comm->LastExterior(), Np);
+            
+            IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]);
+
 
             if (BoundaryConditionSolid == 1) {
                 //TODO IonSolid may also be species-dependent