Feature: enable ABACUS can finish SCF if charge density oscillation i…

…s found (#5421)

* add an input parameter to control the threshold for SCF oscillation

* add a function to tell whether SCF oscillates

* add a member oscillate_esolver in esolver_ks.h, similar to conv_esolver

* add a read-in test for scf_thr_os

* add a read-in test for scf_ene_thr

* add a new parameter scf_os_ndim and corresponding test

* set scf_os_ndim as mixing_ndim if default

* add a UnitTest for new function

* add some docs and comments

* add a new parameter scf_os_stop

* use scf_os_stop to control esovler

* rename scf_thr_os as scf_os_thr

* update the docs of scf_os_thr

* add a value-check for scf_os_thr

docs/advanced/input_files/input-main.md

@@ -88,6 +88,9 @@
     - [scf\_thr](#scf_thr)
     - [scf\_ene\_thr](#scf_ene_thr)
     - [scf\_thr\_type](#scf_thr_type)
+    - [scf\_os\_stop](#scf_os_stop)
+    - [scf\_os\_thr](#scf_os_thr)
+    - [scf\_os\_ndim](#scf_os_ndim)
     - [chg\_extrap](#chg_extrap)
     - [lspinorb](#lspinorb)
     - [noncolin](#noncolin)
@@ -1205,6 +1208,29 @@ Note: In new angle mixing, you should set `mixing_beta_mag >> mixing_beta`. The
 
 - **Default**: 1 (plane-wave basis), or 2 (localized atomic orbital basis).
 
+### scf_os_stop
+
+- **Type**: bool
+- **Description**: For systems that are difficult to converge, the SCF process may exhibit oscillations in charge density, preventing further progress toward the specified convergence criteria and resulting in continuous oscillation until the maximum number of steps is reached; this greatly wastes computational resources. To address this issue, this function allows ABACUS to terminate the SCF process early upon detecting oscillations, thus reducing subsequent meaningless calculations. The detection of oscillations is based on the slope of the logarithm of historical drho values.. To this end, Least Squares Method is used to calculate the slope of the logarithmically taken drho for the previous `scf_os_ndim` iterations. If the calculated slope is larger than `scf_os_thr`, stop the SCF.
+
+  - **0**: The SCF will continue to run regardless of whether there is oscillation or not. 
+  - **1**: If the calculated slope is larger than `scf_os_thr`, stop the SCF.
+
+- **Default**: false
+
+### scf_os_thr
+
+- **Type**: double
+- **Description**: The slope threshold to determine if the SCF is stuck in a charge density oscillation. If the calculated slope is larger than `scf_os_thr`, stop the SCF.
+
+- **Default**: -0.01
+
+### scf_os_ndim
+
+- **Type**: int
+- **Description**: To determine the number of old iterations' `drho` used in slope calculations.
+- **Default**: `mixing_ndim`
+
 ### chg_extrap
 
 - **Type**: String

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -1701,4 +1701,60 @@ void Charge_Mixing::clean_data(std::complex<double>*& data_s, std::complex<doubl
         data_s = nullptr;
         data_hf = nullptr;
     }
+}
+
+bool Charge_Mixing::if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold)
+{
+    ModuleBase::TITLE("Charge_Mixing", "if_scf_oscillate");
+    ModuleBase::timer::tick("Charge_Mixing", "if_scf_oscillate");
+
+    if(threshold >= 0) // close the function
+    {
+        return false;
+    }
+
+    if(this->_drho_history.size() == 0)
+    {
+        this->_drho_history.resize(PARAM.inp.scf_nmax);
+    }
+
+    // add drho into history
+    this->_drho_history[iteration - 1] = drho;
+
+    // check if the history is long enough
+    if(iteration < iternum_used)
+    {
+        return false;
+    }
+
+    // calculate the slope of the last iternum_used iterations' drho
+    double slope = 0.0;
+
+    // Least Squares Method
+    // this part is too short, so I do not design it as a free function in principle
+    double sumX = 0, sumY = 0, sumXY = 0, sumXX = 0;
+    for (int i = iteration - iternum_used; i < iteration; i++)
+    {
+        sumX += i;
+        sumY += std::log10(this->_drho_history[i]);
+        sumXY += i * std::log10(this->_drho_history[i]);
+        sumXX += i * i;
+    }
+    double numerator = iternum_used * sumXY - sumX * sumY;
+    double denominator = iternum_used * sumXX - sumX * sumX;
+    if (denominator == 0) {
+        return false;
+    }
+    slope =  numerator / denominator;
+
+    // if the slope is less than the threshold, return true
+    if(slope > threshold)
+    {
+        return true;
+    }
+
+    return false;
+
+    ModuleBase::timer::tick("Charge_Mixing", "if_scf_oscillate");
+
 }

source/module_elecstate/module_charge/charge_mixing.h

@@ -105,6 +105,9 @@ class Charge_Mixing
     int mixing_restart_step = 0; //which step to restart mixing during SCF
     int mixing_restart_count = 0; // the number of restart mixing during SCF. Do not set mixing_restart_count as bool since I want to keep some flexibility in the future
 
+    // to calculate the slope of drho curve during SCF, which is used to determine if SCF oscillate
+    bool if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold);
+
   private:
 
     // mixing_data
@@ -129,6 +132,8 @@ class Charge_Mixing
     double mixing_angle = 0.0;           ///< mixing angle for nspin=4
     bool mixing_dmr = false;             ///< whether to mixing real space density matrix
 
+    std::vector<double> _drho_history; ///< history of drho used to determine the oscillation, size is scf_nmax
+
     bool new_e_iteration = true;
 
     ModulePW::PW_Basis* rhopw = nullptr;  ///< smooth grid

source/module_elecstate/test/charge_mixing_test.cpp

@@ -1024,3 +1024,41 @@ TEST_F(ChargeMixingTest, MixDivCombTest)
     EXPECT_EQ(datas2, nullptr);
     EXPECT_EQ(datahf2, nullptr);
 }
+
+TEST_F(ChargeMixingTest, SCFOscillationTest)
+{
+    Charge_Mixing CMtest;
+    int scf_nmax = 20;
+    int scf_os_ndim = 3;
+    double scf_os_thr = -0.05;
+    bool scf_oscillate = false;
+    std::vector<double> drho(scf_nmax, 0.0);
+    std::vector<bool> scf_oscillate_ref(scf_nmax, false);
+    drho = {6.83639633652e-05,
+            4.93523029235e-05,
+            3.59230097735e-05,
+            2.68356403913e-05,
+            2.17490806464e-05,
+            2.14231642508e-05,
+            1.67507494811e-05,
+            1.53575889539e-05,
+            1.26504511554e-05,
+            1.04762016224e-05,
+            8.10000162918e-06,
+            7.66427917682e-06,
+            6.70112820094e-06,
+            5.68594436664e-06,
+            4.80120233733e-06,
+            4.86519757184e-06,
+            4.37855804356e-06,
+            4.29922703412e-06,
+            4.36398486331e-06,
+            4.94224615955e-06};
+    scf_oscillate_ref = {false,false,false,false,false,true,false,false,false,false,
+                        false,false,true,false,false,true,true,true,true,true};
+    for (int i = 1; i <= scf_nmax; ++i)
+    {
+        scf_oscillate = CMtest.if_scf_oscillate(i,drho[i-1],scf_os_ndim,scf_os_thr);
+        EXPECT_EQ(scf_oscillate, scf_oscillate_ref[i-1]);
+    } 
+}

source/module_esolver/esolver_ks.cpp

@@ -531,6 +531,11 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
                 this->p_chgmix->mixing_restart_step = iter + 1;
             }
 
+            if (PARAM.inp.scf_os_stop) // if oscillation is detected, SCF will stop
+            {
+                this->oscillate_esolver = this->p_chgmix->if_scf_oscillate(iter, drho, PARAM.inp.scf_os_ndim, PARAM.inp.scf_os_thr);
+            }
+
             // drho will be 0 at this->p_chgmix->mixing_restart step, which is
             // not ground state
             bool not_restart_step = !(iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0);
@@ -639,9 +644,13 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
 #endif //__RAPIDJSON
 
         // 13) check convergence
-        if (this->conv_esolver)
+        if (this->conv_esolver || this->oscillate_esolver)
         {
             this->niter = iter;
+            if (this->oscillate_esolver)
+            {
+                std::cout << " !! Density oscillation is found, STOP HERE !!" << std::endl;
+            }
             break;
         }
 

source/module_esolver/esolver_ks.h

@@ -82,6 +82,7 @@ class ESolver_KS : public ESolver_FP
       protected:
         std::string basisname; // PW or LCAO
         double esolver_KS_ne = 0.0;
+		bool oscillate_esolver = false; // whether esolver is oscillated
 };	
 } // end of namespace
 #endif

source/module_io/read_input_item_elec_stru.cpp

@@ -528,6 +528,36 @@ void ReadInput::item_elec_stru()
         read_sync_double(input.scf_ene_thr);
         this->add_item(item);
     }
+    {
+        Input_Item item("scf_os_stop");
+        item.annotation = "whether to stop scf when oscillation is detected";
+        read_sync_bool(input.scf_os_stop);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("scf_os_thr");
+        item.annotation = "charge density threshold for oscillation";
+        read_sync_double(input.scf_os_thr);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.scf_os_thr >= 0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "scf_os_thr should be negative");
+            }
+        };
+        this->add_item(item);
+    }
+    {
+        Input_Item item("scf_os_ndim");
+        item.annotation = "number of old iterations used for oscillation detection";
+        read_sync_int(input.scf_os_ndim);
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            if (para.input.scf_os_ndim <= 0) // default value
+            {
+                para.input.scf_os_ndim = para.input.mixing_ndim;
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("scf_thr_type");
         item.annotation = "type of the criterion of scf_thr, 1: reci drho for "

source/module_io/test/read_input_ptest.cpp

@@ -164,7 +164,10 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(PARAM.inp.test_force, 0);
     EXPECT_EQ(param.inp.test_stress, 0);
     EXPECT_NEAR(param.inp.scf_thr, 1.0e-8, 1.0e-15);
-    EXPECT_NEAR(param.inp.scf_ene_thr, -1.0, 1.0e-15);
+    EXPECT_EQ(param.inp.scf_os_stop, 1);
+    EXPECT_NEAR(param.inp.scf_os_thr, -0.02, 1.0e-15);
+    EXPECT_EQ(param.inp.scf_os_ndim, 10);
+    EXPECT_NEAR(param.inp.scf_ene_thr, 1.0e-6, 1.0e-15);
     EXPECT_EQ(param.inp.scf_nmax, 50);
     EXPECT_EQ(param.inp.relax_nmax, 1);
     EXPECT_EQ(param.inp.out_stru, 0);

source/module_io/test/support/INPUT

@@ -51,6 +51,10 @@ erf_sigma                      4 #the width of the energy step for reciprocal ve
 fft_mode                       0 #mode of FFTW
 pw_diag_thr                    0.01 #threshold for eigenvalues is cg electron iterations
 scf_thr                        1e-08 #charge density error
+scf_ene_thr                    1e-06 #total energy error threshold
+scf_os_stop                    1 #whether to stop scf when oscillation is detected
+scf_os_thr                     -0.02 #charge density threshold for oscillation
+scf_os_ndim                    10 #number of old iterations used for oscillation detection
 scf_thr_type                   2 #type of the criterion of scf_thr, 1: reci drho for pw, 2: real drho for lcao
 init_wfc                       atomic #start wave functions are from 'atomic', 'atomic+random', 'random' or 'file'
 init_chg                       atomic #start charge is from 'atomic' or file

source/module_parameter/input_parameter.h

@@ -115,6 +115,9 @@ struct Input_para
     double scf_ene_thr = -1.0; ///< energy threshold for scf convergence, in eV
     int scf_thr_type = -1;     ///< type of the criterion of scf_thr, 1: reci drho, 2: real drho
     bool final_scf = false;    ///< whether to do final scf
+    bool scf_os_stop = false;  ///< whether to stop scf when oscillation is detected
+    double scf_os_thr = -0.01;  ///< drho threshold for oscillation
+    int scf_os_ndim = 0;       ///< number of old iterations used for oscillation detection
 
     bool lspinorb = false;   ///< consider the spin-orbit interaction
     bool noncolin = false;   ///< using non-collinear-spin