diff --git a/ABACUS.develop/source/input_conv.cpp b/ABACUS.develop/source/input_conv.cpp
index 1c83346ac0..74fdd9979f 100644
--- a/ABACUS.develop/source/input_conv.cpp
+++ b/ABACUS.develop/source/input_conv.cpp
@@ -37,7 +37,7 @@ void Input_Conv::Convert(void)
     global_pseudo_type = INPUT.pseudo_type;
 	ucell.latName = INPUT.latname; 
 	ucell.ntype = INPUT.ntype;
-	ucell.nelec = INPUT.nelec;
+	CHR.nelec = INPUT.nelec;
 //  ucell.lmaxmax = INPUT.lmaxmax;
 
     NBANDS = INPUT.nbands;
diff --git a/ABACUS.develop/source/module_cell/unitcell_pseudo.cpp b/ABACUS.develop/source/module_cell/unitcell_pseudo.cpp
index b366c47b9c..5dadb02f07 100644
--- a/ABACUS.develop/source/module_cell/unitcell_pseudo.cpp
+++ b/ABACUS.develop/source/module_cell/unitcell_pseudo.cpp
@@ -287,7 +287,7 @@ void UnitCell_pseudo::setup_cell(
 	this->cal_nwfc();
 
 	// setup NBANDS
-	this->cal_nelec();
+	//this->cal_nelec();
 
 	this->cal_meshx();
 
@@ -307,6 +307,7 @@ void UnitCell_pseudo::setup_cell(
 // calculate total number of electrons (nelec) and default
 // number of bands (NBANDS).
 //=========================================================
+/*
 #include "../src_pw/occupy.h"
 void UnitCell_pseudo::cal_nelec(void)
 {
@@ -408,7 +409,7 @@ void UnitCell_pseudo::cal_nelec(void)
 
 	OUT(ofs_running,"NBANDS",NBANDS);
 	return;
-}
+}*/
 
 //===========================================
 // calculate the total number of local basis
diff --git a/ABACUS.develop/source/module_cell/unitcell_pseudo.h b/ABACUS.develop/source/module_cell/unitcell_pseudo.h
index c06de38fe3..6c9ec77d98 100644
--- a/ABACUS.develop/source/module_cell/unitcell_pseudo.h
+++ b/ABACUS.develop/source/module_cell/unitcell_pseudo.h
@@ -23,7 +23,7 @@ class UnitCell_pseudo : public UnitCell
 	int nmax;
 	int nmax_total;//mohan add 2009-09-10
 	int lmax_ppwf;
-	double nelec;
+	//double nelec;
 
 public: // member functions
 	UnitCell_pseudo();
@@ -49,7 +49,7 @@ class UnitCell_pseudo : public UnitCell
 	// cal_meshx	: calculate max number of mesh points in pp file
 	//================================================================
 	void cal_nwfc();
-	void cal_nelec();
+	//void cal_nelec();
 	void cal_meshx();
 	void cal_natomwfc(); 
 	void print_unitcell_pseudo(const string &fn);
diff --git a/ABACUS.develop/source/run_lcao.cpp b/ABACUS.develop/source/run_lcao.cpp
index a04bfb415a..c1fdc6dd1a 100644
--- a/ABACUS.develop/source/run_lcao.cpp
+++ b/ABACUS.develop/source/run_lcao.cpp
@@ -24,6 +24,9 @@ void Run_lcao::lcao_line(void)
     // cell relaxation. b) put NLOCAL and NBANDS as input parameters
     ucell.setup_cell( global_pseudo_dir, out, global_atom_card, ofs_running);
 
+	// setup NBANDS
+	CHR.cal_nelec();  //Yu Liu move here 2021-07-03
+
 	// mohan add 2010-09-06
 	// Yu Liu move here 2021-06-27
 	// because the number of element type
diff --git a/ABACUS.develop/source/run_pw.cpp b/ABACUS.develop/source/run_pw.cpp
index 7742f0157c..7f8fd90c28 100644
--- a/ABACUS.develop/source/run_pw.cpp
+++ b/ABACUS.develop/source/run_pw.cpp
@@ -25,6 +25,9 @@ void Run_pw::plane_wave_line(void)
     ucell.setup_cell( global_pseudo_dir, out, global_atom_card, ofs_running);
     //ucell.setup_cell( global_pseudo_dir , global_atom_card , ofs_running, NLOCAL, NBANDS);
 
+    // setup NBANDS
+	CHR.cal_nelec();  //Yu Liu move here 2021-07-03
+
     // mohan add 2010-09-06
 	// Yu Liu move here 2021-06-27
 	// because the number of element type
diff --git a/ABACUS.develop/source/src_io/berryphase.cpp b/ABACUS.develop/source/src_io/berryphase.cpp
index e2d77f22ba..de0fd6bd88 100644
--- a/ABACUS.develop/source/src_io/berryphase.cpp
+++ b/ABACUS.develop/source/src_io/berryphase.cpp
@@ -14,7 +14,7 @@ berryphase::~berryphase()
 
 void berryphase::get_occupation_bands()
 {
-	double occupied_bands = static_cast<double>(ucell.nelec/DEGSPIN);	
+	double occupied_bands = static_cast<double>(CHR.nelec/DEGSPIN);	
 	if( (occupied_bands - std::floor(occupied_bands)) > 0.0 )
 	{
 		occupied_bands = std::floor(occupied_bands) + 1.0;
diff --git a/ABACUS.develop/source/src_io/epsilon0_vasp.cpp b/ABACUS.develop/source/src_io/epsilon0_vasp.cpp
index ca5edc99ac..95d7940ef2 100644
--- a/ABACUS.develop/source/src_io/epsilon0_vasp.cpp
+++ b/ABACUS.develop/source/src_io/epsilon0_vasp.cpp
@@ -37,7 +37,7 @@ void Epsilon0_vasp::cal_epsilon0()
 	cout << "nomega = "<<nomega<<endl;
 	cout << "eta = "<<eta<<endl;
 	
-	double occupied_bands = static_cast<double>(ucell.nelec/DEGSPIN);
+	double occupied_bands = static_cast<double>(CHR.nelec/DEGSPIN);
 	if( (occupied_bands - std::floor(occupied_bands)) > 0.0 )
 	{
 		occupied_bands = std::floor(occupied_bands) + 1.0;
diff --git a/ABACUS.develop/source/src_io/istate_charge.cpp b/ABACUS.develop/source/src_io/istate_charge.cpp
index 5e6d1c2a80..1fa6af9007 100644
--- a/ABACUS.develop/source/src_io/istate_charge.cpp
+++ b/ABACUS.develop/source/src_io/istate_charge.cpp
@@ -39,13 +39,13 @@ void IState_Charge::begin(void)
 
 	// (1.2) read in LOWF_GAMMA.dat
 	OUT(ofs_running,"LOWF.allocate_flag",LOWF.get_allocate_flag());	
-	cout << " number of electrons = " << ucell.nelec << endl;
+	cout << " number of electrons = " << CHR.nelec << endl;
 
 	// mohan update 2011-03-21
 	// if ucell is odd, it's correct,
 	// if ucell is even, it's also correct.
 	// +1.0e-8 in case like (2.999999999+1)/2
-	fermi_band = static_cast<int>( (ucell.nelec+1)/2 + 1.0e-8 ) ;
+	fermi_band = static_cast<int>( (CHR.nelec+1)/2 + 1.0e-8 ) ;
 	cout << " number of occupied bands = " << fermi_band << endl;
 
 	if(mode == 1)
diff --git a/ABACUS.develop/source/src_io/istate_envelope.cpp b/ABACUS.develop/source/src_io/istate_envelope.cpp
index b156770f04..0dc95d192c 100644
--- a/ABACUS.develop/source/src_io/istate_envelope.cpp
+++ b/ABACUS.develop/source/src_io/istate_envelope.cpp
@@ -31,11 +31,11 @@ void IState_Envelope::begin(void)
 	// if ucell is odd, it's correct,
 	// if ucell is even, it's also correct.
 	// +1.0e-8 in case like (2.999999999+1)/2
-	int fermi_band = static_cast<int>( (ucell.nelec+1)/2 + 1.0e-8 ) ;
+	int fermi_band = static_cast<int>( (CHR.nelec+1)/2 + 1.0e-8 ) ;
 	int bands_below = NBANDS_ISTATE;
 	int bands_above = NBANDS_ISTATE;
 
-	cout << " number of electrons = " << ucell.nelec << endl;
+	cout << " number of electrons = " << CHR.nelec << endl;
 	cout << " number of occupied bands = " << fermi_band << endl;
 	cout << " plot band decomposed charge density below fermi surface with " 
 	<< bands_below << " bands." << endl;
diff --git a/ABACUS.develop/source/src_io/optical.cpp b/ABACUS.develop/source/src_io/optical.cpp
index 810125c6a3..fa59e3c1df 100644
--- a/ABACUS.develop/source/src_io/optical.cpp
+++ b/ABACUS.develop/source/src_io/optical.cpp
@@ -57,7 +57,7 @@ void Optical::cal_epsilon2(const int &nbands)
 
 	double range = maxe - mine;
 	int np = int(range / de) + 1; 
-	int n_occ = static_cast<int>( (ucell.nelec+1)/2 + 1.0e-8 );
+	int n_occ = static_cast<int>( (CHR.nelec+1)/2 + 1.0e-8 );
 
 	cout << " n_occ = " << n_occ << endl;
 	cout << " nbands = " << opt_nbands << endl;	
diff --git a/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp b/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
index c9e5266eed..83585246c1 100644
--- a/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
+++ b/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
@@ -133,7 +133,7 @@ void UnitCell_pseudo::bcast_unitcell_pseudo(void)
 	Parallel_Common::bcast_int( natomwfc );
 	Parallel_Common::bcast_int( lmax );
 	Parallel_Common::bcast_int( lmax_ppwf );
-	Parallel_Common::bcast_double( nelec );
+	Parallel_Common::bcast_double( CHR.nelec );
 
 	bcast_unitcell();
 }
diff --git a/ABACUS.develop/source/src_pw/charge.cpp b/ABACUS.develop/source/src_pw/charge.cpp
index 319c434c7a..dd3aa415d3 100644
--- a/ABACUS.develop/source/src_pw/charge.cpp
+++ b/ABACUS.develop/source/src_pw/charge.cpp
@@ -150,7 +150,7 @@ void Charge::renormalize_rho(void)
     const double sr = this->sum_rho();
 	ofs_warning << setprecision(15);
 	OUT(ofs_warning,"charge before normalized",sr);
-    const double normalize_factor = ucell.nelec / sr;
+    const double normalize_factor = nelec / sr;
 
 	for(int is=0; is<nspin; is++)
 	{
@@ -477,10 +477,10 @@ void Charge::atomic_rho(const int spin_number_need, double** rho_in)const		// Pe
 		ne_tot += ne[is];
 	}
 	OUT(ofs_warning,"total electron number from rho",ne_tot);
-	OUT(ofs_warning,"should be",ucell.nelec);
+	OUT(ofs_warning,"should be",nelec);
 	for(int is=0; is<spin_number_need; ++is)
 		for(int ir=0; ir<pw.nrxx; ++ir)
-			rho_in[is][ir] = rho_in[is][ir] / ne_tot * ucell.nelec;
+			rho_in[is][ir] = rho_in[is][ir] / ne_tot * nelec;
 
 	// if TWO_EFEMI, 
 	// the total magnetism will affect the calculation of
@@ -1062,3 +1062,111 @@ void Charge::init_final_scf()
 	this->allocate_rho_final_scf = true;
     return;
 }
+
+
+//=========================================================
+// calculate total number of electrons (nelec) and default
+// number of bands (NBANDS).
+//=========================================================
+#include "occupy.h"
+void Charge::cal_nelec(void)
+{
+	TITLE("UnitCell_pseudo","cal_nelec");
+	//=======================================================
+	// calculate the total number of electrons in the system
+	// if nelec <>0; use input number (setup.f90)
+	//=======================================================
+
+	ofs_running << "\n SETUP THE ELECTRONS NUMBER" << endl;
+
+	if (nelec == 0)
+	{
+		for (int it = 0; it < ucell.ntype;it++)
+		{
+			stringstream ss1, ss2;
+			ss1 << "electron number of element " << ucell.atoms[it].label;
+			const int nelec_it = ucell.atoms[it].zv * ucell.atoms[it].na;
+			nelec += nelec_it;
+			ss2 << "total electron number of element " << ucell.atoms[it].label; 
+			
+			OUT(ofs_running,ss1.str(),ucell.atoms[it].zv);
+			OUT(ofs_running,ss2.str(),nelec_it);
+		}
+	}
+
+	//OUT(ofs_running,"Total nelec",nelec);
+
+	//=======================================
+	// calculate number of bands (setup.f90)
+	//=======================================
+	double occupied_bands = static_cast<double>(nelec/DEGSPIN);	
+
+	if( (occupied_bands - std::floor(occupied_bands)) > 0.0 )
+	{
+		occupied_bands = std::floor(occupied_bands) + 1.0; //mohan fix 2012-04-16
+	}
+
+	OUT(ofs_running,"occupied bands",occupied_bands);
+	
+	// mohan add 2010-09-04
+    //cout << "nbands(ucell) = " <<NBANDS <<endl;
+	if(NBANDS==occupied_bands)
+	{
+		if( Occupy::gauss() || Occupy::tetra() )
+		{
+			WARNING_QUIT("UnitCell_pseudo::cal_nelec","for smearing, num. of bands > num. of occupied bands");
+		}
+	}
+	
+	if ( CALCULATION!="scf-sto" && CALCULATION!="relax-sto" && CALCULATION!="md-sto" ) //qianrui 2021-2-20
+	{
+	if(NBANDS == 0)
+	{
+		if(NSPIN == 1)
+		{
+			int nbands1 = static_cast<int>(occupied_bands) + 10;
+			int nbands2 = static_cast<int>(1.2 * occupied_bands);
+			NBANDS = max(nbands1, nbands2);
+		}
+		else if (NSPIN ==2 || NSPIN == 4)
+		{
+			int nbands3 = nelec + 20;
+			int nbands4 = 1.2 * nelec;
+			NBANDS = max(nbands3, nbands4);
+		}
+		AUTO_SET("NBANDS",NBANDS);
+	}
+	//else if ( CALCULATION=="scf" || CALCULATION=="md" || CALCULATION=="relax") //pengfei 2014-10-13
+	else
+	{
+		if(NBANDS < occupied_bands) WARNING_QUIT("unitcell","Too few bands!");
+		if(NBANDS < mag.get_nelup() ) 
+		{
+			OUT(ofs_running,"nelup",mag.get_nelup());
+			WARNING_QUIT("unitcell","Too few spin up bands!");
+		}
+		if(NBANDS < mag.get_neldw() )
+        {
+            WARNING_QUIT("unitcell","Too few spin down bands!");
+        }
+	}
+	}
+
+	// mohan update 2021-02-19
+    // mohan add 2011-01-5
+    if(BASIS_TYPE=="lcao" || BASIS_TYPE=="lcao_in_pw")
+    {
+        if( NBANDS > NLOCAL )
+        {
+            WARNING_QUIT("UnitCell_pseudo::cal_nwfc","NLOCAL < NBANDS");
+        }
+        else
+        {
+            OUT(ofs_running,"NLOCAL",NLOCAL);
+            OUT(ofs_running,"NBANDS",NBANDS);
+        }
+    }
+
+	OUT(ofs_running,"NBANDS",NBANDS);
+	return;
+}
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/charge.h b/ABACUS.develop/source/src_pw/charge.h
index d9cf5ded6f..968f0363b8 100644
--- a/ABACUS.develop/source/src_pw/charge.h
+++ b/ABACUS.develop/source/src_pw/charge.h
@@ -17,6 +17,7 @@ class Charge
 
 //==========================================================
 // MEMBER VARIABLES :
+// NAME : total number of electrons
 // NAME : rho (nspin,ncxyz), the charge density in real space
 // NAME : rho_save (nspin,ncxyz), for charge mixing
 // NAME : rhog, charge density in G space
@@ -25,6 +26,7 @@ class Charge
 // NAME : rhog_core [ngm], the core charge in reciprocal space
 //==========================================================
 
+    double nelec;
     double** rho;
     double** rho_save;
 
@@ -47,6 +49,8 @@ class Charge
 
     void set_rho_core(const ComplexMatrix &structure_factor);
 
+    void cal_nelec();  //calculate total number of electrons  YuLiu 2021-07-03
+
     void sum_band(void);
 
     void renormalize_rho(void);
diff --git a/ABACUS.develop/source/src_pw/charge_broyden.cpp b/ABACUS.develop/source/src_pw/charge_broyden.cpp
index f84e6b6ab9..035e4d11b1 100644
--- a/ABACUS.develop/source/src_pw/charge_broyden.cpp
+++ b/ABACUS.develop/source/src_pw/charge_broyden.cpp
@@ -117,11 +117,11 @@ void Charge_Broyden::mix_rho
 	}
 
 	Parallel_Reduce::reduce_double_pool( dr22 );
-	assert( ucell.nelec != 0);
+	assert( nelec != 0);
 	assert( ucell.omega > 0);
 	assert( pw.ncxyz > 0);
 	dr22 *= ucell.omega / static_cast<double>( pw.ncxyz );
-	dr22 /= ucell.nelec;
+	dr22 /= nelec;
 	if(test_charge)ofs_running << " dr2 from real space grid is " << dr22 << endl;
 
 	// mohan add 2011-01-22
diff --git a/ABACUS.develop/source/src_pw/electrons.cpp b/ABACUS.develop/source/src_pw/electrons.cpp
index 38f9143832..32f93d3250 100644
--- a/ABACUS.develop/source/src_pw/electrons.cpp
+++ b/ABACUS.develop/source/src_pw/electrons.cpp
@@ -254,7 +254,7 @@ void Electrons::self_consistent(const int &istep)
                 // if 'dr2 < ETHR * nelec' happen,
                 // in other word, 'dr2 < diago_error'
                 // we update ETHR.
-                diago_error = ETHR*std::max(1.0, ucell.nelec);
+                diago_error = ETHR*std::max(1.0, CHR.nelec);
             }
 
             // if converged is achieved, or the self-consistent error(dr2)
@@ -283,7 +283,7 @@ void Electrons::self_consistent(const int &istep)
 
                     // update ETHR.
                     ofs_running << " Origin ETHR = " << ETHR << endl;
-                    ETHR = 0.1 * dr2 / ucell.nelec;
+                    ETHR = 0.1 * dr2 / CHR.nelec;
                     ofs_running << " New    ETHR = " << ETHR << endl;
                     //goto first_iter_again;
                     goto scf_step;
diff --git a/ABACUS.develop/source/src_pw/magnetism.cpp b/ABACUS.develop/source/src_pw/magnetism.cpp
index d2c464ce33..0881826736 100644
--- a/ABACUS.develop/source/src_pw/magnetism.cpp
+++ b/ABACUS.develop/source/src_pw/magnetism.cpp
@@ -51,7 +51,7 @@ void Magnetism::compute_magnetization()
 		}
 		else
 		{
-			OUT(ofs_running,"nelec",ucell.nelec);
+			OUT(ofs_running,"nelec",CHR.nelec);
 		}
 
 //        cout << "\n tot_mag = " << setprecision(6) << this->tot_magnetization << " Bohr mag/cell" << endl;
