Esolver refactor for LCAO and pw

source/source_esolver/esolver_ks_lcao.cpp

@@ -1,6 +1,8 @@
 #include "esolver_ks_lcao.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
+#include "source_lcao/module_deltaspin/deltaspin_lcao.h"
+#include "source_lcao/dftu_lcao.h"
 #include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
@@ -149,13 +151,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     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->dmat.dm, this->pelec);
-    }
+    init_deltaspin_lcao<TK>(ucell, PARAM.inp, &(this->pv), this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec);
 
     // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
     this->exx_nao.before_scf(ucell, this->kv, orb_, this->p_chgmix, istep, PARAM.inp);
@@ -203,11 +199,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 #endif
 
     // 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);
-    }
+    Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
 
     // 17) update of RDMFT, added by jghan
     if (PARAM.inp.rdmft == true)
@@ -347,15 +339,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     }
 #endif
 
-    if (PARAM.inp.dft_plus_u)
-    {
-        if (istep != 0 || iter != 1)
-        {
-            this->dftu.set_dmr(this->dmat.dm);
-        }
-        // Calculate U and J if Yukawa potential is used
-        this->dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx);
-    }
+    init_dftu_lcao<TK>(istep, iter, PARAM.inp, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);
 
 #ifdef __MLALGO
     // the density matrixes of DeePKS have been updated in each iter
@@ -392,24 +376,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
     bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
     // 2) run the inner lambda loop to contrain atomic moments with the DeltaSpin method
-    bool skip_solve = false;
-    if (PARAM.inp.sc_mag_switch)
-    {
-        spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
-        if (!sc.mag_converged() && this->drho > 0 && this->drho < PARAM.inp.sc_scf_thr)
-        {
-            // optimize lambda to get target magnetic moments, but the lambda is not near target
-            sc.run_lambda_loop(iter - 1);
-            sc.set_mag_converged(true);
-            skip_solve = true;
-        }
-        else if (sc.mag_converged())
-        {
-            // optimize lambda to get target magnetic moments, but the lambda is not near target
-            sc.run_lambda_loop(iter - 1);
-            skip_solve = true;
-        }
-    }
+    bool skip_solve = run_deltaspin_lambda_loop_lcao<TK>(iter - 1, this->drho, PARAM.inp);
 
     // 3) run Hsolver
     if (!skip_solve)
@@ -435,11 +402,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
 #endif
 
     // 5) symmetrize the charge density
-    Symmetry_rho srho;
-    for (int is = 0; is < PARAM.inp.nspin; is++)
-    {
-        srho.begin(is, this->chr, this->pw_rho, ucell.symm);
-    }
+    Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
 
     // 6) calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
@@ -461,36 +424,13 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();
 
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
-    if (PARAM.inp.dft_plus_u)
-    {
-        // old DFT+U method calculates energy correction in esolver,
-        // new DFT+U method calculates energy in Hamiltonian
-        if (PARAM.inp.dft_plus_u == 2)
-        {
-            if (this->dftu.omc != 2)
-            {
-                dftu_cal_occup_m(iter, ucell, dm_vec, this->kv,
-                  this->p_chgmix->get_mixing_beta(), hamilt_lcao, this->dftu);
-            }
-            this->dftu.cal_energy_correction(ucell, istep);
-        }
-		this->dftu.output(ucell);
-		// use the converged occupation matrix for next MD/Relax SCF calculation
-		if (conv_esolver)
-		{
-			this->dftu.initialed_locale = true;
-		}
-	}
+    finish_dftu_lcao<TK>(iter, conv_esolver, PARAM.inp, &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao);
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
     this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp);
 
     // 3) for delta spin
