Update LCAO_matrix and related functions and delete LCAO_gen_fixedH (#…

…4400)

* add set_force and set_stress in LCAO_domain.h

* remove set_force and set_stress in LCAO_matrix

* add new files LCAO_nl_beta and LCAO_nl_mu

* move build_ST_new to LCAO_set_st.cpp

* delete LCAO_gen_fixedH

* delete useless print_Hk and print_Hgamma functions in LCAO_matrix.h and .cpp

* refactor the forces in LCAO

* remove GlobalC in fedm

* update force-related files

* fix a small bug in Makefile

source/Makefile.Objects

@@ -499,23 +499,27 @@ OBJS_LCAO=DM_gamma.o\
       td_velocity.o\
       upsi.o\
       FORCE_STRESS.o\
-      FORCE_gamma.o\
-      FORCE_gamma_edm.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\
-      grid_init.o\
+	  FORCE_gamma.o\
+	  FORCE_k.o\
+	  fvl_dphi_gamma.o\
+	  fvl_dphi_k.o\
+	  fedm_gamma.o\
+	  fedm_k.o\
+	  ftvnl_dphi_gamma.o\
+	  ftvnl_dphi_k.o\
+	  fvnl_dbeta_gamma.o\
+	  fvnl_dbeta_k.o\
+	  grid_init.o\
       spar_dh.o\
       spar_exx.o\
       spar_hsr.o\
       spar_st.o\
       spar_u.o\
       LCAO_matrix.o\
+      LCAO_set_fs.o\
+      LCAO_set_st.o\
+      LCAO_nl_mu.o\
+      LCAO_nl_beta.o\
       LCAO_nnr.o\
       center2_orb-orb11.o\
       center2_orb-orb21.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -163,9 +163,6 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(Input& inp, UnitCell& ucell)
     this->init_basis_lcao(this->orb_con, inp, ucell);
     //------------------init Basis_lcao----------------------
 
