Update DeePKS and do some reconstructions

source/Makefile.Objects

@@ -655,7 +655,7 @@ OBJS_LCAO=evolve_elec.o\
 	  FORCE_k.o\
       stress_tools.o\
       edm.o\
-      pulay_force_stress_center2.o\
+      pulay_fs_center2.o\
 	  grid_init.o\
       spar_dh.o\
       spar_exx.o\

source/source_basis/module_nao/two_center_bundle.h

@@ -1,5 +1,5 @@
-#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H
-#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASIS_MODULE_NAO_TWO_CENTER_BUNDLE_H
+#ifndef TWO_CENTER_BUNDLE_H
+#define TWO_CENTER_BUNDLE_H
 
 #include "source_basis/module_ao/ORB_read.h"
 #include "source_basis/module_nao/two_center_integrator.h"

source/source_esolver/esolver_double_xc.cpp

@@ -6,7 +6,6 @@
 #include "source_lcao/module_deepks/LCAO_deepks.h"
 #include "source_lcao/module_deepks/LCAO_deepks_interface.h"
 #endif
-//-----force& stress-------------------
 #include "source_lcao/FORCE_STRESS.h"
 
 //-----HSolver ElecState Hamilt--------
@@ -16,6 +15,7 @@
 #include "source_hsolver/hsolver_lcao.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_io/write_HS.h" // use ModuleIO::write_hsk()
+#include "source_lcao/setup_deepks.h" // use deepks, mohan add 2025-10-10
 
 namespace ModuleESolver
 {
@@ -388,6 +388,8 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
     // set as base functional Temporarily
     XC_Functional::set_xc_type(PARAM.inp.deepks_out_base);
 
+    this->deepks.dpks_out_type = "base";  // for deepks method
+
     fsl.getForceStress(ucell,
                        PARAM.inp.cal_force,
                        PARAM.inp.cal_stress,
@@ -408,12 +410,10 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->kv,
                        this->pw_rho,
                        this->solvent,
-#ifdef __MLALGO
-                       this->deepks.ld,
-                       "base",
-#endif
+                       this->deepks,
 					   this->exx_nao,
 					   &ucell.symm);
+
     // restore to original xc
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -7,10 +7,6 @@
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
 #include "source_lcao/module_dftu/dftu.h"
-#ifdef __MLALGO
-#include "source_lcao/module_deepks/LCAO_deepks.h"
-#include "source_lcao/module_deepks/LCAO_deepks_interface.h"
-#endif
 #include "source_lcao/FORCE_STRESS.h"
 #include "source_estate/elecstate_lcao.h"
 #include "source_lcao/hamilt_lcao.h"
@@ -203,6 +199,8 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     Force_Stress_LCAO<TK> fsl(this->RA, ucell.nat);
 
+    deepks.dpks_out_type = "tot";  // for deepks method
+
     fsl.getForceStress(ucell,
                        PARAM.inp.cal_force,
                        PARAM.inp.cal_stress,
@@ -223,10 +221,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        this->kv,
                        this->pw_rho,
                        this->solvent,
-#ifdef __MLALGO
-                       this->deepks.ld,
-                       "tot",
-#endif
+                       this->deepks,
                        this->exx_nao,
                        &ucell.symm);
 
@@ -337,43 +332,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
 		}
 #endif
 		elecstate::setup_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step);
-
-/*
-    if (iter == 1 && exx_two_level_step == 0)
-    {
-        std::cout << " WAVEFUN -> CHARGE " << std::endl;
-
-        // calculate the density matrix using read in wave functions
-        // and then calculate the charge density on grid.
-
-        estate->skip_weights = true;
-        elecstate::calculate_weights(estate->ekb,
-                estate->wg,
-                estate->klist,
-                estate->eferm,
-                estate->f_en,
-                estate->nelec_spin,
-                estate->skip_weights);
-
-        elecstate::calEBand(estate->ekb, estate->wg, estate->f_en);
-        elecstate::cal_dm_psi(estate->DM->get_paraV_pointer(), estate->wg, *this->psi, *(estate->DM));
-        estate->DM->cal_DMR();
-
-        estate->psiToRho(*this->psi);
-        estate->skip_weights = false;
-
-        elecstate::cal_ux(ucell);
-
-        //! update the potentials by using new electron charge density
-        estate->pot->update_from_charge(&this->chr, &ucell);
-
-        //! compute the correction energy for metals
-        estate->f_en.descf = estate->cal_delta_escf();
-    }
-
-*/
-
-
 	}
 
 #ifdef __EXX
