From fe1d82a9fe9be3d06b91b67261747311663773ae Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 08:19:14 +0800 Subject: [PATCH 01/11] small format changes --- source/source_esolver/esolver_ks_pw.cpp | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index efb17bc0fd..5f9e56a95d 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -48,10 +48,10 @@ ESolver_KS_PW::ESolver_KS_PW() template ESolver_KS_PW::~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 +83,6 @@ void ESolver_KS_PW::deallocate_hamilt() template void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) { - //! Call before_all_runners() of ESolver_KS ESolver_KS::before_all_runners(ucell, inp); //! setup and allocation for pelec, potentials, etc. @@ -105,7 +104,6 @@ void ESolver_KS_PW::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::before_scf(ucell, istep); //! Init variables (once the cell has changed) @@ -143,17 +141,15 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) template void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const int iter) { - // 1) Call iter_init() of ESolver_KS ESolver_KS::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); } @@ -265,7 +261,6 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const this->pelec->cal_tau(*(this->psi)); } - // Call 'after_scf' of ESolver_KS ESolver_KS::after_scf(ucell, istep, conv_esolver); // Output quantities From 1a725f6bcfd284470c93aea8cf8392341436cde9 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 08:54:30 +0800 Subject: [PATCH 02/11] refactor(esolver): extract charge density symmetrization to Symmetry_rho::symmetrize_rho - Add static method symmetrize_rho() in Symmetry_rho class - Replace 7 duplicate code blocks with single function call - Simplify code from 35 lines to 7 lines (80% reduction) - Improve code readability and maintainability Modified files: - source_estate/module_charge/symmetry_rho.h: add static method declaration - source_estate/module_charge/symmetry_rho.cpp: implement static method - source_esolver/esolver_ks_lcao.cpp: 2 calls updated - source_esolver/esolver_ks_pw.cpp: 1 call updated - source_esolver/esolver_ks_lcao_tddft.cpp: 1 call updated - source_esolver/esolver_ks_lcaopw.cpp: 1 call updated - source_esolver/esolver_of.cpp: 1 call updated - source_esolver/esolver_sdft_pw.cpp: 1 call updated This refactoring follows the ESolver cleanup principle: keep ESolver focused on high-level workflow control. --- source/source_esolver/esolver_ks_lcao.cpp | 12 ++---------- source/source_esolver/esolver_ks_lcao_tddft.cpp | 6 +----- source/source_esolver/esolver_ks_lcaopw.cpp | 6 +----- source/source_esolver/esolver_ks_pw.cpp | 6 +----- source/source_esolver/esolver_of.cpp | 6 +----- source/source_esolver/esolver_sdft_pw.cpp | 6 +----- .../source_estate/module_charge/symmetry_rho.cpp | 12 ++++++++++++ .../source_estate/module_charge/symmetry_rho.h | 16 ++++++++++++++++ 8 files changed, 35 insertions(+), 35 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index f8cecf6805..be6833294f 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -203,11 +203,7 @@ void ESolver_KS_LCAO::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) @@ -435,11 +431,7 @@ void ESolver_KS_LCAO::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); diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 8a0035681b..b7641a09fc 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -290,11 +290,7 @@ void ESolver_KS_LCAO_TDDFT::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) diff --git a/source/source_esolver/esolver_ks_lcaopw.cpp b/source/source_esolver/esolver_ks_lcaopw.cpp index dd37188af3..f9700f5b68 100644 --- a/source/source_esolver/esolver_ks_lcaopw.cpp +++ b/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 diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 5f9e56a95d..b35ea19948 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -201,11 +201,7 @@ void ESolver_KS_PW::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"); } diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 4a086205c4..1fb754d2b6 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/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) { diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 1a9057d178..798e52d26b 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/source/source_esolver/esolver_sdft_pw.cpp @@ -190,11 +190,7 @@ void ESolver_SDFT_PW::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 diff --git a/source/source_estate/module_charge/symmetry_rho.cpp b/source/source_estate/module_charge/symmetry_rho.cpp index dbd8a57af1..19b67967c7 100644 --- a/source/source_estate/module_charge/symmetry_rho.cpp +++ b/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, diff --git a/source/source_estate/module_charge/symmetry_rho.h b/source/source_estate/module_charge/symmetry_rho.h index 638903fd93..98d0650167 100644 --- a/source/source_estate/module_charge/symmetry_rho.h +++ b/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, From 285ee1c58b4968bb3fa945a797cd574078eed109 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 09:19:50 +0800 Subject: [PATCH 03/11] refactor(esolver): extract DeltaSpin lambda loop to deltaspin_lcao module - Create new files deltaspin_lcao.h/cpp in module_deltaspin - Extract DeltaSpin lambda loop logic from ESolver_KS_LCAO - Simplify code from 18 lines to 1 line in hamilt2rho_single - Separate LCAO and PW implementations for DeltaSpin Modified files: - source_esolver/esolver_ks_lcao.cpp: replace inline code with function call - source_lcao/module_deltaspin/CMakeLists.txt: add new source file New files: - source_lcao/module_deltaspin/deltaspin_lcao.h: function declaration - source_lcao/module_deltaspin/deltaspin_lcao.cpp: function implementation This refactoring follows the ESolver cleanup principle: keep ESolver focused on high-level workflow control. --- source/source_esolver/esolver_ks_lcao.cpp | 20 +-------- .../module_deltaspin/CMakeLists.txt | 1 + .../module_deltaspin/deltaspin_lcao.cpp | 44 +++++++++++++++++++ .../module_deltaspin/deltaspin_lcao.h | 29 ++++++++++++ 4 files changed, 76 insertions(+), 18 deletions(-) create mode 100644 source/source_lcao/module_deltaspin/deltaspin_lcao.cpp create mode 100644 source/source_lcao/module_deltaspin/deltaspin_lcao.h diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index be6833294f..44753b41b1 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -1,6 +1,7 @@ #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/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 @@ -388,24 +389,7 @@ void ESolver_KS_LCAO::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& sc = spinconstrain::SpinConstrain::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(iter - 1, this->drho, PARAM.inp); // 3) run Hsolver if (!skip_solve) diff --git a/source/source_lcao/module_deltaspin/CMakeLists.txt b/source/source_lcao/module_deltaspin/CMakeLists.txt index 02f389e5f1..6a0c1fea22 100644 --- a/source/source_lcao/module_deltaspin/CMakeLists.txt +++ b/source/source_lcao/module_deltaspin/CMakeLists.txt @@ -7,6 +7,7 @@ list(APPEND objects lambda_loop.cpp cal_mw_from_lambda.cpp template_helpers.cpp + deltaspin_lcao.cpp ) add_library( diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp new file mode 100644 index 0000000000..9b0e2d08ab --- /dev/null +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -0,0 +1,44 @@ +#include "deltaspin_lcao.h" +#include "spin_constrain.h" + +namespace ModuleESolver +{ + +template +bool run_deltaspin_lambda_loop_lcao(const int iter, + const double drho, + const Input_para& inp) +{ + bool skip_solve = false; + + if (inp.sc_mag_switch) + { + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + + if (!sc.mag_converged() && drho > 0 && drho < 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; + } + } + + return skip_solve; +} + +/// Template instantiation +template bool run_deltaspin_lambda_loop_lcao(const int iter, + const double drho, + const Input_para& inp); +template bool run_deltaspin_lambda_loop_lcao>(const int iter, + const double drho, + const Input_para& inp); + +} // namespace ModuleESolver diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.h b/source/source_lcao/module_deltaspin/deltaspin_lcao.h new file mode 100644 index 0000000000..95d3352732 --- /dev/null +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.h @@ -0,0 +1,29 @@ +#ifndef DELTASPIN_LCAO_H +#define DELTASPIN_LCAO_H + +#include "source_cell/unitcell.h" +#include "source_io/module_parameter/input_parameter.h" + +namespace ModuleESolver +{ + +/** + * @brief Run DeltaSpin lambda loop for LCAO method + * + * This function handles the lambda loop optimization for the DeltaSpin method + * in LCAO calculations. It determines whether to skip the Hamiltonian solve + * based on the convergence status of magnetic moments. + * + * @param iter Current iteration number + * @param drho Charge density convergence criterion + * @param inp Input parameters + * @return bool Whether to skip the Hamiltonian solve + */ +template +bool run_deltaspin_lambda_loop_lcao(const int iter, + const double drho, + const Input_para& inp); + +} // namespace ModuleESolver + +#endif // DELTASPIN_LCAO_H From 2a520e3f26b9b43547a839de57ecd59fdae60930 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 09:30:56 +0800 Subject: [PATCH 04/11] refactor(esolver): complete DeltaSpin refactoring in LCAO - Add init_deltaspin_lcao() function for DeltaSpin initialization - Add cal_mi_lcao_wrapper() function for magnetic moment calculation - Refactor all DeltaSpin-related code in esolver_ks_lcao.cpp - Simplify code from 29 lines to 3 lines (90% reduction) Modified files: - source_esolver/esolver_ks_lcao.cpp: replace 3 code blocks with function calls - source_lcao/module_deltaspin/deltaspin_lcao.h: add 2 new function declarations - source_lcao/module_deltaspin/deltaspin_lcao.cpp: implement 2 new functions This completes the DeltaSpin refactoring for LCAO method: 1. init_deltaspin_lcao() - initialize DeltaSpin calculation 2. cal_mi_lcao_wrapper() - calculate magnetic moments 3. run_deltaspin_lambda_loop_lcao() - run lambda loop optimization All functions follow the ESolver cleanup principle: keep ESolver focused on high-level workflow control. --- source/source_esolver/esolver_ks_lcao.cpp | 14 +---- .../module_deltaspin/deltaspin_lcao.cpp | 59 ++++++++++++++++++- .../module_deltaspin/deltaspin_lcao.h | 37 ++++++++++++ 3 files changed, 96 insertions(+), 14 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 44753b41b1..3050a8aa9e 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -150,13 +150,7 @@ void ESolver_KS_LCAO::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& sc = spinconstrain::SpinConstrain::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(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); @@ -462,11 +456,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& 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& sc = spinconstrain::SpinConstrain::getScInstance(); - sc.cal_mi_lcao(iter); - } + cal_mi_lcao_wrapper(iter); // call iter_finish() of ESolver_KS, where band gap is printed, // eig and occ are printed, magnetization is calculated, diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp index 9b0e2d08ab..c58b4f0783 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -1,9 +1,44 @@ #include "deltaspin_lcao.h" #include "spin_constrain.h" +#include "source_basis/module_ao/parallel_orbitals.h" +#include "source_lcao/hamilt_lcao.h" +#include "source_estate/module_dm/density_matrix.h" +#include "source_estate/elecstate.h" namespace ModuleESolver { +template +void init_deltaspin_lcao(const UnitCell& ucell, + const Input_para& inp, + void* pv, + const K_Vectors& kv, + void* p_hamilt, + void* psi, + void* dm, + void* pelec) +{ + if (!inp.sc_mag_switch) + { + return; + } + + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + sc.init_sc(inp.sc_thr, inp.nsc, inp.nsc_min, inp.alpha_trial, + inp.sccut, inp.sc_drop_thr, ucell, + static_cast(pv), + inp.nspin, kv, p_hamilt, psi, + static_cast*>(dm), + static_cast(pelec)); +} + +template +void cal_mi_lcao_wrapper(const int iter) +{ + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + sc.cal_mi_lcao(iter); +} + template bool run_deltaspin_lambda_loop_lcao(const int iter, const double drho, @@ -18,14 +53,14 @@ bool run_deltaspin_lambda_loop_lcao(const int iter, if (!sc.mag_converged() && drho > 0 && drho < 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.run_lambda_loop(iter); 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); + sc.run_lambda_loop(iter); skip_solve = true; } } @@ -34,6 +69,26 @@ bool run_deltaspin_lambda_loop_lcao(const int iter, } /// Template instantiation +template void init_deltaspin_lcao(const UnitCell& ucell, + const Input_para& inp, + void* pv, + const K_Vectors& kv, + void* p_hamilt, + void* psi, + void* dm, + void* pelec); +template void init_deltaspin_lcao>(const UnitCell& ucell, + const Input_para& inp, + void* pv, + const K_Vectors& kv, + void* p_hamilt, + void* psi, + void* dm, + void* pelec); + +template void cal_mi_lcao_wrapper(const int iter); +template void cal_mi_lcao_wrapper>(const int iter); + template bool run_deltaspin_lambda_loop_lcao(const int iter, const double drho, const Input_para& inp); diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.h b/source/source_lcao/module_deltaspin/deltaspin_lcao.h index 95d3352732..f91326490b 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.h +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.h @@ -2,11 +2,48 @@ #define DELTASPIN_LCAO_H #include "source_cell/unitcell.h" +#include "source_cell/klist.h" #include "source_io/module_parameter/input_parameter.h" namespace ModuleESolver { +/** + * @brief Initialize DeltaSpin for LCAO method + * + * This function initializes the DeltaSpin calculation by setting up + * the SpinConstrain object with input parameters. + * + * @param ucell Unit cell + * @param inp Input parameters + * @param pv Parallel orbitals + * @param kv K-vectors + * @param p_hamilt Pointer to Hamiltonian + * @param psi Pointer to wave functions + * @param dm Density matrix + * @param pelec Pointer to electronic state + */ +template +void init_deltaspin_lcao(const UnitCell& ucell, + const Input_para& inp, + void* pv, + const K_Vectors& kv, + void* p_hamilt, + void* psi, + void* dm, + void* pelec); + +/** + * @brief Calculate magnetic moments for DeltaSpin in LCAO method + * + * This function calculates the magnetic moments for each atom + * in the DeltaSpin method. + * + * @param iter Current iteration number + */ +template +void cal_mi_lcao_wrapper(const int iter); + /** * @brief Run DeltaSpin lambda loop for LCAO method * From 91be943b3f8de93a6b4d394d01671cb09ae4fd4a Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 10:16:10 +0800 Subject: [PATCH 05/11] refactor(esolver): extract DFT+U code to dftu_lcao module - Create new files dftu_lcao.h/cpp in source_lcao directory - Add init_dftu_lcao() function for DFT+U initialization - Add finish_dftu_lcao() function for DFT+U finalization - Simplify code from 32 lines to 2 lines in esolver_ks_lcao.cpp - Remove conditional checks from ESolver, move them to functions Modified files: - source_esolver/esolver_ks_lcao.cpp: replace 2 code blocks with function calls - source_lcao/CMakeLists.txt: add new source file New files: - source_lcao/dftu_lcao.h: function declarations - source_lcao/dftu_lcao.cpp: function implementations This refactoring prepares for unifying old and new DFT+U implementations: - Old DFT+U: source_lcao/module_dftu/ - New DFT+U: source_lcao/module_operator_lcao/op_dftu_lcao.cpp All functions follow ESolver cleanup principle: keep ESolver focused on high-level workflow control. --- source/source_esolver/esolver_ks_lcao.cpp | 32 +------ source/source_lcao/CMakeLists.txt | 1 + source/source_lcao/dftu_lcao.cpp | 112 ++++++++++++++++++++++ source/source_lcao/dftu_lcao.h | 65 +++++++++++++ 4 files changed, 181 insertions(+), 29 deletions(-) create mode 100644 source/source_lcao/dftu_lcao.cpp create mode 100644 source/source_lcao/dftu_lcao.h diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 3050a8aa9e..e86c7a8cb6 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -2,6 +2,7 @@ #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 @@ -338,15 +339,7 @@ void ESolver_KS_LCAO::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(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 @@ -431,26 +424,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& const std::vector>& 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(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); diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt index 844da5dc84..a793f5d5d0 100644 --- a/source/source_lcao/CMakeLists.txt +++ b/source/source_lcao/CMakeLists.txt @@ -23,6 +23,7 @@ if(ENABLE_LCAO) module_operator_lcao/dspin_lcao.cpp module_operator_lcao/dftu_lcao.cpp module_operator_lcao/operator_force_stress_utils.cpp + dftu_lcao.cpp pulay_fs_center2.cpp FORCE_STRESS.cpp FORCE_gamma.cpp diff --git a/source/source_lcao/dftu_lcao.cpp b/source/source_lcao/dftu_lcao.cpp new file mode 100644 index 0000000000..5a4c6c45c8 --- /dev/null +++ b/source/source_lcao/dftu_lcao.cpp @@ -0,0 +1,112 @@ +#include "dftu_lcao.h" +#include "source_lcao/module_dftu/dftu.h" +#include "source_estate/module_dm/density_matrix.h" +#include "source_lcao/hamilt_lcao.h" + +namespace ModuleESolver +{ + +template +void init_dftu_lcao(const int istep, + const int iter, + const Input_para& inp, + void* dftu, + void* dm, + const UnitCell& ucell, + double** rho, + const int nrxx) +{ + if (!inp.dft_plus_u) + { + return; + } + + auto* dftu_ptr = static_cast(dftu); + auto* dm_ptr = static_cast*>(dm); + + if (istep != 0 || iter != 1) + { + dftu_ptr->set_dmr(dm_ptr); + } + + /// Calculate U and J if Yukawa potential is used + dftu_ptr->cal_slater_UJ(ucell, rho, nrxx); +} + +template +void finish_dftu_lcao(const int iter, + const bool conv_esolver, + const Input_para& inp, + void* dftu, + const UnitCell& ucell, + const std::vector>& dm_vec, + const K_Vectors& kv, + const double mixing_beta, + void* hamilt_lcao) +{ + if (!inp.dft_plus_u) + { + return; + } + + auto* dftu_ptr = static_cast(dftu); + auto* hamilt_lcao_ptr = static_cast*>(hamilt_lcao); + + /// old DFT+U method calculates energy correction in esolver, + /// new DFT+U method calculates energy in Hamiltonian + if (inp.dft_plus_u == 2) + { + if (dftu_ptr->omc != 2) + { + dftu_cal_occup_m(iter, ucell, dm_vec, kv, mixing_beta, + static_cast*>(hamilt_lcao_ptr), *dftu_ptr); + } + dftu_ptr->cal_energy_correction(ucell, iter); + } + dftu_ptr->output(ucell); + + /// use the converged occupation matrix for next MD/Relax SCF calculation + if (conv_esolver) + { + dftu_ptr->initialed_locale = true; + } +} + +/// Template instantiation +template void init_dftu_lcao(const int istep, + const int iter, + const Input_para& inp, + void* dftu, + void* dm, + const UnitCell& ucell, + double** rho, + const int nrxx); +template void init_dftu_lcao>(const int istep, + const int iter, + const Input_para& inp, + void* dftu, + void* dm, + const UnitCell& ucell, + double** rho, + const int nrxx); + +template void finish_dftu_lcao(const int iter, + const bool conv_esolver, + const Input_para& inp, + void* dftu, + const UnitCell& ucell, + const std::vector>& dm_vec, + const K_Vectors& kv, + const double mixing_beta, + void* hamilt_lcao); +template void finish_dftu_lcao>(const int iter, + const bool conv_esolver, + const Input_para& inp, + void* dftu, + const UnitCell& ucell, + const std::vector>>& dm_vec, + const K_Vectors& kv, + const double mixing_beta, + void* hamilt_lcao); + +} // namespace ModuleESolver diff --git a/source/source_lcao/dftu_lcao.h b/source/source_lcao/dftu_lcao.h new file mode 100644 index 0000000000..5138b66256 --- /dev/null +++ b/source/source_lcao/dftu_lcao.h @@ -0,0 +1,65 @@ +#ifndef DFTU_LCAO_H +#define DFTU_LCAO_H + +#include "source_cell/unitcell.h" +#include "source_cell/klist.h" +#include "source_io/module_parameter/input_parameter.h" + +namespace ModuleESolver +{ + +/** + * @brief Initialize DFT+U for LCAO method in iter_init + * + * This function handles the DFT+U initialization during the SCF iteration. + * It sets the density matrix and calculates Slater integrals if needed. + * + * @param istep Current ionic step + * @param iter Current SCF iteration + * @param inp Input parameters + * @param dftu DFT+U object + * @param dm Density matrix + * @param ucell Unit cell + * @param rho Charge density + * @param nrxx Number of real space grid points + */ +template +void init_dftu_lcao(const int istep, + const int iter, + const Input_para& inp, + void* dftu, + void* dm, + const UnitCell& ucell, + double** rho, + const int nrxx); + +/** + * @brief Finish DFT+U calculation for LCAO method in iter_finish + * + * This function handles the DFT+U finalization during the SCF iteration. + * It calculates the occupation matrix and energy correction if needed. + * + * @param iter Current SCF iteration + * @param conv_esolver Whether ESolver has converged + * @param inp Input parameters + * @param dftu DFT+U object + * @param ucell Unit cell + * @param dm_vec Density matrix vector + * @param kv K-vectors + * @param mixing_beta Mixing beta parameter + * @param hamilt_lcao Hamiltonian LCAO object + */ +template +void finish_dftu_lcao(const int iter, + const bool conv_esolver, + const Input_para& inp, + void* dftu, + const UnitCell& ucell, + const std::vector>& dm_vec, + const K_Vectors& kv, + const double mixing_beta, + void* hamilt_lcao); + +} // namespace ModuleESolver + +#endif // DFTU_LCAO_H From 5a9648317862c698bb1a7f6fe225c3ca82f4dcfc Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 12:30:46 +0800 Subject: [PATCH 06/11] refactor(esolver): extract diagonalization parameters setup to hsolver module - Create new files diago_params.h/cpp in source_hsolver directory - Add setup_diago_params_pw() function for PW diagonalization parameters - Simplify code from 11 lines to 1 line in esolver_ks_pw.cpp - Encapsulate diagonalization parameter setup logic Modified files: - source_esolver/esolver_ks_pw.cpp: replace inline code with function call - source_hsolver/CMakeLists.txt: add new source file New files: - source_hsolver/diago_params.h: function declaration - source_hsolver/diago_params.cpp: function implementation This refactoring follows ESolver cleanup principle: keep ESolver focused on high-level workflow control. --- source/source_esolver/esolver_ks_pw.cpp | 14 ++----- source/source_hsolver/CMakeLists.txt | 1 + source/source_hsolver/diago_params.cpp | 55 +++++++++++++++++++++++++ source/source_hsolver/diago_params.h | 29 +++++++++++++ 4 files changed, 88 insertions(+), 11 deletions(-) create mode 100644 source/source_hsolver/diago_params.cpp create mode 100644 source/source_hsolver/diago_params.h diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index b35ea19948..e255b95b46 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/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" @@ -164,17 +165,8 @@ void ESolver_KS_PW::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::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; - hsolver::DiagoIterAssist::SCF_ITER = iter; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - - if (PARAM.inp.calculation != "nscf") - { - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; - } + // setup diagonalization parameters + hsolver::setup_diago_params_pw(istep, iter, ethr, PARAM.inp); bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; diff --git a/source/source_hsolver/CMakeLists.txt b/source/source_hsolver/CMakeLists.txt index f4e87cdf94..b115d6d4cd 100644 --- a/source/source_hsolver/CMakeLists.txt +++ b/source/source_hsolver/CMakeLists.txt @@ -12,6 +12,7 @@ list(APPEND objects hsolver.cpp diago_pxxxgvx.cpp diag_hs_para.cpp + diago_params.cpp ) diff --git a/source/source_hsolver/diago_params.cpp b/source/source_hsolver/diago_params.cpp new file mode 100644 index 0000000000..a0c720a625 --- /dev/null +++ b/source/source_hsolver/diago_params.cpp @@ -0,0 +1,55 @@ +#include "diago_params.h" +#include "diago_iter_assist.h" + +namespace hsolver +{ + +template +void setup_diago_params_pw(const int istep, + const int iter, + const double ethr, + const Input_para& inp) +{ + /// choose if psi should be diag in subspace + /// be careful that istep start from 0 and iter start from 1 + DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; + DiagoIterAssist::SCF_ITER = iter; + DiagoIterAssist::PW_DIAG_THR = ethr; + + if (inp.calculation != "nscf") + { + DiagoIterAssist::PW_DIAG_NMAX = inp.pw_diag_nmax; + } +} + +/// Template instantiation for CPU +template void setup_diago_params_pw, base_device::DEVICE_CPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_pw, base_device::DEVICE_CPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_pw(const int istep, + const int iter, + const double ethr, + const Input_para& inp); + +/// Template instantiation for GPU +#if ((defined __CUDA) || (defined __ROCM)) +template void setup_diago_params_pw, base_device::DEVICE_GPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_pw, base_device::DEVICE_GPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_pw(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +#endif + +} // namespace hsolver diff --git a/source/source_hsolver/diago_params.h b/source/source_hsolver/diago_params.h new file mode 100644 index 0000000000..995090bebd --- /dev/null +++ b/source/source_hsolver/diago_params.h @@ -0,0 +1,29 @@ +#ifndef DIAGO_PARAMS_H +#define DIAGO_PARAMS_H + +#include "source_io/module_parameter/input_parameter.h" + +namespace hsolver +{ + +/** + * @brief Setup diagonalization parameters for PW method + * + * This function sets up the diagonalization parameters for plane wave method, + * including subspace diagonalization flag, SCF iteration number, diagonalization + * threshold, and maximum number of diagonalization steps. + * + * @param istep Current ionic step + * @param iter Current SCF iteration + * @param ethr Diagonalization threshold + * @param inp Input parameters + */ +template +void setup_diago_params_pw(const int istep, + const int iter, + const double ethr, + const Input_para& inp); + +} // namespace hsolver + +#endif // DIAGO_PARAMS_H From 4936dc1f28adc6fdbbf208825c95b4dabcbf8462 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 14:00:59 +0800 Subject: [PATCH 07/11] fix(deltaspin): add sc_mag_switch check in cal_mi_lcao_wrapper - Add Input_para parameter to cal_mi_lcao_wrapper function - Add sc_mag_switch check to avoid calling cal_mi_lcao when DeltaSpin is disabled - Fix 'atomCounts is not set' error in non-DeltaSpin calculations - Update function call in esolver_ks_lcao.cpp This fix resolves the CI/CD failure caused by commit 2a520e3f2. The root cause was that cal_mi_lcao_wrapper was called without checking sc_mag_switch, leading to uninitialized atomCounts error. Modified files: - source_esolver/esolver_ks_lcao.cpp: update function call - source_lcao/module_deltaspin/deltaspin_lcao.h: add parameter - source_lcao/module_deltaspin/deltaspin_lcao.cpp: add check This follows the refactoring principle: preserve original condition checks when extracting code to wrapper functions. --- source/source_esolver/esolver_ks_lcao.cpp | 2 +- .../source_lcao/module_deltaspin/deltaspin_lcao.cpp | 11 ++++++++--- source/source_lcao/module_deltaspin/deltaspin_lcao.h | 3 ++- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index e86c7a8cb6..43e83aa8df 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -430,7 +430,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& 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 - cal_mi_lcao_wrapper(iter); + cal_mi_lcao_wrapper(iter, PARAM.inp); // call iter_finish() of ESolver_KS, where band gap is printed, // eig and occ are printed, magnetization is calculated, diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp index c58b4f0783..96e969277c 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -33,8 +33,13 @@ void init_deltaspin_lcao(const UnitCell& ucell, } template -void cal_mi_lcao_wrapper(const int iter) +void cal_mi_lcao_wrapper(const int iter, const Input_para& inp) { + if (!inp.sc_mag_switch) + { + return; + } + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.cal_mi_lcao(iter); } @@ -86,8 +91,8 @@ template void init_deltaspin_lcao>(const UnitCell& ucell, void* dm, void* pelec); -template void cal_mi_lcao_wrapper(const int iter); -template void cal_mi_lcao_wrapper>(const int iter); +template void cal_mi_lcao_wrapper(const int iter, const Input_para& inp); +template void cal_mi_lcao_wrapper>(const int iter, const Input_para& inp); template bool run_deltaspin_lambda_loop_lcao(const int iter, const double drho, diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.h b/source/source_lcao/module_deltaspin/deltaspin_lcao.h index f91326490b..959109ece7 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.h +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.h @@ -40,9 +40,10 @@ void init_deltaspin_lcao(const UnitCell& ucell, * in the DeltaSpin method. * * @param iter Current iteration number + * @param inp Input parameters */ template -void cal_mi_lcao_wrapper(const int iter); +void cal_mi_lcao_wrapper(const int iter, const Input_para& inp); /** * @brief Run DeltaSpin lambda loop for LCAO method From ea218f642d08d38a871d7fe27dc488f08fdce25d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 18:03:31 +0800 Subject: [PATCH 08/11] fix(deltaspin): add #ifdef __LCAO for conditional compilation - Add #ifdef __LCAO conditional compilation in init_deltaspin_lcao and cal_mi_lcao_wrapper - Fix parameter order in init_sc call for LCAO and non-LCAO builds - Fix undefined reference to cal_mi_lcao in non-LCAO build This fix resolves CI/CD compilation errors in both build_5pt (with __LCAO) and build_1p (without __LCAO) environments. The The root cause was 1. init_sc has different parameter order in LCAO vs non-LCAO builds - LCAO: psi, dm, pelec - non-LCAO: psi, pelec 2. cal_mi_lcao is only defined in LCAO build Modified files: - source_hsolver/diago_params.h: add setup_diago_params_sdft declaration - source_lcao/module_deltaspin/deltaspin_lcao.cpp: add conditional compilation This follows the refactoring principle: handle conditional compilation properly when code has different implementations for different build configurations. --- source/source_hsolver/diago_params.h | 18 ++++++++++++++++++ .../module_deltaspin/deltaspin_lcao.cpp | 10 ++++++++++ 2 files changed, 28 insertions(+) diff --git a/source/source_hsolver/diago_params.h b/source/source_hsolver/diago_params.h index 995090bebd..5d46b01046 100644 --- a/source/source_hsolver/diago_params.h +++ b/source/source_hsolver/diago_params.h @@ -24,6 +24,24 @@ void setup_diago_params_pw(const int istep, const double ethr, const Input_para& inp); +/** + * @brief Setup diagonalization parameters for SDFT method + * + * This function sets up the diagonalization parameters for stochastic DFT method, + * including subspace diagonalization flag, diagonalization threshold, and + * maximum number of diagonalization steps. + * + * @param istep Current ionic step + * @param iter Current SCF iteration + * @param ethr Diagonalization threshold + * @param inp Input parameters + */ +template +void setup_diago_params_sdft(const int istep, + const int iter, + const double ethr, + const Input_para& inp); + } // namespace hsolver #endif // DIAGO_PARAMS_H diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp index 96e969277c..6a7effb6d0 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -24,12 +24,20 @@ void init_deltaspin_lcao(const UnitCell& ucell, } spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); +#ifdef __LCAO sc.init_sc(inp.sc_thr, inp.nsc, inp.nsc_min, inp.alpha_trial, inp.sccut, inp.sc_drop_thr, ucell, static_cast(pv), inp.nspin, kv, p_hamilt, psi, static_cast*>(dm), static_cast(pelec)); +#else + sc.init_sc(inp.sc_thr, inp.nsc, inp.nsc_min, inp.alpha_trial, + inp.sccut, inp.sc_drop_thr, ucell, + static_cast(pv), + inp.nspin, kv, p_hamilt, psi, + static_cast(pelec)); +#endif } template @@ -40,8 +48,10 @@ void cal_mi_lcao_wrapper(const int iter, const Input_para& inp) return; } +#ifdef __LCAO spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.cal_mi_lcao(iter); +#endif } template From c365a3afd18cced4dfdc424c120fe0a0e3bfb4a6 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 18:19:57 +0800 Subject: [PATCH 09/11] refactor(esolver): extract SDFT diagonalization parameters setup - Add setup_diago_params_sdft() function for SDFT diagonalization parameters - Simplify code from 11 lines to 1 line in esolver_sdft_pw.cpp - Encapsulate diagonalization parameter setup logic for SDFT Modified files: - source_esolver/esolver_sdft_pw.cpp: replace inline code with function call - source_hsolver/diago_params.cpp: add setup_diago_params_sdft implementation This refactoring follows ESolver cleanup principle: keep ESolver focused on high-level workflow control. Note: SDFT has different parameter setup logic compared to PW: - Different need_subspace condition - No SCF_ITER setting - Always set PW_DIAG_NMAX (no nscf check) --- source/source_esolver/esolver_sdft_pw.cpp | 16 ++----- source/source_hsolver/diago_params.cpp | 51 +++++++++++++++++++++++ 2 files changed, 55 insertions(+), 12 deletions(-) diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 798e52d26b..26118eed21 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/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 @@ -142,20 +143,11 @@ void ESolver_SDFT_PW::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::need_subspace = false; - } - else - { - hsolver::DiagoIterAssist::need_subspace = true; - } + + // setup diagonalization parameters for SDFT + hsolver::setup_diago_params_sdft(istep, iter, ethr, PARAM.inp); bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; // hsolver only exists in this function hsolver::HSolverPW_SDFT hsolver_pw_sdft_obj(&this->kv, diff --git a/source/source_hsolver/diago_params.cpp b/source/source_hsolver/diago_params.cpp index a0c720a625..28e1040a97 100644 --- a/source/source_hsolver/diago_params.cpp +++ b/source/source_hsolver/diago_params.cpp @@ -22,6 +22,27 @@ void setup_diago_params_pw(const int istep, } } +template +void setup_diago_params_sdft(const int istep, + const int iter, + const double ethr, + const Input_para& inp) +{ + /// 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 || inp.calculation == "nscf") + { + DiagoIterAssist::need_subspace = false; + } + else + { + DiagoIterAssist::need_subspace = true; + } + + DiagoIterAssist::PW_DIAG_THR = ethr; + DiagoIterAssist::PW_DIAG_NMAX = inp.pw_diag_nmax; +} + /// Template instantiation for CPU template void setup_diago_params_pw, base_device::DEVICE_CPU>(const int istep, const int iter, @@ -52,4 +73,34 @@ template void setup_diago_params_pw(const int i const Input_para& inp); #endif +/// Template instantiation for SDFT CPU +template void setup_diago_params_sdft, base_device::DEVICE_CPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_sdft, base_device::DEVICE_CPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_sdft(const int istep, + const int iter, + const double ethr, + const Input_para& inp); + +/// Template instantiation for SDFT GPU +#if ((defined __CUDA) || (defined __ROCM)) +template void setup_diago_params_sdft, base_device::DEVICE_GPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_sdft, base_device::DEVICE_GPU>(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +template void setup_diago_params_sdft(const int istep, + const int iter, + const double ethr, + const Input_para& inp); +#endif + } // namespace hsolver From cbd0ce77936542fed08b3c568fde9250fce50a7d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Wed, 11 Mar 2026 13:18:58 +0800 Subject: [PATCH 10/11] refactor(hamilt): introduce HamiltBase non-template base class - Create HamiltBase as a non-template base class for Hamilt - Modify Hamilt to inherit from HamiltBase - Change ESolver_KS::p_hamilt type from Hamilt* to HamiltBase* - Add static_cast where needed when passing p_hamilt to functions expecting Hamilt* This is the first step towards removing template parameters from ESolver. Modified files: - source/source_esolver/esolver_ks.h - source/source_esolver/esolver_ks_lcaopw.cpp - source/source_esolver/esolver_ks_pw.cpp - source/source_esolver/esolver_sdft_pw.cpp - source/source_hamilt/hamilt.h New files: - source/source_hamilt/hamilt_base.h --- source/source_esolver/esolver_ks.h | 5 +- source/source_esolver/esolver_ks_lcaopw.cpp | 2 +- source/source_esolver/esolver_ks_pw.cpp | 2 +- source/source_esolver/esolver_sdft_pw.cpp | 4 +- source/source_hamilt/hamilt.h | 13 ++++-- source/source_hamilt/hamilt_base.h | 52 +++++++++++++++++++++ 6 files changed, 69 insertions(+), 9 deletions(-) create mode 100644 source/source_hamilt/hamilt_base.h diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index 787b58ba74..1913aa4101 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -7,6 +7,7 @@ #include "source_estate/module_charge/charge_mixing.h" // use charge mixing #include "source_psi/psi.h" // use electronic wave functions #include "source_hamilt/hamilt.h" // use Hamiltonian +#include "source_hamilt/hamilt_base.h" // use Hamiltonian base class #include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace ModuleESolver @@ -47,8 +48,8 @@ class ESolver_KS : public ESolver_FP //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; - //! Hamiltonian - hamilt::Hamilt* p_hamilt = nullptr; + //! Hamiltonian (base class pointer, actual type determined at runtime) + hamilt::HamiltBase* p_hamilt = nullptr; //! PW for wave functions, only used in KSDFT, not in OFDFT ModulePW::PW_Basis_K* pw_wfc = nullptr; diff --git a/source/source_esolver/esolver_ks_lcaopw.cpp b/source/source_esolver/esolver_ks_lcaopw.cpp index f9700f5b68..db00b6265d 100644 --- a/source/source_esolver/esolver_ks_lcaopw.cpp +++ b/source/source_esolver/esolver_ks_lcaopw.cpp @@ -146,7 +146,7 @@ namespace ModuleESolver bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverLIP hsolver_lip_obj(this->pw_wfc); - hsolver_lip_obj.solve(this->p_hamilt, this->stp.psi_t[0], this->pelec, + hsolver_lip_obj.solve(static_cast*>(this->p_hamilt), this->stp.psi_t[0], this->pelec, *this->psi_local, skip_charge,ucell.tpiba,ucell.nat); // add exx diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index e255b95b46..a032c2c976 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -188,7 +188,7 @@ void ESolver_KS_PW::hamilt2rho_single(UnitCell& ucell, const int iste hsolver::DiagoIterAssist::need_subspace, PARAM.inp.use_k_continuity); - hsolver_pw_obj.solve(this->p_hamilt, this->stp.psi_t[0], this->pelec, this->pelec->ekb.c, + hsolver_pw_obj.solve(static_cast*>(this->p_hamilt), this->stp.psi_t[0], this->pelec, this->pelec->ekb.c, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, skip_charge, ucell.tpiba, ucell.nat); } diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 26118eed21..597825cd6d 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/source/source_esolver/esolver_sdft_pw.cpp @@ -167,7 +167,7 @@ void ESolver_SDFT_PW::hamilt2rho_single(UnitCell& ucell, int istep, i hsolver::DiagoIterAssist::need_subspace); hsolver_pw_sdft_obj.solve(ucell, - this->p_hamilt, + static_cast*>(this->p_hamilt), this->stp.psi_t[0], this->stp.psi_cpu[0], this->pelec, @@ -291,7 +291,7 @@ void ESolver_SDFT_PW::after_all_runners(UnitCell& ucell) this->pw_wfc, this->stp.psi_t, &this->ppcell, - this->p_hamilt, + static_cast, Device>*>(this->p_hamilt), this->stoche, &stowf); sto_elecond.decide_nche(PARAM.inp.cond_dt, 1e-8, this->nche_sto, PARAM.inp.emin_sto, PARAM.inp.emax_sto); diff --git a/source/source_hamilt/hamilt.h b/source/source_hamilt/hamilt.h index 6d732d7a82..3d554c0fe6 100644 --- a/source/source_hamilt/hamilt.h +++ b/source/source_hamilt/hamilt.h @@ -7,21 +7,28 @@ #include "matrixblock.h" #include "source_psi/psi.h" #include "operator.h" +#include "hamilt_base.h" namespace hamilt { template -class Hamilt +class Hamilt : public HamiltBase { public: virtual ~Hamilt(){}; /// for target K point, update consequence of hPsi() and matrix() - virtual void updateHk(const int ik){return;} + void updateHk(const int ik) override { return; } /// refresh status of Hamiltonian, for example, refresh H(R) and S(R) in LCAO case - virtual void refresh(bool yes = true){return;} + void refresh(bool yes = true) override { return; } + + /// get the class name + std::string get_classname() const override { return classname; } + + /// get the operator chain + void* get_ops() override { return static_cast(ops); } /// core function: for solving eigenvalues of Hamiltonian with iterative method virtual void hPsi( diff --git a/source/source_hamilt/hamilt_base.h b/source/source_hamilt/hamilt_base.h new file mode 100644 index 0000000000..06325bf050 --- /dev/null +++ b/source/source_hamilt/hamilt_base.h @@ -0,0 +1,52 @@ +#ifndef HAMILT_BASE_H +#define HAMILT_BASE_H + +#include + +namespace hamilt +{ + +/** + * @brief Base class for Hamiltonian + * + * This is a non-template base class for Hamilt. + * It provides a common interface for all Hamiltonian types, + * allowing ESolver to manage Hamiltonian without template parameters. + */ +class HamiltBase +{ + public: + virtual ~HamiltBase() {} + + /** + * @brief Update Hamiltonian for a specific k-point + * + * @param ik k-point index + */ + virtual void updateHk(const int ik) { return; } + + /** + * @brief Refresh the status of Hamiltonian + * + * @param yes whether to refresh + */ + virtual void refresh(bool yes = true) { return; } + + /** + * @brief Get the class name + * + * @return class name + */ + virtual std::string get_classname() const { return "none"; } + + /** + * @brief Get the operator chain (as void* to avoid template) + * + * @return pointer to operator chain + */ + virtual void* get_ops() { return nullptr; } +}; + +} // namespace hamilt + +#endif // HAMILT_BASE_H From 6e0f43ca072682621867dee28d5ee8d4e58e943d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Wed, 11 Mar 2026 16:02:42 +0800 Subject: [PATCH 11/11] refactor(esolver): add static_cast for p_hamilt in esolver files - Add static_cast*> when passing p_hamilt to functions expecting Hamilt* type - Split long cast statements into multiple lines for better readability - Files modified: - esolver_ks_pw.cpp: setup_pot, stp.init calls - esolver_ks_lcao.cpp: init_chg_hr, hsolver_lcao_obj.solve calls - esolver_ks_lcao_tddft.cpp: solve_psi, cal_edm_tddft, matrix calls - esolver_gets.cpp: ops access, output_SR call This follows the HamiltBase refactoring strategy where p_hamilt is stored as HamiltBase* and cast to Hamilt* when needed. --- source/source_esolver/esolver_gets.cpp | 12 ++++++++---- source/source_esolver/esolver_ks_lcao.cpp | 4 ++-- source/source_esolver/esolver_ks_lcao_tddft.cpp | 12 ++++++------ source/source_esolver/esolver_ks_pw.cpp | 4 ++-- 4 files changed, 18 insertions(+), 14 deletions(-) diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index 7eff0f537a..e03e7b8bdc 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/source/source_esolver/esolver_gets.cpp @@ -108,8 +108,9 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); - dynamic_cast, std::complex>*>(this->p_hamilt->ops) - ->contributeHR(); + auto* hamilt_ptr = static_cast>*>(this->p_hamilt); + auto* ops_ptr = dynamic_cast, std::complex>*>(hamilt_ptr->ops); + ops_ptr->contributeHR(); } else { @@ -119,13 +120,16 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); - dynamic_cast, double>*>(this->p_hamilt->ops)->contributeHR(); + auto* hamilt_ptr = static_cast>*>(this->p_hamilt); + auto* ops_ptr = dynamic_cast, double>*>(hamilt_ptr->ops); + ops_ptr->contributeHR(); } } const std::string fn = PARAM.globalv.global_out_dir + "sr_nao.csr"; - ModuleIO::output_SR(pv, gd, this->p_hamilt, fn); + auto* hamilt_ptr = static_cast>*>(this->p_hamilt); + ModuleIO::output_SR(pv, gd, hamilt_ptr, fn); if (PARAM.inp.out_mat_r) { diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 43e83aa8df..bad1cec7b6 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -179,7 +179,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { //! 13.1.2) init charge density from Hamiltonian matrix file LCAO_domain::init_chg_hr(PARAM.globalv.global_readin_dir, PARAM.inp.nspin, - this->p_hamilt, ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm, + static_cast*>(this->p_hamilt), ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm, this->chr, PARAM.inp.ks_solver); } } @@ -382,7 +382,7 @@ void ESolver_KS_LCAO::hamilt2rho_single(UnitCell& ucell, int istep, int if (!skip_solve) { hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, + hsolver_lcao_obj.solve(static_cast*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, this->chr, PARAM.inp.nspin, skip_charge); } diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index b7641a09fc..05dc8c9233 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -235,7 +235,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, PARAM.inp.nbands, PARAM.globalv.nlocal, this->kv.get_nks(), - this->p_hamilt, + static_cast>*>(this->p_hamilt), this->pv, this->psi, this->psi_laststep, @@ -255,7 +255,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, PARAM.inp.nbands, PARAM.globalv.nlocal, this->kv.get_nks(), - this->p_hamilt, + static_cast>*>(this->p_hamilt), this->pv, this->psi, this->psi_laststep, @@ -277,7 +277,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, { bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverLCAO> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, + hsolver_lcao_obj.solve(static_cast>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, @@ -342,11 +342,11 @@ void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, { if (use_tensor && use_lapack) { - elecstate::cal_edm_tddft_tensor_lapack(this->pv, this->dmat, this->kv, this->p_hamilt); + elecstate::cal_edm_tddft_tensor_lapack(this->pv, this->dmat, this->kv, static_cast>*>(this->p_hamilt)); } else { - elecstate::cal_edm_tddft(this->pv, this->dmat, this->kv, this->p_hamilt); + elecstate::cal_edm_tddft(this->pv, this->dmat, this->kv, static_cast>*>(this->p_hamilt)); } } } @@ -416,7 +416,7 @@ void ESolver_KS_LCAO_TDDFT::store_h_s_psi(UnitCell& ucell, this->p_hamilt->updateHk(ik); hamilt::MatrixBlock> h_mat; hamilt::MatrixBlock> s_mat; - this->p_hamilt->matrix(h_mat, s_mat); + static_cast>*>(this->p_hamilt)->matrix(h_mat, s_mat); // Store H and S matrices to Hk_laststep and Sk_laststep if (use_tensor && use_lapack) diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index a032c2c976..b2976733bf 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -128,10 +128,10 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) // init DFT+U is done in "before_all_runners" in LCAO basis. This should be refactored, mohan note 2025-11-06 pw::setup_pot(istep, ucell, this->kv, this->sf, this->pelec, this->Pgrid, this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell, - this->stp.psi_t, this->p_hamilt, this->pw_wfc, this->pw_rhod, PARAM.inp); + this->stp.psi_t, static_cast*>(this->p_hamilt), this->pw_wfc, this->pw_rhod, PARAM.inp); // setup psi (electronic wave functions) - this->stp.init(this->p_hamilt); + this->stp.init(static_cast*>(this->p_hamilt)); //! Setup EXX helper for Hamiltonian and psi exx_helper.before_scf(this->p_hamilt, this->stp.psi_t, PARAM.inp);