Reducing the size of esolver_ks_lcao

source/Makefile.Objects

@@ -634,6 +634,7 @@ OBJS_LCAO=evolve_elec.o\
       spar_hsr.o\
       spar_st.o\
       spar_u.o\
+      LCAO_set.o\
       LCAO_set_fs.o\
       LCAO_set_st.o\
       LCAO_nl_mu.o\

source/source_esolver/esolver_ks_lcao.cpp

@@ -17,10 +17,16 @@
 #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" 
 #include "source_io/print_info.h"
 #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
+#include "source_lcao/LCAO_set.h" // mohan add 20251111
+
+
+// tmp
+#include "source_psi/setup_psi.h" // use Setup_Psi
+#include "source_io/read_wfc_nao.h" // use read_wfc_nao
+#include "source_estate/elecstate_tools.h" // use fixed_weights
+
 
 namespace ModuleESolver
 {
@@ -74,62 +80,17 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         return;
     }
 
-    // 5) init electronic wave function psi
-    Setup_Psi<TK>::allocate_psi(this->psi, this->kv, this->pv, inp);
-
-    //! 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->ekb, this->pelec->wg, this->kv.ik2iktot,
-             this->kv.get_nkstot(), inp.nspin))
-        {
-            ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
-        }
-    }
-
+    LCAO_domain::set_psi_occ_dm_chg<TK>(this->kv, this->psi, this->pv, this->pelec,
+      this->dmat, this->chr, inp);
 
-    // 7) init DMK, but DMR is constructed in before_scf()
-    this->dmat.allocate_dm(&this->kv, &this->pv, inp.nspin);
+    LCAO_domain::set_pot<TK>(ucell, this->kv, this->sf, *this->pw_rho, *this->pw_rhod,
+      this->pelec, this->orb_, this->pv, this->locpp, this->dftu,
+      this->solvent, this->exx_nao, this->deepks, inp);
 
-    // 8) init exact exchange calculations
-    this->exx_nao.before_runner(ucell, this->kv, this->orb_, this->pv, inp);
-
-    // 9) initialize DFT+U
-    if (inp.dft_plus_u)
-    {
-        this->dftu.init(ucell, &this->pv, this->kv.get_nks(), &orb_);
-    }
-
-    // 10) init local pseudopotentials
-    this->locpp.init_vloc(ucell, this->pw_rho);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
-
-    // 11) init charge density
-    this->chr.allocate(inp.nspin);
-
-    // 12) init potentials
-    if (this->pelec->pot == nullptr)
-    {
-        this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho,
-          &ucell, &(this->locpp.vloc), &(this->sf), &(this->solvent),
-          &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc));
-    }
-
-    // 13) init deepks
-    this->deepks.before_runner(ucell, this->kv.get_nks(), this->orb_, this->pv, inp);
-
-    // 14) set occupations, tddft does not need to set occupations in the first scf
-    if (inp.ocp && inp.esolver_type != "tddft")
-    {
-        elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec,
-          this->pelec->klist, this->pelec->wg, this->pelec->skip_weights);
-    }
-
-    // 15) if kpar is not divisible by nks, print a warning
+    //! if kpar is not divisible by nks, print a warning
     ModuleIO::print_kpar(this->kv.get_nks(), PARAM.globalv.kpar_lcao);
 
