Feature: Support ultrasoft pseudopotentials (USPP) in PW BPCG (#7904)

* Feature: Enable USPP generalized eigenproblems in PW BPCG

* Reject unsupported USPP SDFT and nspin=4 calculations

* Fix: Synchronize PW projector table caches after cell rescaling

* Fix: Correct float input handling in radial_fft_q

* Fix: Use double becsum for single-precision USPP force and stress

* Reject unsupported USPP calculations before they start

* Fix: Use a dense FFT grid when only the z dimension is increased

* Test: Add USPP BPCG regression coverage with BNDPAR

* Resolve CI/CD failure regarding Makefile

* Update 005_PW_SDFT_MALL_BPCG_GPU reference

* Update BPCG output references after solver correction

Author: AsTonyshment

source/Makefile.Objects

@@ -819,6 +819,7 @@ OBJS_SRCPW=h_ewald_pw.o\
     setup_pot.o\
     setup_pwrho.o\
     setup_pwwfc.o\
+    uspp_support.o\
     update_cell_pw.o\
     dftu_base.o\
     dftu_base_io.o\

source/source_esolver/esolver_fp.cpp

@@ -1,5 +1,6 @@
 #include "esolver_fp.h"
 
+#include "source_base/tool_quit.h"
 #include "source_cell/cal_ux.h"
 #include "source_estate/module_charge/symm_rho.h"
 #include "source_cell/read_pp_ucell.h"
@@ -12,6 +13,7 @@
 #include "source_io/module_parameter/parameter.h"
 
 #include "source_pw/module_pwdft/setup_pwrho.h" // mohan 20251005
+#include "source_pw/module_pwdft/uspp_support.h"
 #include "source_hamilt/module_xc/xc_functional.h" // mohan 20251005
 #include "source_io/module_ctrl/ctrl_output_fp.h"
 #include "source_io/module_chgpot/write_init.h" // write_chg_init, write_pot_init
@@ -69,8 +71,20 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp)
                                             this->inp_->bndpar,
                                             this->inp_->nelec,
                                             this->inp_->nupdown);
+
     elecstate::ParamUpdater::update_from_atoms_info(atoms_info);
 
+    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+    pw::validate_uspp_support(atoms_info.use_uspp,
+                              inp.basis_type,
+                              inp.esolver_type,
+                              inp.nspin,
+                              XC_Functional::get_func_type(),
+                              inp.berry_phase,
+                              inp.towannier90,
+                              inp.cal_cond);
+    GlobalV::ofs_running << XC_Functional::output_info() << std::endl;
+
     //! 2) setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials
     pw::setup_pwrho(ucell, PARAM.globalv.double_grid, this->pw_rho_flag, 
       this->pw_rho, this->pw_rhod, this->pw_big, this->classname, inp);
@@ -112,10 +126,6 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp)
     //! 10) calculate the structure factor
     this->sf.setup(&ucell, Pgrid, this->pw_rhod);
 
-    //! 11) setup the xc functional
-    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
-    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
-
     //! 11) initialize the charge density, we need to first set xc_type,
     // then we can call chr.allocate()
 	this->chr.set_rhopw(this->pw_rhod); // mohan add 20251130

source/source_estate/elecstate_pw.cpp

@@ -49,10 +49,6 @@ ElecStatePW<T, Device>::~ElecStatePW()
             delete[] this->kin_r;
         }
     }
-    if (PARAM.globalv.use_uspp)
-    {
-        delmem_var_h_op()(this->becsum);
-    }
     delmem_complex_op()(this->wfcr);
     delmem_complex_op()(this->wfcr_another_spin);
 }
@@ -291,9 +287,8 @@ void ElecStatePW<T, Device>::cal_becsum(const psi::Psi<T, Device>& psi)
     const int nkb = this->ppcell->nkb;
     this->vkb = this->ppcell->template get_vkb_data<Real>();
     const int nh_tot = this->ppcell->nhm * (this->ppcell->nhm + 1) / 2;
