Fix: instable nonlocal pp in uspp (#5798)

Author: YuLiu98

source/module_esolver/esolver_ks_pw.cpp

@@ -255,7 +255,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
     if (ucell.cell_parameter_updated)
     {
-        this->ppcell.init_vnl(ucell, this->pw_rhod);
+        this->ppcell.rescale_vnl(ucell.omega);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL");
 
         this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz);

source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp

@@ -550,6 +550,8 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_
     ModuleBase::TITLE("pseudopot_cell_vnl", "init_vnl");
     ModuleBase::timer::tick("ppcell_vnl", "init_vnl");
 
+    this->omega_old = cell.omega;
+
     // from init_us_1
     //    a) For each non vanderbilt pseudopotential it computes the D and
     //       the betar in the same form of the Vanderbilt pseudopotential.
@@ -1733,6 +1735,27 @@ void pseudopot_cell_vnl::newd_nc(const int& iat, UnitCell& cell)
     }
 }
 
+// scale the non-local pseudopotential tables
+void pseudopot_cell_vnl::rescale_vnl(const double& omega_in)
+{
+    const double ratio = this->omega_old / omega_in;
+    const double sqrt_ratio = std::sqrt(ratio);
+    this->omega_old = omega_in;
+
+    for (int i = 0; i < this->tab.getSize(); i++)
+    {
+        this->tab.ptr[i] *= sqrt_ratio;
+    }
+    for (int i = 0; i < this->tab_at.getSize(); i++)
+    {
+        this->tab_at.ptr[i] *= sqrt_ratio;
+    }
+    for (int i = 0; i < this->qrad.getSize(); i++)
+    {
+        this->qrad.ptr[i] *= ratio;
+    }
+}
+
 template <>
 float* pseudopot_cell_vnl::get_nhtol_data() const
 {

source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h

@@ -36,6 +36,8 @@ class pseudopot_cell_vnl
 
     void init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_basis);
 
+    void rescale_vnl(const double& omega_in);
+
     template <typename FPTYPE, typename Device>
     void getvnl(Device* ctx, const UnitCell& ucell, const int& ik, std::complex<FPTYPE>* vkb_in) const;
 
@@ -200,6 +202,8 @@ class pseudopot_cell_vnl
 
     Soc soc;
 
+    double omega_old = 0;
+
     /**
      * @brief Compute interpolation table qrad
      *