Feature: add smooth ethr for all iter methods (#5732)

haozhihan · web-flow · commit 72bd505824c0 · 2024-12-16T20:05:48.000+08:00
* update input parameter

* add smooth ethr for all methods

* fix integrated_test

* Update input-main.md
diff --git 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)
+    - [diago\_smooth\_ethr](#diago_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 `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
 
 - **Type**: Integer
diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
@@ -28,8 +28,7 @@ namespace hsolver
 
 #ifdef USE_PAW
 template <typename T, typename Device>
-void HSolverPW<T, Device>::paw_func_in_kloop(const int ik,
-                                             const double tpiba)
+void HSolverPW<T, Device>::paw_func_in_kloop(const int ik, const double tpiba)
 {
     if (this->use_paw)
     {
@@ -97,7 +96,7 @@ void HSolverPW<T, Device>::call_paw_cell_set_currentk(const int ik)
 }
 
 template <typename T, typename Device>
-void HSolverPW<T, Device>::paw_func_after_kloop(psi::Psi<T, Device>& psi, 
+void HSolverPW<T, Device>::paw_func_after_kloop(psi::Psi<T, Device>& psi,
                                                 elecstate::ElecState* pes,
                                                 const double tpiba,
                                                 const int nat)
@@ -211,10 +210,10 @@ void HSolverPW<T, Device>::paw_func_after_kloop(psi::Psi<T, Device>& psi,
 #endif
 
 template <typename T, typename Device>
-void HSolverPW<T, Device>::cal_ethr_band(const double& wk,
-                                         const double* wg,
-                                         const double& ethr,
-                                         std::vector<double>& ethrs)
+void HSolverPW<T, Device>::cal_smooth_ethr(const double& wk,
+                                           const double* wg,
+                                           const double& ethr,
+                                           std::vector<double>& ethrs)
 {
     // threshold for classifying occupied and unoccupied bands
     const double occ_threshold = 1e-2;
@@ -288,7 +287,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* 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<T, Device>::solve(hamilt::Hamilt<T, Device>* 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<T, Device>::PW_DIAG_THR,
-                                ethr_band);
+            this->cal_smooth_ethr(pes->klist->wk[ik],
+                                  &pes->wg(ik, 0),
+                                  DiagoIterAssist<T, Device>::PW_DIAG_THR,
+                                  ethr_band);
         }
 
 #ifdef USE_PAW
@@ -347,7 +345,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
         reinterpret_cast<elecstate::ElecStatePW<T, Device>*>(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<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* 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
@@ -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<double>& ethrs);
+    void cal_smooth_ethr(const double& wk, const double* wg, const double& ethr, std::vector<double>& ethrs);
 
 #ifdef USE_PAW
     void paw_func_in_kloop(const int ik,
diff --git 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
@@ -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
@@ -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
@@ -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
 
diff --git 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
@@ -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
@@ -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

Original file line number	Diff line number	Diff line change
`@@ -314,6 +314,12 @@ void ReadInput::item_elec_stru()`
`314`	`314`	`};`
`315`	`315`	`this->add_item(item);`
`316`	`316`	`}`
	`317`	`+ {`
	`318`	`+ Input_Item item("diago_smooth_ethr");`
	`319`	`+ item.annotation = "smooth ethr for iter methods";`
	`320`	`+ read_sync_bool(input.diago_smooth_ethr);`
	`321`	`+ this->add_item(item);`
	`322`	`+ }`
`317`	`323`	`{`
`318`	`324`	`Input_Item item("pw_diag_ndim");`
`319`	`325`	`item.annotation = "dimension of workspace for Davidson diagonalization";`