-    // becsum on CPU (forces_us / stress_us use CPU dgemm)
-    resmem_var_h_op()(becsum, nh_tot * ucell->nat * PARAM.inp.nspin, "ElecState<PW>::becsum");
-    setmem_var_h_op()(becsum, 0, nh_tot * ucell->nat * PARAM.inp.nspin);
+    const int becsum_size = nh_tot * ucell->nat * PARAM.inp.nspin;
+    this->becsum_.assign(becsum_size, 0.0);
 
     // becp: device buffer for gemm, then D2H for host loops
     T* becp = nullptr;
@@ -429,11 +424,11 @@ void ElecStatePW<T, Device>::cal_becsum(const psi::Psi<T, Device>& psi)
                             {
                                 if (ih == jh)
                                 {
-                                    becsum[index + ijh] += std::real(aux_gk_host[ih * nh_atom + jh]);
+                                    this->becsum_[index + ijh] += static_cast<double>(std::real(aux_gk_host[ih * nh_atom + jh]));
                                 }
                                 else
                                 {
-                                    becsum[index + ijh] += 2.0 * std::real(aux_gk_host[ih * nh_atom + jh]);
+                                    this->becsum_[index + ijh] += 2.0 * static_cast<double>(std::real(aux_gk_host[ih * nh_atom + jh]));
                                 }
                                 ijh++;
                             }
@@ -470,7 +465,7 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
     // add to the charge density in reciprocal space the part which is due to the US augmentation.
     if (PARAM.globalv.use_uspp)
     {
-        this->addusdens_g(becsum, this->charge->rhog);
+        this->addusdens_g(this->charge->rhog);
     }
     // transform back to real space using dense grids
     if (PARAM.globalv.double_grid || PARAM.globalv.use_uspp)
@@ -483,13 +478,14 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
 }
 
 template <typename T, typename Device>
-void ElecStatePW<T, Device>::addusdens_g(const Real* becsum, std::complex<double>** rhog)
+void ElecStatePW<T, Device>::addusdens_g(std::complex<double>** rhog)
 {
     const T one{1, 0};
     const T zero{0, 0};
     const int npw = this->charge->rhopw->npw;
     const int lmaxq = this->ppcell->lmaxq;
     const int nh_tot = this->ppcell->nhm * (this->ppcell->nhm + 1) / 2;
+    const double* becsum = this->becsum_.data();
     Structure_Factor* psf = this->ppcell->psf;
     const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;
 
@@ -577,11 +573,35 @@ void ElecStatePW<T, Device>::addusdens_g(const Real* becsum, std::complex<double
     delmem_var_op()(ylmk0);
 }
 
+// Taoni add 2026-09-02
+// Added to fix USPP single force/stress reading the former float becsum as double.
+// The double-only drivers receive a base ElecState, while becsum belongs to the precision-templated ElecStatePW.
+// Refactor this bridge for true float force/stress.
+template <typename Device>
+const std::vector<double>* get_becsum(const ElecState& elec)
+{
+    const ElecStatePW<std::complex<double>, Device>* double_elec = dynamic_cast<const ElecStatePW<std::complex<double>, Device>*>(&elec);
+    if (double_elec != nullptr)
+    {
+        return &double_elec->get_becsum();
+    }
+
+    const ElecStatePW<std::complex<float>, Device>* single_elec = dynamic_cast<const ElecStatePW<std::complex<float>, Device>*>(&elec);
+    if (single_elec != nullptr)
+    {
+        return &single_elec->get_becsum();
+    }
+
+    return nullptr;
+}
+
 template class ElecStatePW<std::complex<float>, base_device::DEVICE_CPU>;
 template class ElecStatePW<std::complex<double>, base_device::DEVICE_CPU>;
+template const std::vector<double>* get_becsum<base_device::DEVICE_CPU>(const ElecState& elec);
 #if ((defined __CUDA) || (defined __ROCM))
 template class ElecStatePW<std::complex<float>, base_device::DEVICE_GPU>;
 template class ElecStatePW<std::complex<double>, base_device::DEVICE_GPU>;
+template const std::vector<double>* get_becsum<base_device::DEVICE_GPU>(const ElecState& elec);
 #endif 
 
 } // namespace elecstate

