Feature: optimized default settings of mixing_beta and mixing_gg0 (#2036)

dyzheng · dyzheng · web-flow · commit caa372a3927f · 2023-03-10T13:35:51.000+08:00
* Feature: optimized default settings of mixing_beta and mixing_gg0

* Fix: error in CI

---------

Co-authored-by: dyzheng &lt;zhengdy@bjaisi.com&gt;
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -520,8 +520,14 @@ calculations.
 ### mixing_beta
 
 - **Type**: Real
-- **Description**: mixing parameter: 0 means no new charge
-- **Default**: 0.7
+- **Description**: mixing parameter. We recommand the following options:
+  - default: -10.0 means program will autoset **mixing_beta** and **mixing_gg0** before charge mixing method starts. 
+    - Default values of transition metal system are **mixing_beta=0.2** and **mixing_gg0=1.5**; 
+    - Default values of metal system (bandgap <= 1.0 eV) are **mixing_beta=0.2** and **mixing_gg0=0.0**;
+    - Default values of other systems (bandgap > 1.0eV) are **mixing_beta=0.7** and **mixing_gg0=0.0**.
+  - 0: keep charge density unchanged, usually used for restarting with **init_chg=file** or testing.
+  - 0.1 or less: if convergence of SCF calculation is difficult to reach, please try **0 < mixing_beta < 0.1**  
+- **Default**: -10.0
 
 ### mixing_ndim
 
diff --git a/source/module_base/element_elec_config.h b/source/module_base/element_elec_config.h
@@ -251,6 +251,128 @@ const std::map<std::string, int> MinZval
 {"Og", 32}
 };
 
+const std::map<std::string, int> IsTransMetal
+= {
+{"H",  0},   
+{"He", 0},  
+{"Li", 0},  
+{"Be", 0},  
+{"B",  0},   
+{"C",  0},   
+{"N",  0},   
+{"O",  0},   
+{"F",  0},
+{"Ne", 0}, 
+{"Na", 0}, 
+{"Mg", 0}, 
+{"Al", 0}, 
+{"Si", 0}, 
+{"P",  0},  
+{"S",  0},  
+{"Cl", 0}, 
+{"Ar", 0},
+{"K",  0},  
+{"Ca", 0}, 
+{"Sc", 1}, 
+{"Ti", 1}, 
+{"V",  1},  
+{"Cr", 1}, 
+{"Mn", 1}, 
+{"Fe", 1}, 
+{"Co", 1},
+{"Ni", 1}, 
+{"Cu", 1}, 
+{"Zn", 1}, 
+{"Ga", 0}, 
+{"Ge", 0}, 
+{"As", 0}, 
+{"Se", 0}, 
+{"Br", 0}, 
+{"Kr", 0},
+{"Rb", 0}, 
+{"Sr", 0}, 
+{"Y",  1}, 
+{"Zr", 1}, 
+{"Nb", 1}, 
+{"Mo", 1}, 
+{"Tc", 1}, 
+{"Ru", 1}, 
+{"Rh", 1},
+{"Pd", 1}, 
+{"Ag", 1}, 
+{"Cd", 1}, 
+{"In", 0}, 
+{"Sn", 0}, 
+{"Sb", 0}, 
+{"Te", 0}, 
+{"I",  0},  
+{"Xe", 0},
+{"Cs", 0}, 
+{"Ba", 0}, 
+{"La", 1}, 
+{"Ce", 1}, 
+{"Pr", 1}, 
+{"Nd", 1}, 
+{"Pm", 1}, 
+{"Sm", 1}, 
+{"Eu", 1},
+{"Gd", 1}, 
+{"Tb", 1}, 
+{"Dy", 1}, 
+{"Ho", 1}, 
+{"Er", 1}, 
+{"Tm", 1}, 
+{"Yb", 1}, 
+{"Lu", 1}, 
+{"Hf", 1},
+{"Ta", 1}, 
+{"W",  1},  
+{"Re", 1}, 
+{"Os", 1}, 
+{"Ir", 1}, 
+{"Pt", 1}, 
+{"Au", 1}, 
+{"Hg", 1}, 
+{"Tl", 0},
+{"Pb", 0}, 
+{"Bi", 0}, 
+{"Po", 0}, 
+{"At", 0}, 
+{"Rn", 0}, 
+{"Fr", 0},
+{"Ra", 0},
+{"Ac", 1},
+{"Th", 1},
+{"Pa", 1},
+{"U" , 1},
+{"Np", 1},
+{"Pu", 1},
+{"Am", 1},
+{"Cm", 1},
+{"Bk", 1},
+{"Cf", 1},
+{"Es", 1},
+{"Fm", 1},
+{"Md", 1},
+{"No", 1},
+{"Lr", 1},
+{"Rf", 1},
+{"Db", 1},
+{"Sg", 1},
+{"Bh", 1},
+{"Hs", 1},
+{"Mt", 1},
+{"Ds", 1},
+{"Rg", 1},
+{"Cn", 1},
+{"Nh", 0},
+{"Fl", 0},
+{"Mc", 0},
+{"Lv", 0},
+{"Ts", 0},
+{"Og", 0}
+};
+
 }
 
 #endif
