Update ESolver psi and delete some gint codes

source/Makefile.Objects

@@ -277,8 +277,6 @@ OBJS_ESOLVER=esolver.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\
-      lcao_before_scf.o\
-      lcao_after_scf.o\
       esolver_gets.o\
       lcao_others.o\
       esolver_dm2rho.o\
@@ -625,7 +623,6 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_dmk.o\
       unk_overlap_lcao.o\
       read_wfc_nao.o\
-      read_wfc_lcao.o\
       write_wfc_nao.o\
       write_HS_sparse.o\
       single_R_io.o\
@@ -767,6 +764,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     of_stress_pw.o\
     symmetry_rho.o\
     symmetry_rhog.o\
+    setup_psi_pw.o\
     setup_psi.o\
     psi_init.o\
     elecond.o\

source/source_esolver/CMakeLists.txt

@@ -18,8 +18,6 @@ if(ENABLE_LCAO)
   list(APPEND objects
       esolver_ks_lcao.cpp
       esolver_ks_lcao_tddft.cpp
-      lcao_before_scf.cpp
-      lcao_after_scf.cpp
       esolver_gets.cpp
       lcao_others.cpp
       esolver_dm2rho.cpp

source/source_esolver/esolver_double_xc.cpp

@@ -399,8 +399,6 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->pv,
                        this->pelec_base,
                        this->psi,
-                       this->GG, // mohan add 2024-04-01
-                       this->GK, // mohan add 2024-04-01
                        this->two_center_bundle_,
                        this->orb_,
                        force_base,

source/source_esolver/esolver_ks.cpp

@@ -14,6 +14,7 @@
 #include "source_io/output_log.h" // use write_head
 #include "source_estate/elecstate_print.h" // print_etot
 #include "source_io/print_info.h" // print_parameters
+#include "source_psi/setup_psi.h" // mohan add 20251009
 
 namespace ModuleESolver
 {
@@ -28,7 +29,8 @@ ESolver_KS<T, Device>::~ESolver_KS()
 	//****************************************************
 	// do not add any codes in this deconstructor funcion
 	//****************************************************
-	delete this->psi;
+    Setup_Psi<T>::deallocate_psi(this->psi);
+
     delete this->p_hamilt;
     delete this->p_chgmix;
     this->ppcell.release_memory();

source/source_esolver/esolver_ks.h

@@ -71,7 +71,6 @@ class ESolver_KS : public ESolver_FP
     int maxniter;               //! maximum iter steps for scf
     int niter;                  //! iter steps actually used in scf
     bool oscillate_esolver = false; // whether esolver is oscillated
-
     bool scf_nmax_flag = false; // whether scf has reached nmax, mohan add 20250921
 };
 } // namespace ModuleESolver

source/source_esolver/esolver_ks_lcao.cpp

@@ -19,6 +19,9 @@
 #include "source_estate/module_dm/setup_dm.h" // setup dm from electronic wave functions
 #include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
 #include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() 
+#include "source_io/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
+#include "source_psi/setup_psi.h" // mohan add 20251019
+#include "source_io/read_wfc_nao.h" 
 
 namespace ModuleESolver
 {
@@ -73,44 +76,20 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 
     // 5) init electronic wave function psi
-    if (this->psi == nullptr)
-    {
-        int nsk = 0;
-        int ncol = 0;
-        if (PARAM.globalv.gamma_only_local)
-        {
-            nsk = inp.nspin;
-            ncol = this->pv.ncol_bands;
-            if (inp.ks_solver == "genelpa" || inp.ks_solver == "elpa" || inp.ks_solver == "lapack"
-                || inp.ks_solver == "pexsi" || inp.ks_solver == "cusolver"
-                || inp.ks_solver == "cusolvermp")
-            {
-                ncol = this->pv.ncol;
-            }
-        }
-        else
-        {
-            nsk = this->kv.get_nks();
-#ifdef __MPI
-            ncol = this->pv.ncol_bands;
-#else
-            ncol = inp.nbands;
-#endif
-        }
-        this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
-    }
+    Setup_Psi<TK>::allocate_psi(this->psi, this->kv, this->pv, PARAM.inp);
 
-    // 6) read psi from file
+    //! read psi from file
     if (inp.init_wfc == "file" && inp.esolver_type != "tddft")
-	{
-		if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, 
-			 this->pv, *(this->psi), this->pelec, this->pelec->klist->ik2iktot,
-             this->pelec->klist->get_nkstot(), inp.nspin))
+    {
+        if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir,
+             this->pv, *this->psi, this->pelec->ekb, this->pelec->wg, this->kv.ik2iktot,
+             this->kv.get_nkstot(), PARAM.inp.nspin))
         {
             ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
         }
     }
 
