Refactor LCAO step 2

source/Makefile.Objects

@@ -496,6 +496,10 @@ OBJS_LCAO=DM_gamma.o\
       FORCE_gamma_tvnl.o\
       FORCE_gamma_vl.o\
       FORCE_k.o\
+      foverlap_k.o\
+      ftvnl_dphi_k.o\
+      fvl_dphi_k.o\
+      fvnl_dbeta_k.o\
       LCAO_gen_fixedH.o\
       LCAO_hamilt.o\
       LCAO_matrix.o\

source/module_elecstate/elecstate_lcao.cpp

@@ -142,9 +142,9 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
     //------------------------------------------------------------
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    this->uhm->GK.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
+    this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
     Gint_inout inout(this->loc->DM_R, this->charge->rho, Gint_Tools::job_type::rho);
-    this->uhm->GK.cal_gint(&inout);
+    this->gint_k->cal_gint(&inout);
 
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
@@ -153,7 +153,7 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
             ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
         }
         Gint_inout inout1(this->loc->DM_R, this->charge->kin_r, Gint_Tools::job_type::tau);
-        this->uhm->GK.cal_gint(&inout1);
+        this->gint_k->cal_gint(&inout1);
     }
 
     this->charge->renormalize_rho();
@@ -216,17 +216,21 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
     // calculate the charge density on real space grid.
     //------------------------------------------------------------
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    this->uhm->GG.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
+
+    this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
+
     Gint_inout inout(this->loc->DM, this->charge->rho, Gint_Tools::job_type::rho);
-    this->uhm->GG.cal_gint(&inout);
+
+    this->gint_gamma->cal_gint(&inout);
+
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {
             ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
         }
         Gint_inout inout1(this->loc->DM, this->charge->kin_r, Gint_Tools::job_type::tau);
-        this->uhm->GG.cal_gint(&inout1);
+        this->gint_gamma->cal_gint(&inout1);
     }
 
     this->charge->renormalize_rho();

source/module_elecstate/elecstate_lcao.h

@@ -2,7 +2,8 @@
 #define ELECSTATELCAO_H
 
 #include "elecstate.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
+#include "module_hamilt_lcao/module_gint/gint_gamma.h"
+#include "module_hamilt_lcao/module_gint/gint_k.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h"
 #include "module_elecstate/module_dm/density_matrix.h"
@@ -18,22 +19,27 @@ class ElecStateLCAO : public ElecState
                   const K_Vectors* klist_in ,
                   int nks_in,
                   Local_Orbital_Charge* loc_in ,
+                  Gint_Gamma* gint_gamma_in,  //mohan add 2024-04-01
+                  Gint_k* gint_k_in, //mohan add 2024-04-01
                   LCAO_Hamilt* uhm_in ,
                   Local_Orbital_wfc* lowf_in ,
                   ModulePW::PW_Basis* rhopw_in ,
                   ModulePW::PW_Basis_Big* bigpw_in )
     {
         init_ks(chg_in, klist_in, nks_in, rhopw_in, bigpw_in);
         this->loc = loc_in;
-        this->uhm = uhm_in;
+        this->gint_gamma = gint_gamma_in; // mohan add 2024-04-01
+        this->gint_k = gint_k_in; // mohan add 2024-04-01
         this->lowf = lowf_in;
         this->classname = "ElecStateLCAO";
     }
 
     virtual ~ElecStateLCAO()
     {
-        if (this->DM != nullptr)
-            delete this->DM;
+		if (this->DM != nullptr)
+		{
+			delete this->DM;
+		}
     }
 
     // void init(Charge* chg_in):charge(chg_in){} override;
@@ -69,7 +75,9 @@ class ElecStateLCAO : public ElecState
     // void rhoBandK(const psi::Psi<std::complex<double>>& psi);
 
     Local_Orbital_Charge* loc = nullptr;
-    LCAO_Hamilt* uhm = nullptr;
+    Gint_Gamma* gint_gamma = nullptr; // mohan add 2024-04-01
+    Gint_k* gint_k = nullptr; // mohan add 2024-04-01
+    //LCAO_Hamilt* uhm = nullptr; // mohan modify 2024-04-01
     Local_Orbital_wfc* lowf = nullptr;
     DensityMatrix<TK,double>* DM = nullptr;
 
@@ -86,4 +94,4 @@ bool ElecStateLCAO<TK>::need_psi_grid = 1;
 
 } // namespace elecstate
 
-#endif
+#endif

source/module_elecstate/elecstate_lcao_tddft.cpp

@@ -63,9 +63,9 @@ void ElecStateLCAO_TDDFT::psiToRho_td(const psi::Psi<std::complex<double>>& psi)
     //------------------------------------------------------------
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    this->uhm->GK.transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
+    this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
     Gint_inout inout(this->loc->DM_R, this->charge->rho, Gint_Tools::job_type::rho); // rho calculation
-    this->uhm->GK.cal_gint(&inout);
+    this->gint_k->cal_gint(&inout);
 
     this->charge->renormalize_rho();
 
