Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop

 into hotfix

source/module_elecstate/module_charge/charge.cpp

@@ -273,7 +273,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, const ModuleBase::ComplexMa
     {
         for (int is = 0; is < GlobalV::NSPIN; ++is)
         {
-            GlobalC::restart.load_disk("charge", is, rho);
+            GlobalC::restart.load_disk("charge", is, this->nrxx, rho);
         }
         GlobalC::restart.info_load.load_charge_finish = true;
     }

source/module_esolver/esolver_ks_lcao.cpp

@@ -745,7 +745,7 @@ void ESolver_KS_LCAO::eachiterfinish(int iter)
     {
         for (int is = 0; is < GlobalV::NSPIN; ++is)
         {
-            GlobalC::restart.save_disk(*this->UHM.LM, "charge", is, pelec->charge->rho);
+            GlobalC::restart.save_disk(*this->UHM.LM, "charge", is, pelec->charge->nrxx, pelec->charge->rho);
         }
     }
 

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -359,7 +359,7 @@ namespace ModuleESolver
 
         if(GlobalV::CALCULATION == "test_memory")
         {
-            Cal_Test::test_memory();
+            Cal_Test::test_memory(this->pw_rho, this->pw_wfc, GlobalC::CHR_MIX.get_mixing_mode(), GlobalC::CHR_MIX.get_mixing_ndim());
             return;
         }
 
@@ -395,20 +395,15 @@ namespace ModuleESolver
         else if (GlobalV::CALCULATION == "istate")
         {
             IState_Charge ISC(this->psid, this->LOC);
-            ISC.begin(this->UHM.GG, this->pelec);
+            ISC.begin(this->UHM.GG, this->pelec, this->pw_rho, GlobalC::bigpw);
         }
         else if (GlobalV::CALCULATION == "ienvelope")
         {
             IState_Envelope IEP(this->pelec);
             if (GlobalV::GAMMA_ONLY_LOCAL)
-                IEP.begin(this->psid,
-                          this->LOWF,
-                          this->UHM.GG,
-                          INPUT.out_wfc_pw,
-                          GlobalC::wf.out_wfc_r,
-                          GlobalC::kv);
+                IEP.begin(this->psid, this->pw_rho, this->pw_wfc, GlobalC::bigpw, this->LOWF, this->UHM.GG, INPUT.out_wfc_pw, GlobalC::wf.out_wfc_r, GlobalC::kv);
             else
-                IEP.begin(this->psi, this->LOWF, this->UHM.GK, INPUT.out_wfc_pw, GlobalC::wf.out_wfc_r, GlobalC::kv);
+                IEP.begin(this->psi, this->pw_rho, this->pw_wfc, GlobalC::bigpw, this->LOWF, this->UHM.GK, INPUT.out_wfc_pw, GlobalC::wf.out_wfc_r, GlobalC::kv);
         }
         else
         {
@@ -575,14 +570,14 @@ namespace ModuleESolver
         if (GlobalV::CALCULATION == "nscf" && INPUT.towannier90)
         {
             toWannier90 myWannier(GlobalC::kv.nkstot, GlobalC::ucell.G, this->LOWF.wfc_k_grid);
-            myWannier.init_wannier(this->pelec->ekb, GlobalC::kv, nullptr);
+            myWannier.init_wannier(this->pelec->ekb, this->pw_rho, this->pw_wfc, GlobalC::bigpw, GlobalC::kv, nullptr);
         }
 
         // add by jingan
         if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1)
         {
             berryphase bp(this->LOWF);
-            bp.Macroscopic_polarization(this->psi,GlobalC::kv);
+            bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, GlobalC::kv);
         }
 
         //below is for DeePKS NSCF calculation

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -465,7 +465,7 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
         {
             std::stringstream ss_dipole;
             ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DIPOLE";
-            ModuleIO::write_dipole(pelec->charge->rho_save[is], is, istep, ss_dipole.str());
+            ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str());
         }
 
         std::stringstream ssd;

