diff --git a/common/Domain.cpp b/common/Domain.cpp
index d75ec4a5..fc797a8e 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -1286,3 +1286,150 @@ void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
   File.close();
 }
 
+void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData)
+{
+	//........................................................................................
+	// Reading the user-defined input file
+    // NOTE: so far it only supports BC=0 (periodic) and BC=5 (mixed reflection)
+    //       because if checkerboard or inlet/outlet buffer layers are added, the
+    //       value of the void space is undefined. 
+    // NOTE: if BC=5 is used, where the inlet and outlet layers of the domain are modified, 
+    //       user needs to modify the input file accordingly before LBPM simulator read
+    //       the input file.
+	//........................................................................................
+	int rank_offset = 0;
+	int RANK = rank();
+	int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
+	int64_t global_Nx,global_Ny,global_Nz;
+	int64_t i,j,k,n;
+    //TODO These offset we may still need them
+	int64_t xStart,yStart,zStart;
+	xStart=yStart=zStart=0;
+
+	// Read domain parameters
+    // TODO currently the size of the data is still read from Domain{}; 
+    //      but user may have a user-specified size
+	auto size = database->getVector<int>( "n" );
+	auto SIZE = database->getVector<int>( "N" );
+	auto nproc = database->getVector<int>( "nproc" );
+    //TODO currently the funcationality "offset" is disabled as the user-defined input data may have a different size from that of the input domain 
+	if (database->keyExists( "offset" )){
+		auto offset = database->getVector<int>( "offset" );
+		xStart = offset[0];
+		yStart = offset[1];
+		zStart = offset[2];
+	}
+	
+	nx = size[0];
+	ny = size[1];
+	nz = size[2];
+	nprocx = nproc[0];
+	nprocy = nproc[1];
+	nprocz = nproc[2];
+	global_Nx = SIZE[0];
+	global_Ny = SIZE[1];
+	global_Nz = SIZE[2];
+	nprocs=nprocx*nprocy*nprocz;
+
+	double *SegData = NULL;
+	if (RANK==0){
+		printf("User-defined input file: %s (data type: %s)\n",Filename.c_str(),Datatype.c_str());
+        printf("NOTE: currently only BC=0 or 5 supports user-defined input file!\n");
+		// Rank=0 reads the entire segmented data and distributes to worker processes
+		printf("Dimensions of the user-defined input file: %ld x %ld x %ld \n",global_Nx,global_Ny,global_Nz);
+		int64_t SIZE = global_Nx*global_Ny*global_Nz;
+
+		if (Datatype == "double"){
+			printf("Reading input data as double precision floating number\n");
+		    SegData = new double[SIZE];
+			FILE *SEGDAT = fopen(Filename.c_str(),"rb");
+			if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading user-defined file!\n");
+			size_t ReadSeg;
+			ReadSeg=fread(SegData,8,SIZE,SEGDAT);
+			if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading file: %s\n",Filename.c_str());
+			fclose(SEGDAT);
+		}
+        else{
+            ERROR("Error: User-defined input file only supports double-precision floating number!\n"); 
+        }
+		printf("Read file successfully from %s \n",Filename.c_str());
+	}
+	
+	// Get the rank info
+	int64_t N = (nx+2)*(ny+2)*(nz+2);
+
+	// number of sites to use for periodic boundary condition transition zone
+	//int64_t z_transition_size = (nprocz*nz - (global_Nz - zStart))/2;
+	//if (z_transition_size < 0) z_transition_size=0;
+	int64_t z_transition_size = 0;
+
+	//char LocalRankFilename[1000];//just for debug
+	double *loc_id;
+	loc_id = new double [(nx+2)*(ny+2)*(nz+2)];
+
+	// Set up the sub-domains
+	if (RANK==0){
+        printf("Decomposing user-defined input file\n");
+		printf("Distributing subdomains across %i processors \n",nprocs);
+		printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
+		printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
+		printf("Size of transition region: %ld \n", z_transition_size);
+
+		for (int kp=0; kp<nprocz; kp++){
+			for (int jp=0; jp<nprocy; jp++){
+				for (int ip=0; ip<nprocx; ip++){
+					// rank of the process that gets this subdomain
+					int rnk = kp*nprocx*nprocy + jp*nprocx + ip;
+					// Pack and send the subdomain for rnk
+					for (k=0;k<nz+2;k++){
+						for (j=0;j<ny+2;j++){
+							for (i=0;i<nx+2;i++){
+								int64_t x = xStart + ip*nx + i-1;
+								int64_t y = yStart + jp*ny + j-1;
+								// int64_t z = zStart + kp*nz + k-1;
+								int64_t z = zStart + kp*nz + k-1 - z_transition_size;
+								if (x<xStart) 	x=xStart;
+								if (!(x<global_Nx))	x=global_Nx-1;
+								if (y<yStart) 	y=yStart;
+								if (!(y<global_Ny))	y=global_Ny-1;
+								if (z<zStart) 	z=zStart;
+								if (!(z<global_Nz))	z=global_Nz-1;
+								int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
+								int64_t nglobal = z*global_Nx*global_Ny+y*global_Nx+x;
+								loc_id[nlocal] = SegData[nglobal];
+							}
+						}
+					}
+					if (rnk==0){
+						for (k=0;k<nz+2;k++){
+							for (j=0;j<ny+2;j++){
+								for (i=0;i<nx+2;i++){
+									int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
+									UserData[nlocal] = loc_id[nlocal];
+								}
+							}
+						}
+					}
+					else{
+						//printf("Sending data to process %i \n", rnk);
+						MPI_Send(loc_id,N,MPI_DOUBLE,rnk,15,Comm);
+					}
+					// Write the data for this rank data 
+                    // NOTE just for debug
+					//sprintf(LocalRankFilename,"%s.%05i",Filename.c_str(),rnk+rank_offset);
+					//FILE *ID = fopen(LocalRankFilename,"wb");
+					//fwrite(loc_id,8,(nx+2)*(ny+2)*(nz+2),ID);
+					//fclose(ID);
+				}
+			}
+		}
+
+	}
+	else{
+		// Recieve the subdomain from rank = 0
+		//printf("Ready to recieve data %i at process %i \n", N,rank);
+		MPI_Recv(UserData,N,MPI_DOUBLE,0,15,Comm,MPI_STATUS_IGNORE);
+	}
+	//Comm.barrier();
+	MPI_Barrier(Comm);
+}
diff --git a/common/Domain.h b/common/Domain.h
index df2812c1..9aa3e0ad 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -178,6 +178,7 @@ public: // Public variables (need to create accessors instead)
 
     void ReadIDs();
     void Decomp( const std::string& filename );
