fix ig=0

Author: A-006

source/module_elecstate/module_pot/H_Hartree_pw.cpp

@@ -44,18 +44,18 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell,
     double ehart = 0.0;
 
     std::vector<std::complex<double>> vh_g(rho_basis->npw);
+    const int ig0 = rho_basis->ig_gge0;
 #ifdef _OPENMP
 #pragma omp parallel for reduction(+:ehart)
 #endif
     for (int ig = 0; ig < rho_basis->npw; ig++)
     {
-        if (rho_basis->gg[ig] >= 1.0e-8) // LiuXh 20180410
-        {
+        if (ig == ig0) continue; // skip G=0
             const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (cell.tpiba2 * rho_basis->gg[ig]);
 
             ehart += (conj(Porter[ig]) * Porter[ig]).real() * fac;
             vh_g[ig] = fac * Porter[ig];
-        }
+        
     }
 
     Parallel_Reduce::reduce_pool(ehart);

source/module_hamilt_general/module_ewald/H_Ewald_pw.cpp

@@ -124,10 +124,10 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
 	fact = 1.0;
 
     //GlobalV::ofs_running << "\n pwb.gstart = " << pwb.gstart << std::endl;
-
+    const int ig0 = rho_basis->ig_gge0;
     for (int ig = 0; ig < rho_basis->npw; ig++)
     {
-        if(ig == rho_basis->ig_gge0) 
+        if(ig == ig0) 
         { 
             continue;
         }

source/module_hamilt_general/module_surchem/cal_pseudo.cpp

@@ -10,6 +10,7 @@ void surchem::gauss_charge(const UnitCell& cell,
                            Structure_Factor* sf)
 {
     sf->setup_structure_factor(&cell, pgrid, rho_basis); // call strucFac(ntype,ngmc)
+    const int ig0 = rho_basis->ig_gge0; // G=0 index
     for (int it = 0; it < cell.ntype; it++)
     {
         double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];

source/module_hamilt_general/module_surchem/cal_totn.cpp

@@ -12,8 +12,9 @@ void surchem::cal_totn(const UnitCell& cell,
     ModuleBase::GlobalFunc::ZEROS(vloc_g, rho_basis->npw);
 
     rho_basis->real2recip(vlocal, vloc_g);  // now n is vloc in Recispace
+    const int ig0 = rho_basis->ig_gge0; // G=0 index
     for (int ig = 0; ig < rho_basis->npw; ig++) {
-        if(ig==rho_basis->ig_gge0)
+        if(ig==ig0)
         {
             N[ig] = Porter_g[ig];
             continue;
@@ -39,9 +40,10 @@ void surchem::induced_charge(const UnitCell& cell, const ModulePW::PW_Basis* rho
     std::complex<double> *induced_rhog = new complex<double>[rho_basis->npw];
     ModuleBase::GlobalFunc::ZEROS(induced_rhog, rho_basis->npw);
     rho_basis->real2recip(delta_phi, delta_phig);
+    const int ig0 = rho_basis->ig_gge0; // G=0 index
     for (int ig = 0; ig < rho_basis->npw; ig++)
     {   
-        if(rho_basis->ig_gge0 == ig)
+        if(ig==ig0)
         {
             continue;
         }

source/module_hamilt_pw/hamilt_pwdft/forces.cpp

@@ -472,13 +472,17 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
         upperbound = 2.0 * charge * charge * sqrt(2.0 * alpha / ModuleBase::TWO_PI)
                      * erfc(sqrt(ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha));
     } while (upperbound > 1.0e-6);
-
+    const int ig0 = rho_basis->ig_gge0;
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
     for (int ig = 0; ig < rho_basis->npw; ig++)
     {
-        if (rho_basis->gg[ig] !=0)
+        if (ig== ig0)
+        {
+            continue; // skip G=0
+        }
+        else
         {
             aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0)
                     / (rho_basis->gg[ig] * ucell.tpiba2);

source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp

@@ -159,8 +159,7 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 			for(int ig = 0;ig< rho_basis->npw;ig++)
 			{
 				const FPTYPE norm_g = sqrt(rho_basis->gg[ig]);
-				if(norm_g < 1e-4) { 	continue;
-}
+				if(norm_g < 1e-4) { 	continue;}
 				for (int l = 0; l < 3; l++)
 				{
 					for (int m = 0;m< 3;m++)