@@ -410,11 +368,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     {
         this->p_hamilt->refresh();
     }
-    // if (iter == 1 && istep == 0)
-    // {
-    //     // initialize DMR
-    //     this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd);
-    // }
 #endif
 
     if (PARAM.inp.vl_in_h)
@@ -591,11 +544,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     // control the output related to the finished iteration
     ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, estate,
       this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 
-      hamilt_lcao, this->orb_, 
-#ifdef __MLALGO
-      this->deepks.ld,
-#endif
-      exx_nao, iter, istep, conv_esolver, this->scf_ene_thr);
+      hamilt_lcao, this->orb_, this->deepks, 
+      this->exx_nao, iter, istep, conv_esolver, this->scf_ene_thr);
 
 }
 

source/source_esolver/lcao_before_scf.cpp

@@ -6,9 +6,6 @@
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_estate/elecstate_tools.h"
-#ifdef __MLALGO
-#include "source_lcao/module_deepks/LCAO_deepks.h"
-#endif
 #include "source_estate/elecstate_lcao.h"
 #include "source_lcao/LCAO_domain.h"
 #include "source_lcao/module_operator_lcao/op_exx_lcao.h"
@@ -140,25 +137,10 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         );
     }
 
-#ifdef __MLALGO
+
     // 9) 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;
-        // allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
-        DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->deepks.ld.phialpha);
-        // build and save <phi(0)|alpha(R)> at beginning
-        DeePKS_domain::build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd,
-          pv, *(two_center_bundle_.overlap_orb_alpha), this->deepks.ld.phialpha);
-
-        if (PARAM.inp.deepks_out_unittest)
-        {
-            DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_,
-              this->gd, pv, this->deepks.ld.phialpha, GlobalV::MY_RANK);
-        }
-    }
-#endif
+    this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), PARAM.inp);
 
     // 10) prepare sc calculation
     if (PARAM.inp.sc_mag_switch)

source/source_esolver/lcao_others.cpp

@@ -15,9 +15,6 @@
 #include "source_base/timer.h"
 #include "source_cell/module_neighbor/sltk_atom_arrange.h"
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
-#ifdef __MLALGO
-#include "source_lcao/module_deepks/LCAO_deepks.h"
-#endif
 #include "source_lcao/LCAO_domain.h"
 #include "source_lcao/module_operator_lcao/op_exx_lcao.h"
 #include "source_lcao/module_operator_lcao/operator_lcao.h"
@@ -242,35 +239,11 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         );
     }
 
-#ifdef __MLALGO
+
     // 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;
-        // allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
-        DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->deepks.ld.phialpha);
-        // build and save <phi(0)|alpha(R)> at beginning
-        DeePKS_domain::build_phialpha(PARAM.inp.cal_force,
-                                      ucell,
-                                      orb_,
-                                      this->gd,
-                                      pv,
-                                      *(two_center_bundle_.overlap_orb_alpha),
-                                      this->deepks.ld.phialpha);
+    this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), PARAM.inp);
 
