Refactor LCAO code Step 3 (#3914)

* update timer, add std

* update the format of esolver_fp.cpp

* add description in esolver_fp.h

* update a few esolver files

* update esolver description for LCAO

* update some formats of esolver codes

* update esolver_ks_pw.cpp formats

* update formats of esolver_ks_pw.cpp and esolver_ks_pw.h

* update esolver_ks_lcao_tddft.cpp formats

* update format of esolver_sdft_pw.cpp

* update the format of esolver_sdft_pw_tool.cpp

* keep formating esolver_ks_pw.cpp

* formating esolver_of_interface.cpp

* change GlobalC::ucell to ucell in esolver_ks.cpp

* remove some GlobalC::ucell in esolver_sdft_pw.cpp

* refactor the code before getting rid of RA in esolver_lcao

* refactor before getting rid of LOWF

* refactor before getting rid of LCAO_hamilt.h and LCAO_matrix.h

* refactor wavefunc_in_pw

* refactor density matrix

* refactor the format cal_dm_psi.cpp

* format forces.cpp

* refactor esolver_of_tool.cpp

* change member function beforescf in Esolver to before_scf

* change afterscf to after_scf

* change updatepot to update_pot

* change eachiterinit to iter_init, change eachiterfinish to iter_finish

* refactor esolvers, change member function names of most esolvers

* reformat esolver.h

* update tests for esolvers

* add TITLE in esolver_ks_lcao

* update esolver_ks_lcao

* update esolver_lcao

* update timer::tick in esolver_lcao

* try to delete LCAO_Matrix in LCAO_Hamilt, and try to delete Parallel_Orbitals in Force_k

* fix the compiling issue with LCAO_hamilt.hpp

* try to divide the FORCE_k.cpp into several small files

* divide FORCE_k into foverlap_k.cpp ftvnl_dphi_k.cpp fvl_dphi_k.cpp fvnl_dbeta_k.cpp four files

* get rid of UHM in FORCE_k.cpp

* cannot compile, but I have modified some files in order to get rid of Gint_k and Gint_Gamma in UHM

* keep updating, cannot run

* update write_Vxc, cannot run yet

* keep updating gint_gamma and gint_k

* update LCAO_matrix.cpp

* divide force files, and update Makefile.Objects

* update LCAO_hamilt.cpp

* delete genH pointer in UHM

* divide LCAO_hamilt.cpp into small codes, grid_init.cpp is the first one

* update grid_init in esolver_ks_lcao

* add a new namespace named sparse_format, most of the functions that originally belong to LCAO_hamilt should be moved to sparse_format

* cannot find the mismatch of DFTU

* fix the DFTU error

* update

source/Makefile.Objects

@@ -502,6 +502,8 @@ OBJS_LCAO=DM_gamma.o\
       fvnl_dbeta_k.o\
       LCAO_gen_fixedH.o\
       LCAO_hamilt.o\
+      grid_init.o\
+      sparse_format.o\
       LCAO_matrix.o\
       LCAO_nnr.o\
       center2_orb-orb11.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -117,7 +117,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
             &(this->LOC),
             &(this->GG), // mohan add 2024-04-01
             &(this->GK), // mohan add 2024-04-01
-            &(this->UHM),
+            &(this->uhm),
             &(this->LOWF),
             this->pw_rho,
             this->pw_big);
@@ -131,7 +131,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     //------------------init Basis_lcao----------------------
 
     //! pass Hamilt-pointer to Operator
-    this->UHM.genH.LM = this->UHM.LM = &this->LM;
+    this->gen_h.LM = this->uhm.LM = &this->LM;
 
     //! pass basis-pointer to EState and Psi
     this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV = &(this->orb_con.ParaV);
@@ -259,7 +259,7 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
 				&(this->LOC),
                 &(this->GG), // mohan add 2024-04-01
                 &(this->GK), // mohan add 2024-04-01
-				&(this->UHM),
+				&(this->uhm),
 				&(this->LOWF),
 				this->pw_rho,
 				this->pw_big);
