Merge branch 'develop' into cusolver-fix

dzzz2001 · Oct 16, 2024 · bffe5eb · bffe5eb
2 parents 6a378dd + a5c35d9
commit bffe5eb
Show file tree

Hide file tree

Showing 46 changed files with 736 additions and 1,275 deletions.
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -39,7 +39,7 @@ jobs:
             --from-ref ${{ github.event.pull_request.base.sha }}
               --to-ref ${{ github.event.pull_request.head.sha }}
         continue-on-error: true
-      - uses: pre-commit-ci/lite-action@v1.1.0
+      - uses: pre-commit-ci/lite-action@v1.0.3
 
       - name: Build
         run: |

diff --git a/docs/advanced/input_files/input-main.md b/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,
 

diff --git a/examples/lr-tddft/lcao_Si2/INPUT b/examples/lr-tddft/lcao_Si2/INPUT
@@ -37,4 +37,3 @@ out_alllog	1
 
 nvirt 19
 abs_wavelen_range  100 175
-#diago_full_acc 1
diff --git a/python/pyabacus/src/py_diago_dav_subspace.hpp b/python/pyabacus/src/py_diago_dav_subspace.hpp
@@ -106,7 +106,7 @@ class PyDiagoDavSubspace
         double tol,
         int max_iter,
         bool need_subspace,
-        std::vector<bool> is_occupied,
+        std::vector<double> diag_ethr,
         bool scf_type,
         hsolver::diag_comm_info comm_info
     ) {
@@ -141,7 +141,7 @@ class PyDiagoDavSubspace
             comm_info
         );
 
-        return obj->diag(hpsi_func, psi, nbasis, eigenvalue, is_occupied, scf_type);
+        return obj->diag(hpsi_func, psi, nbasis, eigenvalue, diag_ethr.data(), scf_type);
     }
 
 private:

diff --git a/python/pyabacus/src/py_hsolver.cpp b/python/pyabacus/src/py_hsolver.cpp
@@ -59,8 +59,8 @@ void bind_hsolver(py::module& m)
                 The maximum number of iterations.
             need_subspace : bool
                 Whether to use the subspace function.
-            is_occupied : list[bool]
-                A list of boolean values indicating whether the band is occupied,
+            diag_ethr : list[float]
+                A list of float values indicating the thresholds of each band for the diagonalization,
                 meaning that the corresponding eigenvalue is to be calculated.
             scf_type : bool
                 Whether to use the SCF type, which is used to determine the
@@ -76,7 +76,7 @@ void bind_hsolver(py::module& m)
         "tol"_a, 
         "max_iter"_a, 
         "need_subspace"_a, 
