move cel_nelec to charge.cpp

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;

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

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);

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

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

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;

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;

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)

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;

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;	

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();
 }

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;
+}

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);

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

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;

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 ;