Fix: enable force and stress in uspp nscf

YuLiu98 · Dec 23, 2024 · 9c46b77 · 9c46b77
1 parent 178c24c
commit 9c46b77
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 2 deletions.
diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp
@@ -238,7 +238,7 @@ void ElecStatePW<T, Device>::rhoBandK(const psi::Psi<T, Device>& psi)
 }
 
 template <typename T, typename Device>
-void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
+void ElecStatePW<T, Device>::cal_becsum(const psi::Psi<T, Device>& psi)
 {
     const T one{1, 0};
     const T zero{0, 0};
@@ -392,6 +392,12 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
         }
     }
     delmem_complex_op()(this->ctx, becp);
+}
+
+template <typename T, typename Device>
+void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
+{
+    this->cal_becsum(psi);
 
     // transform soft charge to recip space using smooth grids
     T* rhog = nullptr;

diff --git a/source/module_elecstate/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h
@@ -35,6 +35,9 @@ class ElecStatePW : public ElecState
 
     virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
+    //! calculate becsum for uspp
+    void cal_becsum(const psi::Psi<T, Device>& psi);
+
     Real* becsum = nullptr;
 
     //! init rho_data and kin_r_data
@@ -61,7 +64,7 @@ class ElecStatePW : public ElecState
 
     //! calcualte rho for each k
     void rhoBandK(const psi::Psi<T, Device>& psi);
-    
+
     //! add to the charge density in reciprocal space the part which is due to the US augmentation.
     void add_usrho(const psi::Psi<T, Device>& psi);
 

diff --git a/source/module_hsolver/hsolver_lcaopw.cpp b/source/module_hsolver/hsolver_lcaopw.cpp
@@ -281,6 +281,10 @@ void HSolverLIP<T>::solve(hamilt::Hamilt<T>* pHamilt, // ESolver_KS_PW::p_hamilt
     reinterpret_cast<elecstate::ElecStatePW<T>*>(pes)->calEBand();
     if (skip_charge)
     {
+        if (PARAM.globalv.use_uspp)
+        {
+            reinterpret_cast<elecstate::ElecStatePW<T>*>(pes)->cal_becsum(psi);
+        }
         ModuleBase::timer::tick("HSolverLIP", "solve");
         return;
     }

diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
@@ -337,6 +337,10 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
     reinterpret_cast<elecstate::ElecStatePW<T>*>(pes)->calEBand();
     if (skip_charge)
     {
+        if (PARAM.globalv.use_uspp)
+        {
+            reinterpret_cast<elecstate::ElecStatePW<T, Device>*>(pes)->cal_becsum(psi);
+        }
         ModuleBase::timer::tick("HSolverPW", "solve");
         return;
     }