Refactor: remove Hloc, Hloc2, Sloc, Sloc2 in LCAO_Matrix

source/Makefile.Objects

@@ -524,7 +524,6 @@ OBJS_LCAO=DM_gamma.o\
       LCAO_set_fs.o\
       LCAO_set_st.o\
       LCAO_nl_mu.o\
-      LCAO_nl_beta.o\
       LCAO_nnr.o\
       LCAO_set_zero.o\
       LCAO_allocate.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -1010,22 +1010,21 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int iter) {
         && (!GlobalC::exx_info.info_global.separate_loop
             || iter == 1)) // to avoid saving the same value repeatedly
     {
-        std::vector<TK> Hexxk_save(this->orb_con.ParaV.get_local_size());
+        hamilt::HS_Matrix_K<TK> Hexxk_save(&this->orb_con.ParaV, 1);
         for (int ik = 0; ik < this->kv.get_nks(); ++ik) {
-            ModuleBase::GlobalFunc::ZEROS(Hexxk_save.data(), Hexxk_save.size());
+            Hexxk_save.set_zero_hk();
 
-            hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>> opexx_save(
-                &this->LM,
-                nullptr,
-                &Hexxk_save,
-                this->kv);
+            hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>> opexx_save(&Hexxk_save, 
+                                                                         &this->LM, 
+                                                                         nullptr, 
+                                                                         this->kv);
 
             opexx_save.contributeHk(ik);
 
             GlobalC::restart.save_disk("Hexx",
                                        ik,
                                        this->orb_con.ParaV.get_local_size(),
-                                       Hexxk_save.data());
+                                       Hexxk_save.get_hk());
         }
         if (GlobalV::MY_RANK == 0) {
             GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx);

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -170,6 +170,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep) {
             GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr,
             GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK),
             &(this->LM),
+            &this->orb_con.ParaV,
             this->pelec->pot,
             this->kv,
             two_center_bundle_,
@@ -520,11 +521,9 @@ void ESolver_KS_LCAO<std::complex<double>, double>::get_S(void) {
 
     this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL);
 
-    this->LM.ParaV = &this->orb_con.ParaV;
-
     if (this->p_hamilt == nullptr) {
         this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(
-            &this->LM,
+            &this->orb_con.ParaV,
             this->kv,
             *(two_center_bundle_.overlap_orb));
         dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(
@@ -569,7 +568,7 @@ void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S(void) {
     if (this->p_hamilt == nullptr) {
         this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>,
                                                 std::complex<double>>(
-            &this->LM,
+            &this->orb_con.ParaV,
             this->kv,
             *(two_center_bundle_.overlap_orb));
         dynamic_cast<
@@ -593,12 +592,12 @@ void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S(void) {
 }
 
 template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::nscf(void) {
+void ESolver_KS_LCAO<TK, TR>::nscf() {
     ModuleBase::TITLE("ESolver_KS_LCAO", "nscf");
 
     std::cout << " NON-SELF CONSISTENT CALCULATIONS" << std::endl;
 
-    time_t time_start = std::time(NULL);
+    time_t time_start = std::time(nullptr);
 
 #ifdef __EXX
 #ifdef __MPI
@@ -631,7 +630,7 @@ void ESolver_KS_LCAO<TK, TR>::nscf(void) {
                                  "HSolver has not been initialed!");
     }
 
-    time_t time_finish = std::time(NULL);
+    time_t time_finish = std::time(nullptr);
     ModuleBase::GlobalFunc::OUT_TIME("cal_bands", time_start, time_finish);
 
     GlobalV::ofs_running << " end of band structure calculation " << std::endl;

source/module_esolver/esolver_ks_lcao_tmpfunc.cpp

@@ -19,7 +19,8 @@ namespace ModuleESolver
 		GlobalC::dftu.cal_occup_m_gamma(
 				iter, 
 				dm, 
-				this->p_chgmix->get_mixing_beta());
+				this->p_chgmix->get_mixing_beta(),
+                this->p_hamilt);
     }
 
     //! dftu occupation matrix for multiple k-points using dm(complex)

source/module_hamilt_general/matrixblock.h

@@ -12,7 +12,7 @@ template <typename T> struct MatrixBlock
     T* p;
     size_t row;
     size_t col;
-    int* desc;
+    const int* desc;
 };
 
 } // namespace hamilt

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -36,8 +36,7 @@ if(ENABLE_LCAO)
         LCAO_nnr.cpp
 		LCAO_set_fs.cpp
 		LCAO_set_st.cpp
-		LCAO_nl_mu.cpp
-		LCAO_nl_beta.cpp
+		LCAO_nl_mu.cpp	
         LCAO_set_zero.cpp
         LCAO_allocate.cpp
 		record_adj.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -249,15 +249,13 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         }
         else
         {
-            hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(&lm,
+            hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(nullptr, // HK and SK are not used for force&stress
                                                                    kv.kvec_d,
-                                                                   nullptr,
-                                                                   nullptr,
+                                                                   nullptr, // HR are not used for force&stress
                                                                    GlobalC::ucell,
                                                                    &GlobalC::GridD,
                                                                    two_center_bundle.overlap_orb_onsite.get(),
-                                                                   &GlobalC::dftu,
-                                                                   *(lm.ParaV));
+                                                                   &GlobalC::dftu);
 
             tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu);
         }

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -77,7 +77,7 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                               pv,
                               two_center_bundle,
                               &GlobalC::GridD,