@@ -90,11 +90,11 @@ double Magnetism::get_nelup(void)
 //===============================================================
 //  this type of electrons are used as "fixed" magnetization.
 //===============================================================
-		nelup = 0.5 * ucell.nelec + 0.5 * tot_magnetization;
+		nelup = 0.5 * CHR.nelec + 0.5 * tot_magnetization;
 	}
 	else
 	{
-		nelup = 0.5 * ucell.nelec;
+		nelup = 0.5 * CHR.nelec;
 	}
     return nelup;
 
@@ -104,7 +104,7 @@ double Magnetism::get_nelup(void)
 //	double nelup = 0.0;
 //	for(int i=0; i<pw.ntype; i++)
 //	{
-//		nelup += ucell.nelec * (1.0+start_magnetization[i])/2.0/pw.ntype;
+//		nelup += CHR.nelec * (1.0+start_magnetization[i])/2.0/pw.ntype;
 //	}
 //	return nelup;
 }
@@ -118,11 +118,11 @@ double Magnetism::get_neldw(void)
 //===============================================================
 //  this type of electrons are used as "fixed" magnetization.
 //===============================================================
-		neldw = 0.5 * ucell.nelec - 0.5 * tot_magnetization;
+		neldw = 0.5 * CHR.nelec - 0.5 * tot_magnetization;
 	}
 	else
 	{
-		neldw = 0.5 * ucell.nelec;
+		neldw = 0.5 * CHR.nelec;
 	}
     return neldw ;
 