-        "is_occupied"_a, 
+        "diag_ethr"_a, 
         "scf_type"_a, 
         "comm_info"_a)
         .def("set_psi", &py_hsolver::PyDiagoDavSubspace::set_psi, R"pbdoc(

diff --git a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py
@@ -25,7 +25,7 @@ def dav_subspace(
     tol: float = 1e-2,
     max_iter: int = 1000,
     need_subspace: bool = False,
-    is_occupied: Union[List[bool], None] = None,
+    diag_ethr: Union[List[float], None] = None,
     scf_type: bool = False
 ) -> Tuple[NDArray[np.float64], NDArray[np.complex128]]:
     """ A function to diagonalize a matrix using the Davidson-Subspace method.
@@ -52,10 +52,8 @@ def dav_subspace(
         The maximum number of iterations, by default 1000.
     need_subspace : bool, optional
         Whether to use subspace function, by default False.
-    is_occupied : List[bool] | None, optional
-        The list of occupied bands, by default None. This indicates how many eigenvalues 
-        need to be calculated, starting from the smallest eigenvalue. Only the energy levels 
-        occupied by electrons (occupied) need to be calculated.
+    diag_ethr : List[float] | None, optional
+        The list of thresholds of bands, by default None.
     scf_type : bool, optional
         Indicates whether the calculation is a self-consistent field (SCF) calculation. 
         If True, the initial precision of eigenvalue calculation can be coarse. 
@@ -72,8 +70,8 @@ def dav_subspace(
     if not callable(mvv_op):
         raise TypeError("mvv_op must be a callable object.")
 
-    if is_occupied is None:
-        is_occupied = [True] * num_eigs
+    if diag_ethr is None:
+        diag_ethr = [tol] * num_eigs
 
     if init_v.ndim != 1 or init_v.dtype != np.complex128:
         init_v = init_v.flatten().astype(np.complex128, order='C')
@@ -93,7 +91,7 @@ def dav_subspace(
         tol,
         max_iter,
         need_subspace,
-        is_occupied,
+        diag_ethr,
         scf_type,
         comm_info
     )
@@ -113,7 +111,6 @@ def davidson(
     tol: float = 1e-2,
     max_iter: int = 1000,
     use_paw: bool = False,
-    # is_occupied: Union[List[bool], None] = None,
     # scf_type: bool = False
 ) -> Tuple[NDArray[np.float64], NDArray[np.complex128]]:
     """ A function to diagonalize a matrix using the Davidson-Subspace method.

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -568,12 +568,9 @@ OBJS_LCAO=evolve_elec.o\
       FORCE_STRESS.o\
 	  FORCE_gamma.o\
 	  FORCE_k.o\
-	  fvl_dphi_gamma.o\
-	  fvl_dphi_k.o\
-	  fedm_gamma.o\
-	  fedm_k.o\
-	  ftvnl_dphi_gamma.o\
-	  ftvnl_dphi_k.o\
+      stress_tools.o\
+      edm.o\
+      pulay_force_stress_center2.o\
 	  fvnl_dbeta_gamma.o\
 	  fvnl_dbeta_k.o\
 	  grid_init.o\

diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp
@@ -252,31 +252,6 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge
     this->wg.create(nk_in, PARAM.inp.nbands);
 }
 
-void set_is_occupied(std::vector<bool>& is_occupied,
-                     elecstate::ElecState* pes,
-                     const int i_scf,
-                     const int nk,
-                     const int nband,
-                     const bool diago_full_acc)
-{
-    if (i_scf != 0 && diago_full_acc == false)
-    {
-        for (int i = 0; i < nk; i++)
-        {
-            if (pes->klist->wk[i] > 0.0)
-            {
-                for (int j = 0; j < nband; j++)
-                {
-                    if (pes->wg(i, j) / pes->klist->wk[i] < 0.01)
-                    {
-                        is_occupied[i * nband + j] = false;
-                    }
-                }
-            }
-        }
-    }
-};
-
 
 
 } // namespace elecstate
diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h
@@ -194,13 +194,5 @@ class ElecState
     bool skip_weights = false;
 };
 
-// This is an independent function under the elecstate namespace and does not depend on any class.
-void set_is_occupied(std::vector<bool>& is_occupied,
-                     elecstate::ElecState* pes,
-                     const int i_scf,
-                     const int nk,
-                     const int nband,
-                     const bool diago_full_acc);
-
 } // namespace elecstate
 #endif
diff --git a/source/module_elecstate/potentials/pot_local.cpp b/source/module_elecstate/potentials/pot_local.cpp
@@ -29,6 +29,12 @@ void PotLocal::cal_fixed_v(double *vl_pseudo // store the local pseudopotential
         }
     }
 
+    /// 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.
     const_cast<ModulePW::PW_Basis *>(this->rho_basis_)->recip2real(vg, vl_pseudo);

diff --git a/source/module_elecstate/potentials/pot_local.h b/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,6 +25,11 @@ class PotLocal : public PotBase
 
     void cal_fixed_v(double* vl_pseudo) override;
 
+    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
+    double* vl_of_0_ = nullptr;
+
     // std::vector<double> vltot;
 
     const ModuleBase::matrix* vloc_ = nullptr; // local pseduopotentials

diff --git a/source/module_elecstate/potentials/potential_new.h b/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;

