Refactor: Combine gamma-only and multi-k versions of some functions in DeePKS. (#5717)

* Add support for INPUT deepks_v_delta>0 in multi-k points DeePKS calculations

* Refactor: Change LCAO_Deepks_Interface to template class.

* Remove the h_mat and h_mat_k variables in LCAO_Deepks and change H_V_delta to form consistent with H_V_delta_k.

* Change functions in deepks_hmat to template.

* Combine gamma-only and multi-k for v_delta_precalc.

* Change functions about v_delta_precalc and psialpha in deepks_v_delta calculations to templates.

* Change save_npy_h to template.

* Change some functions in LCAO_deepks_io to templates.

* Remove ld.V_deltaR.

* Change cal_orbital_precalc to template.

* Remove orbital_precalc_k.cpp.

* Change cal_gdmx into template function.

* [pre-commit.ci lite] apply automatic fixes

* Update LCAO_deepks_interface.cpp

* Update FORCE_STRESS.cpp

* Update FORCE_gamma.cpp

* Update deepks_lcao.cpp

* Update LCAO_deepks.cpp

* Update LCAO_deepks.cpp

* Update LCAO_deepks.h

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

Author: ErjieWu

source/Makefile.Objects

@@ -201,14 +201,11 @@ OBJS_DEEPKS=LCAO_deepks.o\
         deepks_hmat.o\
         LCAO_deepks_interface.o\
         orbital_precalc.o\
-        orbital_precalc_k.o\
         cal_gdmx.o\
-        cal_gdmx_k.o\
         cal_gedm.o\
         cal_gvx.o\
         cal_descriptor.o\
         v_delta_precalc.o\
-        v_delta_precalc_k.o\
 
 
 OBJS_ELECSTAT=elecstate.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -961,7 +961,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     // 6) write Hamiltonian and Overlap matrix
     for (int ik = 0; ik < this->kv.get_nks(); ++ik)
     {
-        if (PARAM.inp.out_mat_hs[0] || PARAM.inp.deepks_v_delta)
+        if (PARAM.inp.out_mat_hs[0])
         {
             this->p_hamilt->updateHk(ik);
         }
@@ -1000,12 +1000,6 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
                                     this->pv,
                                     GlobalV::DRANK);
             }
-#ifdef __DEEPKS
-            if (PARAM.inp.deepks_out_labels && PARAM.inp.deepks_v_delta)
-            {
-                DeePKS_domain::save_h_mat(h_mat.p, this->pv.nloc, ik);
-            }
-#endif
         }
     }
 
@@ -1023,24 +1017,30 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
 
     //! 8) Write DeePKS information
 #ifdef __DEEPKS
-    std::shared_ptr<LCAO_Deepks> ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {});
-    LCAO_Deepks_Interface LDI = LCAO_Deepks_Interface(ld_shared_ptr);
-    ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
-    LDI.out_deepks_labels(this->pelec->f_en.etot,
-                          this->pelec->klist->get_nks(),
-                          ucell.nat,
-                          PARAM.globalv.nlocal,
-                          this->pelec->ekb,
-                          this->pelec->klist->kvec_d,
-                          ucell,
-                          orb_,
-                          GlobalC::GridD,
-                          &(this->pv),
-                          *(this->psi),
-                          dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                          PARAM.inp.deepks_v_delta);
-
-    ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
+    if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0))
+    {
+        hamilt::HamiltLCAO<TK, TR>* p_ham_deepks
+            = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+        std::shared_ptr<LCAO_Deepks> ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {});
+        LCAO_Deepks_Interface<TK, TR> LDI(ld_shared_ptr);
+
+        ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
+        LDI.out_deepks_labels(this->pelec->f_en.etot,
+                            this->pelec->klist->get_nks(),
+                            ucell.nat,
+                            PARAM.globalv.nlocal,
+                            this->pelec->ekb,
+                            this->pelec->klist->kvec_d,
+                            ucell,
+                            orb_,
+                            GlobalC::GridD,
+                            &(this->pv),
+                            *(this->psi),
+                            dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+                            p_ham_deepks);
+
+        ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
+    }
 #endif
 
     //! 9) Perform RDMFT calculations

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -522,7 +522,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                     {
                         const std::vector<std::vector<double>>& dm_gamma
                             = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec)->get_DM()->get_DMK_vector();