diff --git a/ABACUS.develop/source/src_pw/occupy.cpp b/ABACUS.develop/source/src_pw/occupy.cpp
index 7195528b21..ac536d1559 100644
--- a/ABACUS.develop/source/src_pw/occupy.cpp
+++ b/ABACUS.develop/source/src_pw/occupy.cpp
@@ -48,7 +48,7 @@ void Occupy::calculate_weights(void)
         else
         {
 			// -1 means don't need to consider spin.
-            iweights(kv.nks, kv.wk, NBANDS, ucell.nelec, wf.ekb, en.ef, wf.wg, -1, kv.isk);
+            iweights(kv.nks, kv.wk, NBANDS, CHR.nelec, wf.ekb, en.ef, wf.wg, -1, kv.isk);
         }
     }
     else if (use_tetrahedron_method)
@@ -56,7 +56,7 @@ void Occupy::calculate_weights(void)
         WARNING_QUIT("calculate_weights","not implemented yet,coming soon!");
 //		if(my_rank == 0)
 //		{
-//			tweights(kv.nkstot, nspin, NBANDS, ucell.nelec, ntetra,tetra, wf.et, en.ef, wf.wg);
+//			tweights(kv.nkstot, nspin, NBANDS, CHR.nelec, ntetra,tetra, wf.et, en.ef, wf.wg);
 //		}
     }
     else if (use_gaussian_broadening)