-        if (PARAM.inp.deepks_out_unittest)
-        {
-            DeePKS_domain::check_phialpha(PARAM.inp.cal_force,
-                                          ucell,
-                                          orb_,
-                                          this->gd,
-                                          pv,
-                                          this->deepks.ld.phialpha,
-                                          GlobalV::MY_RANK);
-        }
-    }
-#endif
     if (PARAM.inp.sc_mag_switch)
     {
         spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();

source/source_io/berryphase.h

@@ -4,7 +4,6 @@
 #ifdef __LCAO
 #include "unk_overlap_lcao.h"
 #endif
-//#include "source_basis/module_pw/pw_basis.h"
 #include "source_basis/module_pw/pw_basis_k.h"
 #include "source_cell/klist.h"
 #include "source_psi/psi.h"

source/source_io/ctrl_iter_lcao.cpp

@@ -21,9 +21,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
         Charge_Mixing* p_chgmix, // charge mixing *
 		hamilt::HamiltLCAO<TK, TR>* p_hamilt, // hamiltonian *
 		LCAO_Orbitals &orb, // orbital info *
-#ifdef __MLALGO
-		LCAO_Deepks<TK>& ld,
-#endif
+        Setup_DeePKS<TK> &deepks,
         Exx_NAO<TK> &exx_nao,
         int &iter,
         const int istep,
@@ -65,7 +63,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
     {
         if (iter % inp.deepks_out_freq_elec == 0 )
         {
-            std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&ld, [](LCAO_Deepks<TK>*) {});
+            std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&deepks.ld, [](LCAO_Deepks<TK>*) {});
             LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
             deepks_interface.out_deepks_labels(pelec->f_en.etot, kv.get_nks(),
@@ -91,9 +89,7 @@ template void ctrl_iter_lcao<double, double>(UnitCell& ucell, // unit cell *
         Charge_Mixing* p_chgmix, // charge mixing *
 		hamilt::HamiltLCAO<double, double>* p_hamilt, // hamiltonian *
 		LCAO_Orbitals &orb, // orbital info *
-#ifdef __MLALGO
-		LCAO_Deepks<double>& ld,
-#endif
+        Setup_DeePKS<double> &deepks,
         Exx_NAO<double> &exx_nao,
         int &iter,
         const int istep,
@@ -112,9 +108,7 @@ template void ctrl_iter_lcao<std::complex<double>, double>(UnitCell& ucell, // u
         Charge_Mixing* p_chgmix, // charge mixing *
 		hamilt::HamiltLCAO<std::complex<double>, double>* p_hamilt, // hamiltonian *
 		LCAO_Orbitals &orb, // orbital info *
-#ifdef __MLALGO
-		LCAO_Deepks<std::complex<double>>& ld,
-#endif
+        Setup_DeePKS<std::complex<double>> &deepks,
         Exx_NAO<std::complex<double>> &exx_nao,
         int &iter,
         const int istep,
@@ -133,9 +127,7 @@ template void ctrl_iter_lcao<std::complex<double>, std::complex<double>>(UnitCel
         Charge_Mixing* p_chgmix, // charge mixing *
 		hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>* p_hamilt, // hamiltonian *
 		LCAO_Orbitals &orb, // orbital info *
-#ifdef __MLALGO
-		LCAO_Deepks<std::complex<double>>& ld,
-#endif
+        Setup_DeePKS<std::complex<double>> &deepks,
         Exx_NAO<std::complex<double>> &exx_nao,
         int &iter,
         const int istep,

source/source_io/ctrl_iter_lcao.h

@@ -9,6 +9,7 @@
 #include "source_estate/module_charge/charge_mixing.h" // use charge mixing
 #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO<TK, TR>
 #include "source_lcao/setup_exx.h" // mohan add 20251008
+#include "source_lcao/setup_deepks.h" // mohan add 20251010
 
 namespace ModuleIO
 {
@@ -25,9 +26,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
         Charge_Mixing* p_chgmix, // charge mixing *
 		hamilt::HamiltLCAO<TK, TR>* p_hamilt, // hamiltonian *
 		LCAO_Orbitals &orb, // orbital info *
-#ifdef __MLALGO
-		LCAO_Deepks<TK>& ld,
-#endif
+        Setup_DeePKS<TK> &deepks,
         Exx_NAO<TK> &exx_nao,
         int &iter,
         const int istep,

source/source_lcao/CMakeLists.txt

@@ -22,7 +22,7 @@ if(ENABLE_LCAO)
         module_operator_lcao/td_pot_hybrid.cpp
         module_operator_lcao/dspin_lcao.cpp
         module_operator_lcao/dftu_lcao.cpp
-        pulay_force_stress_center2.cpp
+        pulay_fs_center2.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_k.cpp

source/source_lcao/FORCE.h

@@ -13,6 +13,7 @@
 #include "source_lcao/module_gint/gint_gamma.h"
 #include "source_lcao/module_gint/gint_k.h"
 #include "source_psi/psi.h"
+#include "source_lcao/setup_deepks.h"
 
 #ifndef TGINT_H
 #define TGINT_H
@@ -63,11 +64,9 @@ class Force_LCAO
                 ModuleBase::matrix& stvnl_dphi,
                 ModuleBase::matrix& svnl_dbeta,
                 ModuleBase::matrix& svl_dphi,
-#ifdef __MLALGO
                 ModuleBase::matrix& fvnl_dalpha,
                 ModuleBase::matrix& svnl_dalpha,
-                LCAO_Deepks<T>& ld,
-#endif
+                Setup_DeePKS<T>& deepks,
                 typename TGint<T>::type& gint,
                 const TwoCenterBundle& two_center_bundle,
                 const LCAO_Orbitals& orb,