Refactor: Replace nlm_save in DeePKS by HContainer object phialpha.

docs/advanced/input_files/input-main.md

@@ -2060,7 +2060,7 @@ Warning: this function is not robust enough for the current version. Please try
 - **Type**: int
 - **Availability**: numerical atomic orbital basis
 - **Description**: Include V_delta label for DeePKS training. When `deepks_out_labels` is true and `deepks_v_delta` > 0, ABACUS will output h_base.npy, v_delta.npy and h_tot.npy(h_tot=h_base+v_delta). 
-  Meanwhile, when `deepks_v_delta` equals 1, ABACUS will also output v_delta_precalc.npy, which is used to calculate V_delta during DeePKS training. However, when the number of atoms grows, the size of v_delta_precalc.npy will be very large. In this case, it's recommended to set `deepks_v_delta` as 2, and ABACUS will output psialpha.npy and grad_evdm.npy but not v_delta_precalc.npy. These two files are small and can be used to calculate v_delta_precalc in the procedure of training DeePKS.
+  Meanwhile, when `deepks_v_delta` equals 1, ABACUS will also output v_delta_precalc.npy, which is used to calculate V_delta during DeePKS training. However, when the number of atoms grows, the size of v_delta_precalc.npy will be very large. In this case, it's recommended to set `deepks_v_delta` as 2, and ABACUS will output phialpha.npy and grad_evdm.npy but not v_delta_precalc.npy. These two files are small and can be used to calculate v_delta_precalc in the procedure of training DeePKS.
 - **Default**: 0
 
 ### deepks_out_unittest

source/Makefile.Objects

@@ -189,13 +189,12 @@ OBJS_CELL=atom_pseudo.o\
     check_atomic_stru.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
-        deepks_fgamma.o\
-        deepks_fk.o\
+        deepks_force.o\
         LCAO_deepks_odelta.o\
         LCAO_deepks_io.o\
         LCAO_deepks_mpi.o\
         LCAO_deepks_pdm.o\
-        LCAO_deepks_psialpha.o\
+        LCAO_deepks_phialpha.o\
         LCAO_deepks_torch.o\
         LCAO_deepks_vdelta.o\
         deepks_hmat.o\

source/module_esolver/lcao_before_scf.cpp

@@ -205,17 +205,19 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 
 #ifdef __DEEPKS
-    // for each ionic step, the overlap <psi|alpha> must be rebuilt
+    // for each ionic step, the overlap <phi|alpha> must be rebuilt
     // since it depends on ionic positions
     if (PARAM.globalv.deepks_setorb)
     {
         const Parallel_Orbitals* pv = &this->pv;
-        // build and save <psi(0)|alpha(R)> at beginning
-        GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
+        // allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
+        GlobalC::ld.allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
+        // build and save <phi(0)|alpha(R)> at beginning
+        GlobalC::ld.build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
 
         if (PARAM.inp.deepks_out_unittest)
         {
-            GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
+            GlobalC::ld.check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
         }
     }
 #endif

source/module_esolver/lcao_others.cpp

@@ -211,17 +211,19 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     }
 
 #ifdef __DEEPKS
-    // for each ionic step, the overlap <psi|alpha> must be rebuilt
+    // for each ionic step, the overlap <phi|alpha> must be rebuilt
     // since it depends on ionic positions
     if (PARAM.globalv.deepks_setorb)
     {
         const Parallel_Orbitals* pv = &this->pv;
-        // build and save <psi(0)|alpha(R)> at beginning
-        GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
+        // allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
+        GlobalC::ld.allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
+        // build and save <phi(0)|alpha(R)> at beginning
+        GlobalC::ld.build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
 
         if (PARAM.inp.deepks_out_unittest)
         {
-            GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
+            GlobalC::ld.check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
         }
     }
 #endif

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -12,12 +12,11 @@
 #include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10
 #include "module_hamilt_general/module_vdw/vdw.h"
 #include "module_parameter/parameter.h"
-#ifdef __DEEPKS
 #include "module_elecstate/elecstate_lcao.h"
+#ifdef __DEEPKS
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"    //caoyu add for deepks 2021-06-03
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h" // mohan add 2024-07-22
 #endif
-#include "module_elecstate/elecstate_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dspin_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h"
@@ -511,7 +510,14 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                     {
                         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, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress);