@@ -313,10 +313,10 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 			GlobalV::TEST_STRESS,
 			this->LOC,
             this->orb_con.ParaV, 
-            this->LM,
 			this->pelec,
 			this->psi,
-			this->UHM,
+            this->LM,
+            this->gen_h, // mohan add 2024-04-02
             this->GG, // mohan add 2024-04-01
             this->GK, // mohan add 2024-04-01
 			force,
@@ -426,7 +426,7 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
         ModuleIO::write_proj_band_lcao(
             this->psi,
-            this->UHM,
+            this->uhm,
             this->pelec,
             this->kv,
             GlobalC::ucell,
@@ -437,7 +437,7 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
 		ModuleIO::out_dos_nao(
 				this->psi,
-				this->UHM,
+				this->uhm,
 				this->pelec->ekb,
 				this->pelec->wg,
 				INPUT.dos_edelta_ev,
@@ -1204,7 +1204,8 @@ ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Spars
 			istep,
 			this->pelec->pot->get_effective_v(),
 			*this->LOWF.ParaV,
-			this->UHM,
+			this->uhm,
+            this->gen_h, // mohan add 2024-04-02
             this->GK, // mohan add 2024-04-01
 			this->LM,
 			this->kv,

source/module_esolver/esolver_ks_lcao.h

@@ -77,7 +77,9 @@ namespace ModuleESolver
         Local_Orbital_Charge LOC;
 
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
-        LCAO_Hamilt UHM;
+        LCAO_Hamilt uhm;
+
+        LCAO_gen_fixedH gen_h; // mohan add 2024-04-02
 
         // used for k-dependent grid integration.
         Gint_k GK;

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -23,6 +23,7 @@
 #include "module_io/dm_io.h"
 
 #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
 
 namespace ModuleESolver
 {
@@ -69,7 +70,7 @@ void ESolver_KS_LCAO<TK, TR>::set_matrix_grid(Record_adj& ra)
     // (2)For each atom, calculate the adjacent atoms in different cells
     // and allocate the space for H(R) and S(R).
     // If k point is used here, allocate HlocR after atom_arrange.
-    Parallel_Orbitals* pv = this->UHM.LM->ParaV;
+    Parallel_Orbitals* pv = this->LM.ParaV;
     ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL);
     if (!GlobalV::GAMMA_ONLY_LOCAL)
     {
@@ -126,7 +127,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     }
 
     // prepare grid in Gint
-	this->UHM.grid_prepare(
+	LCAO_domain::grid_prepare(
 			this->GridT, 
 			this->GG,
 			this->GK,
@@ -144,7 +145,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
         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->GG) : nullptr,
             GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK),
-            &(this->UHM.genH),
+            &(this->gen_h),
             &(this->LM),
             &(this->LOC),
             this->pelec->pot,
@@ -234,7 +235,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     // since it depends on ionic positions
     if (GlobalV::deepks_setorb)
     {
-        const Parallel_Orbitals* pv = this->UHM.LM->ParaV;
+        const Parallel_Orbitals* pv = this->LM.ParaV;
         // build and save <psi(0)|alpha(R)> at beginning
         GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD, GlobalC::UOT);
 
@@ -476,14 +477,15 @@ void ESolver_KS_LCAO<std::complex<double>, double>::get_S()
                          GlobalV::test_atom_input);
 
     this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL);
-    this->UHM.genH.LM->ParaV = &this->orb_con.ParaV;
+    this->LM.ParaV = &this->orb_con.ParaV;
     if (this->p_hamilt == nullptr)
     {
-        this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(this->UHM.genH.LM, this->kv);
+        this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(&this->LM, this->kv);
         dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(this->p_hamilt->ops)->contributeHR();
     }
-    ModuleIO::output_S_R(this->UHM, this->p_hamilt, "SR.csr");
+    ModuleIO::output_S_R(this->uhm, this->p_hamilt, "SR.csr");
 }
+
 template <>
 void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S()
 {
@@ -503,15 +505,15 @@ void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S()
                          GlobalV::test_atom_input);
 
     this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL);