@@ -110,5 +110,4 @@ void ElecStateLCAO_TDDFT::calculate_weights_td()
 }
 
 
-
-} // namespace elecstate
+} // namespace elecstate

source/module_elecstate/elecstate_lcao_tddft.h

@@ -3,7 +3,7 @@
 
 #include "elecstate.h"
 #include "elecstate_lcao.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
+#include "module_hamilt_lcao/module_gint/gint_k.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h"
 
@@ -16,14 +16,14 @@ class ElecStateLCAO_TDDFT : public ElecStateLCAO<std::complex<double>>
                         const K_Vectors* klist_in ,
                         int nks_in,
                         Local_Orbital_Charge* loc_in ,
-                        LCAO_Hamilt* uhm_in ,
+                        Gint_k* gint_k_in, //mohan add 2024-04-01
                         Local_Orbital_wfc* lowf_in ,
                         ModulePW::PW_Basis* rhopw_in ,
                         ModulePW::PW_Basis_Big* bigpw_in )
     {
         init_ks(chg_in, klist_in, nks_in, rhopw_in, bigpw_in);        
         this->loc = loc_in;
-        this->uhm = uhm_in;
+        this->gint_k = gint_k_in;
         this->lowf = lowf_in;
         this->classname = "ElecStateLCAO_TDDFT";
     }
@@ -33,4 +33,4 @@ class ElecStateLCAO_TDDFT : public ElecStateLCAO<std::complex<double>>
 
 } // namespace elecstate
 
-#endif
+#endif

source/module_esolver/esolver_ks_lcao.cpp