source/module_esolver/esolver_ks_pw.cpp

@@ -329,7 +329,10 @@ void ESolver_KS_PW<FPTYPE, Device>::othercalculation(const int istep)
     ModuleBase::timer::tick("ESolver_KS_PW", "othercalculation");
     if (GlobalV::CALCULATION == "test_memory")
     {
-        Cal_Test::test_memory();
+        Cal_Test::test_memory(this->pw_rho,
+                              this->pw_wfc,
+                              GlobalC::CHR_MIX.get_mixing_mode(),
+                              GlobalC::CHR_MIX.get_mixing_ndim());
         return;
     }
 
@@ -854,14 +857,14 @@ void ESolver_KS_PW<FPTYPE, Device>::postprocess()
         if (winput::out_spillage <= 2)
         {
             Numerical_Basis numerical_basis;
-            numerical_basis.output_overlap(this->psi[0], GlobalC::sf, GlobalC::kv);
+            numerical_basis.output_overlap(this->psi[0], GlobalC::sf, GlobalC::kv, GlobalC::wfcpw);
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "BASIS OVERLAP (Q and S) GENERATION.");
         }
     }
 
     if (GlobalC::wf.out_wfc_r == 1) // Peize Lin add 2021.11.21
     {
-        ModuleIO::write_psi_r_1(this->psi[0], "wfc_realspace", true, GlobalC::kv);
+        ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, GlobalC::kv);
     }
 
     if (INPUT.cal_cond)
@@ -949,7 +952,7 @@ void ESolver_KS_PW<FPTYPE, Device>::nscf()
     if (INPUT.towannier90)
     {
         toWannier90 myWannier(GlobalC::kv.nkstot, GlobalC::ucell.G);
-        myWannier.init_wannier(this->pelec->ekb, GlobalC::kv, this->psi);
+        myWannier.init_wannier(this->pelec->ekb, this->pw_rho, this->pw_wfc, GlobalC::bigpw, GlobalC::kv, this->psi);
     }
 
     //=======================================================
@@ -959,7 +962,7 @@ void ESolver_KS_PW<FPTYPE, Device>::nscf()
     if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1)
     {
         berryphase bp;
-        bp.Macroscopic_polarization(this->psi, GlobalC::kv);
+        bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, GlobalC::kv);
     }
 
     ModuleBase::timer::tick("ESolver_KS_PW", "nscf");

source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp

@@ -248,7 +248,7 @@ const double *rab, const int &l, double* table)
 }
 
 void Wavefunc_in_pw::produce_local_basis_in_pw(const int &ik,
-                                               ModulePW::PW_Basis_K *wfc_basis,
+                                               const ModulePW::PW_Basis_K *wfc_basis,
                                                ModuleBase::ComplexMatrix &psi,
                                                const ModuleBase::realArray &table_local)
 {

source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h

@@ -34,7 +34,7 @@ namespace Wavefunc_in_pw
 		const double &beta);
 
     void produce_local_basis_in_pw(const int &ik,
-                                   ModulePW::PW_Basis_K *wfc_basis,
+                                   const ModulePW::PW_Basis_K *wfc_basis,
                                    ModuleBase::ComplexMatrix &psi,
                                    const ModuleBase::realArray &table_local);
 

source/module_io/berryphase.cpp

@@ -228,7 +228,13 @@ void berryphase::set_kpoints(const K_Vectors& kv, const int direction)
 }
 
 #include "../module_base/complexmatrix.h"
-double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv)
+double berryphase::stringPhase(int index_str,
+                               int nbands,
+                               const int npwx,
+                               const psi::Psi<std::complex<double>>* psi_in,
+                               const ModulePW::PW_Basis* rhopw,
+                               const ModulePW::PW_Basis_K* wfcpw,
+                               const K_Vectors& kv)
 {
 	std::complex<double> zeta(1.0, 0.0);
 	ModuleBase::ComplexMatrix mat(nbands,nbands);
@@ -271,13 +277,13 @@ double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::co
 								ModuleBase::Vector3<double> tem_G(0.0,0.0,1.0);
 								G = tem_G;
 							}