+    void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
     void CommunicateMeshHalo(DoubleArray &Mesh);
     void CommInit(); 
     int PoreCount();
diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 9a711330..edd601e8 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -85,6 +85,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
                                        double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
 
 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);
 
diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp
index 98d35142..5d7ef87a 100644
--- a/cpu/Ion.cpp
+++ b/cpu/Ion.cpp
@@ -220,6 +220,22 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
 	}
 }
 
+extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
+{
+	int n;
+    double DenInit;
+	for (n=0; n<Np; n++){
+        DenInit = Den[n];
+		dist[0*Np+n] = 0.25*DenInit;
+		dist[1*Np+n] = 0.125*DenInit;		
+		dist[2*Np+n] = 0.125*DenInit;	
+		dist[3*Np+n] = 0.125*DenInit;	
+		dist[4*Np+n] = 0.125*DenInit;	
+		dist[5*Np+n] = 0.125*DenInit;	
+		dist[6*Np+n] = 0.125*DenInit;	
+	}
+}
+
 extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
     
     int n;
diff --git a/gpu/Ion.cu b/gpu/Ion.cu
index 877e5591..601d5310 100644
--- a/gpu/Ion.cu
+++ b/gpu/Ion.cu
@@ -263,6 +263,27 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenI
 	}
 }
 
