Fix the incorrect calculation of partial charge for multi-k-points under both PW and LCAO when symmetry is turned on

source/Makefile.Objects

@@ -517,8 +517,8 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_proj_band_lcao.o\
       write_istate_info.o\
       nscf_fermi_surf.o\
-      istate_charge.o\
-      istate_envelope.o\
+      get_pchg.o\
+      get_wf.o\
       io_dmk.o\
       unk_overlap_lcao.o\
       read_wfc_nao.o\

source/module_elecstate/module_charge/symmetry_rho.cpp

@@ -1,14 +1,13 @@
 #include "symmetry_rho.h"
+
 #include "module_hamilt_general/module_xc/xc_functional.h"
 
 Symmetry_rho::Symmetry_rho()
 {
-
 }
 
 Symmetry_rho::~Symmetry_rho()
 {
-
 }
 
 void Symmetry_rho::begin(const int& spin_now,
@@ -17,96 +16,137 @@ void Symmetry_rho::begin(const int& spin_now,
                          Parallel_Grid& Pgrid,
                          ModuleSymmetry::Symmetry& symm) const
 {
-	assert(spin_now < 4);//added by zhengdy-soc
+    assert(spin_now < 4); // added by zhengdy-soc
+
+    if (ModuleSymmetry::Symmetry::symm_flag != 1) {
+        return;
+}
+    // both parallel and serial
+    // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
+    // {
+    // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
+    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now],
+    // rho_basis,Pgrid,symm);
+    // }
+    // else	//space group, do rho_symm in reciprocal space
+    {
+        rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]);
+        psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm); // need to modify
+        rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]);
+
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // Use std::vector to manage kin_g instead of raw pointer
+            std::vector<std::complex<double>> kin_g(CHR.ngmc);
+            rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data());
+            psymmg(kin_g.data(), rho_basis, Pgrid, symm);
+            rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]);
+        }
+    }
+    return;
+}
+
+void Symmetry_rho::begin(const int& spin_now,
+                         double** rho,
+                         std::complex<double>** rhog,
+                         int ngmc,
+                         double** kin_r,
+                         const ModulePW::PW_Basis* rho_basis,
+                         Parallel_Grid& Pgrid,
+                         ModuleSymmetry::Symmetry& symm) const
+{
+    assert(spin_now < 4); // added by zhengdy-soc
 
-	if(ModuleSymmetry::Symmetry::symm_flag != 1) return;
-	// both parallel and serial
+    if (ModuleSymmetry::Symmetry::symm_flag != 1)
+    {
+        return;
+    }
+    // both parallel and serial
     // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
     // {
     // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
-    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], rho_basis,Pgrid,symm);
+    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now],
+    // rho_basis,Pgrid,symm);
     // }
     // else	//space group, do rho_symm in reciprocal space
-	{
-		rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]);
-		psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm);	//need to modify
-		rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]);
+    {
+        rho_basis->real2recip(rho[spin_now], rhog[spin_now]);
+        psymmg(rhog[spin_now], rho_basis, Pgrid, symm);
+        rho_basis->recip2real(rhog[spin_now], rho[spin_now]);
 
-		if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) 
-		{
-			std::complex<double>* kin_g = new std::complex<double>[CHR.ngmc];
-			rho_basis->real2recip(CHR.kin_r[spin_now], kin_g);
-			psymmg(kin_g,  rho_basis,Pgrid,symm);
-			rho_basis->recip2real(kin_g, CHR.kin_r[spin_now]);
-			delete[] kin_g;
-		}	
-	}
-	return;
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // Use std::vector to manage kin_g instead of raw pointer
+            std::vector<std::complex<double>> kin_g(ngmc);
+            rho_basis->real2recip(kin_r[spin_now], kin_g.data());
+            psymmg(kin_g.data(), rho_basis, Pgrid, symm);
+            rho_basis->recip2real(kin_g.data(), kin_r[spin_now]);
+        }
+    }
+    return;
 }
 