diff --git a/source/module_elecstate/energy.cpp b/source/module_elecstate/energy.cpp
@@ -511,6 +511,11 @@ void energy::cal_converged(elecstate::ElecState* pelec)
 
 void energy::cal_bandgap(const elecstate::ElecState* pelec)
 {
+	if(pelec->ekb.nr == 0 || pelec->ekb.nc == 0)
+	{//which means no homo and no lumo
+		this->bandgap = 0.0;
+		return;
+	}
 	int nbands = GlobalV::NBANDS;
 	int nks = GlobalC::kv.nks;
 	double homo = pelec->ekb(0,0);
@@ -534,6 +539,12 @@ void energy::cal_bandgap(const elecstate::ElecState* pelec)
 
 void energy::cal_bandgap_updw(const elecstate::ElecState* pelec)
 {
+	if(pelec->ekb.nr == 0 || pelec->ekb.nc == 0)
+	{//which means no homo and no lumo
+		this->bandgap_up = 0.0;
+		this->bandgap_dw = 0.0;
+		return;
+	}
 	int nbands = GlobalV::NBANDS;
 	int nks = GlobalC::kv.nks;
 	double homo_up = pelec->ekb(0,0);
diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -3,6 +3,7 @@
 #include "module_base/inverse_matrix.h"
 #include "module_base/parallel_reduce.h"
 #include "module_base/timer.h"
+#include "module_base/element_elec_config.h"
 
 Charge_Mixing::Charge_Mixing()
 {
@@ -47,6 +48,60 @@ void Charge_Mixing::set_mixing
     return;
 }
 
+void Charge_Mixing::need_auto_set()
+{
+	this->autoset = true;
+}
+
+void Charge_Mixing::auto_set(const double& bandgap_in, const UnitCell& ucell_)
+{
+	//auto set parameters once
+	if(!this->autoset)
+	{
+		return;
+	}
+	else {
+		this->autoset = false;
+	}
+	GlobalV::ofs_running<<"--------------AUTO-SET---------------"<<std::endl;
+	//0.2 for metal and 0.7 for others
+	if(bandgap_in * ModuleBase::Ry_to_eV < 1.0)
+	{
+		this->mixing_beta = 0.2;
+	}
+	else
+	{
+		this->mixing_beta = 0.7;
+	}
+	GlobalV::ofs_running<<"      Autoset mixing_beta to "<<this->mixing_beta<<std::endl;
+
+	bool has_trans_metal = false;
+	// find elements of cell
+	for(int it=0;it<ucell_.ntype;it++)
+	{
+		if(ModuleBase::IsTransMetal.find(ucell_.atoms[it].ncpp.psd) != ModuleBase::IsTransMetal.end())
+		{
+			if(ModuleBase::IsTransMetal.at(ucell_.atoms[it].ncpp.psd)) 
+			{
+				has_trans_metal = true;
+			}
+		}
+	}
+	// auto set kerker mixing for trans metal system
+	if(has_trans_metal)
+	{
+		this->mixing_gg0 = 1.5;
+	}
+	else
+	{
+		this->mixing_gg0 = 0.0;
+	}
+	GlobalV::ofs_running<<"      Autoset mixing_gg0 to "<<this->mixing_gg0<<std::endl;
+	GlobalV::ofs_running<<"-------------------------------------"<<std::endl;
+	//auto set for inhomogeneous system
+
+}
+
 double Charge_Mixing::get_drho(Charge* chr, const double nelec)
 {
 	for (int is=0; is<GlobalV::NSPIN; is++)
diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h
@@ -13,6 +13,7 @@
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
 #include "module_base/matrix.h"
+#include "module_cell/unitcell.h"
 #include "charge.h"
 class Charge_Mixing
 {
@@ -32,6 +33,9 @@ class Charge_Mixing
 		const bool &mixing_tau_in
     );//mohan add mixing_gg0_in 2014-09-27
 
+	void need_auto_set();
+	void auto_set(const double& bandgap_in, const UnitCell& ucell_);
+
 	double get_drho(Charge* chr, const double nelec);
 
     void mix_rho(const int &iter, Charge* chr);// mix rho
@@ -131,6 +135,9 @@ class Charge_Mixing
 
 	std::complex<double>*** dF; // dF(i) = rhog(i) - rhog_save(i), (GlobalV::NSPIN, rstep, rhopw->npw)
 	std::complex<double>*** dn; // dn(i) = rhog(i+1) - rhog(i), (GlobalV::NSPIN, rstep, rhopw->npw)
+
+	private: 
+	bool autoset = false;
 };
 
 #endif
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
@@ -219,9 +219,26 @@ namespace ModuleESolver
                     }
                     else
                     {
-                        //charge mixing
+                        //----------charge mixing---------------
+                        //before first calling mix_rho(), bandgap and cell structure can be analyzed to get better default parameters
+                        if(iter == 1)
+                        {
+                            double bandgap_for_autoset = 0.0;
+                            if (!GlobalV::TWO_EFERMI)
+                            {
+                                GlobalC::en.cal_bandgap(this->pelec);
+                                bandgap_for_autoset = GlobalC::en.bandgap;
+                            }
+                            else
+                            {
+                                GlobalC::en.cal_bandgap_updw(this->pelec);
+                                bandgap_for_autoset = std::min(GlobalC::en.bandgap_up, GlobalC::en.bandgap_dw);
+                            }
+                            GlobalC::CHR_MIX.auto_set(bandgap_for_autoset, GlobalC::ucell);
+                        }
                         //conv_elec = this->estate.mix_rho();
                         GlobalC::CHR_MIX.mix_rho(iter, pelec->charge);
+                        //----------charge mixing done-----------
                     }
                 }
 #ifdef __MPI
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
@@ -264,7 +264,7 @@ void Input::Default(void)
     //  charge mixing
     //----------------------------------------------------------
     mixing_mode = "pulay";