+__global__  void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np){
+
+	int n;
+    double DenInit;
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
+		if (n<Np) {
+            DenInit = Den[n];
+            dist[0*Np+n] = 0.25*DenInit;
+            dist[1*Np+n] = 0.125*DenInit;		
+            dist[2*Np+n] = 0.125*DenInit;	
+            dist[3*Np+n] = 0.125*DenInit;	
+            dist[4*Np+n] = 0.125*DenInit;	
+            dist[5*Np+n] = 0.125*DenInit;	
+            dist[6*Np+n] = 0.125*DenInit;	
+		}
+	}
+}
+
 __global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
 
     int n;
@@ -345,6 +366,18 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
 	//cudaProfilerStop();
 }
 
+extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np){
+
+	//cudaProfilerStart();
+	dvc_ScaLBL_D3Q7_Ion_Init_FromFile<<<NBLOCKS,NTHREADS >>>(dist,Den,Np);
+
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_Ion_Init_FromFile: %s \n",cudaGetErrorString(err));
+	}
+	//cudaProfilerStop();
+}
+
 extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
 
 	//cudaProfilerStart();
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 7366c0b4..ac53d5cc 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -136,7 +136,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
     }
     //read initial ion concentration list; INPUT unit [mol/m^3]
     //it must be converted to LB unit [mol/lu^3]
-	if (ion_db->keyExists("IonConcentrationList")){
+    if (ion_db->keyExists("IonConcentrationList")){
         IonConcentration.clear();
 	    IonConcentration = ion_db->getVector<double>( "IonConcentrationList" );
         if (IonConcentration.size()!=number_ion_species){
@@ -544,6 +544,35 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
 	}
 }
 
+void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion)
+{
+    double *Ci_host;
+    Ci_host = new double[N];
+	double VALUE=0.f;
+
+    Mask->ReadFromFile(File_ion[0],File_ion[1],Ci_host);
+
+	for (int k=0;k<Nz;k++){
+		for (int j=0;j<Ny;j++){
+			for (int i=0;i<Nx;i++){
+				int idx = Map(i,j,k);
+				if (!(idx < 0)){
+				    int n = k*Nx*Ny+j*Nx+i;
+                    //NOTE: default user-input unit: mol/m^3
+                    VALUE = Ci_host[n]*(h*h*h*1.0e-18);
+                    if (VALUE<0.0){
+                        ERROR("Error: Ion concentration value must be a positive number! \n");
+                    }
+                    else{
+                        Ci[idx] = VALUE;
+                    }
+                }
+			}
+		}
+	}
+    delete [] Ci_host;
+}
+
 void ScaLBL_IonModel::Create(){
 	/*
 	 *  This function creates the variables needed to run a LBM 
@@ -601,8 +630,6 @@ void ScaLBL_IonModel::Create(){
         ScaLBL_DeviceBarrier();
         delete [] IonSolid_host;
     }
-
-
 }        
 
 void ScaLBL_IonModel::Initialize(){
@@ -610,8 +637,26 @@ void ScaLBL_IonModel::Initialize(){
 	 * This function initializes model
 	 */
     if (rank==0)    printf ("LB Ion Solver: initializing D3Q7 distributions\n");
-	for (int ic=0; ic<number_ion_species; ic++){
-        ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np); 
+	if (ion_db->keyExists("IonConcentrationFile")){
+        //TODO: Need to figure out how to deal with multi-species concentration initialization
+        //NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
+		auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
+        double *Ci_host;
+        Ci_host = new double[number_ion_species*Np];
+        for (int ic=0; ic<number_ion_species; ic++){
+            AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
+        }
+	    ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
+	    MPI_Barrier(comm);
+        for (int ic=0; ic<number_ion_species; ic++){
+            ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np); 
+        }
+        delete [] Ci_host;
+    }
+    else{
+        for (int ic=0; ic<number_ion_species; ic++){
+            ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np); 
+        }
     }
     if (rank==0)    printf ("LB Ion Solver: initializing charge density\n");
 	for (int ic=0; ic<number_ion_species; ic++){
diff --git a/models/IonModel.h b/models/IonModel.h
index 0bc881af..995898bd 100644
--- a/models/IonModel.h
+++ b/models/IonModel.h
@@ -89,5 +89,6 @@ private:
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);    	
     void AssignSolidBoundary(double *ion_solid);
+    void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion);
     void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
 };