-							
-							mat(nb,mb) = pw_method.unkdotp_G0(ik_1, ik_2, nb, mb, psi_in, G);
-						}
+
+                            mat(nb, mb) = pw_method.unkdotp_G0(rhopw, wfcpw, ik_1, ik_2, nb, mb, psi_in, G);
+                        }
 						else 
 						{
-							mat(nb, mb) = pw_method.unkdotp_G(ik_1, ik_2, nb, mb, psi_in);
-						}
+                            mat(nb, mb) = pw_method.unkdotp_G(wfcpw, ik_1, ik_2, nb, mb, psi_in);
+                        }
 					}
 					else
 					{
@@ -298,11 +304,12 @@ double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::co
 								ModuleBase::Vector3<double> tem_G(0.0,0.0,1.0);
 								G = tem_G;
 							}
-
-							mat(nb,mb) = pw_method.unkdotp_soc_G0(ik_1, ik_2, nb, mb, psi_in, G);							
-						}
-						else  mat(nb, mb) = pw_method.unkdotp_soc_G(ik_1, ik_2, nb, mb, psi_in);
-					}
+
+                            mat(nb, mb) = pw_method.unkdotp_soc_G0(rhopw, wfcpw, ik_1, ik_2, nb, mb, psi_in, G);
+                        }
+                        else
+                            mat(nb, mb) = pw_method.unkdotp_soc_G(wfcpw, ik_1, ik_2, nb, mb, npwx, psi_in);
+                    }
 
 				} // nb
 
@@ -377,7 +384,14 @@ double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::co
 	return log(zeta).imag();
 }
 
-void berryphase::Berry_Phase(int nbands, double &pdl_elec_tot, int &mod_elec_tot, const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv)
+void berryphase::Berry_Phase(int nbands,
+                             double& pdl_elec_tot,
+                             int& mod_elec_tot,
+                             const int npwx,
+                             const psi::Psi<std::complex<double>>* psi_in,
+                             const ModulePW::PW_Basis* rhopw,
+                             const ModulePW::PW_Basis_K* wfcpw,
+                             const K_Vectors& kv)
 {		
 	std::complex<double> cave = 0.0;
 	double *phik = new double[total_string];
@@ -397,8 +411,8 @@ void berryphase::Berry_Phase(int nbands, double &pdl_elec_tot, int &mod_elec_tot
 
 	for(int istring = 0; istring < total_string; istring++)
 	{
-		phik[istring] = stringPhase(istring,nbands, psi_in, kv);
-		// transfer phase to complex number
+        phik[istring] = stringPhase(istring, nbands, npwx, psi_in, rhopw, wfcpw, kv);
+        // transfer phase to complex number
 		cphik[istring] = std::complex<double>(cos(phik[istring]),sin(phik[istring]));	
 		cave = cave + std::complex<double>(wistring[istring],0.0) * cphik[istring];
 
@@ -448,8 +462,11 @@ void berryphase::Berry_Phase(int nbands, double &pdl_elec_tot, int &mod_elec_tot
 
 }
 
-
-void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv)
+void berryphase::Macroscopic_polarization(const int npwx,
+                                          const psi::Psi<std::complex<double>>* psi_in,
+                                          const ModulePW::PW_Basis* rhopw,
+                                          const ModulePW::PW_Basis_K* wfcpw,
+                                          const K_Vectors& kv)
 {	
 	get_occupation_bands();
 #ifdef __LCAO	
@@ -577,9 +594,9 @@ void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>*
 			set_kpoints(kv, direction);
 			double pdl_elec_tot = 0.0;
 			int mod_elec_tot = 0;
-			Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, psi_in, kv);
-		
-			const double rmod = GlobalC::ucell.a1.norm() * GlobalC::ucell.lat0;
+            Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv);
+
+            const double rmod = GlobalC::ucell.a1.norm() * GlobalC::ucell.lat0;
 			const double unit1 = rmod;
 			const double unit2 = rmod / GlobalC::ucell.omega;
 			const double unit3 = ( rmod / GlobalC::ucell.omega ) * ( 1.60097e-19/pow(5.29177e-11,2) );
@@ -612,9 +629,9 @@ void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>*
 			set_kpoints(kv, direction);
 			double pdl_elec_tot = 0.0;
 			int mod_elec_tot = 0;
-			Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, psi_in, kv);
-		
-			const double rmod = GlobalC::ucell.a2.norm() * GlobalC::ucell.lat0;
+            Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv);
+
+            const double rmod = GlobalC::ucell.a2.norm() * GlobalC::ucell.lat0;
 			const double unit1 = rmod;
 			const double unit2 = rmod / GlobalC::ucell.omega;
 			const double unit3 = ( rmod / GlobalC::ucell.omega ) * ( 1.60097e-19/pow(5.29177e-11,2) );