-    this->UHM.genH.LM->ParaV = &this->orb_con.ParaV;
+    this->LM.ParaV = &this->orb_con.ParaV;
     if (this->p_hamilt == nullptr)
     {
         this->p_hamilt
-            = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(this->UHM.genH.LM, this->kv);
+            = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(&this->LM, this->kv);
         dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>*>(this->p_hamilt->ops)
             ->contributeHR();
     }
-    ModuleIO::output_S_R(this->UHM, this->p_hamilt, "SR.csr");
+    ModuleIO::output_S_R(this->uhm, this->p_hamilt, "SR.csr");
 }
 
 template <typename TK, typename TR>

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -100,7 +100,7 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
     //------------------init Hamilt_lcao----------------------
 
     // pass Hamilt-pointer to Operator
-    this->UHM.genH.LM = this->UHM.LM = &this->LM;
+    this->gen_h.LM = this->uhm.LM = &this->LM;
     // pass basis-pointer to EState and Psi
     this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV;
 
@@ -455,8 +455,9 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
                         kv,
                         tmp_DM->get_paraV_pointer(),
                         this->RA,
-                        this->UHM);
-    }
+						this->LM, // mohan add 2024-04-02
+						this->gen_h); // mohan add 2024-02
+	}
     ESolver_KS_LCAO<std::complex<double>, double>::after_scf(istep);
 }
 

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -27,6 +27,8 @@ if(ENABLE_LCAO)
         fvnl_dbeta_k.cpp
         LCAO_gen_fixedH.cpp
         LCAO_hamilt.cpp
+        grid_init.cpp
+        sparse_format.cpp
         LCAO_matrix.cpp
         LCAO_nnr.cpp
         record_adj.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -31,10 +31,10 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const bool istests,
 										  Local_Orbital_Charge &loc,
 										  Parallel_Orbitals &pv,
-										  LCAO_Matrix &lm,
                                           const elecstate::ElecState* pelec,
                                           const psi::Psi<T>* psi,
-										  LCAO_Hamilt& uhm,
+										  LCAO_Matrix& lm,
+                                          LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 										  Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 										  Gint_k &gint_k, // mohan add 2024-04-01
 										  ModuleBase::matrix& fcs,
@@ -159,7 +159,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 #else
                                      svl_dphi,
 #endif
-									 uhm,
+                                     gen_h, // mohan add 2024-04-02
 									 gint_gamma,
 									 gint_k,
 									 pv,
@@ -236,18 +236,19 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         }
         if(GlobalV::dft_plus_u == 2)
         {
-            GlobalC::dftu.force_stress(pelec, *uhm.LM, force_dftu, stress_dftu, kv);
+            GlobalC::dftu.force_stress(pelec, lm, force_dftu, stress_dftu, kv);
         }
         else
         {
-            hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(uhm.LM,
-                                                                 kv.kvec_d,
-                                                                 nullptr,
-                                                                 nullptr,
-                                                                 GlobalC::ucell,
-                                                                 &GlobalC::GridD,
-                                                                 &GlobalC::dftu,
-                                                                 *(uhm.LM->ParaV));
+			hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(
+					&lm,
+					kv.kvec_d,
+					nullptr,
+					nullptr,
+					GlobalC::ucell,
+					&GlobalC::GridD,
+					&GlobalC::dftu,
+					*(lm.ParaV));
             tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu);
         }
     }