-                        GlobalC::ld.cal_gdmx(dm_gamma[0], ucell, orb, GlobalC::GridD, isstress);
+                        GlobalC::ld.cal_gdmx(dm_gamma, ucell, orb, GlobalC::GridD, kv.get_nks(), kv.kvec_d, isstress);
                     }
                     else
                     {
@@ -531,13 +531,13 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                   ->get_DM()
                                   ->get_DMK_vector();
 
-                        GlobalC::ld.cal_gdmx_k(dm_k,
-                                               ucell,
-                                               orb,
-                                               GlobalC::GridD,
-                                               kv.get_nks(),
-                                               kv.kvec_d,
-                                               isstress);
+                        GlobalC::ld.cal_gdmx(dm_k,
+                                             ucell,
+                                             orb,
+                                             GlobalC::GridD,
+                                             kv.get_nks(),
+                                             kv.kvec_d,
+                                             isstress);
                     }
                     if (PARAM.inp.deepks_out_unittest)
                     {

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -260,15 +260,16 @@ void Force_LCAO<double>::ftable(const bool isforce,
 
         if (PARAM.inp.deepks_out_unittest)
         {
-            LCAO_deepks_io::print_dm(dm_gamma[0], PARAM.globalv.nlocal, this->ParaV->nrow);
+            const int nks = 1; // 1 for gamma-only
+            LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma); 
 
             GlobalC::ld.check_projected_dm();
 
             GlobalC::ld.check_descriptor(ucell, PARAM.globalv.global_out_dir);
 
             GlobalC::ld.check_gedm();
 
-            GlobalC::ld.cal_e_delta_band(dm_gamma);
+            GlobalC::ld.cal_e_delta_band(dm_gamma,nks);
 
             std::ofstream ofs("E_delta_bands.dat");
             ofs << std::setprecision(10) << GlobalC::ld.e_delta_band;

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp

@@ -186,7 +186,7 @@ void DeePKS<OperatorLCAO<std::complex<double>, double>>::contributeHR()
     {
         ModuleBase::timer::tick("DeePKS", "contributeHR");
 
-        GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD);
+        GlobalC::ld.cal_projected_DM(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD);
         GlobalC::ld.cal_descriptor(this->ucell->nat);
         // calculate dE/dD
         GlobalC::ld.cal_gedm(this->ucell->nat);
@@ -219,7 +219,7 @@ void DeePKS<OperatorLCAO<std::complex<double>, std::complex<double>>>::contribut
     {
         ModuleBase::timer::tick("DeePKS", "contributeHR");
 
-        GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD);
+        GlobalC::ld.cal_projected_DM(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD);
         GlobalC::ld.cal_descriptor(this->ucell->nat);
         // calculate dE/dD
         GlobalC::ld.cal_gedm(this->ucell->nat);
@@ -497,7 +497,7 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::cal_HR_IJR(const double* hr_i
 
 inline void get_h_delta_k(int ik, double*& h_delta_k)
 {
-    h_delta_k = GlobalC::ld.H_V_delta.data();
+    h_delta_k = GlobalC::ld.H_V_delta[ik].data();
     return;
 }
 inline void get_h_delta_k(int ik, std::complex<double>*& h_delta_k)

source/module_hamilt_lcao/module_deepks/CMakeLists.txt

@@ -13,14 +13,11 @@ if(ENABLE_DEEPKS)
       deepks_hmat.cpp
       LCAO_deepks_interface.cpp
       orbital_precalc.cpp
-      orbital_precalc_k.cpp
       cal_gdmx.cpp
-      cal_gdmx_k.cpp
       cal_gedm.cpp
       cal_gvx.cpp      
       cal_descriptor.cpp
       v_delta_precalc.cpp
-      v_delta_precalc_k.cpp
   )
 
   add_library(

source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp

@@ -17,7 +17,6 @@
 //4. subroutines that are related to V_delta:
 //  - allocate_V_delta : allocates H_V_delta; if calculating force, it also calls
 //      init_gdmx, as well as allocating F_delta
-//  - allocate_V_deltaR : allcoates H_V_deltaR, for multi-k calculations
 
 #ifdef __DEEPKS
 
@@ -35,7 +34,6 @@ LCAO_Deepks::LCAO_Deepks()
     alpha_index = new ModuleBase::IntArray[1];
     inl_index = new ModuleBase::IntArray[1];
     inl_l = nullptr;
-    H_V_deltaR = nullptr;
     gedm = nullptr;
 }
 
@@ -45,7 +43,6 @@ LCAO_Deepks::~LCAO_Deepks()
     delete[] alpha_index;
     delete[] inl_index;
     delete[] inl_l;
-    delete[] H_V_deltaR;
 
     //=======1. to use deepks, pdm is required==========
     //delete pdm**
@@ -92,7 +89,10 @@ void LCAO_Deepks::init(
 
     int tot_inl = tot_inl_per_atom * nat;
 
-    if(PARAM.inp.deepks_equiv) tot_inl = nat;
+    if(PARAM.inp.deepks_equiv) 
+    {
+        tot_inl = nat;
+    }
 
     this->lmaxd = lm;
     this->nmaxd = nm;
@@ -143,25 +143,6 @@ void LCAO_Deepks::init(
 
     this->pv = &pv_in;
 
-    if(PARAM.inp.deepks_v_delta)
-    {
-        //allocate and init h_mat
-        if(PARAM.globalv.gamma_only_local)
-        {
-            int nloc=this->pv->nloc;
-            this->h_mat.resize(nloc,0.0);
-        }
-        else
-        {
-            int nloc=this->pv->nloc;
-            this->h_mat_k.resize(nks);
-            for (int ik = 0; ik < nks; ik++)
-            {
-                this->h_mat_k[ik].resize(nloc,std::complex<double>(0.0,0.0));
-            }
-        }
-    }
-
     return;
 }
 
@@ -335,8 +316,9 @@ void LCAO_Deepks::allocate_V_delta(const int nat, const int nks)
     //initialize the H matrix H_V_delta
     if(PARAM.globalv.gamma_only_local)
     {
-        this->H_V_delta.resize(pv->nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->H_V_delta.data(), pv->nloc);
+        H_V_delta.resize(1); // the first dimension is for the consistence with H_V_delta_k
+        this->H_V_delta[0].resize(pv->nloc);
+        ModuleBase::GlobalFunc::ZEROS(this->H_V_delta[0].data(), pv->nloc);
     }
     else
     {
@@ -387,15 +369,6 @@ void LCAO_Deepks::allocate_V_delta(const int nat, const int nks)
     return;
 }
 
-void LCAO_Deepks::allocate_V_deltaR(const int nnr)
-{
-    ModuleBase::TITLE("LCAO_Deepks", "allocate_V_deltaR");
-    GlobalV::ofs_running << nnr << std::endl;
-    delete[] H_V_deltaR;
-    H_V_deltaR = new double[nnr];
-    ModuleBase::GlobalFunc::ZEROS(H_V_deltaR, nnr);
-}
-
 void LCAO_Deepks::init_orbital_pdm_shell(const int nks)
 {
 
@@ -541,12 +514,12 @@ void LCAO_Deepks::del_v_delta_pdm_shell(const int nks,const int nlocal)
 
 void LCAO_Deepks::dpks_cal_e_delta_band(const std::vector<std::vector<double>>& dm, const int nks)
 {
-    this->cal_e_delta_band(dm);
+    this->cal_e_delta_band(dm, nks);
 }
 
 void LCAO_Deepks::dpks_cal_e_delta_band(const std::vector<std::vector<std::complex<double>>>& dm, const int nks)
 {
-    this->cal_e_delta_band_k(dm, nks);
+    this->cal_e_delta_band(dm, nks);
 }
 
 #endif