@@ -647,9 +664,9 @@ void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>*
 			set_kpoints(kv, direction);
 			double pdl_elec_tot = 0.0;
 			int mod_elec_tot = 0;
-			Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, psi_in, kv);
-		
-			const double rmod = GlobalC::ucell.a3.norm() * GlobalC::ucell.lat0;
+            Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv);
+
+            const double rmod = GlobalC::ucell.a3.norm() * GlobalC::ucell.lat0;
 			const double unit1 = rmod;
 			const double unit2 = rmod / GlobalC::ucell.omega;
 			const double unit3 = ( rmod / GlobalC::ucell.omega ) * ( 1.60097e-19/pow(5.29177e-11,2) );

source/module_io/berryphase.h

@@ -4,6 +4,8 @@
 #ifdef __LCAO
 #include "unk_overlap_lcao.h"
 #endif
+#include "module_basis/module_pw/pw_basis.h"
+#include "module_basis/module_pw/pw_basis_k.h"
 #include "module_cell/klist.h"
 #include "module_psi/psi.h"
 
@@ -39,13 +41,30 @@ class berryphase
 
 	void set_kpoints(const K_Vectors& kv, const int direction);
 
-	double stringPhase(int index_str, int nbands, const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv);
+    double stringPhase(int index_str,
+                       int nbands,
+                       const int npwx,
+                       const psi::Psi<std::complex<double>>* psi_in,
+                       const ModulePW::PW_Basis* rhopw,
+                       const ModulePW::PW_Basis_K* wfcpw,
+                       const K_Vectors& kv);
 
-	void Berry_Phase(int nbands, double &pdl_elec_tot, int &mod_elec_tot, const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv);
+    void Berry_Phase(int nbands,
+                     double& pdl_elec_tot,
+                     int& mod_elec_tot,
+                     const int npwx,
+                     const psi::Psi<std::complex<double>>* psi_in,
+                     const ModulePW::PW_Basis* rhopw,
+                     const ModulePW::PW_Basis_K* wfcpw,
+                     const K_Vectors& kv);
 
-	void Macroscopic_polarization(const psi::Psi<std::complex<double>>* psi_in, const K_Vectors& kv);
+    void Macroscopic_polarization(const int npwx,
+                                  const psi::Psi<std::complex<double>>* psi_in,
+                                  const ModulePW::PW_Basis* rhopw,
+                                  const ModulePW::PW_Basis_K* wfcpw,
+                                  const K_Vectors& kv);
 
-	std::string outFormat(const double polarization, const double modulus, const ModuleBase::Vector3<double> project);
+    std::string outFormat(const double polarization, const double modulus, const ModuleBase::Vector3<double> project);
 
 };
 

source/module_io/cal_test.cpp

@@ -45,31 +45,31 @@ double Cal_Test::meigts123;
 
 double Cal_Test::mtot;
 