-    // 4) redundant ParaV and LM pointers
-    this->gen_h.LM = &this->LM;
-
     //! pass basis-pointer to EState and Psi
     /*
     Inform: on getting rid of ORB_control and Parallel_Orbitals
@@ -381,7 +378,6 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 			this->pelec,
 			this->psi,
             this->LM,
-            this->gen_h, // mohan add 2024-04-02
             this->GG, // mohan add 2024-04-01
             this->GK, // mohan add 2024-04-01
             uot_,
@@ -1433,7 +1429,6 @@ ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Spars
 			istep,
 			this->pelec->pot->get_effective_v(),
 			this->orb_con.ParaV,
-            this->gen_h, // mohan add 2024-04-06
             this->GK, // mohan add 2024-04-01
             uot_,
 			this->LM,

source/module_esolver/esolver_ks_lcao.h

@@ -75,8 +75,6 @@ namespace ModuleESolver
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
         Local_Orbital_Charge LOC;
 
-        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

@@ -145,7 +145,6 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr,
             GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK),
-            &(this->gen_h),
             &(this->LM),
             &(this->LOC),
             this->pelec->pot,

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -91,7 +91,6 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell)
 
     // this part will be updated soon 
     // pass Hamilt-pointer to Operator
-    this->gen_h.LM = &this->LM;
     this->LOC.ParaV = this->LM.ParaV;;
     this->LOWF.ParaV = this->LM.ParaV;
 
@@ -435,8 +434,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
                                 uot_,
                                 tmp_DM->get_paraV_pointer(),
                                 this->RA,
-                                this->LM,     // mohan add 2024-04-02
-                                this->gen_h); // mohan add 2024-02
+                                this->LM);     // mohan add 2024-04-02
     }
     ESolver_KS_LCAO<std::complex<double>, double>::after_scf(istep);
 }

source/module_esolver/io_npz.cpp

@@ -117,7 +117,7 @@ void ESolver_KS_LCAO<TK, TR>::read_mat_npz(std::string& zipname, hamilt::HContai
                 assert(orbital_info[iw*3] == GlobalC::ucell.atoms[it].iw2n[iw]);
                 assert(orbital_info[iw*3+1] == GlobalC::ucell.atoms[it].iw2l[iw]);
                 const int im = GlobalC::ucell.atoms[it].iw2m[iw];
-                const int m = (im % 2 == 0) ? -im/2 : (im+1)/2; // copied from LCAO_gen_fixedH.cpp
+                const int m = (im % 2 == 0) ? -im/2 : (im+1)/2; 
                 assert(orbital_info[iw*3+2] == m);
             }
         }
@@ -413,7 +413,7 @@ void ESolver_KS_LCAO<TK, TR>::output_mat_npz(std::string& zipname, const hamilt:
                 orbital_info[iw*3] = GlobalC::ucell.atoms[it].iw2n[iw];
                 orbital_info[iw*3+1] = GlobalC::ucell.atoms[it].iw2l[iw];
                 const int im = GlobalC::ucell.atoms[it].iw2m[iw];
-                const int m = (im % 2 == 0) ? -im/2 : (im+1)/2; // copied from LCAO_gen_fixedH.cpp
+                const int m = (im % 2 == 0) ? -im/2 : (im+1)/2; 
                 orbital_info[iw*3+2] = m;
             }
             shape={(size_t)GlobalC::ucell.atoms[it].nw,3};

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -17,15 +17,15 @@ if(ENABLE_LCAO)
         operator_lcao/dftu_lcao.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
-        FORCE_gamma_edm.cpp
-        FORCE_gamma_tvnl.cpp
-        FORCE_gamma_vl.cpp
         FORCE_k.cpp
-        foverlap_k.cpp
-        ftvnl_dphi_k.cpp
+        fvl_dphi_gamma.cpp
         fvl_dphi_k.cpp
+        fedm_gamma.cpp
+        fedm_k.cpp
+        ftvnl_dphi_gamma.cpp
+        ftvnl_dphi_k.cpp
+        fvnl_dbeta_gamma.cpp
         fvnl_dbeta_k.cpp
-        LCAO_gen_fixedH.cpp
         grid_init.cpp
         spar_dh.cpp
         spar_exx.cpp
@@ -34,7 +34,11 @@ if(ENABLE_LCAO)
         spar_u.cpp
         LCAO_matrix.cpp
         LCAO_nnr.cpp
-        record_adj.cpp
+		LCAO_set_fs.cpp
+		LCAO_set_st.cpp
+		LCAO_nl_mu.cpp
+		LCAO_nl_beta.cpp
+		record_adj.cpp
         center2_orb-orb11.cpp
         center2_orb-orb21.cpp
         center2_orb-orb22.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h

@@ -6,7 +6,6 @@
 #include "module_base/matrix.h"
 #include "module_elecstate/module_dm/density_matrix.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_psi/psi.h"
 #include "module_hamilt_lcao/module_gint/gint_gamma.h"
@@ -44,8 +43,10 @@ class Force_LCAO
     elecstate::Potential* pot;
 
     // orthonormal force + contribution from T and VNL
-    void ftable(const bool isforce,
+    void ftable(
+        const bool isforce,
         const bool isstress,
+        const UnitCell& ucell,
         const psi::Psi<T>* psi,
         const elecstate::ElecState* pelec,
         ModuleBase::matrix& foverlap,
@@ -59,7 +60,6 @@ class Force_LCAO
 #ifdef __DEEPKS
         ModuleBase::matrix& svnl_dalpha,
 #endif
-        LCAO_gen_fixedH& gen_h, // mohan add 2024-04-02
         typename TGint<T>::type& gint,
         const ORB_gen_tables* uot,
         const Parallel_Orbitals& pv,
@@ -71,7 +71,6 @@ class Force_LCAO
     // get the ds, dt, dvnl.
     void allocate(const Parallel_Orbitals& pv,
         LCAO_Matrix& lm,
-        LCAO_gen_fixedH& gen_h,
         const ORB_gen_tables* uot,
         const int& nks = 0,
         const std::vector<ModuleBase::Vector3<double>>& kvec_d = {});
@@ -88,22 +87,24 @@ class Force_LCAO
     // forces related to energy density matrix
     //-------------------------------------------------------------
 
-    void cal_foverlap(const bool isforce,
+    void cal_fedm(
+        const bool isforce,
         const bool isstress,
+        const UnitCell& ucell,
+        const elecstate::DensityMatrix<T, double>* dm,
         const psi::Psi<T>* psi,
         const Parallel_Orbitals& pv,
         const elecstate::ElecState* pelec,
         LCAO_Matrix& lm,
         ModuleBase::matrix& foverlap,
         ModuleBase::matrix& soverlap,
         const K_Vectors* kv = nullptr,
-        Record_adj* ra = nullptr,
-        const elecstate::DensityMatrix<T, double>* DM = nullptr);
+        Record_adj* ra = nullptr);
 
     //-------------------------------------------------------------
     // forces related to kinetic and non-local pseudopotentials
     //--------------------------------------------------------------
-    void cal_ftvnl_dphi(const elecstate::DensityMatrix<T, double>* DM,
+    void cal_ftvnl_dphi(const elecstate::DensityMatrix<T, double>* dm,
         const Parallel_Orbitals& pv,
         const UnitCell& ucell,
         LCAO_Matrix& lm,
@@ -113,7 +114,7 @@ class Force_LCAO
         ModuleBase::matrix& stvnl_dphi,
         Record_adj* ra = nullptr);
 
-    void cal_fvnl_dbeta(const elecstate::DensityMatrix<T, double>* DM,
+    void cal_fvnl_dbeta(const elecstate::DensityMatrix<T, double>* dm,
         const Parallel_Orbitals& pv,
         const UnitCell& ucell,
         const LCAO_Orbitals& orb,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -33,7 +33,6 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const elecstate::ElecState* pelec,
                                           const psi::Psi<T>* psi,
 										  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
                                           const ORB_gen_tables* uot,
@@ -156,7 +155,6 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 #ifdef __DEEPKS
         svnl_dalpha,
 #endif
-        gen_h, // mohan add 2024-04-02
         gint_gamma,
         gint_k,
         uot,
@@ -744,7 +742,6 @@ void Force_Stress_LCAO<double>::integral_part(
 #if __DEEPKS
     ModuleBase::matrix& svnl_dalpha,
 #endif
-    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
     const ORB_gen_tables* uot,
@@ -755,6 +752,7 @@ void Force_Stress_LCAO<double>::integral_part(
 
     flk.ftable(isforce,
         isstress,
+        GlobalC::ucell,
         psi,
         pelec,
         foverlap,
@@ -768,7 +766,6 @@ void Force_Stress_LCAO<double>::integral_part(
 #if __DEEPKS
         svnl_dalpha,
 #endif
-        gen_h,
         gint_gamma,
         uot,
         pv,
@@ -795,7 +792,6 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #if __DEEPKS
     ModuleBase::matrix& svnl_dalpha,
 #endif
-    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma,
 	Gint_k &gint_k,
     const ORB_gen_tables* uot,
@@ -805,6 +801,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 {
     flk.ftable(isforce,
         isstress,
+        GlobalC::ucell,
         psi,
         pelec,
         foverlap,
@@ -818,7 +815,6 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #if __DEEPKS
         svnl_dalpha,
 #endif
-        gen_h,
         gint_k,
         uot,
         pv,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h

@@ -36,7 +36,6 @@ class Force_Stress_LCAO
 		const elecstate::ElecState* pelec,
         const psi::Psi<T>* psi,
 		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
         const ORB_gen_tables* uot,
@@ -88,7 +87,6 @@ class Force_Stress_LCAO
 #if __DEEPKS
         ModuleBase::matrix& svnl_dalpha,
 #endif
-        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma,
 		Gint_k &gint_k,
         const ORB_gen_tables* uot,