-    if (PARAM.inp.sc_mag_switch)
-    {
-        spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
-        sc.cal_mi_lcao(iter);
-    }
+    cal_mi_lcao_wrapper<TK>(iter, PARAM.inp);
 
     // call iter_finish() of ESolver_KS, where band gap is printed,
     // eig and occ are printed, magnetization is calculated,

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -290,11 +290,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
     // Symmetrize the charge density only for ground state
     if (istep <= 1)
     {
-        Symmetry_rho srho;
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            srho.begin(is, this->chr, this->pw_rho, ucell.symm);
-        }
+        Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
     }
 #ifdef __EXX
     if (GlobalC::exx_info.info_ri.real_number)

source/source_esolver/esolver_ks_lcaopw.cpp

@@ -157,11 +157,7 @@ namespace ModuleESolver
         }
 #endif
 
-        Symmetry_rho srho;
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
-        }
+        Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rhod, ucell.symm);
 
         // deband is calculated from "output" charge density calculated
         // in sum_band

source/source_esolver/esolver_ks_pw.cpp

@@ -6,6 +6,7 @@
 
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_hsolver/hsolver_pw.h"
+#include "source_hsolver/diago_params.h"
 
 #include "source_hsolver/kernels/hegvd_op.h"
 #include "source_io/module_parameter/parameter.h"
@@ -48,10 +49,10 @@ ESolver_KS_PW<T, Device>::ESolver_KS_PW()
 template <typename T, typename Device>
 ESolver_KS_PW<T, Device>::~ESolver_KS_PW()
 {
-	//****************************************************
-	// do not add any codes in this deconstructor funcion
-	//****************************************************
-	// delete Hamilt
+    //****************************************************
+    // do not add any codes in this deconstructor funcion
+    //****************************************************
+    // delete Hamilt
     this->deallocate_hamilt();
 
     // mohan add 2025-10-12
@@ -83,7 +84,6 @@ void ESolver_KS_PW<T, Device>::deallocate_hamilt()
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
-    //! Call before_all_runners() of ESolver_KS
     ESolver_KS<T, Device>::before_all_runners(ucell, inp);
 
     //! setup and allocation for pelec, potentials, etc. 
@@ -105,7 +105,6 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     ModuleBase::TITLE("ESolver_KS_PW", "before_scf");
     ModuleBase::timer::tick("ESolver_KS_PW", "before_scf");
 
-    //! Call before_scf() of ESolver_KS
     ESolver_KS<T, Device>::before_scf(ucell, istep);
 
     //! Init variables (once the cell has changed)
@@ -143,17 +142,15 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const int iter)
 {
-    // 1) Call iter_init() of ESolver_KS
     ESolver_KS<T, Device>::iter_init(ucell, istep, iter);
 
-    // 2) perform charge mixing for KSDFT using pw basis
     module_charge::chgmixing_ks_pw(iter, this->p_chgmix, this->dftu, PARAM.inp);
 
-    // 3) mohan move harris functional here, 2012-06-05
+    // mohan move harris functional here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
 
-    // 4) update local occupations for DFT+U
+    // update local occupations for DFT+U
     // should before lambda loop in DeltaSpin
     pw::iter_init_dftu_pw(iter, istep, this->dftu, this->stp.psi_t, this->pelec->wg, ucell, PARAM.inp);
 }
@@ -168,17 +165,8 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
     this->pelec->f_en.eband = 0.0;
     this->pelec->f_en.demet = 0.0;
 
-    // choose if psi should be diag in subspace
-    // be careful that istep start from 0 and iter start from 1
-    // if (iter == 1)
-    hsolver::DiagoIterAssist<T, Device>::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true;
-    hsolver::DiagoIterAssist<T, Device>::SCF_ITER = iter;
-    hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
-
-    if (PARAM.inp.calculation != "nscf")
-    {
-        hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax;
-    }
+    // setup diagonalization parameters
+    hsolver::setup_diago_params_pw<T, Device>(istep, iter, ethr, PARAM.inp);
 
     bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
@@ -205,11 +193,7 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
     }
 
     // symmetrize the charge density
-    Symmetry_rho srho;
-    for (int is = 0; is < PARAM.inp.nspin; is++)
-    {
-        srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
-    }
+    Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rhod, ucell.symm);
 
     ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2rho_single");
 }