-void Cal_Test::test_memory(void)
+void Cal_Test::test_memory(const ModulePW::PW_Basis* rhopw, const ModulePW::PW_Basis_K* wfcpw, const std::string chr_mixing_mode, const int chr_mixing_ndim)
 {
 	ModuleBase::TITLE("Cal_Test","test_memory");
 
-	const int ngmw = Cal_Test::cal_np(GlobalC::wfcpw->ggecut, GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz);
-	const int ngmc = Cal_Test::cal_np(GlobalC::rhopw->ggecut, GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz);
+	const int ngmw = Cal_Test::cal_np(wfcpw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz);
+	const int ngmc = Cal_Test::cal_np(rhopw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz);
 
 	std::cout << " number of atoms = " << GlobalC::ucell.nat << std::endl;
 	std::cout << " plane wave number for wave functions = " << ngmw << std::endl;
 	std::cout << " plane wave number for chage density  = " << ngmc << std::endl;
 
-	mporter = ModuleBase::Memory::calculate_mem ( GlobalC::rhopw->nxyz, "double");
+	mporter = ModuleBase::Memory::calculate_mem ( rhopw->nxyz, "double");
 
 	mrho = mporter;
 	mrho_save = mrho;
 	mrho_core = mrho;
 
 	// (2) memory for charge mixing
-	std::cout << " Mixing mode = " << GlobalC::CHR_MIX.get_mixing_mode() << std::endl;
-	if(GlobalC::CHR_MIX.get_mixing_mode() == "pulay")
+	std::cout << " Mixing mode = " << chr_mixing_mode << std::endl;
+	if(chr_mixing_mode == "pulay")
 	{
-		std::cout << " Mixing dimension = " << GlobalC::CHR_MIX.get_mixing_ndim() << std::endl;
-		mRrho = GlobalC::CHR_MIX.get_mixing_ndim() * mrho;
-		mdRrho = (GlobalC::CHR_MIX.get_mixing_ndim()-1) * mrho;
-		mdrho = (GlobalC::CHR_MIX.get_mixing_ndim()-1) * mrho;
+		std::cout << " Mixing dimension = " << chr_mixing_ndim << std::endl;
+		mRrho = chr_mixing_ndim * mrho;
+		mdRrho = (chr_mixing_ndim-1) * mrho;
+		mdrho = (chr_mixing_ndim-1) * mrho;
 		mrho_save2 = mrho;
 //		std::cout << " Memory for pulay mixing: " << mrho << " MB" << std::endl;	
 	}
@@ -97,7 +97,7 @@ void Cal_Test::test_memory(void)
 	mgg = ModuleBase::Memory::calculate_mem( ngmc, "double");
 	mig123 = ModuleBase::Memory::calculate_mem( ngmc*3, "int");
 	mstrucFac = ModuleBase::Memory::calculate_mem( GlobalC::ucell.ntype*ngmc, "cdouble");
-	meigts123 = ModuleBase::Memory::calculate_mem( GlobalC::ucell.nat * (2*GlobalC::rhopw->nx+1+2*GlobalC::rhopw->ny+1+2*GlobalC::rhopw->nz+1), "cdouble");
+	meigts123 = ModuleBase::Memory::calculate_mem( GlobalC::ucell.nat * (2*rhopw->nx+1+2*rhopw->ny+1+2*rhopw->nz+1), "cdouble");
 
 //	std::cout << " Memory for "
 

source/module_io/cal_test.h

@@ -1,9 +1,14 @@
 #ifndef CAL_TEST_H
 #define CAL_TEST_H
+#include "module_basis/module_pw/pw_basis.h"
+#include "module_basis/module_pw/pw_basis_k.h"
 
 namespace Cal_Test
 {
-	void test_memory(void);
+	void test_memory(const ModulePW::PW_Basis* rhopw,
+		const ModulePW::PW_Basis_K* wfcpw,
+		const std::string chr_mixing_mode,
+		const int chr_mixing_ndim);
 	int cal_np(const double &ggcut, const int &n1, const int &n2, const int &n3);
 	void print_mem(const int &nproc);