-                              lm.Sloc.data());
+                              nullptr);
 
     // calculate dT in LCAP
     // allocation dt
@@ -101,11 +101,11 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                               pv,
                               two_center_bundle,
                               &GlobalC::GridD,
-                              lm.Hloc_fixed.data());
+                              nullptr);
 
     LCAO_domain::build_Nonlocal_mu_new(lm,
                                        fsr,
-                                       lm.Hloc_fixed.data(),
+                                       nullptr,
                                        cal_deri,
                                        GlobalC::ucell,
                                        GlobalC::ORB,
@@ -116,8 +116,9 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
     if (INPUT.cal_syns)
     {
         cal_deri = false;
-
-        lm.zeros_HSgamma('S');
+        ModuleBase::WARNING_QUIT("cal_syns", "this function has been broken and will be fixed later.");
+        /*
+        std::vector<double> Sloc(pv.nloc, 0.0);
 
         LCAO_domain::build_ST_new(lm,
                                   fsr,
@@ -128,7 +129,7 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                                   pv,
                                   two_center_bundle,
                                   &GlobalC::GridD,
-                                  lm.Sloc.data(),
+                                  Sloc.data(),
                                   INPUT.cal_syns,
                                   INPUT.dmax);
 
@@ -147,7 +148,7 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                            GlobalV::DRANK);
 
         ModuleIO::save_mat(0,
-                           lm.Sloc.data(),
+                           Sloc.data(),
                            GlobalV::NLOCAL,
                            bit,
                            GlobalV::out_ndigits,
@@ -156,7 +157,7 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                            "S",
                            "data-" + std::to_string(0),
                            pv,
-                           GlobalV::DRANK);
+                           GlobalV::DRANK);*/
     }
 
     ModuleBase::timer::tick("Force_LCAO", "allocate");

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -131,7 +131,7 @@ void Force_LCAO<std::complex<double>>::allocate(const Parallel_Orbitals& pv,
     // calculate dVnl=<phi|dVnl|dphi> in LCAO
     LCAO_domain::build_Nonlocal_mu_new(lm,
                                        fsr,
-                                       lm.Hloc_fixed.data(),
+                                       nullptr,
                                        cal_deri,
                                        GlobalC::ucell,
                                        GlobalC::ORB,
@@ -143,6 +143,9 @@ void Force_LCAO<std::complex<double>>::allocate(const Parallel_Orbitals& pv,
     {
         cal_deri = false;
 
+        ModuleBase::WARNING_QUIT("cal_syns",
+                                 "This function has been broken and will be fixed later.");
+
         LCAO_domain::build_ST_new(lm,
                                   fsr,
                                   'S',
@@ -158,10 +161,9 @@ void Force_LCAO<std::complex<double>>::allocate(const Parallel_Orbitals& pv,
 
         for (int ik = 0; ik < nks; ik++)
         {
-            lm.zeros_HSk('S');
-            lm.folding_fixedH(ik, kvec_d, 1);
+
             bool bit = false; // LiuXh, 2017-03-21
-            ModuleIO::save_mat(0,
+            /*ModuleIO::save_mat(0,
                                lm.Hloc2.data(),
                                GlobalV::NLOCAL,
                                bit,
@@ -183,7 +185,7 @@ void Force_LCAO<std::complex<double>>::allocate(const Parallel_Orbitals& pv,
                                "S",
                                "data-" + std::to_string(ik),
                                pv,
-                               GlobalV::DRANK);
+                               GlobalV::DRANK);*/
         }
     }
 

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp

@@ -12,12 +12,6 @@ void divide_HS_in_frag(LCAO_Matrix& lm, const bool isGamma, Parallel_Orbitals& p
     //(1), (2): set up matrix division have been moved into ORB_control
     // just pass `ParaV` as pointer is enough
     lm.ParaV = &pv;
-    // (3) allocate for S, H_fixed, H, and S_diag
-    if (isGamma) {
-        LCAO_domain::allocate_HS_gamma(lm, lm.ParaV->nloc);
-    } else {
-        LCAO_domain::allocate_HS_k(lm, lm.ParaV->nloc);
-    }
 #ifdef __DEEPKS
     // wenfei 2021-12-19
     // preparation for DeePKS
@@ -48,51 +42,4 @@ void divide_HS_in_frag(LCAO_Matrix& lm, const bool isGamma, Parallel_Orbitals& p
     return;
 }
 
-void allocate_HS_gamma(LCAO_Matrix& lm, const long& nloc) {
-    ModuleBase::TITLE("LCAO_Matrix", "allocate_HS_gamma");
-
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nloc", nloc);
-
-    if (nloc == 0) {
-        return;
-    }
-
-    // because we initilize in the constructor function
-    // with dimension '1', so here we reconstruct these
-    // matrices
-
-    lm.Sloc.resize(nloc);
-    lm.Hloc_fixed.resize(nloc);
-    lm.Hloc.resize(nloc);
-
-    ModuleBase::GlobalFunc::ZEROS(lm.Sloc.data(), nloc);
-    ModuleBase::GlobalFunc::ZEROS(lm.Hloc_fixed.data(), nloc);
-    ModuleBase::GlobalFunc::ZEROS(lm.Hloc.data(), nloc);
-
-    return;
-}
-
-void allocate_HS_k(LCAO_Matrix& lm, const long& nloc) {
-    ModuleBase::TITLE("LCAO_Matrix", "allocate_HS_k");
-
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nloc", nloc);
-
-    if (nloc == 0) {
-        return; // mohan fix bug 2012-05-25
-    }
-
-    // because we initilize in the constructor function
-    // with dimension '1', so here we reconstruct these
-    // matrices
-    lm.Sloc2.resize(nloc);
-    lm.Hloc_fixed2.resize(nloc);
-    lm.Hloc2.resize(nloc);
-
-    ModuleBase::GlobalFunc::ZEROS(lm.Sloc2.data(), nloc);
-    ModuleBase::GlobalFunc::ZEROS(lm.Hloc_fixed2.data(), nloc);
-    ModuleBase::GlobalFunc::ZEROS(lm.Hloc2.data(), nloc);
-
-    return;
-}
-
 }

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h

@@ -16,14 +16,6 @@
 namespace LCAO_domain
 {
 
-// can used in gamma algorithm.
-void build_Nonlocal_beta_new(LCAO_Matrix& lm,
-                             double* Hloc,
-                             const UnitCell& ucell,
-                             const LCAO_Orbitals& orb,
-                             const TwoCenterIntegrator& intor_orb_beta,
-                             Grid_Driver* GridD);
-
 void build_Nonlocal_mu_new(LCAO_Matrix& lm,
                            ForceStressArrays& fsr, // mohan 2024-06-16
                            double* HlocR,
@@ -173,10 +165,6 @@ void zeros_HSR(const char& mtype, LCAO_HS_Arrays& HS_arrays);
 
 void divide_HS_in_frag(LCAO_Matrix& lm, const bool isGamma, Parallel_Orbitals& pv, const int& nks);
 
-void allocate_HS_k(LCAO_Matrix& lm, const long& nloc);
-
-void allocate_HS_gamma(LCAO_Matrix& lm, const long& nloc);
-
 } // namespace LCAO_domain
 
 #endif