@@ -75,7 +75,7 @@ void Occupy::calculate_weights(void)
 		else
 		{
 			// -1 means is no related to spin.
-        	Occupy::gweights(kv.nks, kv.wk, NBANDS, ucell.nelec, gaussian_parameter, gaussian_type,
+        	Occupy::gweights(kv.nks, kv.wk, NBANDS, CHR.nelec, gaussian_parameter, gaussian_type,
                  wf.ekb, en.ef, en.demet, wf.wg, -1, kv.isk);
 
 		}
diff --git a/ABACUS.develop/source/src_pw/sto_elec.cpp b/ABACUS.develop/source/src_pw/sto_elec.cpp
index ad20a82c99..da4e416bc9 100644
--- a/ABACUS.develop/source/src_pw/sto_elec.cpp
+++ b/ABACUS.develop/source/src_pw/sto_elec.cpp
@@ -92,15 +92,15 @@ void Stochastic_Elec::scf_stochastic(const int &istep)
         //this->update_ethr(iter);
         if(FINAL_SCF && iter==1) 
 		{
-			ETHR = 1.0e-4/ucell.nelec; //smaller ETHR than KS-DFT
+			ETHR = 1.0e-4/CHR.nelec; //smaller ETHR than KS-DFT
 		}
         else 
 		{
 			if (iter == 2)
         	{
-            	ETHR = 1.0e-4/ucell.nelec;
+            	ETHR = 1.0e-4/CHR.nelec;
         	}
-			ETHR = std::min( ETHR, 0.1*dr2/ std::max(1.0, ucell.nelec));
+			ETHR = std::min( ETHR, 0.1*dr2/ std::max(1.0, CHR.nelec));
         }
 
         
@@ -229,7 +229,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep)
 			// if 'dr2 < ETHR * nelec' happen,
 			// in other word, 'dr2 < diago_error'
 			// we update ETHR.