@@ -743,7 +744,7 @@ void Force_Stress_LCAO<double>::integral_part(
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-	LCAO_Hamilt &uhm,
+    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 	Gint_k &gint_k, // mohan add 2024-04-01
 	Parallel_Orbitals &pv,
@@ -769,7 +770,7 @@ void Force_Stress_LCAO<double>::integral_part(
 #else
         svl_dphi,
 #endif
-        uhm,
+        gen_h,
         gint_gamma,
         lm);
     return;
@@ -797,7 +798,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-	LCAO_Hamilt &uhm,
+    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma,
 	Gint_k &gint_k,
 	Parallel_Orbitals &pv,
@@ -823,7 +824,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #else
                      svl_dphi,
 #endif
-					 uhm,
+					 gen_h,
                      gint_k,
 					 pv,
 					 lm,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h

@@ -34,10 +34,10 @@ class Force_Stress_LCAO
         const bool istests,
 		Local_Orbital_Charge& loc,
 		Parallel_Orbitals &pv,
-		LCAO_Matrix &lm,
 		const elecstate::ElecState* pelec,
         const psi::Psi<T>* psi,
-		LCAO_Hamilt& uhm,
+		LCAO_Matrix &lm,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 		Gint_k &gint_k, // mohan add 2024-04-01
         ModuleBase::matrix& fcs,
@@ -92,7 +92,7 @@ class Force_Stress_LCAO
 #else
         ModuleBase::matrix& svl_dphi,
 #endif
-		LCAO_Hamilt &uhm,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma,
 		Gint_k &gint_k,
 	    Parallel_Orbitals &pv,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -38,7 +38,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 #else
                                     ModuleBase::matrix& svl_dphi,
 #endif
-									LCAO_Hamilt &uhm,
+		                            LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
                                     Gint_Gamma &gint_gamma,
 									LCAO_Matrix &lm)
 {
@@ -54,7 +54,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 
     // allocate DSloc_x, DSloc_y, DSloc_z
     // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
-    this->allocate_gamma(*this->ParaV, uhm.genH, lm);
+    this->allocate_gamma(*this->ParaV, gen_h, lm);
 
     // calculate the 'energy density matrix' here.
     this->cal_foverlap(isforce, isstress, psid, pelec, lm, foverlap, soverlap);

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h

@@ -45,7 +45,7 @@ class Force_LCAO_gamma
 #else
                       ModuleBase::matrix& svl_dphi,
 #endif
-					  LCAO_Hamilt& uhm,
+                      LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 					  Gint_Gamma &gint_gamma,
 					  LCAO_Matrix& lm);
 

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -50,7 +50,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 #else
                             ModuleBase::matrix& svl_dphi,
 #endif
-                            LCAO_Hamilt &uhm,
+                            LCAO_gen_fixedH &gen_h,
                             Gint_k &gint_k,
                             Parallel_Orbitals &pv,
                             LCAO_Matrix &lm,
@@ -65,7 +65,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 	this->allocate_k(
 			pv, 
 			lm, 
-            uhm.genH,
+            gen_h,
 			kv.nks, 
 			kv.kvec_d);
 
@@ -148,7 +148,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
             GlobalC::ld.check_v_delta_k(pv->nnr);
             for (int ik = 0; ik < kv.nks; ik++)
             {
-                uhm.LM->folding_fixedH(ik, kv.kvec_d);
+                LM->folding_fixedH(ik, kv.kvec_d);
             }
             GlobalC::ld.cal_e_delta_band_k(dm_k, kv.nks);
             std::ofstream ofs("E_delta_bands.dat");

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h

@@ -42,7 +42,7 @@ class Force_LCAO_k : public Force_LCAO_gamma
 #else
                   ModuleBase::matrix& svl_dphi,
 #endif
-				  LCAO_Hamilt &uhm,
+                  LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
                   Gint_k &gint_k,
 				  Parallel_Orbitals &pv,
 				  LCAO_Matrix &lm,

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h

@@ -0,0 +1,24 @@
+#ifndef LCAO_DOMAIN_H
+#define LCAO_DOMAIN_H
+
+#include "module_cell/module_neighbor/sltk_atom_arrange.h"
+#include "module_cell/module_neighbor/sltk_grid_driver.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
+
+
+namespace LCAO_domain
+{
+
+//! prepare grid integration 
+void grid_prepare(
+		const Grid_Technique& gt, 
+        Gint_Gamma &gint_gamma,
+        Gint_k &gint_k,
+		const ModulePW::PW_Basis& rhopw, 
+		const ModulePW::PW_Basis_Big& bigpw);
+
+
+}
+
+#endif