-void Symmetry_rho::psymm(double* rho_part, const ModulePW::PW_Basis *rho_basis, Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const
+void Symmetry_rho::psymm(double* rho_part,
+                         const ModulePW::PW_Basis* rho_basis,
+                         Parallel_Grid& Pgrid,
+                         ModuleSymmetry::Symmetry& symm) const
 {
 #ifdef __MPI
-	// (1) reduce all rho from the first pool.
-	double* rhotot;
-	if(GlobalV::MY_RANK == 0)
-	{
-		rhotot = new double[rho_basis->nxyz];
-		ModuleBase::GlobalFunc::ZEROS(rhotot, rho_basis->nxyz);
-	}
-	Pgrid.reduce_to_fullrho(rhotot, rho_part);
+    // (1) reduce all rho from the first pool.
+    std::vector<double> rhotot;
+    if (GlobalV::MY_RANK == 0)
+    {
+        rhotot.resize(rho_basis->nxyz);
+        ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz);
+    }
+    Pgrid.reduce_to_fullrho(rhotot.data(), rho_part);
 
-	// (2)
-	if(GlobalV::MY_RANK==0)
-	{
-		symm.rho_symmetry(rhotot, rho_basis->nx, rho_basis->ny, rho_basis->nz);
+    // (2)
+    if (GlobalV::MY_RANK == 0)
+    {
+        symm.rho_symmetry(rhotot.data(), rho_basis->nx, rho_basis->ny, rho_basis->nz);
 #else
-		symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz);
+    symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz);
 #endif
-		/*
-		int count = 0;
-		GlobalV::ofs_running << scientific;
-		for(int iz=0; iz<rho_basis->nz; iz++)
-		{
-			GlobalV::ofs_running << "\n iz=" << iz;
-			for(int iy=0; iy<rho_basis->ny; iy++)
-			{
-				for(int ix=0; ix<rho_basis->nx; ix++)
-				{
-					if(count%5==0) GlobalV::ofs_running << "\n";
-					++count;
-					GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz];
-				}
-			}
-		}
-		*/
-	#ifdef __MPI
-	}
-
-	// (3)
-	const int ncxy = rho_basis->nx * rho_basis->ny;
-	double* zpiece = new double[ncxy];
-	for(int iz=0; iz<rho_basis->nz; iz++)
-	{
-		//GlobalV::ofs_running << "\n iz=" << iz;
-		ModuleBase::GlobalFunc::ZEROS(zpiece, ncxy);
-		if(GlobalV::MY_RANK==0)
-		{
-			for(int ix=0; ix<rho_basis->nx; ix++)
-			{
-				for(int iy=0; iy<rho_basis->ny; iy++)
-				{
-					const int ir = ix * rho_basis->ny + iy;
-					zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz];
-					//rho[ir*nczp+znow] = zpiece[ir];
-				}
-			}
-		}
-		Pgrid.zpiece_to_all(zpiece,iz, rho_part);
-	}
+        /*
+        int count = 0;
+        GlobalV::ofs_running << scientific;
+        for(int iz=0; iz<rho_basis->nz; iz++)
+        {
+            GlobalV::ofs_running << "\n iz=" << iz;
+            for(int iy=0; iy<rho_basis->ny; iy++)
+            {
+                for(int ix=0; ix<rho_basis->nx; ix++)
+                {
+                    if(count%5==0) GlobalV::ofs_running << "\n";
+                    ++count;
+                    GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz];
+                }
+            }
+        }
+        */
+#ifdef __MPI
+    }
 
-	if(GlobalV::MY_RANK==0)		delete[] rhotot;
-	delete[] zpiece;
+    // (3)
+    const int ncxy = rho_basis->nx * rho_basis->ny;
+    std::vector<double> zpiece(ncxy);
+    for (int iz = 0; iz < rho_basis->nz; iz++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(zpiece.data(), ncxy);
+        if (GlobalV::MY_RANK == 0)
+        {
+            for (int ix = 0; ix < rho_basis->nx; ix++)
+            {
+                for (int iy = 0; iy < rho_basis->ny; iy++)
+                {
+                    const int ir = ix * rho_basis->ny + iy;
+                    zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz];
+                }
+            }
+        }
+        Pgrid.zpiece_to_all(zpiece.data(), iz, rho_part);
+    }
 #endif
-	return;
-}
+    return;
+}

source/module_elecstate/module_charge/symmetry_rho.h

@@ -1,37 +1,57 @@
 #ifndef SYMMETRY_RHO_H
 #define SYMMETRY_RHO_H