-			diago_error = ETHR*std::max(1.0, ucell.nelec);
+			diago_error = ETHR*std::max(1.0, CHR.nelec);
 		}
 
 		// if converged is achieved, or the self-consistent error(dr2)
@@ -266,7 +266,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep)
 
 				// update ETHR.
 				ofs_running << " Origin ETHR = " << ETHR << endl;
-				ETHR = dr2 / ucell.nelec;
+				ETHR = dr2 / CHR.nelec;
 				ofs_running << " New    ETHR = " << ETHR << endl;
 				//                  goto first_iter_again;
 			}
diff --git a/ABACUS.develop/source/src_pw/sto_iter.cpp b/ABACUS.develop/source/src_pw/sto_iter.cpp
index 750425d529..b4f2962c2b 100644
--- a/ABACUS.develop/source/src_pw/sto_iter.cpp
+++ b/ABACUS.develop/source/src_pw/sto_iter.cpp
@@ -25,7 +25,7 @@ void Stochastic_Iter::init(int &dim, int& chetype)
     nchip = STO_WF.nchip;
     stotype = STO_WF.stotype;
     //wait for init
-    targetne = ucell.nelec;
+    targetne = CHR.nelec;
     stoche.init( dim, chetype );
     stohchi.init();
     stohchi.get_GRA_index();
