deepmodeling · cheerly-pku · Mar 28, 2026 · Mar 28, 2026 · Apr 7, 2026 · Apr 8, 2026
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -66,6 +66,7 @@
     - [use\_k\_continuity](#use_k_continuity)
     - [pw\_diag\_nmax](#pw_diag_nmax)
     - [pw\_diag\_ndim](#pw_diag_ndim)
+    - [pw\_diag\_rr\_step](#pw_diag_rr_step)
     - [diago\_cg\_prec](#diago_cg_prec)
   - [Numerical atomic orbitals related variables](#numerical-atomic-orbitals-related-variables)
     - [lmaxmax](#lmaxmax)
@@ -1058,7 +1059,7 @@
 ### pw_diag_thr
 
 - **Type**: Real
-- **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 &lt;= 1e-3.
+- **Description**: Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. 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 &lt;= 1e-3.
 - **Default**: 0.01
 
 ### diago_smooth_ethr
@@ -1077,16 +1078,24 @@
 ### pw_diag_nmax
 
 - **Type**: Integer
-- **Availability**: *[`basis_type`](#basis_type)==pw and [`ks_solver`](#ks_solver) in [cg, dav, dav_subspace, bpcg]*
-- **Description**: Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method.
+- **Availability**: *[`basis_type`](#basis_type)==pw and [`ks_solver`](#ks_solver) in [cg, dav, dav_subspace, bpcg, ppcg]*
+- **Description**: Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg/ppcg method.
 - **Default**: 50
 
 ### pw_diag_ndim
 
 - **Type**: Integer
-- **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.
+- **Availability**: *[`basis_type`](#basis_type)==pw and [`ks_solver`](#ks_solver) in [dav, dav_subspace, ppcg]*
+- **Description**: Only useful when you use ks_solver = dav, dav_subspace, or ppcg. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method, and the block size for the PPCG 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
 
+### pw_diag_rr_step
+
+- **Type**: Integer
+- **Availability**: *[`basis_type`](#basis_type)==pw and [`ks_solver`](#ks_solver)==ppcg*
+- **Description**: Only useful when you use ks_solver = ppcg. It controls how often (in subspace iterations) H and S are re-applied to reset the accumulated rounding drift after the Rayleigh-Ritz rotation. A larger value reduces the number of H/S applications and thus the wall time without changing the iteration count in well-conditioned cases; a smaller value is more robust against rounding drift in ill-conditioned problems.
+- **Default**: 16
+
 ### diago_cg_prec
 
 - **Type**: Integer
@@ -1192,6 +1201,7 @@
   - cg: The conjugate-gradient (CG) method.
   - dav: The Davidson algorithm.
   - dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method.
+  - ppcg: The projection preconditioned conjugate-gradient method. It is optimized and validated for CPU plane-wave calculations; non-CPU devices use a transitional host/device bridge.
   - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use.
 
   For numerical atomic orbitals basis,

diff --git a/docs/advanced/scf/hsolver.md b/docs/advanced/scf/hsolver.md
@@ -4,7 +4,7 @@
 
 Method of explicit solving KS-equation can be chosen by variable "ks_solver" in INPUT file.
 
-When "basis_type = pw", `ks_solver` can be `cg`, `bpcg` or `dav`. The default setting `cg` is recommended, which is band-by-band conjugate gradient diagonalization method. There is a large probability that the use of setting of `dav` , which is block Davidson diagonalization method, can be tried to improve performance.  
+When "basis_type = pw", `ks_solver` can be `cg`, `bpcg`, `dav`, `dav_subspace`, or `ppcg`. The default setting `cg` is recommended, which is a band-by-band conjugate-gradient diagonalization method. The `dav` and `dav_subspace` settings use Davidson-style subspace diagonalization and can be tried to improve performance. The `ppcg` setting uses the projection preconditioned conjugate-gradient method (a restarted block method that keeps a bounded `2*nband` subspace). It targets the many-eigenpair regime where the bounded memory and block operations pay off; it is optimized and validated for CPU plane-wave calculations, and non-CPU devices use a transitional host/device bridge. The PPCG block size / Rayleigh-Ritz interval is controlled by `pw_diag_ndim`.
 
 When "basis_type = lcao", `ks_solver` can be `genelpa` or `scalapack_gvx`. The default setting `genelpa` is recommended, which is based on ELPA (EIGENVALUE SOLVERS FOR PETAFLOP APPLICATIONS) (https://elpa.mpcdf.mpg.de/) and the kernel is auto choosed by GENELPA(https://github.com/pplab/GenELPA), usually faster than the setting of "scalapack_gvx", which is based on ScaLAPACK(Scalable Linear Algebra PACKage)  
 

diff --git a/docs/parameters.yaml b/docs/parameters.yaml
@@ -543,6 +543,7 @@ parameters:
       * cg: The conjugate-gradient (CG) method.
       * dav: The Davidson algorithm.
       * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method.
+      * ppcg: The projection preconditioned conjugate-gradient method. It is optimized and validated for CPU plane-wave calculations; non-CPU devices use a transitional host/device bridge.
       * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use.
 
       For numerical atomic orbitals basis,
@@ -976,7 +977,7 @@ parameters:
     category: Plane wave related variables
     type: Real
     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.
+      Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. 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_value: "0.01"
     unit: ""
     availability: ""
@@ -1000,18 +1001,26 @@ parameters:
     category: Plane wave related variables
     type: Integer
     description: |
-      Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method.
+      Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg/ppcg method.
     default_value: "50"
     unit: ""
-    availability: "basis_type==pw and ks_solver in [cg, dav, dav_subspace, bpcg]"
+    availability: "basis_type==pw and ks_solver in [cg, dav, dav_subspace, bpcg, ppcg]"
   - name: pw_diag_ndim
     category: Plane wave related variables
     type: Integer
     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.
+      Only useful when you use ks_solver = dav, dav_subspace, or ppcg. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method, and the block size for the PPCG 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_value: "4"
     unit: ""
-    availability: ""
+    availability: "basis_type==pw and ks_solver in [dav, dav_subspace, ppcg]"
+  - name: pw_diag_rr_step
+    category: Plane wave related variables
+    type: Integer
+    description: |
+      Only useful when you use ks_solver = ppcg. It controls how often (in subspace iterations) H and S are re-applied to reset the accumulated rounding drift after the Rayleigh-Ritz rotation. A larger value reduces the number of H/S applications and thus the wall time without changing the iteration count in well-conditioned cases; a smaller value is more robust against rounding drift in ill-conditioned problems.
+    default_value: "16"
+    unit: ""
+    availability: basis_type==pw and ks_solver==ppcg
   - name: diago_cg_prec
     category: Plane wave related variables
     type: Integer

diff --git a/source/source_base/parallel_reduce.cpp b/source/source_base/parallel_reduce.cpp
@@ -104,6 +104,7 @@ template void Parallel_Reduce::reduce_pool<double>(double&);
 template void Parallel_Reduce::reduce_pool<std::complex<double>>(std::complex<double>&);
 
 template void Parallel_Reduce::reduce_pool<int>(int*, const int);
+template void Parallel_Reduce::reduce_pool<float>(float*, const int);
 template void Parallel_Reduce::reduce_pool<double>(double*, const int);
 template void Parallel_Reduce::reduce_pool<std::complex<float>>(std::complex<float>*, const int);
 template void Parallel_Reduce::reduce_pool<std::complex<double>>(std::complex<double>*, const int);

diff --git a/source/source_hsolver/CMakeLists.txt b/source/source_hsolver/CMakeLists.txt
@@ -4,6 +4,7 @@ list(APPEND objects
     diago_david.cpp
     diago_dav_subspace.cpp
     diago_bpcg.cpp
+    diago_ppcg.cpp
     para_lin_tf.cpp
     hsolver_pw.cpp
     hsolver_lcaopw.cpp

diff --git a/source/source_hsolver/diago_iter_assist.h b/source/source_hsolver/diago_iter_assist.h
@@ -21,6 +21,8 @@ class DiagoIterAssist
   public:
     static Real PW_DIAG_THR;
     static int PW_DIAG_NMAX;
+    static int PW_DIAG_NDIM;
+    static int PW_DIAG_RR_STEP;
 
     static Real LCAO_DIAG_THR;
     static int LCAO_DIAG_NMAX;
@@ -157,6 +159,12 @@ typename DiagoIterAssist<T, Device>::Real DiagoIterAssist<T, Device>::avg_iter =
 template <typename T, typename Device>
 int DiagoIterAssist<T, Device>::PW_DIAG_NMAX = 30;
 
+template <typename T, typename Device>
+int DiagoIterAssist<T, Device>::PW_DIAG_NDIM = 4;
+
+template <typename T, typename Device>
+int DiagoIterAssist<T, Device>::PW_DIAG_RR_STEP = 16;
+
 template <typename T, typename Device>
 typename DiagoIterAssist<T, Device>::Real DiagoIterAssist<T, Device>::PW_DIAG_THR = 1.0e-2;
 
@@ -179,4 +187,4 @@ template <typename T, typename Device>
 T DiagoIterAssist<T, Device>::zero = static_cast<T>(0.0);
 } // namespace hsolver
 
-#endif
+#endif
diff --git a/source/source_hsolver/diago_params.cpp b/source/source_hsolver/diago_params.cpp
@@ -15,6 +15,8 @@ void setup_diago_params_pw(const int istep,
     DiagoIterAssist<T, Device>::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true;
     DiagoIterAssist<T, Device>::SCF_ITER = iter;
     DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
+    DiagoIterAssist<T, Device>::PW_DIAG_NDIM = inp.pw_diag_ndim;
+    DiagoIterAssist<T, Device>::PW_DIAG_RR_STEP = inp.pw_diag_rr_step;
 
     if (inp.calculation != "nscf")
     {
@@ -41,6 +43,8 @@ void setup_diago_params_sdft(const int istep,
 
     DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
     DiagoIterAssist<T, Device>::PW_DIAG_NMAX = inp.pw_diag_nmax;
+    DiagoIterAssist<T, Device>::PW_DIAG_NDIM = inp.pw_diag_ndim;
+    DiagoIterAssist<T, Device>::PW_DIAG_RR_STEP = inp.pw_diag_rr_step;
 }
 
 /// Template instantiation for CPU