-#include "module_elecstate/module_charge/charge.h"
 #include "module_basis/module_pw/pw_basis.h"
-#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h"
-
 #include "module_cell/module_symmetry/symmetry.h"
+#include "module_elecstate/module_charge/charge.h"
+#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h"
 
 class Symmetry_rho
 {
-	public:
-	Symmetry_rho();
-	~Symmetry_rho();
+  public:
+    Symmetry_rho();
+    ~Symmetry_rho();
 
     void begin(const int& spin_now,
                const Charge& CHR,
                const ModulePW::PW_Basis* pw,
                Parallel_Grid& Pgrid,
                ModuleSymmetry::Symmetry& symm) const;
 
+    void begin(const int& spin_now,
+               double** rho,
+               std::complex<double>** rhog,
+               int ngmc,
+               double** kin_r,
+               const ModulePW::PW_Basis* pw,
+               Parallel_Grid& Pgrid,
+               ModuleSymmetry::Symmetry& symm) const;
+
   private:
-	// in real space:
-	void psymm(double *rho_part, const ModulePW::PW_Basis *pw, Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const;
-	// in reciprocal space:
-	void psymmg(std::complex<double>* rhog_part, const ModulePW::PW_Basis *rho_basis, 
-			Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const;
+    // in real space:
+    void psymm(double* rho_part,
+               const ModulePW::PW_Basis* pw,
+               Parallel_Grid& Pgrid,
+               ModuleSymmetry::Symmetry& symm) const;
+    // in reciprocal space:
+    void psymmg(std::complex<double>* rhog_part,
+                const ModulePW::PW_Basis* rho_basis,
+                Parallel_Grid& Pgrid,
+                ModuleSymmetry::Symmetry& symm) const;
 #ifdef __MPI
-	void reduce_to_fullrhog(const ModulePW::PW_Basis *rho_basis, std::complex<double>* rhogtot, 
-			std::complex<double>* rhogin, int* ig2isztot, const int* ig2iszin, int max_npw) const;
-	void rhog_piece_to_all(const ModulePW::PW_Basis *rho_basis,
-			std::complex<double>* rhogtot, std::complex<double>* rhog_part) const;
+    void reduce_to_fullrhog(const ModulePW::PW_Basis* rho_basis,
+                            std::complex<double>* rhogtot,
+                            std::complex<double>* rhogin,
+                            int* ig2isztot,
+                            const int* ig2iszin,
+                            int max_npw) const;
+    void rhog_piece_to_all(const ModulePW::PW_Basis* rho_basis,
+                           std::complex<double>* rhogtot,
+                           std::complex<double>* rhog_part) const;
 #endif
-	void get_ixyz2ipw(const ModulePW::PW_Basis *rho_basis, const int* ig2isztot, 
-			const int* fftixy2is, int* ixyz2ipw) const;	//(ix, iy, iz) -> (ip, ig)
+    void get_ixyz2ipw(const ModulePW::PW_Basis* rho_basis,
+                      const int* ig2isztot,
+                      const int* fftixy2is,
+                      int* ixyz2ipw) const; //(ix, iy, iz) -> (ip, ig)
 };
 
 #endif

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -8,8 +8,8 @@
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_io/berryphase.h"
-#include "module_io/istate_charge.h"
-#include "module_io/istate_envelope.h"
+#include "module_io/get_pchg.h"
+#include "module_io/get_wf.h"
 #include "module_io/to_wannier90_lcao.h"
 #include "module_io/to_wannier90_lcao_in_pw.h"
 #include "module_io/write_HS_R.h"
@@ -479,8 +479,10 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
         } else {
             ISC.begin(this->GK,
                     this->pelec->charge->rho,
+                    this->pelec->charge->rhog,
                     this->pelec->wg,
                     this->pelec->eferm.get_all_ef(),
+                    this->pw_rho,
                     this->pw_rho->nrxx,
                     this->pw_rho->nplane,
                     this->pw_rho->startz_current,
@@ -502,7 +504,9 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
                     &GlobalC::ucell,
                     &GlobalC::GridD,
                     this->kv,
-                    PARAM.inp.if_separate_k);
+                    PARAM.inp.if_separate_k,
+                    &GlobalC::Pgrid,
+                    this->pelec->charge->ngmc);
         }
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge");
     } else if (cal_type == "get_wf") {