diff --git a/source/module_elecstate/potentials/potential_types.cpp b/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
         {

diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -352,15 +352,6 @@ void ESolver_KS_PW<T, Device>::hamilt2density(const int istep, const int iter, c
         hsolver::DiagoIterAssist<T, Device>::SCF_ITER = iter;
         hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
         hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax;
-
-        std::vector<bool> is_occupied(this->kspw_psi->get_nk() * this->kspw_psi->get_nbands(), true);
-
-        elecstate::set_is_occupied(is_occupied,
-                                   this->pelec,
-                                   hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
-                                   this->kspw_psi->get_nk(),
-                                   this->kspw_psi->get_nbands(),
-                                   PARAM.inp.diago_full_acc);
 
         hsolver::HSolverPW<T, Device> hsolver_pw_obj(this->pw_wfc, 
                                                      &this->wf, 
@@ -383,7 +374,6 @@ void ESolver_KS_PW<T, Device>::hamilt2density(const int istep, const int iter, c
                            this->kspw_psi[0],
                            this->pelec,
                            this->pelec->ekb.c,
-                           is_occupied,
                            GlobalV::RANK_IN_POOL,
                            GlobalV::NPROC_IN_POOL,
                            false);

diff --git a/source/module_esolver/pw_fun.cpp b/source/module_esolver/pw_fun.cpp
@@ -71,15 +71,6 @@ void ESolver_KS_PW<T, Device>::hamilt2estates(const double ethr)
     hsolver::DiagoIterAssist<T, Device>::need_subspace = false;
     hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
 
-    std::vector<bool> is_occupied(this->kspw_psi->get_nk() * this->kspw_psi->get_nbands(), true);
-
-    elecstate::set_is_occupied(is_occupied,
-                                this->pelec,
-                                hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
-                                this->kspw_psi->get_nk(),
-                                this->kspw_psi->get_nbands(),
-                                PARAM.inp.diago_full_acc);
-
     hsolver::HSolverPW<T, Device> hsolver_pw_obj(this->pw_wfc, 
                                                  &this->wf, 
 
@@ -101,7 +92,6 @@ void ESolver_KS_PW<T, Device>::hamilt2estates(const double ethr)
                          this->kspw_psi[0],
                          this->pelec,
                          this->pelec->ekb.c,
-                         is_occupied,
                          GlobalV::RANK_IN_POOL,
                          GlobalV::NPROC_IN_POOL,
                          true);

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
@@ -15,15 +15,12 @@ if(ENABLE_LCAO)
         operator_lcao/td_nonlocal_lcao.cpp
         operator_lcao/sc_lambda_lcao.cpp
         operator_lcao/dftu_lcao.cpp
+        pulay_force_stress_center2.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp
-        fvl_dphi_gamma.cpp
-        fvl_dphi_k.cpp
-        fedm_gamma.cpp
-        fedm_k.cpp
-        ftvnl_dphi_gamma.cpp
-        ftvnl_dphi_k.cpp
+        stress_tools.cpp
+        edm.cpp
         fvnl_dbeta_gamma.cpp
         fvnl_dbeta_k.cpp
         grid_init.cpp

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h
@@ -94,7 +94,7 @@ class Force_LCAO
                   const bool isstress,
                   ForceStressArrays& fsr,
                   const UnitCell& ucell,
-                  const elecstate::DensityMatrix<T, double>* dm,
+                  const elecstate::DensityMatrix<T, double>& dm,
                   const psi::Psi<T>* psi,
                   const Parallel_Orbitals& pv,
                   const elecstate::ElecState* pelec,
@@ -136,12 +136,16 @@ class Force_LCAO
                       typename TGint<T>::type& gint,
                       ModuleBase::matrix& fvl_dphi,
                       ModuleBase::matrix& svl_dphi);
+
+    elecstate::DensityMatrix<T, double> cal_edm(const elecstate::ElecState* pelec,
+        const psi::Psi<T>& psi,
+        const elecstate::DensityMatrix<T, double>& dm,
+        const K_Vectors& kv,
+        const Parallel_Orbitals& pv,
+        const int& nspin, 
+        const int& nbands,
+        const UnitCell& ucell,
+        Record_adj& ra) const;
 };
 
-// this namespace used to store global function for some stress operation
-namespace StressTools
-{
-// set upper matrix to whole matrix
-void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix);
-} // namespace StressTools
 #endif