@@ -111,10 +111,13 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     // autoset nbands in ElecState, it should before basis_init (for Psi 2d divid)
     if (this->pelec == nullptr)
     {
-        this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr),
+        this->pelec = new elecstate::ElecStateLCAO<TK>(
+            &(this->chr),
             &(this->kv),
             this->kv.nks,
             &(this->LOC),
+            &(this->GG), // mohan add 2024-04-01
+            &(this->GK), // mohan add 2024-04-01
             &(this->UHM),
             &(this->LOWF),
             this->pw_rho,
@@ -255,6 +258,8 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
 				&(this->kv),
 				this->kv.nks,
 				&(this->LOC),
+                &(this->GG), // mohan add 2024-04-01
+                &(this->GK), // mohan add 2024-04-01
 				&(this->UHM),
 				&(this->LOWF),
 				this->pw_rho,
@@ -313,6 +318,8 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 			this->pelec,
 			this->psi,
 			this->UHM,
+            this->GG, // mohan add 2024-04-01
+            this->GK, // mohan add 2024-04-01
 			force,
 			this->scs,
 			this->sf,
@@ -665,7 +672,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
         // update Gint_K
         if (!GlobalV::GAMMA_ONLY_LOCAL)
         {
-            this->UHM.GK.renew();
+            this->GK.renew();
         }
         // update real space Hamiltonian
         this->p_hamilt->refresh();
@@ -803,7 +810,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
     {
         if (!GlobalV::GAMMA_ONLY_LOCAL && hsolver::HSolverLCAO<TK>::out_mat_hs[0])
         {
-            this->UHM.GK.renew(true);
+            this->GK.renew(true);
         }
         for (int ik = 0; ik < this->kv.nks; ++ik)
         {
@@ -1016,9 +1023,24 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     bool out_exc = true;    // tmp, add parameter!
     if (GlobalV::out_mat_xc)
     {
-        ModuleIO::write_Vxc<TK, TR>(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK,
-            *this->psi, GlobalC::ucell, this->sf, *this->pw_rho, *this->pw_rhod, GlobalC::ppcell.vloc,
-            *this->pelec->charge, this->UHM, this->LM, this->LOC, this->kv, this->pelec->wg, GlobalC::GridD);
+		ModuleIO::write_Vxc<TK, TR>(
+				GlobalV::NSPIN, 
+				GlobalV::NLOCAL, 
+				GlobalV::DRANK,
+				*this->psi, 
+				GlobalC::ucell, 
+				this->sf, 
+				*this->pw_rho, 
+				*this->pw_rhod, 
+				GlobalC::ppcell.vloc,
+				*this->pelec->charge, 
+				this->GG,
+				this->GK,
+				this->LM, 
+				this->LOC, 
+				this->kv, 
+				this->pelec->wg, 
+				GlobalC::GridD);
     }
 
 #ifdef __EXX
@@ -1181,14 +1203,16 @@ ModuleIO::Output_DM1 ESolver_KS_LCAO<TK, TR>::create_Output_DM1(int istep)
 template <typename TK, typename TR>
 ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Sparse(int istep)
 {
-	return ModuleIO::Output_Mat_Sparse<TK>(hsolver::HSolverLCAO<TK>::out_mat_hsR,
+	return ModuleIO::Output_Mat_Sparse<TK>(
+            hsolver::HSolverLCAO<TK>::out_mat_hsR,
 			hsolver::HSolverLCAO<TK>::out_mat_dh,
 			hsolver::HSolverLCAO<TK>::out_mat_t,
 			INPUT.out_mat_r,
 			istep,
 			this->pelec->pot->get_effective_v(),
 			*this->LOWF.ParaV,
 			this->UHM,
+            this->GK, // mohan add 2024-04-01
 			this->LM,
 			this->kv,
 			this->p_hamilt);

source/module_esolver/esolver_ks_lcao.h

@@ -6,6 +6,10 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
+// for grid integration
+#include "module_hamilt_lcao/module_gint/gint_gamma.h"
+#include "module_hamilt_lcao/module_gint/gint_k.h"
+
 #include "module_basis/module_ao/ORB_control.h"
 #ifdef __EXX
 #include "module_ri/Mix_DMk_2D.h"
@@ -75,6 +79,12 @@ namespace ModuleESolver
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
         LCAO_Hamilt UHM;
 
+        // used for k-dependent grid integration.
+        Gint_k GK;
+
+		// used for gamma only algorithms.
+		Gint_Gamma GG;
+
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
         LCAO_Matrix LM;
 

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -126,7 +126,12 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     }
 
     // prepare grid in Gint
-    this->UHM.grid_prepare(this->GridT, *this->pw_rho, *this->pw_big);
+	this->UHM.grid_prepare(
+			this->GridT, 
+			this->GG,
+			this->GK,
+			*this->pw_rho, 
+			*this->pw_big);
 
     // init Hamiltonian
     if (this->p_hamilt != nullptr)
@@ -137,8 +142,8 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     if (this->p_hamilt == nullptr)
     {
         elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
-        this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(GlobalV::GAMMA_ONLY_LOCAL ? &(this->UHM.GG) : nullptr,
-            GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->UHM.GK),
+        this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr,
+            GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK),
             &(this->UHM.genH),
             &(this->LM),
             &(this->LOC),
@@ -194,29 +199,29 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
         if (GlobalV::GAMMA_ONLY_LOCAL)
         {
             Gint_inout inout(this->LOC.DM, this->pelec->charge->rho, Gint_Tools::job_type::rho);
-            this->UHM.GG.cal_gint(&inout);
+            this->GG.cal_gint(&inout);
             if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
             {
                 for (int is = 0; is < GlobalV::NSPIN; is++)
                 {
                     ModuleBase::GlobalFunc::ZEROS(this->pelec->charge->kin_r[0], this->pw_rho->nrxx);
                 }
                 Gint_inout inout1(this->LOC.DM, this->pelec->charge->kin_r, Gint_Tools::job_type::tau);
-                this->UHM.GG.cal_gint(&inout1);
+                this->GG.cal_gint(&inout1);
             }
         }
         else
         {
             Gint_inout inout(this->LOC.DM_R, this->pelec->charge->rho, Gint_Tools::job_type::rho);
-            this->UHM.GK.cal_gint(&inout);
+            this->GK.cal_gint(&inout);
             if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
             {
                 for (int is = 0; is < GlobalV::NSPIN; is++)
                 {
                     ModuleBase::GlobalFunc::ZEROS(this->pelec->charge->kin_r[0], this->pw_rho->nrxx);
                 }
                 Gint_inout inout1(this->LOC.DM_R, this->pelec->charge->kin_r, Gint_Tools::job_type::tau);
-                this->UHM.GK.cal_gint(&inout1);
+                this->GK.cal_gint(&inout1);
             }
         }
 
@@ -389,7 +394,7 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
     else if (GlobalV::CALCULATION == "get_pchg")
     {
         IState_Charge ISC(this->psi, this->LOC);
-        ISC.begin(this->UHM.GG,
+        ISC.begin(this->GG,
                   this->pelec,
                   this->pw_rho,
                   this->pw_big,
@@ -411,7 +416,7 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
                       this->pw_wfc,
                       this->pw_big,
                       this->LOWF,
-                      this->UHM.GG,
+                      this->GG,
                       INPUT.out_wfc_pw,
                       this->wf.out_wfc_r,
                       this->kv,
@@ -427,7 +432,7 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
                       this->pw_wfc,
                       this->pw_big,
                       this->LOWF,
-                      this->UHM.GK,
+                      this->GK,
                       INPUT.out_wfc_pw,
                       this->wf.out_wfc_r,
                       this->kv,

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -79,8 +79,8 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
         this->pelec = new elecstate::ElecStateLCAO_TDDFT(&(this->chr),
                                                          &(kv),
                                                          kv.nks,
-                                                         &(this->LOC),
-                                                         &(this->UHM),
+														 &(this->LOC),
+														 &(this->GK), // mohan add 2024-04-01
                                                          &(this->LOWF),
                                                          this->pw_rho,
                                                          pw_big);
@@ -264,7 +264,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
     {
         if (!GlobalV::GAMMA_ONLY_LOCAL)
         {
-            this->UHM.GK.renew(true);
+            this->GK.renew(true);
         }
         for (int ik = 0; ik < kv.nks; ++ik)
         {