source/source_estate/elecstate_pw.h

@@ -1,6 +1,8 @@
 #ifndef ELECSTATEPW_H
 #define ELECSTATEPW_H
 
+#include <vector>
+
 #include <source_base/macros.h>
 
 #include "elecstate.h"
@@ -41,7 +43,15 @@ class ElecStatePW : public ElecState
     //! calculate becsum for uspp
     void cal_becsum(const psi::Psi<T, Device>& psi);
 
-    Real* becsum = nullptr;
+    /**
+     * @brief Return the USPP projector occupancy coefficients.
+     *
+     * @return Read-only canonical double-precision coefficients.
+     */
+    const std::vector<double>& get_becsum() const
+    {
+        return becsum_;
+    }
 
     //! init rho_data and kin_r_data
     void init_rho_data();
@@ -74,7 +84,7 @@ class ElecStatePW : public ElecState
 
     //! Non-local pseudopotentials
     //! \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
-    void addusdens_g(const Real* becsum, std::complex<double>** rhog);
+    void addusdens_g(std::complex<double>** rhog);
 
     Device * ctx = {};
 
@@ -89,6 +99,8 @@ class ElecStatePW : public ElecState
     T* wfcr_another_spin = nullptr;
 
   private:
+    std::vector<double> becsum_;
+
     using meta_op = hamilt::meta_pw_op<Real, Device>;
     using elecstate_pw_op = elecstate::elecstate_pw_op<Real, Device>;
 
@@ -106,16 +118,22 @@ class ElecStatePW : public ElecState
     using syncmem_complex_d2h_op = base_device::memory::synchronize_memory_op<T, base_device::DEVICE_CPU, Device>;
     using syncmem_complex_h2d_op = base_device::memory::synchronize_memory_op<T, Device, base_device::DEVICE_CPU>;
 
-    using resmem_var_h_op = base_device::memory::resize_memory_op<Real, base_device::DEVICE_CPU>;
-    using delmem_var_h_op = base_device::memory::delete_memory_op<Real, base_device::DEVICE_CPU>;
-    using setmem_var_h_op = base_device::memory::set_memory_op<Real, base_device::DEVICE_CPU>;
     using syncmem_var_h2d_op = base_device::memory::synchronize_memory_op<Real, Device, base_device::DEVICE_CPU>;
     using syncmem_var_d2h_op = base_device::memory::synchronize_memory_op<Real, base_device::DEVICE_CPU, Device>;
 
     using gemv_op = ModuleBase::gemv_op<T, Device>;
     using gemm_op = ModuleBase::gemm_op<T, Device>;
 };
 
+/**
+ * @brief Return the USPP projector occupancy coefficients of a PW electronic state.
+ *
+ * @param elec Electronic state to inspect.
+ * @return Read-only coefficients, or nullptr when the state is not a supported PW state for Device.
+ */
+template <typename Device>
+const std::vector<double>* get_becsum(const ElecState& elec);
+
 } // namespace elecstate
 
 #endif

source/source_estate/test/elecstate_pw_test.cpp

@@ -242,6 +242,7 @@ TEST_F(ElecStatePWTest, ConstructorDouble)
     EXPECT_EQ(elecstate_pw_d->charge, chg);
     EXPECT_EQ(elecstate_pw_d->klist, klist);
     EXPECT_EQ(elecstate_pw_d->bigpw, bigpw);
+    EXPECT_TRUE(elecstate_pw_d->get_becsum().empty());
 }
 
 TEST_F(ElecStatePWTest, ConstructorSingle)
@@ -257,6 +258,7 @@ TEST_F(ElecStatePWTest, ConstructorSingle)
     EXPECT_EQ(elecstate_pw_s->charge, chg);
     EXPECT_EQ(elecstate_pw_s->klist, klist);
     EXPECT_EQ(elecstate_pw_s->bigpw, bigpw);
+    EXPECT_TRUE(elecstate_pw_s->get_becsum().empty());
 }
 
 TEST_F(ElecStatePWTest, InitRhoDataDouble)