Fix: comments in PR

Author: dyzheng

docs/advanced/input_files/input-main.md

@@ -39,7 +39,6 @@
     - [pw\_diag\_thr](#pw_diag_thr)
     - [pw\_diag\_nmax](#pw_diag_nmax)
     - [pw\_diag\_ndim](#pw_diag_ndim)
-    - [diago\_full\_acc](#diago_full_acc)
     - [erf\_ecut](#erf_ecut)
     - [fft\_mode](#fft_mode)
     - [erf\_height](#erf_height)
@@ -779,12 +778,6 @@ These variables are used to control the plane wave related parameters.
 - **Description**: Only useful when you use `ks_solver = dav` or `ks_solver = dav_subspace`. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method. A larger value may yield a smaller number of iterations in the algorithm but uses more memory and more CPU time in subspace diagonalization.
 - **Default**: 4
 
-### diago_full_acc
-
-- **Type**: bool
-- **Description**: Only useful when you use `ks_solver = dav_subspace`. If `TRUE`, all the empty states are diagonalized at the same level of accuracy of the occupied ones. Otherwise the empty states are diagonalized using a larger threshold (10-5) (this should not affect total energy, forces, and other ground-state properties).
-- **Default**: false
-
 ### erf_ecut
 
 - **Type**: Real
@@ -925,7 +918,7 @@ calculations.
   - **cg**: cg method.
   - **bpcg**: bpcg method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment.
   - **dav**: the Davidson algorithm.
-  - **dav_subspace**: subspace Davidson algorithm
+  - **dav_subspace**: Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method.
 
   For atomic orbitals basis,
 

examples/lr-tddft/lcao_Si2/INPUT

@@ -37,4 +37,3 @@ out_alllog	1
 
 nvirt 19
 abs_wavelen_range  100 175
-#diago_full_acc 1

source/module_elecstate/potentials/pot_local.cpp

@@ -8,8 +8,6 @@
 namespace elecstate
 {
 
-double PotLocal::vl_of_0 = 0.0;
-
 //==========================================================
 // This routine computes the local potential in real space
 //==========================================================
@@ -31,7 +29,11 @@ void PotLocal::cal_fixed_v(double *vl_pseudo // store the local pseudopotential
         }
     }
 
-    PotLocal::vl_of_0 = vg[0].real();
+    /// save the V_local at G=0
+    if(this->rho_basis_->npw > 0)
+    {
+        *vl_of_0_ = vg[0].real();
+    }
 
     // recip2real should be a const function, but now it isn't
     // a dangerous usage appears here, which should be fix in the future.

source/module_elecstate/potentials/pot_local.h

@@ -12,8 +12,9 @@ class PotLocal : public PotBase
   public:
     PotLocal(const ModuleBase::matrix* vloc_in, // local pseduopotentials
              const ModuleBase::ComplexMatrix* sf_in,
-             const ModulePW::PW_Basis* rho_basis_in)
-        : vloc_(vloc_in), sf_(sf_in)
+             const ModulePW::PW_Basis* rho_basis_in,
+             double& vl_of_0)
+        : vloc_(vloc_in), sf_(sf_in), vl_of_0_(&vl_of_0)
     {
         assert(this->vloc_->nr == this->sf_->nr);
         this->rho_basis_ = rho_basis_in;
@@ -24,14 +25,10 @@ class PotLocal : public PotBase
 
     void cal_fixed_v(double* vl_pseudo) override;
 
-    /// @brief get the value of vloc at G=0;
-    /// @return vl(0)
-    static double get_vl_of_0() { return vl_of_0; }
-
     private:
 
     /// @brief save the value of vloc at G=0; this is a static member because there is only one vl(0) for all instances
-    static double vl_of_0;
+    double* vl_of_0_ = nullptr;
 
     // std::vector<double> vltot;
 

source/module_elecstate/potentials/potential_new.h

@@ -169,6 +169,14 @@ class Potential : public PotBase
         return this->v_effective_fixed.data();
     }
 
+
+    /// @brief get the value of vloc at G=0;
+    /// @return vl(0)
+    double get_vl_of_0() const
+    {
+        return this->vl_of_0;
+    }
+
   private:
     void cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff) override;
     void cal_fixed_v(double* vl_pseudo) override;
@@ -196,6 +204,8 @@ class Potential : public PotBase
     double* etxc_ = nullptr;
     double* vtxc_ = nullptr;
 
+    double vl_of_0 = 0.0;
+
     std::vector<PotBase*> components;
 
     const UnitCell* ucell_ = nullptr;

source/module_elecstate/potentials/potential_types.cpp

@@ -27,7 +27,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type)
     {
         if(!PARAM.inp.use_paw)
         {
-            return new PotLocal(this->vloc_, &(this->structure_factors_->strucFac), this->rho_basis_);
+            return new PotLocal(this->vloc_, &(this->structure_factors_->strucFac), this->rho_basis_, this->vl_of_0);
         }
         else
         {

source/module_hsolver/hsolver_pw.cpp

@@ -15,7 +15,6 @@
 #include "module_hsolver/diago_dav_subspace.h"
 #include "module_hsolver/diago_david.h"
 #include "module_hsolver/diago_iter_assist.h"
-#include "module_elecstate/potentials/pot_local.h"
 
 #include <algorithm>
 #include <vector>
@@ -212,19 +211,23 @@ void HSolverPW<T, Device>::paw_func_after_kloop(psi::Psi<T, Device>& psi, elecst
 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)
 {
+    // threshold for classifying occupied and unoccupied bands
+    const double occ_threshold = 1e-2;
+    // diagonalization threshold limitation for unoccupied bands
+    const double ethr_limit = 1e-5;
     if(wk > 0.0)
     {
         // Note: the idea of threshold for unoccupied bands (1e-5) comes from QE
         // In ABACUS, We applied a smoothing process to this truncation to avoid abrupt changes in energy errors between different bands.
-        const double ethr_unocc = std::max(1e-5, ethr);
+        const double ethr_unocc = std::max(ethr_limit, ethr);
         for (int i = 0; i < ethrs.size(); i++)
         {
             double band_weight = wg[i] / wk;
-            if (band_weight > 1e-2)
+            if (band_weight > occ_threshold)
             {
                 ethrs[i] = ethr; 
             }
-            else if(band_weight > 1e-5)
+            else if(band_weight > ethr_limit)
             {// similar energy difference for different bands when band_weight in range [1e-5, 1e-2]
                 ethrs[i] = std::min(ethr_unocc, ethr / band_weight);
             }
@@ -283,7 +286,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
         this->updatePsiK(pHamilt, psi, ik);
 
         // template add precondition calculating here
-        update_precondition(precondition, ik, this->wfc_basis->npwk[ik]);
+        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.
@@ -589,7 +592,7 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
 }
 
 template <typename T, typename Device>
-void HSolverPW<T, Device>::update_precondition(std::vector<Real>& h_diag, const int ik, const int npw)
+void HSolverPW<T, Device>::update_precondition(std::vector<Real>& h_diag, const int ik, const int npw, const Real vl_of_0)
 {
     h_diag.assign(h_diag.size(), 1.0);
     int precondition_type = 2;
@@ -616,7 +619,7 @@ void HSolverPW<T, Device>::update_precondition(std::vector<Real>& h_diag, const
 
             if (this->method == "dav_subspace")
             {
-                h_diag[ig] = g2kin + Real(elecstate::PotLocal::get_vl_of_0());
+                h_diag[ig] = g2kin + vl_of_0;
             }
             else
             {

source/module_hsolver/hsolver_pw.h

@@ -63,7 +63,7 @@ class HSolverPW
     void updatePsiK(hamilt::Hamilt<T, Device>* pHamilt, psi::Psi<T, Device>& psi, const int ik);
 
     // calculate the precondition array for diagonalization in PW base
-    void update_precondition(std::vector<Real>& h_diag, const int ik, const int npw);
+    void update_precondition(std::vector<Real>& h_diag, const int ik, const int npw, const Real vl_of_0);
 
     void output_iterInfo();
 

source/module_hsolver/hsolver_pw_sdft.cpp

@@ -44,7 +44,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt<std::complex<double>>* pHamilt,
         {
             this->updatePsiK(pHamilt, psi, ik);
             // template add precondition calculating here
-            update_precondition(precondition, ik, this->wfc_basis->npwk[ik]);
+            update_precondition(precondition, ik, this->wfc_basis->npwk[ik], pes->pot->get_vl_of_0());
             /// solve eigenvector and eigenvalue for H(k)
             double* p_eigenvalues = &(pes->ekb(ik, 0));
             this->hamiltSolvePsiK(pHamilt, psi, precondition, p_eigenvalues);

source/module_hsolver/test/hsolver_supplementary_mock.h

@@ -1,6 +1,5 @@
 #pragma once
 #include "module_elecstate/elecstate.h"
-#include "module_elecstate/potentials/pot_local.h"
 
 namespace elecstate
 {
@@ -58,8 +57,6 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge
     return;
 }
 
-double PotLocal::vl_of_0 = 0.0;
-
 } // namespace elecstate