From fe1d82a9fe9be3d06b91b67261747311663773ae Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 10 Mar 2026 08:19:14 +0800 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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