Feature: enable restart charge density mixing during SCF

docs/advanced/input_files/input-main.md

@@ -1003,6 +1003,13 @@ We recommend the following options:
   For systems that are difficult to converge, one could try increasing the value of 'mixing_ndim' to enhance the stability of the self-consistent field (SCF) calculation.
 - **Default**: 8
 
+### mixing_restart
+
+- **Type**: Integer
+- **Description**: At `mixing_restart`-th iteration, SCF will restart by using output charge density from perivos iteration as input charge density directly, and start a new mixing. `mixing_restart=0|1` means SCF starts from scratch.
+
+- **Default**: 0
+
 ### mixing_gg0
 
 - **Type**: Real

source/module_base/global_variable.cpp

@@ -248,6 +248,7 @@ std::string of_kernel_file = "WTkernel.txt";
 std::string MIXING_MODE = "broyden";
 double MIXING_BETA = 0.7;
 int MIXING_NDIM = 8;
+int MIXING_RESTART = 0;
 double MIXING_GG0 = 1.00;
 double MIXING_BETA_MAG = 1.6;
 double MIXING_GG0_MAG = 1.00;

source/module_base/global_variable.h

@@ -277,6 +277,7 @@ extern std::string of_kernel_file; // The name of WT kernel file.
 extern std::string MIXING_MODE;
 extern double MIXING_BETA;
 extern int MIXING_NDIM;
+extern int MIXING_RESTART;
 extern double MIXING_GG0;
 extern bool MIXING_TAU;
 extern double MIXING_BETA_MAG;

source/module_esolver/esolver_ks.cpp

@@ -410,7 +410,7 @@ namespace ModuleESolver
                         }
                     }
 
-                    this->conv_elec = (drho < this->scf_thr);
+                    this->conv_elec = (drho < this->scf_thr && iter!=GlobalV::MIXING_RESTART);
 
                     // If drho < hsolver_error in the first iter or drho < scf_thr, we do not change rho.
                     if (drho < hsolver_error || this->conv_elec)
@@ -436,8 +436,16 @@ namespace ModuleESolver
                         //     }
                         //     p_chgmix->auto_set(bandgap_for_autoset, GlobalC::ucell);
                         // }
-
-                        p_chgmix->mix_rho(pelec->charge);
+                        // mixing will restart after GlobalV::MIXING_RESTART steps
+                        // So, GlobalV::MIXING_RESTART=1 means mix from scratch
+                        if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1)
+                        {
+                            // do not mix charge density
+                        }
+                        else
+                        {
+                            p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing
+                        }
                         if (GlobalV::SCF_THR_TYPE == 2) pelec->charge->renormalize_rho(); // renormalize rho in R-space would induce a error in K-space
                         //----------charge mixing done-----------
                     }
@@ -467,6 +475,11 @@ namespace ModuleESolver
                     bool stop = this->do_after_converge(iter);
                     if(stop) break;
                 }