@@ -153,13 +153,13 @@ void Stochastic_Iter::itermu(int &iter)
     if(iter == 1)
     {
         dmu = 2;
-        th_ne = 0.1 * DRHO2 * ucell.nelec;
+        th_ne = 0.1 * DRHO2 * CHR.nelec;
         cout<<"th_ne "<<th_ne<<endl;
     }
     else
     {
         dmu = 0.1;
-        th_ne = DRHO2 * 1e-2 * ucell.nelec;
+        th_ne = DRHO2 * 1e-2 * CHR.nelec;
     }
     sumpolyval();
     mu = mu0 - dmu;
diff --git a/ABACUS.develop/source/src_pw/threshold_elec.cpp b/ABACUS.develop/source/src_pw/threshold_elec.cpp
index 88c3e836a0..680e96546d 100644
--- a/ABACUS.develop/source/src_pw/threshold_elec.cpp
+++ b/ABACUS.develop/source/src_pw/threshold_elec.cpp
@@ -23,7 +23,7 @@ void Threshold_Elec::set_ethr(void) const
     {
         if (ETHR == 0.0)
         {
-            ETHR = 0.1 * std::min(1.0e-2, DRHO2 / ucell.nelec);
+            ETHR = 0.1 * std::min(1.0e-2, DRHO2 / CHR.nelec);
         }
     }
     //=================