+                        GlobalC::ld.cal_gdmx(dm_gamma,
+                                             ucell,
+                                             orb,
+                                             gd,
+                                             kv.get_nks(),
+                                             kv.kvec_d,
+                                             GlobalC::ld.phialpha,
+                                             isstress);
                     }
                     else
                     {
@@ -520,7 +526,8 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                                   ->get_DM()
                                   ->get_DMK_vector();
 
-                        GlobalC::ld.cal_gdmx(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress);
+                        GlobalC::ld
+                            .cal_gdmx(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, GlobalC::ld.phialpha, isstress);
                     }
                     if (PARAM.inp.deepks_out_unittest)
                     {

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -1,19 +1,19 @@
 #include "FORCE.h"
 #include "module_base/memory.h"
-#include "module_parameter/parameter.h"
 #include "module_base/parallel_reduce.h"
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_parameter/parameter.h"
 #ifdef __DEEPKS
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks on 20210813
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h"
 #endif
 #include "module_cell/module_neighbor/sltk_grid_driver.h" //GridD
 #include "module_elecstate/elecstate_lcao.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
-#include "module_io/write_HS.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"
+#include "module_io/write_HS.h"
 
 template <>
 void Force_LCAO<double>::allocate(const UnitCell& ucell,
@@ -147,28 +147,28 @@ void Force_LCAO<double>::finish_ftable(ForceStressArrays& fsr)
     return;
 }
 
-//template <>
-//void Force_LCAO<double>::test(Parallel_Orbitals& pv, double* mm, const std::string& name)
+// template <>
+// void Force_LCAO<double>::test(Parallel_Orbitals& pv, double* mm, const std::string& name)
 //{
-//    std::cout << "\n PRINT " << name << std::endl;
-//    std::cout << std::setprecision(6) << std::endl;
-//    for (int i = 0; i < PARAM.globalv.nlocal; i++)
-//    {
-//        for (int j = 0; j < PARAM.globalv.nlocal; j++)
-//        {
-//            if (std::abs(mm[i * PARAM.globalv.nlocal + j]) > 1.0e-5)
-//            {
-//                std::cout << std::setw(12) << mm[i * PARAM.globalv.nlocal + j];
-//            }
-//            else
-//            {
-//                std::cout << std::setw(12) << "0";
-//            }
-//        }
-//        std::cout << std::endl;
-//    }
-//    return;
-//}
+//     std::cout << "\n PRINT " << name << std::endl;
+//     std::cout << std::setprecision(6) << std::endl;
+//     for (int i = 0; i < PARAM.globalv.nlocal; i++)
+//     {
+//         for (int j = 0; j < PARAM.globalv.nlocal; j++)
+//         {
+//             if (std::abs(mm[i * PARAM.globalv.nlocal + j]) > 1.0e-5)
+//             {
+//                 std::cout << std::setw(12) << mm[i * PARAM.globalv.nlocal + j];
+//             }
+//             else
+//             {
+//                 std::cout << std::setw(12) << "0";
+//             }
+//         }
+//         std::cout << std::endl;
+//     }
+//     return;
+// }
 
 // be called in force_lo.cpp
 template <>
@@ -210,20 +210,40 @@ void Force_LCAO<double>::ftable(const bool isforce,
     // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
     this->allocate(ucell, gd, pv, fsr, two_center_bundle, orb);
 
-    const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z };
-    const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 };
+    const double* dSx[3] = {fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z};
+    const double* dSxy[6] = {fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33};
     // calculate the force related to 'energy density matrix'.
-    PulayForceStress::cal_pulay_fs(foverlap, soverlap,
+    PulayForceStress::cal_pulay_fs(
+        foverlap,
+        soverlap,
         this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
-        ucell, pv, dSx, dSxy, isforce, isstress);
-
-    const double* dHx[3] = { fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z };
-    const double* dHxy[6] = { fsr.DHloc_fixed_11, fsr.DHloc_fixed_12, fsr.DHloc_fixed_13, fsr.DHloc_fixed_22, fsr.DHloc_fixed_23, fsr.DHloc_fixed_33 };
-    //tvnl_dphi
+        ucell,
+        pv,
+        dSx,
+        dSxy,
+        isforce,
+        isstress);
+
+    const double* dHx[3] = {fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z};
+    const double* dHxy[6] = {fsr.DHloc_fixed_11,
+                             fsr.DHloc_fixed_12,
+                             fsr.DHloc_fixed_13,
+                             fsr.DHloc_fixed_22,
+                             fsr.DHloc_fixed_23,
+                             fsr.DHloc_fixed_33};
+    // tvnl_dphi
     PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);
 
     // vl_dphi
-    PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);
+    PulayForceStress::cal_pulay_fs(fvl_dphi,
+                                   svl_dphi,
+                                   *dm,
+                                   ucell,
+                                   pelec->pot,
+                                   gint,
+                                   isforce,
+                                   isstress,
+                                   false /*reset dm to gint*/);
 
 #ifdef __DEEPKS
     if (PARAM.inp.deepks_scf)
@@ -237,21 +257,21 @@ void Force_LCAO<double>::ftable(const bool isforce,
 
         GlobalC::ld.cal_gedm(ucell.nat);
 
-        const int nks=1;
-		DeePKS_domain::cal_f_delta_gamma(dm_gamma, 
-				                         ucell, 
-				                         orb, 
-				                         gd, 
-                                         *this->ParaV,
-                                         GlobalC::ld.lmaxd,
-                                         nks, 
-				                         kv->kvec_d, 
-                                         GlobalC::ld.nlm_save,
-                                         GlobalC::ld.gedm,
-                                         GlobalC::ld.inl_index,
-                                         GlobalC::ld.F_delta,
-				                         isstress, 
-				                         svnl_dalpha);
+        const int nks = 1;
+        DeePKS_domain::cal_f_delta<double>(dm_gamma,
+                                           ucell,
+                                           orb,
+                                           gd,
+                                           *this->ParaV,
+                                           GlobalC::ld.lmaxd,
+                                           nks,
+                                           kv->kvec_d,
+                                           GlobalC::ld.phialpha,
+                                           GlobalC::ld.gedm,
+                                           GlobalC::ld.inl_index,
+                                           GlobalC::ld.F_delta,
+                                           isstress,
+                                           svnl_dalpha);
 
 #ifdef __MPI
         Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc);
@@ -265,15 +285,15 @@ void Force_LCAO<double>::ftable(const bool isforce,
         if (PARAM.inp.deepks_out_unittest)
         {
             const int nks = 1; // 1 for gamma-only
-            LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma); 
+            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,nks);
+            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;