+                // notice for restart
+                if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1)
+                {
+                    std::cout<<"SCF restart after this step!"<<std::endl;
+                }
             }
             afterscf(istep);
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -491,8 +491,19 @@ namespace ModuleESolver
     template <typename TK, typename TR>
     void ESolver_KS_LCAO<TK, TR>::eachiterinit(const int istep, const int iter)
 {
-    if (iter == 1)
+    if (iter == 1 || iter == GlobalV::MIXING_RESTART)
+    {
+        if (iter == GlobalV::MIXING_RESTART) // delete mixing and re-construct it to restart 
+        {
+            this->p_chgmix->set_mixing(GlobalV::MIXING_MODE,
+                                GlobalV::MIXING_BETA,
+                                GlobalV::MIXING_NDIM,
+                                GlobalV::MIXING_GG0,
+                                GlobalV::MIXING_TAU,
+                                GlobalV::MIXING_BETA_MAG);
+        }
         this->p_chgmix->mix_reset();
+    }
 
     // mohan update 2012-06-05
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();

source/module_esolver/esolver_ks_pw.cpp

@@ -492,9 +492,19 @@ void ESolver_KS_PW<T, Device>::othercalculation(const int istep)
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::eachiterinit(const int istep, const int iter)
 {
-    if (iter == 1)
+    if (iter == 1 || iter == GlobalV::MIXING_RESTART)
+    {
+        if (iter == GlobalV::MIXING_RESTART) // delete mixing and re-construct it to restart 
+        {
+            this->p_chgmix->set_mixing(GlobalV::MIXING_MODE,
+                                GlobalV::MIXING_BETA,
+                                GlobalV::MIXING_NDIM,
+                                GlobalV::MIXING_GG0,
+                                GlobalV::MIXING_TAU,
+                                GlobalV::MIXING_BETA_MAG);
+        }
         this->p_chgmix->mix_reset();
-
+    }
     // mohan move harris functional to here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();

source/module_io/input.cpp

@@ -303,6 +303,7 @@ void Input::Default(void)
     mixing_mode = "broyden";
     mixing_beta = -10;
     mixing_ndim = 8;
+    mixing_restart = 0;
     mixing_gg0 = 1.00;       // use Kerker defaultly
     mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4
     mixing_gg0_mag = 0.0;    // defaultly exclude Kerker from mixing magnetic density
@@ -1258,6 +1259,10 @@ bool Input::Read(const std::string& fn)
         {
             read_value(ifs, mixing_ndim);
         }
+        else if (strcmp("mixing_restart", word) == 0)
+        {
+            read_value(ifs, mixing_restart);
+        }
         else if (strcmp("mixing_gg0", word) == 0) // mohan add 2014-09-27
         {
             read_value(ifs, mixing_gg0);
@@ -3325,6 +3330,7 @@ void Input::Bcast()
     Parallel_Common::bcast_string(mixing_mode);
     Parallel_Common::bcast_double(mixing_beta);
     Parallel_Common::bcast_int(mixing_ndim);
+    Parallel_Common::bcast_int(mixing_restart);
     Parallel_Common::bcast_double(mixing_gg0); // mohan add 2014-09-27
     Parallel_Common::bcast_double(mixing_beta_mag);
     Parallel_Common::bcast_double(mixing_gg0_mag);

source/module_io/input.h

@@ -232,6 +232,7 @@ class Input
     std::string mixing_mode; // "plain","broyden",...
     double mixing_beta; // 0 : no_mixing
     int mixing_ndim; // used in Broyden method
+    int mixing_restart;
     double mixing_gg0; // used in kerker method. mohan add 2014-09-27
     double mixing_beta_mag;
     double mixing_gg0_mag;

source/module_io/input_conv.cpp

@@ -750,6 +750,7 @@ void Input_Conv::Convert(void)
     GlobalV::MIXING_MODE = INPUT.mixing_mode;
     GlobalV::MIXING_BETA = INPUT.mixing_beta;
     GlobalV::MIXING_NDIM = INPUT.mixing_ndim;
+    GlobalV::MIXING_RESTART = INPUT.mixing_restart;
     GlobalV::MIXING_GG0 = INPUT.mixing_gg0;
     GlobalV::MIXING_BETA_MAG = INPUT.mixing_beta_mag;
     GlobalV::MIXING_GG0_MAG = INPUT.mixing_gg0_mag;

source/module_io/parameter_pool.cpp

@@ -818,6 +818,10 @@ void input_parameters_set(std::map<std::string, InputParameter> input_parameters
     {
         INPUT.mixing_ndim = *static_cast<int*>(input_parameters["mixing_ndim"].get());
     }
+    else if (input_parameters.count("mixing_restart") != 0)
+    {
+        INPUT.mixing_restart = *static_cast<int*>(input_parameters["mixing_restart"].get());
+    }
     else if (input_parameters.count("mixing_gg0") != 0)
     {
         INPUT.mixing_gg0 = *static_cast<double*>(input_parameters["mixing_gg0"].get());

source/module_io/test/input_conv_test.cpp

@@ -183,6 +183,7 @@ TEST_F(InputConvTest, Conv)
 	EXPECT_EQ(GlobalV::sc_mag_switch,0);
     EXPECT_TRUE(GlobalV::decay_grad_switch);
     EXPECT_EQ(GlobalV::sc_file, "sc.json");
+	EXPECT_EQ(GlobalV::MIXING_RESTART,0);
 }
 
 TEST_F(InputConvTest, ConvRelax)

source/module_io/test/input_test_para.cpp

@@ -381,6 +381,7 @@ TEST_F(InputParaTest, Bcast)
     EXPECT_TRUE(INPUT.mdp.dump_virial);
     EXPECT_FALSE(INPUT.mixing_tau);
     EXPECT_FALSE(INPUT.mixing_dftu);
+    EXPECT_EQ(INPUT.mixing_restart,0);
     EXPECT_EQ(INPUT.out_bandgap, 0);
     EXPECT_EQ(INPUT.out_mat_t, 0);
 

source/module_io/test/write_input_test.cpp

@@ -384,13 +384,16 @@ TEST_F(write_input, Mixing7)
     std::string output((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());
     EXPECT_THAT(output, testing::HasSubstr("#Parameters (7.Charge Mixing)"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_type                    broyden #plain; pulay; broyden"));
-    EXPECT_THAT(output,
-                testing::HasSubstr("mixing_beta                    0.7 #mixing parameter: 0 means no new charge"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_beta                    0.7 #mixing parameter: 0 means no new charge"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_ndim                    8 #mixing dimension in pulay or broyden"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_gg0                     0 #mixing parameter in kerker"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_beta_mag                -10 #mixing parameter for magnetic density"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_gg0_mag                 0 #mixing parameter in kerker"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_gg0_min                 0.1 #the minimum kerker coefficient"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_angle                   -10 #angle mixing parameter for non-colinear calculations"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_tau                     0 #whether to mix tau in mGGA calculation"));
-    EXPECT_THAT(output,
-                testing::HasSubstr("mixing_dftu                    0 #whether to mix locale in DFT+U calculation"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_dftu                    0 #whether to mix locale in DFT+U calculation"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_restart                 0 #which step to restart mixing during SCF"));
     EXPECT_THAT(output, testing::HasSubstr(""));
     ifs.close();
     remove("write_input_test.log");

source/module_io/write_input.cpp

@@ -248,6 +248,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_type", mixing_mode, "plain; pulay; broyden");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta", mixing_beta, "mixing parameter: 0 means no new charge");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_ndim", mixing_ndim, "mixing dimension in pulay or broyden");
+    ModuleBase::GlobalFunc::OUTP(ofs, "mixing_restart", mixing_restart, "which step to restart mixing during SCF");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0", mixing_gg0, "mixing parameter in kerker");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta_mag", mixing_beta_mag, "mixing parameter for magnetic density");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0_mag", mixing_gg0_mag, "mixing parameter in kerker");