remove PAW-related codes (#7277)

mohanchen · abacus_fixer · web-flow · commit 5271ce2af1f5 · 2026-04-24T12:33:10.000+08:00
* remove PAW-related codes

* update paw

* fix bug

* fix bugs

* fix bug

---------

Co-authored-by: abacus_fixer &lt;mohanchen@pku.eud.cn&gt;
diff --git a/cmake/FindMKL.cmake b/cmake/FindMKL.cmake
@@ -112,10 +112,7 @@ endif()
 
 endif() # MKL::MKL
 
-# For compatibility with legacy libpaw_interface CMakeLists.txt
-if(TARGET MKL::MKL AND NOT TARGET IntelMKL::MKL)
-  add_library(IntelMKL::MKL ALIAS MKL::MKL)
-endif()
+
 
 # In oneAPI 2022, MKL_SCALAPACK might not be linked properly
 if(NOT TARGET MKL::MKL_SCALAPACK)
diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md
@@ -53,7 +53,6 @@ For those who are interested in the source code, the following figure shows the
 |   `-- module_pw               Data structures and relevant methods for planewave involved calculations
 |-- source_cell                 The module for defining the unit cell and its operations, and reading pseudopotentials.
 |   |-- module_neighbor         The module for finding the neighbors of each atom in the unit cell.
-|   |-- module_paw              The module for performing PAW calculations.
 |   |-- module_symmetry         The module for finding the symmetry operations of the unit cell.
 |-- source_estate               The module for defining the electronic state and its operations.
 |   |-- module_charge           The module for calculating the charge density, charge mixing
diff --git a/docs/advanced/pp_orb.md b/docs/advanced/pp_orb.md
@@ -127,7 +127,7 @@ Users can find pseudopotentials in the following links:
 - [Quantum ESPRESSO](https://www.quantum-espresso.org/pseudopotentials): the official website of Quantum ESPRESSO, where you can find a large number of pseudopotential files.
 - [Stantard Solid State Pseudopotential library](https://www.materialscloud.org/sssp): a library of **high-quality** pseudopotentials for solid-state calculations, with **a large number of tests on efficiency and precison**.
 - [PWmat](http://www.pwmat.com/potential-download): a website that provides a large number of pseudopotential files, various kinds of semi-core constructed pseudopotentials are included. **Several sets (with or without f-electrons/noncolinear core correction) of Lanthanide pseudopotentials are also available**.
-- [THEOS](http://theossrv1.epfl.ch/Main/Pseudopotentials): PSlibrary 0.3.1, a library of pseudopotentials for DFT calculations, including ultrasoft, paw, norm-conserving both full-relativistic and scalar-relativistic pseudopotentials.
+- [THEOS](http://theossrv1.epfl.ch/Main/Pseudopotentials): PSlibrary 0.3.1, a library of pseudopotentials for DFT calculations, including ultrasoft, norm-conserving both full-relativistic and scalar-relativistic pseudopotentials.
 - [ABACUS@USTC](https://abacus.ustc.edu.cn/pseudo/list.htm): **ABACUS official website** where you can find a large number of pseudopotential files and numerical atomic orbital files.
 - [BLPS](https://github.com/PrincetonUniversity/BLPSLibrary): BLPS format pseudopotential library
 
diff --git a/python/pyabacus/src/ModuleESolver/components/diagonalizer_wrapper.hpp b/python/pyabacus/src/ModuleESolver/components/diagonalizer_wrapper.hpp
@@ -156,7 +156,6 @@ class DiagonalizerWrapper : public IDiagonalizer<TK>
             config_.tolerance,
             diag_ethr,
             config_.max_iterations,
-            config_.use_paw,
             comm_info
         );
 
diff --git a/python/pyabacus/src/ModuleESolver/interfaces/i_diagonalizer.hpp b/python/pyabacus/src/ModuleESolver/interfaces/i_diagonalizer.hpp
@@ -58,7 +58,6 @@ struct DiagConfig
     double tolerance = 1e-6;     ///< Convergence tolerance
     int max_iterations = 100;    ///< Maximum iterations
     int dav_ndim = 4;            ///< Davidson subspace dimension multiplier
-    bool use_paw = false;        ///< Use PAW method
     int nproc_in_pool = 1;       ///< Number of processes in pool
 };
 
diff --git a/python/pyabacus/src/hsolver/diago_adapter.hpp b/python/pyabacus/src/hsolver/diago_adapter.hpp
@@ -84,7 +84,6 @@ class PyDiagoDavidAdapter
         double tol,
         std::vector<double>& diag_ethr,
         int max_iter,
-        bool use_paw,
         ::hsolver::diag_comm_info comm_info)
     {
         auto hpsi_func = make_hpsi_func_fstyle<T>(mm_op);
@@ -95,7 +94,6 @@ class PyDiagoDavidAdapter
             nband_,
             nbasis_,
             dav_ndim,
-            use_paw,
             comm_info
         );
 
diff --git a/python/pyabacus/src/hsolver/py_diago_david.hpp b/python/pyabacus/src/hsolver/py_diago_david.hpp
@@ -106,7 +106,6 @@ class PyDiagoDavid
         double tol,
         std::vector<double>& diag_ethr,
         int max_iter,
-        bool use_paw,
         hsolver::diag_comm_info comm_info
     ) {
         auto hpsi_func = [mm_op] (
@@ -143,7 +142,6 @@ class PyDiagoDavid
             nband, 
             nbasis, 
             dav_ndim, 
-            use_paw,
             comm_info
         );
 
diff --git a/python/pyabacus/src/hsolver/py_hsolver.cpp b/python/pyabacus/src/hsolver/py_hsolver.cpp
@@ -145,16 +145,13 @@ void bind_hsolver(py::module& m)
                 The tolerance vector.
             max_iter : int
                 The maximum number of iterations.
-            use_paw : bool
-                Whether to use the projector augmented wave method.
         )pbdoc",
         "mm_op"_a,
         "precond_vec"_a,
         "dav_ndim"_a,
         "tol"_a,
         "diag_ethr"_a,
         "max_iter"_a,
-        "use_paw"_a,
         "comm_info"_a)
         .def("set_psi", &PyDiagoDavidAdapter::set_psi, R"pbdoc(
             Set the initial guess of the eigenvectors, i.e. the wave functions.
diff --git a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py
@@ -136,7 +136,6 @@ def davidson(
     tol: float = 1e-2,
     max_iter: int = 1000,
     diag_ethr: Union[List[float], None] = None,
-    use_paw: bool = False,
     # scf_type: bool = False
 ) -> Tuple[NDArray[np.float64], NDArray[np.complex128]]:
     """ A function to diagonalize a matrix using the Davidson-Subspace method.
@@ -163,8 +162,6 @@ def davidson(
         The maximum number of iterations, by default 1000.
     diag_ethr : List[float] | None, optional
         The list of thresholds of bands, by default None.    
-    use_paw : bool, optional
-        Whether to use projector augmented wave (PAW) method, by default False.
     
     Returns
     -------
@@ -195,7 +192,6 @@ def davidson(
         tol,
         diag_ethr,
         max_iter,
-        use_paw,
         comm_info
     )
     
diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp
@@ -216,7 +216,6 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
                                                      PARAM.inp.calculation,
                                                      PARAM.inp.basis_type,
                                                      PARAM.inp.ks_solver,
-                                                     false,
                                                      PARAM.globalv.use_uspp,
                                                      PARAM.inp.nspin,
                                                      hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp
@@ -158,7 +158,6 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
                                                            PARAM.inp.calculation,
                                                            PARAM.inp.basis_type,
                                                            PARAM.inp.ks_solver,
-                                                           false,
                                                            PARAM.globalv.use_uspp,
                                                            PARAM.inp.nspin,
                                                            hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
diff --git a/source/source_hsolver/diago_david.cpp b/source/source_hsolver/diago_david.cpp
@@ -16,9 +16,8 @@ DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in,
                                   const int nband_in,
                                   const int dim_in,
                                   const int david_ndim_in,
-                                  const bool use_paw_in,
                                   const diag_comm_info& diag_comm_in)
-    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in)
+    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), diag_comm(diag_comm_in)
 {
     this->device = base_device::get_device_type(this->ctx);
     this->precondition = precondition_in;
diff --git a/source/source_hsolver/diago_david.h b/source/source_hsolver/diago_david.h
@@ -45,7 +45,6 @@ class DiagoDavid
      * @param[in] david_ndim_in Dimension of the reduced basis set of Davidson.
      *                      `david_ndim_in` * `nband_in` is the maximum allowed size of
      *                      the reduced basis set before \b restart of Davidson.
-     * @param[in] use_paw_in Flag indicating whether to use PAW.
      * @param[in] diag_comm_in Communication information for diagonalization.
      *
      * @tparam T The data type of the matrices and arrays.
@@ -57,7 +56,6 @@ class DiagoDavid
                const int nband_in,
                const int dim_in,
                const int david_ndim_in,
-               const bool use_paw_in,
                const diag_comm_info& diag_comm_in);
 
     /**
@@ -142,7 +140,6 @@ class DiagoDavid
       const int notconv_max = 0); // Maximum number of allowed non-converged eigenvectors
 
   private:
-    bool use_paw = false;
     int test_david = 0;
 
     diag_comm_info diag_comm;
diff --git a/source/source_hsolver/hsolver_pw.cpp b/source/source_hsolver/hsolver_pw.cpp
@@ -109,9 +109,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
             // update H(k) for each k point
             pHamilt->updateHk(ik);
 
-#ifdef USE_PAW
-            this->paw_func_in_kloop(ik, tpiba);
-#endif
+
 
             // update psi pointer for each k point
             psi.fix_k(ik);
@@ -133,9 +131,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
                                     ethr_band);
             }
 
-#ifdef USE_PAW
-            this->call_paw_cell_set_currentk(ik);
-#endif
+
 
             // solve eigenvector and eigenvalue for H(k)
             this->hamiltSolvePsiK(pHamilt, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks);
@@ -156,9 +152,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
             // update H(k) for each k point
             pHamilt->updateHk(ik);
 
-#ifdef USE_PAW
-            this->paw_func_in_kloop(ik, tpiba);
-#endif
+
 
             // update psi pointer for each k point
             psi.fix_k(ik);
@@ -175,9 +169,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
                                     ethr_band);
             }
 
-#ifdef USE_PAW
-            this->call_paw_cell_set_currentk(ik);
-#endif
+
 
             // solve eigenvector and eigenvalue for H(k)
             this->hamiltSolvePsiK(pHamilt, psi, precondition, eigenvalues.data() + ik * psi.get_nbands(), this->wfc_basis->nks);
@@ -350,7 +342,7 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
         const int nband = psi.get_nbands();            /// number of eigenpairs sought
         const int ld_psi = psi.get_nbasis();           /// leading dimension of psi
 
-        DiagoDavid<T, Device> david(pre_condition.data(), nband, dim, PARAM.inp.pw_diag_ndim, this->use_paw, comm_info);
+        DiagoDavid<T, Device> david(pre_condition.data(), nband, dim, PARAM.inp.pw_diag_ndim, comm_info);
         // do diag and add davidson iteration counts up to avg_iter
         DiagoIterAssist<T, Device>::avg_iter += static_cast<double>(
              david.diag(hpsi_func,
diff --git a/source/source_hsolver/hsolver_pw.h b/source/source_hsolver/hsolver_pw.h
@@ -28,7 +28,6 @@ class HSolverPW
               const std::string calculation_type_in,
               const std::string basis_type_in,
               const std::string method_in,
-              const bool use_paw_in,
               const bool use_uspp_in,
               const int nspin_in,
               const int scf_iter_in,
@@ -37,7 +36,7 @@ class HSolverPW
               const bool need_subspace_in,
               const bool use_k_continuity_in = false)
         : wfc_basis(wfc_basis_in), calculation_type(calculation_type_in), basis_type(basis_type_in), method(method_in),
-          use_paw(use_paw_in), use_uspp(use_uspp_in), nspin(nspin_in), scf_iter(scf_iter_in),
+          use_uspp(use_uspp_in), nspin(nspin_in), scf_iter(scf_iter_in),
           diag_iter_max(diag_iter_max_in), diag_thr(diag_thr_in), need_subspace(need_subspace_in),
           use_k_continuity(use_k_continuity_in) {};
 
@@ -76,7 +75,6 @@ class HSolverPW
     const std::string calculation_type;
     const std::string basis_type;
     const std::string method;
-    const bool use_paw;
     const bool use_uspp;
     const int nspin;
 
@@ -100,14 +98,7 @@ class HSolverPW
     /// @brief calculate the threshold for iterative-diagonalization for each band
     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,
-                           const double tpiba);
 
-    void call_paw_cell_set_currentk(const int ik);
-
-    void paw_func_after_kloop(psi::Psi<T, Device>& psi, elecstate::ElecState* pes,const double tpiba,const int nat);
-#endif
 
     // K-point continuity related members
     std::vector<int> k_order;
diff --git a/source/source_hsolver/hsolver_pw_sdft.h b/source/source_hsolver/hsolver_pw_sdft.h
@@ -20,7 +20,6 @@ class HSolverPW_SDFT : public HSolverPW<T, Device>
                    const std::string calculation_type_in,
                    const std::string basis_type_in,
                    const std::string method_in,
-                   const bool use_paw_in,
                    const bool use_uspp_in,
                    const int nspin_in,
                    const int scf_iter_in,
@@ -31,7 +30,6 @@ class HSolverPW_SDFT : public HSolverPW<T, Device>
                                calculation_type_in,
                                basis_type_in,
                                method_in,
-                               use_paw_in,
                                use_uspp_in,
                                nspin_in,
                                scf_iter_in,
diff --git a/source/source_hsolver/test/diago_david_float_test.cpp b/source/source_hsolver/test/diago_david_float_test.cpp
@@ -93,7 +93,7 @@ class DiagoDavPrepare
 		const int dim = phi.get_current_ngk() ;
 		const int nband = phi.get_nbands();
 		const int ld_psi =phi.get_nbasis();
-		hsolver::DiagoDavid<std::complex<float>> dav(precondition, nband, dim, order, false, comm_info);
+		hsolver::DiagoDavid<std::complex<float>> dav(precondition, nband, dim, order, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<float>>::PW_DIAG_NMAX = maxiter;
 		hsolver::DiagoIterAssist<std::complex<float>>::PW_DIAG_THR = eps;
diff --git a/source/source_hsolver/test/diago_david_real_test.cpp b/source/source_hsolver/test/diago_david_real_test.cpp
@@ -92,7 +92,7 @@ class DiagoDavPrepare
 		const int dim = phi.get_current_ngk();
         const int nband = phi.get_nbands();
         const int ld_psi = phi.get_nbasis();
-        hsolver::DiagoDavid<double> dav(precondition, nband, dim, order, false, comm_info);
+        hsolver::DiagoDavid<double> dav(precondition, nband, dim, order, comm_info);
 
         hsolver::DiagoIterAssist<double>::PW_DIAG_NMAX = maxiter;
         hsolver::DiagoIterAssist<double>::PW_DIAG_THR = eps;
diff --git a/source/source_hsolver/test/diago_david_test.cpp b/source/source_hsolver/test/diago_david_test.cpp
@@ -97,7 +97,7 @@ class DiagoDavPrepare
 		const int dim = phi.get_current_ngk();
 		const int nband = phi.get_nbands();
 		const int ld_psi = phi.get_nbasis();
-		hsolver::DiagoDavid<std::complex<double>> dav(precondition, nband, dim, order, false, comm_info);
+		hsolver::DiagoDavid<std::complex<double>> dav(precondition, nband, dim, order, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<double>>::PW_DIAG_NMAX = maxiter;
 		hsolver::DiagoIterAssist<std::complex<double>>::PW_DIAG_THR = eps;
diff --git a/source/source_hsolver/test/hsolver_pw_sup.h b/source/source_hsolver/test/hsolver_pw_sup.h
@@ -123,9 +123,8 @@ DiagoDavid<T, Device>::DiagoDavid(const Real* precondition_in,
                                   const int nband_in,
                                   const int dim_in,
                                   const int david_ndim_in,
-                                  const bool use_paw_in,
                                   const diag_comm_info& diag_comm_in)
-    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), use_paw(use_paw_in), diag_comm(diag_comm_in) {
+    : nband(nband_in), dim(dim_in), nbase_x(david_ndim_in * nband_in), david_ndim(david_ndim_in), diag_comm(diag_comm_in) {
     this->device = base_device::get_device_type(this->ctx);
     this->precondition = precondition_in;
 
diff --git a/source/source_hsolver/test/test_hsolver_sdft.cpp b/source/source_hsolver/test/test_hsolver_sdft.cpp
@@ -280,7 +280,6 @@ class TestHSolverPW_SDFT : public ::testing::Test
             "scf",
             "pw",
             "cg",
-            false,
             PARAM.sys.use_uspp,
             PARAM.input.nspin,
             hsolver::DiagoIterAssist<std::complex<double>>::SCF_ITER,
diff --git a/source/source_lcao/module_lr/hsolver_lrtd.hpp b/source/source_lcao/module_lr/hsolver_lrtd.hpp
@@ -86,7 +86,6 @@ namespace LR
                                                  nband,
                                                  dim,
                                                  PARAM.inp.pw_diag_ndim,
-                                                 false,
                                                  comm_info);
                     std::vector<double> ethr_band(nband, diag_ethr);
                     hsolver::DiagoIterAssist<T>::avg_iter += static_cast<double>(david.diag(hpsi_func, spsi_func,
diff --git a/source/source_pw/module_pwdft/forces.cpp b/source/source_pw/module_pwdft/forces.cpp
@@ -49,7 +49,6 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
     ModuleBase::matrix forcecc(nat, 3);
     ModuleBase::matrix forcenl(nat, 3);
     ModuleBase::matrix forcescc(nat, 3);
-    ModuleBase::matrix forcepaw(nat, 3);
     ModuleBase::matrix forceonsite(nat, 3);
 
     // Force due to local ionic potential