+
     // 7) init DMK, but DMR is constructed in before_scf()
     dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin);
 
@@ -183,6 +162,147 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     return;
 }
 
+
+template <typename TK, typename TR>
+void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
+{
+    ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf");
+
+    //! 1) call before_scf() of ESolver_KS.
+    ESolver_KS<TK>::before_scf(ucell, istep);
+
+    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
+    if(!estate)
+    {
+        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","pelec does not exist");
+    }
+
+    //! 2) find search radius
+    double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
+      PARAM.inp.out_level, orb_.get_rcutmax_Phi(), ucell.infoNL.get_rcutmax_Beta(),
+      PARAM.globalv.gamma_only_local);
+
+    //! 3) use search_radius to search adj atoms
+    atom_arrange::search(PARAM.globalv.search_pbc, GlobalV::ofs_running,
+      this->gd, ucell, search_radius, PARAM.inp.test_atom_input);
+
+    //! 4) initialize NAO basis set
+    // here new is a unique pointer, which will be deleted automatically
+    gint_info_.reset(
+        new ModuleGint::GintInfo(
+        this->pw_big->nbx, this->pw_big->nby, this->pw_big->nbz,
+        this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz,
+        0, 0, this->pw_big->nbzp_start,
+        this->pw_big->nbx, this->pw_big->nby, this->pw_big->nbzp,
+        orb_.Phi, ucell, this->gd));
+    ModuleGint::Gint::set_gint_info(gint_info_.get());
+
+    // 7) For each atom, calculate the adjacent atoms in different cells
+    // and allocate the space for H(R) and S(R).
+    // If k point is used here, allocate HlocR after atom_arrange.
+    this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
+
+    // 8) initialize the Hamiltonian operators
+    // if atom moves, then delete old pointer and add a new one
+    if (this->p_hamilt != nullptr)
+    {
+        delete this->p_hamilt;
+        this->p_hamilt = nullptr;
+    }
+    if (this->p_hamilt == nullptr)
+    {
+        elecstate::DensityMatrix<TK, double>* DM = estate->get_DM();
+
+        this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
+            PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
+            PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
+            ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
+            two_center_bundle_, orb_, DM, this->deepks
+#ifdef __EXX
+            ,
+            istep,
+            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->two_level_step : &this->exx_nao.exc->two_level_step,
+            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->get_Hexxs() : nullptr,
+            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exx_nao.exc->get_Hexxs()
+#endif
+        );
+    }
+
+    // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
+    // since it depends on ionic positions
+    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)
+    {
+        spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
+        sc.init_sc(PARAM.inp.sc_thr, PARAM.inp.nsc, PARAM.inp.nsc_min, PARAM.inp.alpha_trial,
+                   PARAM.inp.sccut, PARAM.inp.sc_drop_thr, ucell, &(this->pv),
+                   PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->pelec);
+    }
+
+    // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
+#ifdef __EXX
+    if (PARAM.inp.calculation != "nscf")
+    {
+        if (GlobalC::exx_info.info_ri.real_number)
+        {
+            this->exx_nao.exd->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_);
+        }
+        else
+        {
+            this->exx_nao.exc->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_);
+        }
+    }
+#endif
+
+    // 12) init_scf, should be before_scf? mohan add 2025-03-10
+    this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+
+    // 13) initalize DMR
+    // DMR should be same size with Hamiltonian(R)
+    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+    if(!hamilt_lcao)
+    {
+        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
+    }
+    estate->get_DM()->init_DMR(*hamilt_lcao->getHR());
+
+#ifdef __MLALGO
+    // 14) initialize DMR of DeePKS
+    this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd);
+#endif
+
+    // 15) two cases are considered:
+    // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
+    // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
+    if (istep > 0)
+    {
+        estate->get_DM()->cal_DMR();
+    }
+
+    // 16) the electron charge density should be symmetrized,
+    Symmetry_rho srho;
+    for (int is = 0; is < PARAM.inp.nspin; is++)
+    {
+        srho.begin(is, this->chr, this->pw_rho, ucell.symm);
+    }
+
+    // 17) why we need to set this sentence? mohan add 2025-03-10
+    this->p_hamilt->non_first_scf = istep;
+
+    // 18) update of RDMFT, added by jghan
+    if (PARAM.inp.rdmft == true)
+    {
+        rdmft_solver.update_ion(ucell, *(this->pw_rho), this->locpp.vloc, this->sf.strucFac);
+    }
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf");
+    return;
+}
+
+
 template <typename TK, typename TR>
 double ESolver_KS_LCAO<TK, TR>::cal_energy()
 {
@@ -201,29 +321,13 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     deepks.dpks_out_type = "tot";  // for deepks method
 
-    fsl.getForceStress(ucell,
-                       PARAM.inp.cal_force,
-                       PARAM.inp.cal_stress,
-                       PARAM.inp.test_force,
-                       PARAM.inp.test_stress,
-                       this->gd,
-                       this->pv,
-                       this->pelec,
-                       this->psi,
-                       this->GG, // mohan add 2024-04-01
-                       this->GK, // mohan add 2024-04-01
-                       two_center_bundle_,
-                       orb_,
-                       force,
-                       this->scs,
-                       this->locpp,
-                       this->sf,
-                       this->kv,
-                       this->pw_rho,
-                       this->solvent,
-                       this->deepks,
-                       this->exx_nao,
-                       &ucell.symm);
+    fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, 
+                       PARAM.inp.test_force, PARAM.inp.test_stress,
+                       this->gd, this->pv, this->pelec, this->psi,
+                       two_center_bundle_, orb_, force, this->scs,
+                       this->locpp, this->sf, this->kv,
+                       this->pw_rho, this->solvent, this->deepks,
+                       this->exx_nao, &ucell.symm);
 
     // delete RA after cal_force
     this->RA.delete_grid();
