diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 4d88b728..5c3deba1 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -301,7 +301,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDi
 extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
 extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np);
 
-extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);
+extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np);
 
 
 // LBM Poisson solver
diff --git a/cuda/Ion.cu b/cuda/Ion.cu
index cdae76e8..1650449a 100644
--- a/cuda/Ion.cu
+++ b/cuda/Ion.cu
@@ -566,7 +566,7 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, in
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
 
     int n;
     double Ci;//ion concentration of species i
@@ -579,10 +579,17 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDe
 		//........Get 1-D index for this thread....................
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
+
             Ci = Den[n+ion_component*Np];
             CD = ChargeDensity[n];
+            if (ion_component == 0) CD=0.0;
             CD_tmp = F*IonValence*Ci;
-            ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
+            ChargeDensity[n] = CD + CD_tmp;
+
+ //           Ci = Den[n+ion_component*Np];
+   //         CD = ChargeDensity[n];
+     //       CD_tmp = F*IonValence*Ci;
+       //     ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
 		}
 	}
 }
@@ -660,7 +667,7 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
 
 	//cudaProfilerStart();
 	dvc_ScaLBL_D3Q7_Ion_ChargeDensity<<<NBLOCKS,NTHREADS >>>(Den,ChargeDensity,IonValence,ion_component,start,finish,Np);
diff --git a/hip/Ion.hip b/hip/Ion.hip
index 18be1aee..ea30063f 100644
--- a/hip/Ion.hip
+++ b/hip/Ion.hip
@@ -564,7 +564,7 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, in
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
 
     int n;
     double Ci;//ion concentration of species i
@@ -579,8 +579,9 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDe
 		if (n<finish) {
             Ci = Den[n+ion_component*Np];
             CD = ChargeDensity[n];
+	    if (ion_component == 0) CD=0.0;
             CD_tmp = F*IonValence*Ci;
-            ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
+            ChargeDensity[n] = CD + CD_tmp;
 		}
 	}
 }
@@ -658,7 +659,7 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
 
 	//cudaProfilerStart();
 	dvc_ScaLBL_D3Q7_Ion_ChargeDensity<<<NBLOCKS,NTHREADS >>>(Den,ChargeDensity,IonValence,ion_component,start,finish,Np);
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 2ac4133e..3479bb95 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -1443,17 +1443,35 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
 
     //Compute charge density for Poisson equation
     for (size_t ic = 0; ic < number_ion_species; ic++) {
-        ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic,
+      int Valence = IonValence[ic];
+      if (rank==0) printf("compute charge density for ion %i, Valence =%i \n", ic,Valence);
+
+      ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic,
                                       ScaLBL_Comm->FirstInterior(),
                                       ScaLBL_Comm->LastInterior(), Np);
-        ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0,
+        ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, 0,
                                       ScaLBL_Comm->LastExterior(), Np);
     }
 
+    DoubleArray Charge(Nx,Ny,Nz);
+    ScaLBL_Comm->RegularLayout(Map, ChargeDensity, Charge);
+    double charge_sum=0.0;
+    double charge_sum_total=0.0;
+    for (int k=1; k<Nz-1; k++){
+      for (int j=1; j<Ny-1; j++){
+	for (int i=1; i<Nx-1; i++){
+	  charge_sum += Charge(i,j,k);
+	}
+      }
+    }
+    printf("  Local charge value = %.8g (rank=%i)\n",charge_sum, rank);
     ScaLBL_Comm->Barrier();
     comm.barrier();
-    if (rank==0) printf(" IonMembrane: completeted full step  \n");
-    fflush(stdout);
+    
+    ScaLBL_Comm->Barrier();
+    comm.barrier();
+    //if (rank==0) printf(" IonMembrane: completeted full step  \n");
+    //fflush(stdout);
     //************************************************************************/
     //if (rank==0) printf("-------------------------------------------------------------------\n");
     //// Compute the walltime per timestep
diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp
index 39e0803f..4db1f3a2 100644
--- a/models/PoissonSolver.cpp
+++ b/models/PoissonSolver.cpp
@@ -557,6 +557,7 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times
 		ScaLBL_Comm->Barrier(); comm.barrier();
 		//************************************************************************/
 
+		
         // Check convergence of steady-state solution
         if (timestep==2){
             //save electric potential for convergence check