diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index f8cecf6805..43e83aa8df 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -1,6 +1,8 @@ #include "esolver_ks_lcao.h" #include "source_estate/elecstate_tools.h" #include "source_lcao/module_deltaspin/spin_constrain.h" +#include "source_lcao/module_deltaspin/deltaspin_lcao.h" +#include "source_lcao/dftu_lcao.h" #include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp #include "source_estate/module_charge/symmetry_rho.h" #include "source_lcao/LCAO_domain.h" // need DeePKS_init @@ -149,13 +151,7 @@ void ESolver_KS_LCAO::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); @@ -203,11 +199,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) @@ -347,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 @@ -392,24 +376,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) @@ -435,11 +402,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); @@ -461,36 +424,13 @@ 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); // 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, 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_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 efb17bc0fd..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" @@ -48,10 +49,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 +84,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 +105,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 +142,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); } @@ -168,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; @@ -205,11 +193,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"); } @@ -265,7 +249,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 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..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, @@ -190,11 +182,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, 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..28e1040a97 --- /dev/null +++ b/source/source_hsolver/diago_params.cpp @@ -0,0 +1,106 @@ +#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 +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, + 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 + +/// 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 diff --git a/source/source_hsolver/diago_params.h b/source/source_hsolver/diago_params.h new file mode 100644 index 0000000000..5d46b01046 --- /dev/null +++ b/source/source_hsolver/diago_params.h @@ -0,0 +1,47 @@ +#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); + +/** + * @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/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 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..6a7effb6d0 --- /dev/null +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -0,0 +1,114 @@ +#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(); +#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 +void cal_mi_lcao_wrapper(const int iter, const Input_para& inp) +{ + if (!inp.sc_mag_switch) + { + return; + } + +#ifdef __LCAO + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + sc.cal_mi_lcao(iter); +#endif +} + +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); + 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); + skip_solve = true; + } + } + + return skip_solve; +} + +/// 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, 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, + 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..959109ece7 --- /dev/null +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.h @@ -0,0 +1,67 @@ +#ifndef DELTASPIN_LCAO_H +#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 + * @param inp Input parameters + */ +template +void cal_mi_lcao_wrapper(const int iter, const Input_para& inp); + +/** + * @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