Refactor : decouple charge and charge_mixing (#1488)

* fix : latname in unit test elecstate_lcao_test

* refactor : removed glbalc::hm and globalc::ufft

* fix : removed 'use_fft.h' in headers

* fix : unit tests in module_elecstate

* fix : comment out lobalc::hm from cuda codes

* differentiate different gint in title and timer

* checks relax params only for relax & cell-relax

* keep the total gint timer

* refactor : separate charge_mixing from charge

* remove 'kerker' and 'pulay-kerker', and use gg0
to control whether kerker is used

* clean up lcao_hamilt

* modify treatment of new_e_iteration

* fix : electrons.cu and electrons_hip.cpp

* set mixing parameters to be private

* fix test cases

* fix : inclusion of parallel_common in input_update.cpp

* removed some unused code in charge_broyden

* remove obsolete electrons.cpp; make pulay-related variables private

* fix makefile

* merge charge_pulay and charge_broyden into charge_mixing

Co-authored-by: wenfei-li <liwenfei@gmail.com>

docs/advanced/input_files/input-main.md

@@ -494,11 +494,9 @@ calculations.
 ### mixing_type
 
 - **Type**: String
-- **Description**: Charge mixing methods.
+- **Description**: Charge mixing methods. We offer the following 3 options:
   - plain: Just simple mixing.
-  - kerker: Use kerker method, which is the mixing method in G space.
   - pulay: Standard Pulay method.
-  - pulay-kerker:
   - broyden: Broyden method.
 - **Default**: pulay
 
@@ -517,8 +515,8 @@ calculations.
 ### mixing_gg0
 
 - **Type**: Real
-- **Description**: used in pulay-kerker mixing method
-- **Default**: 1.5
+- **Description**: When set to a positive number, the high frequency wave vectors will be suppressed by multiplying a scaling factor $\frac{k^2}{k^2+gg0^2}$; if set to 0, then no Kerker scaling is performed.
+- **Default**: 0.0
 
 ### gamma_only
 

docs/advanced/scf/converge.md

@@ -8,7 +8,7 @@ By mixing the electron density with that obtained from previous steps, numerical
 
 For each of the mixing types, we also provide variables for controlling relevant parameters, including `mixing_beta`, `mixing_ndim`, and `mixing_gg0`.
 
-The default choice is `pulay`, which should work fine in most cases. If convergence issue arises in metallic systems, inclusion of Kerker preconditioning may be helpful, namely, to use `kerker` or `pulay-kerker` mixing methods.
+The default choice is `pulay`, which should work fine in most cases. If convergence issue arises in metallic systems, inclusion of Kerker preconditioning may be helpful, which can be achieved by setting [mixing_gg0](../input_files/input-main.md#mixing_gg0) to be a positive number. For the default pulay method, a choice of 1.5 might be a good start.
 
 A large `mixing_beta` means a larger change in electron density for each SCF step. For well-behaved systems, a larger `mixing_beta` leads to faster convergence. However, for some difficult cases, a smaller `mixing_beta` is preferred to avoid numerical instabilities.
 

source/Makefile.Objects

@@ -552,8 +552,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     wf_atomic.o\
     wf_igk.o\
     xc_3.o\
-    hamilt.o\
-    electrons.o
+    hamilt.o\  
 
 OBJS_VDW=vdw.o\
     vdwd2_parameters.o\

source/input.cpp

@@ -1014,11 +1014,6 @@ bool Input::Read(const std::string &fn)
         else if (strcmp("mixing_type", word) == 0)
         {
             read_value(ifs, mixing_mode);
-            // 2015-06-15, xiaohui
-            if (mixing_mode == "pulay-kerker")
-            {
-                mixing_gg0 = 1.5;
-            }
         }
         else if (strcmp("mixing_beta", word) == 0)
         {
@@ -3065,12 +3060,6 @@ void Input::Check(void)
         }
     }
 
-    // 2015-06-15, xiaohui
-    if (mixing_mode == "pulay" && mixing_gg0 > 0.0)
-    {
-        ModuleBase::WARNING("Input", "To use pulay-kerker mixing method, please set mixing_type=pulay-kerker");
-    }
-
     if (berry_phase)
     {
         if (basis_type == "pw")

source/input_conv.cpp

@@ -471,7 +471,7 @@ void Input_Conv::Convert(void)
     //----------------------------------------------------------
     // charge mixing(3/3)
     //----------------------------------------------------------
-    GlobalC::CHR.set_mixing(INPUT.mixing_mode,
+    GlobalC::CHR_MIX.set_mixing(INPUT.mixing_mode,
                             INPUT.mixing_beta,
                             INPUT.mixing_ndim,
                             INPUT.mixing_gg0); // mohan modify 2014-09-27, add mixing_gg0

source/input_update.cpp

@@ -9,6 +9,7 @@
 #include "input.h"
 #include "module_relax/relax_old/ions_move_basic.h"
 #include "src_io/optical.h"
+#include "src_parallel/parallel_common.h"
 #ifdef __LCAO
 #include "src_lcao/FORCE_STRESS.h"
 #include "src_lcao/local_orbital_charge.h"
@@ -160,10 +161,18 @@ bool Update_input::Read(const std::string &fn)
         else if (strcmp("mixing_beta", word) == 0)
         {
             read_value(ifs, mixing_beta);
-			if(mixing_beta!=GlobalC::CHR.mixing_beta)
+#ifdef __MPI
+            Parallel_Common::bcast_double( mixing_beta );
+#endif
+			if(mixing_beta!=GlobalC::CHR_MIX.get_mixing_beta())
 			{
-				this->change(GlobalV::ofs_warning,"mixing_beta",GlobalC::CHR.mixing_beta,mixing_beta);
-    			GlobalC::CHR.mixing_beta = mixing_beta;
+                double mixing_beta_old = GlobalC::CHR_MIX.get_mixing_beta();
+				this->change(GlobalV::ofs_warning,"mixing_beta",mixing_beta_old,mixing_beta);
+    			GlobalC::CHR_MIX.set_mixing(
+                    GlobalC::CHR_MIX.get_mixing_mode(),
+                    mixing_beta,
+                    GlobalC::CHR_MIX.get_mixing_ndim(),
+                    GlobalC::CHR_MIX.get_mixing_gg0());
 			}
         }
 		// 8
@@ -274,7 +283,6 @@ bool Update_input::Read(const std::string &fn)
     return true;
 }//end read_parameters
 
-#include "src_parallel/parallel_common.h"
 #ifdef __MPI
 void Update_input::Bcast()
 {
@@ -288,7 +296,6 @@ void Update_input::Bcast()
     Parallel_Common::bcast_double( GlobalV::SCF_THR );
     Parallel_Common::bcast_int( GlobalV::SCF_NMAX );
     Parallel_Common::bcast_int( GlobalV::RELAX_NMAX );
-    Parallel_Common::bcast_double( GlobalC::CHR.mixing_beta );
     Parallel_Common::bcast_int( GlobalC::en.printe );
     Parallel_Common::bcast_string( GlobalC::pot.chg_extrap );//xiaohui modify 2015-02-01
     Parallel_Common::bcast_int( GlobalC::CHR.out_chg );

source/module_base/test/memory_test.cpp

@@ -113,8 +113,8 @@ TEST_F(MemoryTest, printall)
 {
 	ofs.open("tmp");
 	// total memory is an internal parameter and added inside the class Memory
-	ModuleBase::Memory::record("Charge_Pulay","Rrho",1024*1024,"ModuleBase::Vector3<double>");
-	ModuleBase::Memory::record("Charge_pulay","drho",1024*1024,"AtomLink");
+	ModuleBase::Memory::record("Charge_Mixing","Rrho",1024*1024,"ModuleBase::Vector3<double>");
+	ModuleBase::Memory::record("Charge_Mixing","drho",1024*1024,"AtomLink");
 	ModuleBase::Memory::record("wavefunc","evc",1024*1024,"float");
 	ModuleBase::Memory::print_all(ofs);
 	ofs.close();

source/module_elecstate/test/updaterhok_pw_test.h

@@ -5,7 +5,7 @@
 #include "src_pw/magnetism.h"
 #include "src_pw/wf_atomic.h"
 #include "src_pw/wavefunc.h"
-#include "src_pw/charge_broyden.h"
+#include "src_pw/charge_mixing.h"
 #include "src_pw/potential.h"
 #include "module_cell/atom_pseudo.h"
 #include "module_cell/atom_spec.h"
@@ -44,10 +44,6 @@ Atom_pseudo::Atom_pseudo(){}
 Atom_pseudo::~Atom_pseudo(){}
 Charge_Mixing::Charge_Mixing(){}
 Charge_Mixing::~Charge_Mixing(){}
-Charge_Pulay::Charge_Pulay(){}
-Charge_Pulay::~Charge_Pulay(){}
-Charge_Broyden::Charge_Broyden(){}
-Charge_Broyden::~Charge_Broyden(){}
 Potential::Potential(){}
 Potential::~Potential(){}
 InfoNonlocal::InfoNonlocal(){}

source/module_esolver/esolver_ks.cpp

@@ -11,7 +11,7 @@
 //--------------Temporary----------------
 #include "../module_base/global_variable.h"
 #include "../src_pw/global.h"
-#include "../src_pw/charge_broyden.h"
+#include "../src_pw/charge_mixing.h"
 #include "../module_base/timer.h"
 //---------------------------------------
 
@@ -206,7 +206,8 @@ namespace ModuleESolver
                 {
                     // double drho = this->estate.caldr2(); 
                     // EState should be used after it is constructed.
-                    drho = GlobalC::CHR.get_drho();
+                    drho = GlobalC::CHR_MIX.get_drho(GlobalC::CHR.rho, GlobalC::CHR.rho_save,
+                        GlobalC::CHR.rhog, GlobalC::CHR.rhog_save, GlobalC::CHR.nelec);
                     double hsolver_error = 0.0;
                     if (firstscf)
                     {
@@ -217,7 +218,8 @@ namespace ModuleESolver
                         {
                             diag_ethr = this->phsol->reset_diagethr(GlobalV::ofs_running, hsolver_error, drho);
                             this->hamilt2density(istep, iter, diag_ethr);
-                            drho = GlobalC::CHR.get_drho();
+                            drho = GlobalC::CHR_MIX.get_drho(GlobalC::CHR.rho, GlobalC::CHR.rho_save,
+                                GlobalC::CHR.rhog, GlobalC::CHR.rhog_save, GlobalC::CHR.nelec);
                             hsolver_error = this->phsol->cal_hsolerror();
                         }
                     }
@@ -233,7 +235,7 @@ namespace ModuleESolver
                     {
                         //charge mixing
                         //conv_elec = this->estate.mix_rho();
-                        GlobalC::CHR.mix_rho(iter);
+                        GlobalC::CHR_MIX.mix_rho(iter, GlobalC::CHR.rho, GlobalC::CHR.rho_save, GlobalC::CHR.rhog, GlobalC::CHR.rhog_save);
                     }
                 }
 #ifdef __MPI

source/module_esolver/esolver_ks_lcao.cpp

@@ -161,7 +161,7 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell)
     }
     else
     {
-        this->pelec = new elecstate::ElecStateLCAO((Charge*)(&(GlobalC::CHR)),
+        this->pelec = new elecstate::ElecStateLCAO(&(GlobalC::CHR),
                                                    &(GlobalC::kv),
                                                    GlobalC::kv.nks,
                                                    GlobalV::NBANDS,
@@ -320,21 +320,7 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter)
 
     // mohan add 2010-07-16
     // used for pulay mixing.
-    if (iter == 1)
-    {
-        GlobalC::CHR.set_new_e_iteration(true);
-    }
-    else
-    {
-        GlobalC::CHR.set_new_e_iteration(false);
-    }
-
-    if (GlobalV::FINAL_SCF && iter == 1)
-    {
-        GlobalC::CHR.irstep = 0;
-        GlobalC::CHR.idstep = 0;
-        GlobalC::CHR.totstep = 0;
-    }
+    if (iter == 1) GlobalC::CHR_MIX.reset(GlobalV::FINAL_SCF);
 
     // mohan update 2012-06-05
     GlobalC::en.calculate_harris(1);

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -108,7 +108,7 @@ void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell)
     }
     else
     {
-        this->pelec = new elecstate::ElecStateLCAO_TDDFT((Charge*)(&(GlobalC::CHR)),
+        this->pelec = new elecstate::ElecStateLCAO_TDDFT(   &(GlobalC::CHR),
                                                             &(GlobalC::kv),
                                                             GlobalC::kv.nks,
                                                             GlobalV::NBANDS,
@@ -127,21 +127,7 @@ void ESolver_KS_LCAO_TDDFT::eachiterinit(const int istep, const int iter)
 
     // mohan add 2010-07-16
     // used for pulay mixing.
-    if (iter == 1)
-    {
-        GlobalC::CHR.set_new_e_iteration(true);
-    }
-    else
-    {
-        GlobalC::CHR.set_new_e_iteration(false);
-    }
-
-    if (GlobalV::FINAL_SCF && iter == 1)
-    {
-        GlobalC::CHR.irstep = 0;
-        GlobalC::CHR.idstep = 0;
-        GlobalC::CHR.totstep = 0;
-    }
+    if (iter == 1) GlobalC::CHR_MIX.reset(GlobalV::FINAL_SCF);
 
     // mohan update 2012-06-05
     GlobalC::en.calculate_harris(1);

source/module_esolver/esolver_ks_pw.cpp

@@ -11,7 +11,6 @@
 #include "../src_pw/symmetry_rho.h"
 #include "../src_io/print_info.h"
 #include "../src_pw/H_Ewald_pw.h"
-#include "../src_pw/electrons.h"
 #include "../src_pw/occupy.h"
 #include "../src_io/chi0_standard.h"
 #include "../src_io/chi0_hilbert.h"
@@ -131,7 +130,7 @@ namespace ModuleESolver
         //init ElecState,
         if(this->pelec == nullptr)
         {
-            this->pelec = new elecstate::ElecStatePW( GlobalC::wfcpw, (Charge*)(&(GlobalC::CHR)), (K_Vectors*)(&(GlobalC::kv)), GlobalV::NBANDS);
+            this->pelec = new elecstate::ElecStatePW( GlobalC::wfcpw, &(GlobalC::CHR), (K_Vectors*)(&(GlobalC::kv)), GlobalV::NBANDS);
         }
         //init HSolver
         if(this->phsol == nullptr)
@@ -266,15 +265,7 @@ namespace ModuleESolver
     void ESolver_KS_PW::eachiterinit(const int istep, const int iter)
     {
         // mohan add 2010-07-16
-        if (iter == 1) GlobalC::CHR.set_new_e_iteration(true);
-        else GlobalC::CHR.set_new_e_iteration(false);
-
-        if (GlobalV::FINAL_SCF && iter == 1)
-        {
-            GlobalC::CHR.irstep = 0;
-            GlobalC::CHR.idstep = 0;
-            GlobalC::CHR.totstep = 0;
-        }
+        if (iter == 1) GlobalC::CHR_MIX.reset(GlobalV::FINAL_SCF);
 
         // mohan move harris functional to here, 2012-06-05
         // use 'rho(in)' and 'v_h and v_xc'(in)

source/module_esolver/esolver_of.cpp

@@ -13,7 +13,6 @@
 #include "../src_pw/symmetry_rho.h"
 #include "../src_io/print_info.h"
 #include "../src_pw/H_Ewald_pw.h"
-#include "../src_pw/electrons.h"
 //-----force-------------------
 #include "../src_pw/forces.h"
 //-----stress------------------

source/module_esolver/esolver_sdft_pw.cpp

@@ -37,7 +37,7 @@ void ESolver_SDFT_PW::Init(Input &inp, UnitCell &ucell)
     ESolver_KS::Init(inp,ucell);
 
 
-    this->pelec = new elecstate::ElecStatePW_SDFT( GlobalC::wfcpw, (Charge*)(&(GlobalC::CHR)), (K_Vectors*)(&(GlobalC::kv)), GlobalV::NBANDS);
+    this->pelec = new elecstate::ElecStatePW_SDFT( GlobalC::wfcpw, &(GlobalC::CHR), (K_Vectors*)(&(GlobalC::kv)), GlobalV::NBANDS);
 
     // Inititlize the charge density.
     this->pelec->charge->allocate(GlobalV::NSPIN, GlobalC::rhopw->nrxx, GlobalC::rhopw->npw);

source/module_gint/gint_gamma.h

@@ -31,7 +31,7 @@ class Gint_Gamma : public Gint
     // there is an additional step in calculating vlocal for gamma point
     // namely the redistribution of Hamiltonian from grid to 2D block format
     // hence we have an additional layer outside the unified interface
-	void cal_vlocal(Gint_inout *inout);
+	void cal_vlocal(Gint_inout *inout, const bool new_e_iteration);
 
     //------------------------------------------------------
     // in gint_gamma_env.cpp 
@@ -54,7 +54,7 @@ class Gint_Gamma : public Gint
     //------------------------------------------------------
     // method for redistributing the Hamiltonian
     // from grid to 2D format
-    void vl_grid_to_2D(const int lgd, LCAO_Matrix& lm);
+    void vl_grid_to_2D(const int lgd, LCAO_Matrix& lm, const bool new_e_iteration);
 
     ///===============================
     /// Use MPI_Alltoallv to convert a grid distributed matrix

source/module_gint/gint_gamma_vl.cpp

@@ -27,7 +27,7 @@ extern "C"
     void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol);
 }
 
-void Gint_Gamma::cal_vlocal(Gint_inout *inout)
+void Gint_Gamma::cal_vlocal(Gint_inout *inout, const bool new_e_iteration)
 {
 	const int max_size = GlobalC::GridT.max_atom;
 	const int lgd = GlobalC::GridT.lgd;
@@ -42,7 +42,7 @@ void Gint_Gamma::cal_vlocal(Gint_inout *inout)
 
         this->cal_gint(inout);
 
-		this->vl_grid_to_2D(lgd,inout->lm[0]);
+		this->vl_grid_to_2D(lgd,inout->lm[0], new_e_iteration);
 
 		if (max_size >0 && lgd > 0) delete[] pvpR_grid;
 	}
@@ -192,11 +192,11 @@ inline int setBufferParameter(
 }
 #endif
 
-void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm)
+void Gint_Gamma::vl_grid_to_2D(const int lgd_now, LCAO_Matrix &lm, const bool new_e_iteration)
 {
     // setup send buffer and receive buffer size
     // OUT(GlobalV::ofs_running, "Start transforming vlocal from grid distribute to 2D block");
-    if(GlobalC::CHR.get_new_e_iteration())
+    if(new_e_iteration)
     {
         ModuleBase::timer::tick("Gint_Gamma","distri_vl_index");
         #ifdef __MPI

source/module_hamilt/ks_lcao/operator_lcao.h

@@ -76,6 +76,9 @@ class OperatorLCAO : public Operator<T>
     // Hamiltonian matrix which are stored in LCAO_Matrix and calculated in OperatorLCAO
     LCAO_Matrix* LM = nullptr;
 
+    protected:
+    bool new_e_iteration = true;
+
     private:
     void get_hs_pointers();
 
@@ -86,7 +89,6 @@ class OperatorLCAO : public Operator<T>
 
     //only used for Gamma_only case
     bool allocated_smatrix = false;
-    bool new_e_iteration = true;
 
     //fixed HR matrix folding to HK
     void folding_fixed(const int ik);

source/module_hamilt/ks_lcao/veff_lcao.cpp

@@ -67,12 +67,12 @@ void Veff<OperatorLCAO<double>>::contributeHk(int ik)
     if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5)
     {
         Gint_inout inout(GlobalC::pot.vr_eff1, GlobalC::pot.vofk_eff1, this->LM, Gint_Tools::job_type::vlocal_meta);
-        this->GG->cal_vlocal(&inout);
+        this->GG->cal_vlocal(&inout, new_e_iteration);
     }
     else
     {
         Gint_inout inout(GlobalC::pot.vr_eff1, this->LM, Gint_Tools::job_type::vlocal);
-        this->GG->cal_vlocal(&inout);
+        this->GG->cal_vlocal(&inout, new_e_iteration);
     }
 
     ModuleBase::timer::tick("Veff", "contributeHk");

source/module_hsolver/hsolver_pw_sdft.cpp

@@ -113,7 +113,7 @@ namespace hsolver
 			Symmetry_rho srho;
 			for(int is=0; is < GlobalV::NSPIN; is++)
 			{
-				srho.begin(is, GlobalC::CHR,GlobalC::rhopw, GlobalC::Pgrid, GlobalC::symm);
+				srho.begin(is, GlobalC::CHR, GlobalC::rhopw, GlobalC::Pgrid, GlobalC::symm);
 			}
 		}
 		else