-    // 16) init rdmft, added by jghan
+    //! init rdmft, added by jghan
     if (inp.rdmft == true)
     {
         rdmft_solver.init(this->pv, ucell,
@@ -310,8 +271,6 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
 
     ESolver_KS<TK>::after_all_runners(ucell);
 
-    const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1;
-
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
 	if(!hamilt_lcao)
 	{

source/source_lcao/CMakeLists.txt

@@ -33,6 +33,7 @@ if(ENABLE_LCAO)
         spar_hsr.cpp
         spar_st.cpp
         spar_u.cpp
+        LCAO_set.cpp
 		LCAO_set_fs.cpp
 		LCAO_set_st.cpp
 		LCAO_nl_mu.cpp	

source/source_lcao/LCAO_set.cpp

@@ -0,0 +1,147 @@
+#include "source_lcao/LCAO_set.h"
+#include "source_io/module_parameter/parameter.h"
+#include "source_psi/setup_psi.h" // use Setup_Psi
+#include "source_io/read_wfc_nao.h" // use read_wfc_nao
+#include "source_estate/elecstate_tools.h" // use fixed_weights
+
+template <typename TK>
+void LCAO_domain::set_psi_occ_dm_chg(
+		const K_Vectors &kv, // k-points
+		psi::Psi<TK>* &psi, // coefficients of NAO basis
+		const Parallel_Orbitals &pv, // parallel scheme of NAO basis
+		elecstate::ElecState* pelec, // eigen values and weights
+		LCAO_domain::Setup_DM<TK> &dmat, // density matrix 
+		Charge &chr, // charge density 
+		const Input_para &inp) // input parameters
+{
+
+    //! 1) init electronic wave function psi
+    Setup_Psi<TK>::allocate_psi(psi, kv, pv, inp);
+
+    //! 2) read psi from file
+    if (inp.init_wfc == "file" && inp.esolver_type != "tddft")
+    {
+        if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir,
+             pv, *psi, pelec->ekb, pelec->wg, kv.ik2iktot,
+             kv.get_nkstot(), inp.nspin))
+        {
+            ModuleBase::WARNING_QUIT("set_psi_occ_dm_chg", "read electronic wave functions failed");
+        }
+    }
+
+    //! 3) set occupations, tddft does not need to set occupations in the first scf
+    if (inp.ocp && inp.esolver_type != "tddft")
+    {
+        elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec,
+          &kv, pelec->wg, pelec->skip_weights);
+    }
+
+    //! 4) init DMK, but DMR is constructed in before_scf()
+    dmat.allocate_dm(&kv, &pv, inp.nspin);
+
+    //! 5) init charge density
+    chr.allocate(inp.nspin);
+
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "CHARGE");
+
+    return;
+}
+
+
+template <typename TK>
+void LCAO_domain::set_pot(
+        UnitCell &ucell, // not const because of dftu
+		K_Vectors &kv, // not const due to exx 
+	    Structure_Factor& sf, // will be modified in potential	
+		const ModulePW::PW_Basis &pw_rho, 
+		const ModulePW::PW_Basis &pw_rhod, 
+		elecstate::ElecState* pelec,
+		const LCAO_Orbitals& orb,
+		Parallel_Orbitals &pv, // not const due to deepks 
+		pseudopot_cell_vl &locpp, 
+        Plus_U &dftu,
+        surchem& solvent,
+        Exx_NAO<TK> &exx_nao,
+        Setup_DeePKS<TK> &deepks,
+        const Input_para &inp)
+{
+    //! 1) init local pseudopotentials
+    locpp.init_vloc(ucell, &pw_rho);
+
+    //! 2) init potentials
+    if (pelec->pot == nullptr)
+    {
+        // where is the pot deleted?
+        pelec->pot = new elecstate::Potential(&pw_rhod, &pw_rho,
+          &ucell, &locpp.vloc, &sf, &solvent,
+          &(pelec->f_en.etxc), &(pelec->f_en.vtxc));
+    }
+
+    //! 3) initialize DFT+U
+    if (inp.dft_plus_u)
+    {
+        dftu.init(ucell, &pv, kv.get_nks(), &orb);
+    }
+
+    //! 4) init exact exchange calculations
+    exx_nao.before_runner(ucell, kv, orb, pv, inp);
+
+    //! 5) init deepks
+    deepks.before_runner(ucell, kv.get_nks(), orb, pv, inp);
+
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "POTENTIALS");
+
+    return;
+}
+
+
+
+template void LCAO_domain::set_psi_occ_dm_chg<double>(
+		const K_Vectors &kv, // k-points
+		psi::Psi<double>* &psi, // coefficients of NAO basis
+		const Parallel_Orbitals &pv, // parallel scheme of NAO basis
+		elecstate::ElecState* pelec, // eigen values and weights
+		LCAO_domain::Setup_DM<double> &dmat, // density matrix 
+		Charge &chr, // charge density 
+		const Input_para &inp);
+
+template void LCAO_domain::set_psi_occ_dm_chg<std::complex<double>>(
+		const K_Vectors &kv, // k-points
+		psi::Psi<std::complex<double>>* &psi, // coefficients of NAO basis
+		const Parallel_Orbitals &pv, // parallel scheme of NAO basis
+		elecstate::ElecState* pelec, // eigen values and weights
+		LCAO_domain::Setup_DM<std::complex<double>> &dmat, // density matrix 
+		Charge &chr, // charge density 
+		const Input_para &inp);
+
+template void LCAO_domain::set_pot<double>(
+        UnitCell &ucell,
+		K_Vectors &kv, 
+	    Structure_Factor& sf,	
+		const ModulePW::PW_Basis &pw_rho, 
+		const ModulePW::PW_Basis &pw_rhod, 
+		elecstate::ElecState* pelec,
+		const LCAO_Orbitals& orb,
+		Parallel_Orbitals &pv, 
+		pseudopot_cell_vl &locpp, 
+        Plus_U &dftu,
+        surchem& solvent,
+        Exx_NAO<double> &exx_nao,
+        Setup_DeePKS<double> &deepks,
+        const Input_para &inp);
+
+template void LCAO_domain::set_pot<std::complex<double>>(
+        UnitCell &ucell,
+	    K_Vectors &kv, 
+	    Structure_Factor& sf,	
+		const ModulePW::PW_Basis &pw_rho, 
+		const ModulePW::PW_Basis &pw_rhod, 
+		elecstate::ElecState* pelec,
+		const LCAO_Orbitals& orb,
+		Parallel_Orbitals &pv, 
+		pseudopot_cell_vl &locpp, 
+        Plus_U &dftu,
+        surchem& solvent,
+        Exx_NAO<std::complex<double>> &exx_nao,
+        Setup_DeePKS<std::complex<double>> &deepks,
+        const Input_para &inp);