@@ -526,6 +630,50 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 
 }
 
+template <typename TK, typename TR>
+void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)
+{
+    ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");
+
+    //! 1) call after_scf() of ESolver_KS
+    ESolver_KS<TK>::after_scf(ucell, istep, conv_esolver);
+
+    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
+    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+
+    if(!estate)
+    {
+        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","pelec does not exist");
+    }
+
+    if(!hamilt_lcao)
+    {
+        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","p_hamilt does not exist");
+    }
+
+    //! 2) output of lcao every few ionic steps
+    ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
+            PARAM.inp, this->kv, estate, this->pv,
+            this->gd, this->psi, hamilt_lcao,
+            this->two_center_bundle_, this->GK,
+            this->orb_, this->pw_wfc, this->pw_rho,
+            this->GridT, this->pw_big, this->sf,
+            this->rdmft_solver, this->deepks, this->exx_nao, 
+            this->conv_esolver, this->scf_nmax_flag,
+            istep);
+
+
+    //! 3) Clean up RA, which is used to serach for adjacent atoms
+    if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress)
+    {
+        RA.delete_grid();
+    }
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");
+}
+
+
 template class ESolver_KS_LCAO<double, double>;
 template class ESolver_KS_LCAO<std::complex<double>, double>;
 template class ESolver_KS_LCAO<std::complex<double>, std::complex<double>>;

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -102,9 +102,10 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::before_all_runners(UnitCell& ucell, cons
 		if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, 
 					this->pv, 
 					*(this->psi), 
-					this->pelec, 
-                    this->pelec->klist->ik2iktot,
-                    this->pelec->klist->get_nkstot(),
+					this->pelec->ekb,
+                    this->pelec->wg, 
+                    this->kv.ik2iktot,
+                    this->kv.get_nkstot(),
 					PARAM.inp.nspin,
                     0,
                     TD_info::estep_shift))

source/source_esolver/esolver_ks_pw.h

@@ -1,7 +1,7 @@
 #ifndef ESOLVER_KS_PW_H
 #define ESOLVER_KS_PW_H
 #include "./esolver_ks.h"
-#include "source_psi/setup_psi.h" // mohan add 20251012
+#include "source_psi/setup_psi_pw.h" // mohan add 20251012
 #include "source_pw/module_pwdft/VSep_in_pw.h"
 #include "source_pw/module_pwdft/global.h"
 #include "source_pw/module_pwdft/module_exx_helper/exx_helper.h"
@@ -53,7 +53,7 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
     virtual void deallocate_hamilt();
 
     // Electronic wave function psi
-    Setup_Psi<T, Device> stp;
+    Setup_Psi_pw<T, Device> stp;
 
     // DFT-1/2 method
     VSep* vsep_cell = nullptr;