diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index ffc39d5fd0..92f8e11f1f 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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 diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index a139d5dfd1..8f346d3b3d 100644 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -1701,4 +1701,60 @@ void Charge_Mixing::clean_data(std::complex*& data_s, std::complex= 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"); + } \ No newline at end of file diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 8317b2fb17..a863865987 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/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 _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 diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index da199d8d10..907596808c 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/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 drho(scf_nmax, 0.0); + std::vector 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]); + } +} diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 76e9d49ba3..bcc5361a7a 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -531,6 +531,11 @@ void ESolver_KS::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::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; } diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 26a2535885..ad01e0e47a 100644 --- a/source/module_esolver/esolver_ks.h +++ b/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 diff --git a/source/module_io/read_input_item_elec_stru.cpp b/source/module_io/read_input_item_elec_stru.cpp index 196ab540ee..e039b7f043 100644 --- a/source/module_io/read_input_item_elec_stru.cpp +++ b/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 " diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index a36cb5cd7d..69d6d13a1c 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/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); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 95a987e8e7..5fec3b43b6 100644 --- a/source/module_io/test/support/INPUT +++ b/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 diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index e637357c06..43b4de56da 100644 --- a/source/module_parameter/input_parameter.h +++ b/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