source/source_lcao/LCAO_set.h

@@ -0,0 +1,58 @@
+#ifndef LCAO_SET_H
+#define LCAO_SET_H
+
+#include "source_cell/unitcell.h"
+#include "source_cell/klist.h"
+#include "source_psi/psi.h"
+#include "source_estate/elecstate.h"
+#include "source_lcao/setup_dm.h"
+#include "source_pw/module_pwdft/structure_factor.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_hamilt/module_surchem/surchem.h"
+#include "source_pw/module_pwdft/VL_in_pw.h"
+#include "source_lcao/module_deepks/LCAO_deepks.h"
+#include "source_lcao/module_dftu/dftu.h"
+#include "source_lcao/setup_exx.h"
+#include "source_lcao/setup_deepks.h"
+
+namespace LCAO_domain
+{
+
+/**
+ * @brief set up wave functions, occupation numbers,
+ * density matrix and charge density 
+ */
+template <typename TK>
+void set_psi_occ_dm_chg(
+		const K_Vectors &kv, // k-points
+		psi::Psi<TK>* &psi, // coefficients of NAO basis
+		const Parallel_Orbitals &pv, // parallel scheme of NAO basis
+		elecstate::ElecState* pelec, // eigen values and weights
+		LCAO_domain::Setup_DM<TK> &dmat, // density matrix 
+		Charge &chr, // charge density 
+		const Input_para& inp); // input parameters
+
+/**
+ * @brief set up potentials, including local pseudopotentials,
+ * +U potential, solvent potential, exx potential and deepks potential 
+ */
+template <typename TK>
+void set_pot(
+        UnitCell &ucell,
+		K_Vectors &kv, 
+	    Structure_Factor& sf,	
+		const ModulePW::PW_Basis &pw_rho, 
+		const ModulePW::PW_Basis &pw_rhod, 
+		elecstate::ElecState* pelec,
+		const LCAO_Orbitals& orb,
+		Parallel_Orbitals &pv, 
+		pseudopot_cell_vl &locpp, 
+        Plus_U &dftu,
+        surchem& solvent,
+        Exx_NAO<TK> &exx_nao,
+        Setup_DeePKS<TK> &deepks,
+        const Input_para &inp);
+
+} // end namespace
+
+#endif