From c7e7c02b134d7a7467d92a7c14d8d92253413110 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Sun, 15 Dec 2024 06:52:31 +0000 Subject: [PATCH 1/4] update input parameter --- docs/advanced/input_files/input-main.md | 7 +++++++ source/module_io/read_input_item_elec_stru.cpp | 6 ++++++ source/module_io/test/read_input_ptest.cpp | 1 + source/module_io/test/support/INPUT | 1 + source/module_parameter/input_parameter.h | 1 + 5 files changed, 16 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 2b5d1009af..dc13a45c80 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -37,6 +37,7 @@ - [ndx, ndy, ndz](#ndx-ndy-ndz) - [pw\_seed](#pw_seed) - [pw\_diag\_thr](#pw_diag_thr) + - [diag\_smooth\_ethr](#diag_smooth_ethr) - [pw\_diag\_nmax](#pw_diag_nmax) - [pw\_diag\_ndim](#pw_diag_ndim) - [erf\_ecut](#erf_ecut) @@ -777,6 +778,12 @@ These variables are used to control the plane wave related parameters. - **Description**: Only used when you use `ks_solver = cg/dav/dav_subspace/bpcg`. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. **For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.** - **Default**: 0.01 +### diago_smooth_ethr + +- **Type**: bool +- **Description**: If `FALSE`, all the empty states are diagonalized at the same level of accuracy of the occupied ones. If `TRUE` the empty states are diagonalized using a larger threshold (10-5) (this should not affect total energy, forces, and other ground-state properties, but computational efficiency will be improved.). +- **Default**: false + ### pw_diag_nmax - **Type**: Integer diff --git a/source/module_io/read_input_item_elec_stru.cpp b/source/module_io/read_input_item_elec_stru.cpp index 75026089c8..089f842de5 100644 --- a/source/module_io/read_input_item_elec_stru.cpp +++ b/source/module_io/read_input_item_elec_stru.cpp @@ -314,6 +314,12 @@ void ReadInput::item_elec_stru() }; this->add_item(item); } + { + Input_Item item("diago_smooth_ethr"); + item.annotation = "smooth ethr for iter methods"; + read_sync_bool(input.diago_smooth_ethr); + this->add_item(item); + } { Input_Item item("pw_diag_ndim"); item.annotation = "dimension of workspace for Davidson diagonalization"; diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 33608f6569..11acb41f43 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -154,6 +154,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.diago_cg_prec, 1); EXPECT_EQ(param.inp.pw_diag_ndim, 4); EXPECT_DOUBLE_EQ(param.inp.pw_diag_thr, 1.0e-2); + EXPECT_FALSE(param.inp.diago_smooth_ethr); EXPECT_EQ(param.inp.nb2d, 0); EXPECT_EQ(param.inp.nurse, 0); EXPECT_EQ(param.inp.t_in_h, 1); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 580a6ae803..944e46ca2f 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -50,6 +50,7 @@ erf_height 20 #the height of the energy step for reciprocal erf_sigma 4 #the width of the energy step for reciprocal vectors fft_mode 0 #mode of FFTW pw_diag_thr 0.01 #threshold for eigenvalues is cg electron iterations +diago_smooth_ethr false #smooth ethr for iter methods 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 diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index fe86fbfefc..9b2151c945 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -86,6 +86,7 @@ struct Input_para int nspin = 1; ///< LDA ; LSDA ; non-linear spin int pw_diag_nmax = 50; double pw_diag_thr = 0.01; ///< used in cg method + bool diago_smooth_ethr = false; ///< smooth ethr for iter methods int pw_diag_ndim = 4; ///< dimension of workspace for Davidson diagonalization int diago_cg_prec = 1; ///< mohan add 2012-03-31 From 0dc51f7b1eaa462619bdf4ccc29ef125c118a35e Mon Sep 17 00:00:00 2001 From: haozhihan Date: Sun, 15 Dec 2024 07:12:08 +0000 Subject: [PATCH 2/4] add smooth ethr for all methods --- docs/advanced/input_files/input-main.md | 2 +- source/module_hsolver/hsolver_pw.cpp | 32 ++++++++++++------------- source/module_hsolver/hsolver_pw.h | 2 +- 3 files changed, 17 insertions(+), 19 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index dc13a45c80..17319e4d07 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -37,7 +37,7 @@ - [ndx, ndy, ndz](#ndx-ndy-ndz) - [pw\_seed](#pw_seed) - [pw\_diag\_thr](#pw_diag_thr) - - [diag\_smooth\_ethr](#diag_smooth_ethr) + - [diago\_smooth\_ethr](#diago_smooth_ethr) - [pw\_diag\_nmax](#pw_diag_nmax) - [pw\_diag\_ndim](#pw_diag_ndim) - [erf\_ecut](#erf_ecut) diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 83ef0ed833..1a8e7292f6 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -28,8 +28,7 @@ namespace hsolver #ifdef USE_PAW template -void HSolverPW::paw_func_in_kloop(const int ik, - const double tpiba) +void HSolverPW::paw_func_in_kloop(const int ik, const double tpiba) { if (this->use_paw) { @@ -97,7 +96,7 @@ void HSolverPW::call_paw_cell_set_currentk(const int ik) } template -void HSolverPW::paw_func_after_kloop(psi::Psi& psi, +void HSolverPW::paw_func_after_kloop(psi::Psi& psi, elecstate::ElecState* pes, const double tpiba, const int nat) @@ -211,10 +210,10 @@ void HSolverPW::paw_func_after_kloop(psi::Psi& psi, #endif template -void HSolverPW::cal_ethr_band(const double& wk, - const double* wg, - const double& ethr, - std::vector& ethrs) +void HSolverPW::cal_smooth_ethr(const double& wk, + const double* wg, + const double& ethr, + std::vector& ethrs) { // threshold for classifying occupied and unoccupied bands const double occ_threshold = 1e-2; @@ -288,7 +287,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, pHamilt->updateHk(ik); #ifdef USE_PAW - this->paw_func_in_kloop(ik,tpiba); + this->paw_func_in_kloop(ik, tpiba); #endif /// update psi pointer for each k point @@ -297,14 +296,13 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // template add precondition calculating here update_precondition(precondition, ik, this->wfc_basis->npwk[ik], Real(pes->pot->get_vl_of_0())); - // only dav_subspace method used smooth threshold for all bands now, - // for other methods, this trick can be added in the future to accelerate calculation without accuracy loss. - if (this->method == "dav_subspace") + // use smooth threshold for all iter methods + if (PARAM.inp.diago_smooth_ethr == true) { - this->cal_ethr_band(pes->klist->wk[ik], - &pes->wg(ik, 0), - DiagoIterAssist::PW_DIAG_THR, - ethr_band); + this->cal_smooth_ethr(pes->klist->wk[ik], + &pes->wg(ik, 0), + DiagoIterAssist::PW_DIAG_THR, + ethr_band); } #ifdef USE_PAW @@ -347,7 +345,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, reinterpret_cast*>(pes)->psiToRho(psi); #ifdef USE_PAW - this->paw_func_after_kloop(psi, pes,tpiba,nat); + this->paw_func_after_kloop(psi, pes, tpiba, nat); #endif ModuleBase::timer::tick("HSolverPW", "solve"); @@ -473,7 +471,7 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, } else if (this->method == "bpcg") { - const int nband = psi.get_nbands(); + const int nband = psi.get_nbands(); const int nbasis = psi.get_nbasis(); auto ngk_pointer = psi.get_ngk_pointer(); // hpsi_func (X, HX, ld, nvec) -> HX = H(X), X and HX blockvectors of size ld x nvec diff --git a/source/module_hsolver/hsolver_pw.h b/source/module_hsolver/hsolver_pw.h index 49d8b2806e..97b8143eba 100644 --- a/source/module_hsolver/hsolver_pw.h +++ b/source/module_hsolver/hsolver_pw.h @@ -88,7 +88,7 @@ class HSolverPW private: /// @brief calculate the threshold for iterative-diagonalization for each band - void cal_ethr_band(const double& wk, const double* wg, const double& ethr, std::vector& ethrs); + void cal_smooth_ethr(const double& wk, const double* wg, const double& ethr, std::vector& ethrs); #ifdef USE_PAW void paw_func_in_kloop(const int ik, From 5fd1398539184058c30e22370371ddfedb9b4e7a Mon Sep 17 00:00:00 2001 From: haozhihan Date: Sun, 15 Dec 2024 08:27:31 +0000 Subject: [PATCH 3/4] fix integrated_test --- tests/integrate/104_PW_NC_magnetic/INPUT | 1 + tests/integrate/140_PW_15_SO/INPUT | 1 + tests/integrate/160_PW_DJ_PK_PU_SO/INPUT | 1 + 3 files changed, 3 insertions(+) diff --git a/tests/integrate/104_PW_NC_magnetic/INPUT b/tests/integrate/104_PW_NC_magnetic/INPUT index d761feb2c1..c01847313f 100644 --- a/tests/integrate/104_PW_NC_magnetic/INPUT +++ b/tests/integrate/104_PW_NC_magnetic/INPUT @@ -22,6 +22,7 @@ mixing_beta 0.2 mixing_ndim 10 ks_solver dav_subspace +diago_smooth_ethr true pw_diag_ndim 2 basis_type pw symmetry 0 diff --git a/tests/integrate/140_PW_15_SO/INPUT b/tests/integrate/140_PW_15_SO/INPUT index 1c58840435..5df7c5728d 100644 --- a/tests/integrate/140_PW_15_SO/INPUT +++ b/tests/integrate/140_PW_15_SO/INPUT @@ -34,6 +34,7 @@ lspinorb 1 basis_type pw ks_solver dav_subspace +diago_smooth_ethr true pw_diag_ndim 2 chg_extrap second-order out_dm 0 diff --git a/tests/integrate/160_PW_DJ_PK_PU_SO/INPUT b/tests/integrate/160_PW_DJ_PK_PU_SO/INPUT index 8e4f45ab0e..3b28a05ff4 100644 --- a/tests/integrate/160_PW_DJ_PK_PU_SO/INPUT +++ b/tests/integrate/160_PW_DJ_PK_PU_SO/INPUT @@ -27,6 +27,7 @@ mixing_dmr 1 mixing_gg0 1.1 ks_solver dav_subspace +diago_smooth_ethr true pw_diag_ndim 2 basis_type pw gamma_only 0 From e53f43fec0dd1548f21ec7408f89f893b90c169d Mon Sep 17 00:00:00 2001 From: Haozhi Han Date: Mon, 16 Dec 2024 16:59:48 +0800 Subject: [PATCH 4/4] Update input-main.md --- docs/advanced/input_files/input-main.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 17319e4d07..a2b137e8a4 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -781,7 +781,7 @@ These variables are used to control the plane wave related parameters. ### diago_smooth_ethr - **Type**: bool -- **Description**: If `FALSE`, all the empty states are diagonalized at the same level of accuracy of the occupied ones. If `TRUE` the empty states are diagonalized using a larger threshold (10-5) (this should not affect total energy, forces, and other ground-state properties, but computational efficiency will be improved.). +- **Description**: If `TRUE`, the smooth threshold strategy, which applies a larger threshold (10e-5) for the empty states, will be implemented in the diagonalization methods. (This strategy should not affect total energy, forces, and other ground-state properties, but computational efficiency will be improved.) If `FALSE`, the smooth threshold strategy will not be applied. - **Default**: false ### pw_diag_nmax