fix bug of stress and force

Author: Qianruipku

source/module_hamilt_pw/hamilt_pwdft/forces_nl.cpp

@@ -33,6 +33,7 @@ void Forces<FPTYPE, Device>::cal_force_nl(ModuleBase::matrix& forcenl,
     hamilt::FS_Nonlocal_tools<FPTYPE, Device> nl_tools(nlpp_in, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb);
 
     const int nks = wfc_basis->nks;
+    const int max_nbands = wg.nc;
     for (int ik = 0; ik < nks; ik++) // loop k points
     {
         // skip zero weights to speed up
@@ -46,19 +47,19 @@ void Forces<FPTYPE, Device>::cal_force_nl(ModuleBase::matrix& forcenl,
             }
         }
         const int npm = nbands_occ;
-        nl_tools.cal_vkb(ik, npm);
+        nl_tools.cal_vkb(ik, max_nbands);
         // calculate becp = <psi|beta> for all beta functions
         nl_tools.cal_becp(ik, npm, &psi_in[0](ik,0,0));
-        nl_tools.reduce_pool_becp(npm);
+        nl_tools.reduce_pool_becp(max_nbands);
         for (int ipol = 0; ipol < 3; ipol++)
         {
-            nl_tools.cal_vkb_deri_f(ik, npm, ipol);
+            nl_tools.cal_vkb_deri_f(ik, max_nbands, ipol);
             // calculate dbecp = <psi|\nabla beta> for all beta functions
-            nl_tools.cal_dbecp_f(ik, npm, npm, ipol, &psi_in[0](ik,0,0));
+            nl_tools.cal_dbecp_f(ik, max_nbands, npm, ipol, &psi_in[0](ik,0,0));
             nl_tools.revert_vkb(ik, ipol);
         }
         // calculate the force_i = \sum_{n,k}f_{nk}\sum_I \sum_{lm,l'm'}D_{l,l'}^{I} becp * dbecp_i
-        nl_tools.cal_force(ik, npm, npm, true, force);
+        nl_tools.cal_force(ik, max_nbands, npm, true, force);
     } // end ik
 
     syncmem_var_d2h_op()(this->cpu_ctx, this->ctx, forcenl.c, force, forcenl.nr * forcenl.nc);

source/module_hamilt_pw/hamilt_pwdft/stress_func_nl.cpp

@@ -37,6 +37,7 @@ void Stress_Func<FPTYPE, Device>::stress_nl(ModuleBase::matrix& sigma,
     hamilt::FS_Nonlocal_tools<FPTYPE, Device> nl_tools(nlpp_in, &ucell_in, p_kv, wfc_basis, p_sf, wg, &ekb);
 
     const int nks = p_kv->get_nks();
+    const int max_nbands = wg.nc;
     for (int ik = 0; ik < nks; ik++) // loop k points
     {
         // skip zero weights to speed up
@@ -51,17 +52,17 @@ void Stress_Func<FPTYPE, Device>::stress_nl(ModuleBase::matrix& sigma,
         }
         const int npm = nbands_occ;
 
-        nl_tools.cal_vkb(ik, npm);
+        nl_tools.cal_vkb(ik, max_nbands);
         // calculate becp = <psi|beta> for all beta functions
         nl_tools.cal_becp(ik, npm, &psi_in[0](ik,0,0));
-        nl_tools.reduce_pool_becp(npm);
+        nl_tools.reduce_pool_becp(max_nbands);
         // calculate dbecp = <psi|d(beta)/dR> for all beta functions
         // calculate stress = \sum <psi|d(beta_j)/dR> * <psi|beta_i> * D_{ij}
         for (int ipol = 0; ipol < 3; ipol++)
         {
             for (int jpol = 0; jpol <= ipol; jpol++)
             {
-                nl_tools.cal_vkb_deri_s(ik, npm, ipol, jpol);
+                nl_tools.cal_vkb_deri_s(ik, max_nbands, ipol, jpol);
                 nl_tools.cal_dbecp_s(ik, npm, &psi_in[0](ik,0,0));
                 nl_tools.cal_stress(ik, npm, true, ipol, jpol, stress_device);
             }

tests/integrate/187_PW_SDFT_ALL_GPU/INPUT

@@ -9,7 +9,7 @@ symmetry		0
 pseudo_dir	../../PP_ORB
 device          gpu
 
-kpar            1
+kpar            2
 bndpar          1
 
 nbands			0
@@ -21,7 +21,7 @@ cal_force       1
 cal_stress      1
 
 #Parameters (2.Iteration)
-ecutwfc			30
+ecutwfc			20
 scf_thr			1e-6
 scf_nmax		30
 

tests/integrate/187_PW_SDFT_ALL_GPU/result.ref

@@ -1,5 +1,5 @@
-etotref -108.9478430350318803
-etotperatomref -54.4739215175
-totalforceref 187.066308
-totalstressref 248810.173686
-totaltimeref 40.81
+etotref -105.2612355454176907
+etotperatomref -52.6306177727
+totalforceref 197.906706
+totalstressref 254537.682905
+totaltimeref 20.82