-    mixing_beta = 0.7;
+    mixing_beta = -10.0;
     mixing_ndim = 8;
     mixing_gg0 = 0.00; // used in kerker method. mohan add 2014-09-27
     mixing_tau = false;
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
@@ -493,6 +493,15 @@ void Input_Conv::Convert(void)
                                 INPUT.mixing_ndim,
                                 INPUT.mixing_gg0,
                                 INPUT.mixing_tau); // mohan modify 2014-09-27, add mixing_gg0
+    //using bandgap to auto set mixing_beta
+    if(std::abs(INPUT.mixing_beta + 10.0) < 1e-6)
+    {
+        GlobalC::CHR_MIX.need_auto_set();
+    }
+    else if(INPUT.mixing_beta > 1.0 || INPUT.mixing_beta<0.0)
+    {
+        ModuleBase::WARNING("INPUT", "You'd better set mixing_beta to [0.0, 1.0]!");
+    }
 
     //----------------------------------------------------------
     // iteration
diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp
@@ -137,7 +137,7 @@ TEST_F(InputTest, Default)
         EXPECT_EQ(INPUT.smearing_method,"fixed");
         EXPECT_DOUBLE_EQ(INPUT.smearing_sigma,0.01);
         EXPECT_EQ(INPUT.mixing_mode,"pulay");
-        EXPECT_DOUBLE_EQ(INPUT.mixing_beta,0.7);
+        EXPECT_DOUBLE_EQ(INPUT.mixing_beta,-10.0);
         EXPECT_EQ(INPUT.mixing_ndim,8);
         EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,0.00);
         EXPECT_EQ(INPUT.init_wfc,"atomic");