@@ -265,7 +249,6 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const
         this->pelec->cal_tau(*(this->psi));
     }
 
-    // Call 'after_scf' of ESolver_KS
     ESolver_KS<T, Device>::after_scf(ucell, istep, conv_esolver);
 
     // Output quantities

source/source_esolver/esolver_of.cpp

@@ -234,11 +234,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
 
     this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
 
-    Symmetry_rho srho;
-    for (int is = 0; is < PARAM.inp.nspin; is++)
-    {
-        srho.begin(is, this->chr, this->pw_rho, ucell.symm);
-    }
+    Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
 
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {

source/source_esolver/esolver_sdft_pw.cpp

@@ -8,6 +8,7 @@
 #include "source_pw/module_stodft/sto_forces.h"
 #include "source_pw/module_stodft/sto_stress_pw.h"
 #include "source_hsolver/diago_iter_assist.h"
+#include "source_hsolver/diago_params.h"
 #include "source_io/module_parameter/parameter.h"
 
 #include <algorithm>
@@ -142,20 +143,11 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
     // reset energy
     this->pelec->f_en.eband = 0.0;
     this->pelec->f_en.demet = 0.0;
-    // choose if psi should be diag in subspace
-    // be careful that istep start from 0 and iter start from 1
-    if (istep == 0 && iter == 1 || PARAM.inp.calculation == "nscf")
-    {
-        hsolver::DiagoIterAssist<T, Device>::need_subspace = false;
-    }
-    else
-    {
-        hsolver::DiagoIterAssist<T, Device>::need_subspace = true;
-    }
+
+    // setup diagonalization parameters for SDFT
+    hsolver::setup_diago_params_sdft<T, Device>(istep, iter, ethr, PARAM.inp);
 
     bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
-    hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR = ethr;
-    hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax;
 
     // hsolver only exists in this function
     hsolver::HSolverPW_SDFT<T, Device> hsolver_pw_sdft_obj(&this->kv,
@@ -190,11 +182,7 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
 
     if (PARAM.globalv.ks_run)
     {
-        Symmetry_rho srho;
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            srho.begin(is, this->chr, this->pw_rho, ucell.symm);
-        }
+        Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
         this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
     }
     else

source/source_estate/module_charge/symmetry_rho.cpp

@@ -10,6 +10,18 @@ Symmetry_rho::~Symmetry_rho()
 {
 }
 
+void Symmetry_rho::symmetrize_rho(const int nspin,
+                                   const Charge& chr,
+                                   const ModulePW::PW_Basis* pw,
+                                   ModuleSymmetry::Symmetry& symm)
+{
+    Symmetry_rho srho;
+    for (int is = 0; is < nspin; is++)
+    {
+        srho.begin(is, chr, pw, symm);
+    }
+}
+
 void Symmetry_rho::begin(const int& spin_now,
                          const Charge& chr,
                          const ModulePW::PW_Basis* rho_basis,

source/source_estate/module_charge/symmetry_rho.h

@@ -11,6 +11,22 @@ class Symmetry_rho
     Symmetry_rho();
     ~Symmetry_rho();
 
+    /**
+     * @brief Symmetrize charge density for all spin channels
+     * 
+     * This is a static helper function that symmetrizes the charge density
+     * for all spin channels by calling begin() for each spin.
+     * 
+     * @param nspin Number of spin channels
+     * @param chr Charge object containing the density
+     * @param pw Plane wave basis
+     * @param symm Symmetry object
+     */
+    static void symmetrize_rho(const int nspin,
+                               const Charge& chr,
+                               const ModulePW::PW_Basis* pw,
+                               ModuleSymmetry::Symmetry& symm);
+
     void begin(const int& spin_now,
                const Charge& CHR,
                const ModulePW::PW_Basis* pw,

source/source_hsolver/CMakeLists.txt

@@ -12,6 +12,7 @@ list(APPEND objects
     hsolver.cpp
     diago_pxxxgvx.cpp
     diag_hs_para.cpp
+    diago_params.cpp
 
 )