@@ -76,12 +76,12 @@ void Threshold_Elec::update_ethr(const int &iter)
 		//----------------------------
 		if(BASIS_TYPE=="lcao")
 		{
-			ETHR = std::min( ETHR, 0.01*dr2/ std::max(1.0, ucell.nelec));
+			ETHR = std::min( ETHR, 0.01*dr2/ std::max(1.0, CHR.nelec));
 		}
 		// mohan update 2009-09-04
 		else
 		{
-			ETHR = std::min( ETHR, 0.1*dr2/ std::max(1.0, ucell.nelec));
+			ETHR = std::min( ETHR, 0.1*dr2/ std::max(1.0, CHR.nelec));
 			//cout << " new ethr = " << ETHR << endl;
 		}
 
diff --git a/ABACUS.develop/tests/127_PW_15_PK_AF/result.ref b/ABACUS.develop/tests/127_PW_15_PK_AF/result.ref
index 6f7d9a47d1..c9452b587f 100644
--- a/ABACUS.develop/tests/127_PW_15_PK_AF/result.ref
+++ b/ABACUS.develop/tests/127_PW_15_PK_AF/result.ref
@@ -1,3 +1,3 @@
-etotref -6141.0777981366600216
-etotperatomref -3070.5388990683
-totaltimeref 6.160
+etotref -6141.0777980752291114
+etotperatomref -3070.5388990376
+totaltimeref 6.457
diff --git a/ABACUS.develop/tests/150_PW_15_CR_VDW3/result.ref b/ABACUS.develop/tests/150_PW_15_CR_VDW3/result.ref
index 5e1a4b70bd..f612ece177 100644
--- a/ABACUS.develop/tests/150_PW_15_CR_VDW3/result.ref
+++ b/ABACUS.develop/tests/150_PW_15_CR_VDW3/result.ref
@@ -2,4 +2,4 @@ etotref -4009.5574820816177635
 etotperatomref -2004.7787410408
 totalforceref 0.412136
 totalstressref 5111.513024
-totaltimeref +2.398
+totaltimeref +2.443