diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index d5dc5f18ddf..181635234dc 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -190,6 +190,10 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter) iter_time = ModuleBase::get_time(); + // Reset DeltaSpin RMS for this SCF iteration. If DeltaSpin is enabled, use 0 + // (meaning "lambda loop not yet run this iteration"); otherwise -1 (no RMS column). + this->ds_rms_ = this->inp_->sc_mag_switch ? 0.0 : -1.0; + if (this->inp_->esolver_type == "ksdft") { diag_ethr = hsolver::set_diagethr_ks(this->inp_->basis_type, this->inp_->esolver_type, @@ -279,8 +283,8 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool & double duration = ModuleBase::get_duration(iter_time, ModuleBase::get_time()); // print energies - elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, drho, - dkin, duration, diag_ethr); + elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, drho, + dkin, duration, diag_ethr, 0, true, this->ds_rms_); #ifdef __RAPIDJSON diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index c0b69a87984..8f7119cee71 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -70,6 +70,9 @@ class ESolver_KS : public ESolver_FP double scf_ene_thr; //! scf energy threshold double drho; //! the difference between rho_in (before HSolver) and rho_out (After HSolver) double hsolver_error; //! the error of HSolver + /// DeltaSpin RMS from the most recent lambda optimization loop; -1.0 means no DeltaSpin. + /// Set by ESolver_KS_LCAO after run_lambda_loop, read by ESolver_KS::iter_finish when calling print_etot. + double ds_rms_ = -1.0; int maxniter; //! maximum iter steps for scf int niter; //! iter steps actually used in scf bool oscillate_esolver = false; // whether esolver is oscillated diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 9f958d68a79..fb2646d222a 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -430,18 +430,20 @@ void ESolver_KS_LCAO::hamilt2rho_single(UnitCell& ucell, int istep, int spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); if (this->inp_->sc_lambda_strategy == "linear_scan") { - sc.run_lambda_linear_scan(iter - 1); + sc.run_lambda_linear_scan(iter - 1, GlobalV::ofs_running); skip_solve = true; } else if (!sc.mag_converged() && this->drho > 0 && this->drho < this->inp_->sc_scf_thr) { - sc.run_lambda_loop(iter - 1); + sc.run_lambda_loop(iter - 1, true, GlobalV::ofs_running); + this->ds_rms_ = sc.get_last_rms_error(); sc.set_mag_converged(true); skip_solve = true; } else if (sc.mag_converged()) { - sc.run_lambda_loop(iter - 1); + sc.run_lambda_loop(iter - 1, true, GlobalV::ofs_running); + this->ds_rms_ = sc.get_last_rms_error(); skip_solve = true; } } diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index f2fccf7e3c3..e46bf3ea60d 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -19,6 +19,7 @@ #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type used via general_exx_info_ #include "source_io/module_ctrl/ctrl_output_pw.h" // mohan add 20250927 #include "source_pw/module_pwdft/deltaspin_pw.h" // mohan add 20250309 +#include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/setup_pot.h" // mohan add 20250929 #include "source_pw/module_pwdft/update_cell_pw.h" // mohan add 20250309 #include "source_pw/module_pwdft/setup_dftu_pw.h" // mohan add 20250309 @@ -227,7 +228,14 @@ void ESolver_KS_PW::hamilt2rho_single(UnitCell& ucell, const int iste bool skip_charge = this->inp_->calculation == "nscf" ? true : false; // run the inner lambda loop to contrain atomic moments with the DeltaSpin method - bool skip_solve = pw::run_deltaspin_lambda_loop(iter - 1, this->drho, *this->inp_); + bool skip_solve = pw::run_deltaspin_lambda_loop(iter - 1, this->drho, *this->inp_, GlobalV::ofs_running); + if (skip_solve) + { + // Fetch the most recent DeltaSpin RMS for display in the SCF iteration table. + spinconstrain::SpinConstrain>& sc + = spinconstrain::SpinConstrain>::getScInstance(); + this->ds_rms_ = sc.get_last_rms_error(); + } if (!skip_solve) { diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 961c0302e5c..12fb5259974 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -38,6 +38,8 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) ModuleBase::timer::start("ESolver_KS_LCAO", "others"); const std::string cal_type = this->inp_->calculation; + const std::string global_out_dir = PARAM.globalv.global_out_dir; + const bool gamma_only_local = PARAM.globalv.gamma_only_local; if (cal_type == "test_memory") { @@ -78,7 +80,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) this->inp_->out_level, orb_.get_rcutmax_Phi(), ucell.infoNL->get_rcutmax_Beta(), - PARAM.globalv.gamma_only_local); + gamma_only_local); atom_arrange::search(PARAM.globalv.search_pbc, GlobalV::ofs_running, @@ -108,7 +110,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell, this->gd, this->pv, gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -175,14 +177,14 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) // pelec should be initialized before these calculations elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, - istep, PARAM.globalv.global_out_dir, *this->inp_, this->pelec); + istep, global_out_dir, *this->inp_, this->pelec); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); Get_pchg_lcao get_pchg(this->psi, &(this->pv)); - if (PARAM.globalv.gamma_only_local) + if (gamma_only_local) { get_pchg.begin(this->chr.rho, this->pelec->wg, @@ -196,7 +198,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) this->Pgrid, &this->gd, this->kv, - PARAM.globalv.global_out_dir, + global_out_dir, GlobalV::ofs_running); } else @@ -215,7 +217,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) this->Pgrid, &this->gd, this->kv, - PARAM.globalv.global_out_dir, + global_out_dir, GlobalV::ofs_running, this->inp_->if_separate_k, this->chr.ngmc); @@ -226,7 +228,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); Get_wf_lcao get_wf(this->pelec); - if (PARAM.globalv.gamma_only_local) + if (gamma_only_local) { get_wf.begin(ucell, this->psi, @@ -241,7 +243,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) this->inp_->nbands, this->inp_->nspin, PARAM.globalv.nlocal, - PARAM.globalv.global_out_dir, + global_out_dir, GlobalV::ofs_running); } else @@ -259,7 +261,7 @@ void ESolver_KS_LCAO::others(BaseCell& basecell, const int istep) this->inp_->nbands, this->inp_->nspin, PARAM.globalv.nlocal, - PARAM.globalv.global_out_dir, + global_out_dir, GlobalV::ofs_running); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting wave function"); diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index efa1dd8c62f..5ce1dfca450 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -44,7 +44,8 @@ void print_scf_iterinfo(const std::string& ks_solver, const std::vector& drho, const int& wrho, const double& time, - const int& wtime) + const int& wtime, + const double& ds_rms) { std::map iter_header_dict = {{"cg", "CG"}, @@ -91,6 +92,12 @@ void print_scf_iterinfo(const std::string& ks_solver, { td_fmt.emplace_back(" %" + std::to_string(wrho) + ".4e"); } + // DeltaSpin RMS column (optional, same width/format as DRHO/DKIN) + if (ds_rms >= 0) + { + th_fmt.emplace_back(" %" + std::to_string(wrho) + "s"); + td_fmt.emplace_back(" %" + std::to_string(wrho) + ".4e"); + } // time column, trivial th_fmt.emplace_back(" %" + std::to_string(wtime) + "s\n"); td_fmt.emplace_back(" %" + std::to_string(wtime) + ".2f\n"); @@ -132,6 +139,13 @@ void print_scf_iterinfo(const std::string& ks_solver, titles.push_back(FmtCore::center("DKIN", wrho)); values.push_back(drho[1]); } + // DeltaSpin RMS column: shown only when a valid RMS value is provided (>= 0). + // Placed after DKIN (if any) and before TIME. + if (ds_rms >= 0) + { + titles.push_back(FmtCore::center("RMS", wrho)); + values.push_back(ds_rms); + } titles.push_back(FmtCore::center("TIME/s", wtime)); values.push_back(time); std::string buf; @@ -168,7 +182,8 @@ void print_etot(const Magnetism& magnet, const double& duration, const double& pw_diag_thr, const double& avg_iter, - const bool print) + const bool print, + const double& ds_rms) { ModuleBase::TITLE("energy", "print_etot"); const int iter = iter_in; @@ -395,7 +410,8 @@ void print_etot(const Magnetism& magnet, drho, 12, duration, - 6); + 6, + ds_rms); } return; } diff --git a/source/source_estate/elecstate_print.h b/source/source_estate/elecstate_print.h index 8b54c4fb643..65383e21447 100644 --- a/source/source_estate/elecstate_print.h +++ b/source/source_estate/elecstate_print.h @@ -17,6 +17,7 @@ namespace elecstate const double& duration, const double& pw_diag_thr = 0, const double& avg_iter = 0, - bool print = true); + bool print = true, + const double& ds_rms = -1.0); } #endif diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index d07098bcd8c..e6f615b7b96 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -10,6 +10,7 @@ #include "../module_wf/get_wf_pw.h" #include "../module_wf/write_wfc_pw.h" // use write_wfc_pw #include "source_base/formatter.h" +#include "source_lcao/module_deltaspin/lambda_loop_helper.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/elecond.h" #include "source_pw/module_pwdft/onsite_proj.h" // use projector @@ -214,7 +215,7 @@ void ModuleIO::ctrl_scf_pw(const int istep, { spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); sc.cal_mi_pw(); - sc.print_Mag_Force(GlobalV::ofs_running); + spinconstrain::print_Mag_Force(sc, GlobalV::ofs_running); } //------------------------------------------------------------------ diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index e6bb5e7f2e9..201ed8273c6 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -25,6 +25,7 @@ #include "../module_dm/write_dmr.h" // use ModuleIO::write_dmr() #include "../module_dos/write_dos_lcao.h" // use ModuleIO::write_dos_lcao() #include "../module_wf/write_wfc_nao.h" // use ModuleIO::write_wfc_nao() +#include "source_lcao/module_deltaspin/lambda_loop_helper.h" // print_Mi / print_Mag_Force free functions #include "source_lcao/module_deltaspin/spin_constrain.h" // use spinconstrain::SpinConstrain #include "source_lcao/module_operator_lcao/ekinetic.h" // use hamilt::EKinetic #ifdef __MLALGO @@ -554,8 +555,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, { spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.cal_mi_lcao(istep); - sc.print_Mi(GlobalV::ofs_running); - sc.print_Mag_Force(GlobalV::ofs_running); + spinconstrain::print_Mi(sc, GlobalV::ofs_running); + spinconstrain::print_Mag_Force(sc, GlobalV::ofs_running); } //------------------------------------------------------------------ diff --git a/source/source_lcao/module_deltaspin/CMakeLists.txt b/source/source_lcao/module_deltaspin/CMakeLists.txt index 33e667a7e2a..32a91868403 100644 --- a/source/source_lcao/module_deltaspin/CMakeLists.txt +++ b/source/source_lcao/module_deltaspin/CMakeLists.txt @@ -8,8 +8,8 @@ list(APPEND objects cal_mw_from_lambda.cpp template_helpers.cpp deltaspin_lcao.cpp - sc_parse_json.cpp cal_mw_helper.cpp + mi_tools.cpp ) add_library( diff --git a/source/source_lcao/module_deltaspin/basic_funcs.cpp b/source/source_lcao/module_deltaspin/basic_funcs.cpp index 83b101de641..e044994c1e5 100644 --- a/source/source_lcao/module_deltaspin/basic_funcs.cpp +++ b/source/source_lcao/module_deltaspin/basic_funcs.cpp @@ -1,8 +1,8 @@ #include "basic_funcs.h" #include -#include "source_base/formatter.h" #include "source_base/constants.h" +#include "source_base/formatter.h" double maxval_abs_2d(const std::vector>& array) { @@ -142,8 +142,13 @@ void print_2d(const std::string info, const std::vector - -#include "source_base/matrix.h" -#include "source_base/name_angular.h" -#include "source_base/parallel_reduce.h" #include "source_base/tool_title.h" #include "source_base/timer.h" -#include "source_pw/module_pwdft/onsite_proj.h" #include "spin_constrain.h" -#include "source_io/module_parameter/parameter.h" #ifdef __LCAO #include "source_estate/elecstate_lcao.h" #include "source_lcao/hamilt_lcao.h" @@ -90,76 +83,8 @@ void spinconstrain::SpinConstrain>::cal_mi_lcao(const int& #endif -/** - * @brief Calculate atomic magnetic moments using projector overlap (PW basis). - * - * @details For each k-point: - * 1. Tabulate atomic projectors: set up |alpha_{l,m}> for each atom - * 2. Compute becp = via overlap_proj_psi - * 3. Decompose becp into magnetic moments via accumulate_Mi_from_becp - * - * The magnetic moment is computed as: - * Mi = sum_{k,i} w_{k,i} * - * where P_at is the atomic projector and sigma are the Pauli matrices. - * - * Finally, Mi is summed across all MPI k-pool ranks since each pool only - * has a subset of k-points. - */ -template <> -void spinconstrain::SpinConstrain>::cal_mi_pw() -{ - ModuleBase::TITLE("module_deltaspin", "cal_mi_pw"); - ModuleBase::timer::start("spinconstrain::SpinConstrain", "cal_mi_pw"); - - this->zero_Mi(); - if(PARAM.inp.device == "cpu") - { - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - // Loop over k-points to calculate Mi of sum_{k,i,l,m} - std::complex* psi_pointer = nullptr; - psi::Psi, base_device::DEVICE_CPU>* psi_t = static_cast, base_device::DEVICE_CPU>*>(this->psi); - const int nbands = psi_t->get_nbands(); - const int nks = psi_t->get_nk(); - const int npol = psi_t->get_npol(); - for(int ik = 0; ik < nks; ik++) - { - psi_t->fix_k(ik); - psi_pointer = psi_t->get_pointer(); - onsite_p->tabulate_atomic(ik); // Set up atomic projectors for this k-point - onsite_p->overlap_proj_psi(nbands * npol, psi_pointer); // Compute becp = - const std::complex* becp = onsite_p->get_h_becp(); - int nkb = onsite_p->get_tot_nproj(); - this->accumulate_Mi_from_becp(becp, nkb, nbands, npol, ik, - &this->pelec->wg(ik, 0), &onsite_p->get_nh(0)); - } - } -#if ((defined __CUDA) || (defined __ROCM)) - else - { - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - std::complex* psi_pointer = nullptr; - psi::Psi, base_device::DEVICE_GPU>* psi_t = static_cast, base_device::DEVICE_GPU>*>(this->psi); - const int nbands = psi_t->get_nbands(); - const int nks = psi_t->get_nk(); - const int npol = psi_t->get_npol(); - for(int ik = 0; ik < nks; ik++) - { - psi_t->fix_k(ik); - psi_pointer = psi_t->get_pointer(); - onsite_p->tabulate_atomic(ik); - onsite_p->overlap_proj_psi(nbands * npol, psi_pointer); - const std::complex* becp = onsite_p->get_h_becp(); - int nkb = onsite_p->get_size_becp() / nbands / npol; - this->accumulate_Mi_from_becp(becp, nkb, nbands, npol, ik, - &this->pelec->wg(ik, 0), &onsite_p->get_nh(0)); - } - } -#endif - // MPI reduction: sum Mi across all k-pool ranks - Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, GlobalV::NPROC_IN_POOL, &(this->Mi_[0][0]), 3 * this->Mi_.size()); - - ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mi_pw"); -} +// cal_mi_pw() has been moved to source/source_pw/module_pwdft/deltaspin_pw_impl.cpp +// because it depends on PW-specific OnsiteProjector. /// @brief Set the DeltaSpin operator pointer for LCAO magnetic moment calculation template <> diff --git a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp index b0019697a16..a77630a1119 100644 --- a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp +++ b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp @@ -4,17 +4,14 @@ #include "source_hsolver/diago_iter_assist.h" #include "source_io/module_parameter/parameter.h" #include "spin_constrain.h" +#include "mi_tools.h" #include "source_pw/module_pwdft/onsite_proj.h" #include "source_base/parallel_reduce.h" -#include "source_base/kernels/math_kernel_op.h" #include "source_hsolver/hsolver_lcao.h" -#include "source_hsolver/hsolver_pw.h" -#include "source_estate/elecstate_pw.h" #include "source_estate/elecstate_tools.h" #ifdef __LCAO #include "source_estate/elecstate_lcao.h" -#include "source_estate/elecstate_tools.h" #include "source_estate/module_dm/cal_dm_psi.h" #include "source_lcao/module_operator_lcao/dspin_lcao.h" #endif @@ -52,437 +49,15 @@ * Solution: Ensure cal_mw_from_lambda() is called at the start of each SCF step. */ -/** - * @brief Compute DeltaSpin correction to the subspace Hamiltonian. - * - * @details Adds the constraint term to H in the projector subspace: - * H += becp^† * ps, where ps = delta_lambda * becp - * - * For non-collinear (npol=2), this implements the full 2x2 Pauli matrix: - * H_delta = | lambda_z lambda_x + i*lambda_y | - * | lambda_x - i*lambda_y -lambda_z | - * - * For collinear (npol=1), only the diagonal z-component with spin_sign: - * H_delta = lambda_z * spin_sign - * - * @param h_tmp Subspace Hamiltonian (nbands x nbands, modified in place) - * @param becp_k Projector coefficients for k-point ik - * @param delta_lambda Lambda change per atom (or full lambda if full_update) - * @param nbands Number of bands - * @param nkb Total number of projectors - * @param nh_iat Number of projectors per atom - * @param ik K-point index (for spin_sign lookup in collinear mode) - * @param full_update If true, compute delta = lambda_current - lambda_at_save - */ -template <> -void spinconstrain::SpinConstrain>::calculate_delta_hcc(std::complex* h_tmp, const std::complex* becp_k, const ModuleBase::Vector3* delta_lambda, const int nbands, const int nkb, const int* nh_iat, const int ik, bool full_update) -{ - ModuleBase::TITLE("spinconstrain::SpinConstrain", "calculate_delta_hcc"); - ModuleBase::timer::start("spinconstrain::SpinConstrain", "calculate_delta_hcc"); - - // If full_update, compute actual delta = lambda_current - lambda_at_save - // This applies only the CHANGE in lambda, not the full lambda value - std::vector> actual_delta; - const ModuleBase::Vector3* effective_lambda = delta_lambda; - if (full_update) - { - int nat = this->get_nat(); - actual_delta.resize(nat); - for (int iat = 0; iat < nat; iat++) - { - actual_delta[iat] = delta_lambda[iat] - this->lambda_in_sub_[iat]; - } - effective_lambda = actual_delta.data(); - } - - int sum = 0; // Running sum of projectors across atoms - int size_ps = nkb * this->npol_ * nbands; // Total size of ps array - std::complex* becp_cpu = nullptr; - - // Handle GPU/CPU memory for becp - if(PARAM.inp.device == "gpu") - { -#if ((defined __CUDA) || (defined __ROCM)) - base_device::memory::resize_memory_op, base_device::DEVICE_CPU>()(becp_cpu, size_ps); - base_device::memory::synchronize_memory_op, base_device::DEVICE_CPU, base_device::DEVICE_GPU>()(becp_cpu, becp_k, size_ps); -#endif - } - else if (PARAM.inp.device == "cpu") - { - becp_cpu = const_cast*>(becp_k); - } - - // Compute modified projector coefficients: ps = delta_lambda * becp - std::vector> ps(size_ps, 0.0); - if(this->npol_ == 2) - { - // ============================================================= - // nspin=4 (non-collinear): full Pauli matrix treatment - // ============================================================= - // For each atom, construct 2x2 coefficients: - // | lambda_z lambda_x + i*lambda_y | - // | lambda_x - i*lambda_y -lambda_z | - // Then: ps_up = coeff0 * becp_up + coeff2 * becp_dn - // ps_dn = coeff1 * becp_up + coeff3 * becp_dn - for (int iat = 0; iat < this->Mi_.size(); iat++) - { - const int nproj = nh_iat[iat]; - const std::complex coefficients0(effective_lambda[iat][2], 0.0); - const std::complex coefficients1(effective_lambda[iat][0] , effective_lambda[iat][1]); - const std::complex coefficients2(effective_lambda[iat][0] , -1 * effective_lambda[iat][1]); - const std::complex coefficients3(-1 * effective_lambda[iat][2], 0.0); - for (int ib = 0; ib < nbands * this->npol_; ib += this->npol_) - { - for (int ip = 0; ip < nproj; ip++) - { - const int becpind = ib * nkb + sum + ip; - const std::complex becp1 = becp_cpu[becpind]; - const std::complex becp2 = becp_cpu[becpind + nkb]; - ps[becpind] += coefficients0 * becp1 - + coefficients2 * becp2; - ps[becpind + nkb] += coefficients1 * becp1 - + coefficients3 * becp2; - } - } - sum += nproj; - } - } - else if(this->npol_ == 1) - { - // ============================================================= - // nspin=2 (collinear): only z-component with spin_sign - // ============================================================= - // ps = lambda_z * spin_sign * becp - // spin_sign = +1 for spin-up k-points, -1 for spin-down - for (int iat = 0; iat < this->Mi_.size(); iat++) - { - const int nproj = nh_iat[iat]; - double coefficients0 = effective_lambda[iat][2] * this->get_spin_sign(ik); - for (int ib = 0; ib < nbands; ib++) - { - for (int ip = 0; ip < nproj; ip++) - { - const int becpind = ib * nkb + sum + ip; - const std::complex becp1 = becp_cpu[becpind]; - ps[becpind] += coefficients0 * becp1; - } - } - sum += nproj; - } - } - - // Copy ps to GPU if needed - std::complex* ps_pointer = nullptr; - if(PARAM.inp.device == "gpu") - { -#if ((defined __CUDA) || (defined __ROCM)) - base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(ps_pointer, size_ps); - base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_CPU>()(ps_pointer, ps.data(), size_ps); -#endif - } - else if (PARAM.inp.device == "cpu") - { - ps_pointer = ps.data(); - } - - // ============================================================= - // H += becp^† * ps (GEMM: C = alpha * A^† * B + beta * C) - // A = becp_k (npm x nbands), B = ps (npm x nbands), C = h_tmp (nbands x nbands) - // ============================================================= - char transa = 'C'; // Conjugate transpose of becp - char transb = 'N'; // Normal ps - const int npm = nkb * this->npol_; - if (PARAM.inp.device == "gpu") - { -#if ((defined __CUDA) || (defined __ROCM)) - ModuleBase::gemm_op, base_device::DEVICE_GPU>()( - transa, - transb, - nbands, - nbands, - npm, - &ModuleBase::ONE, - becp_k, - npm, - ps_pointer, - npm, - &ModuleBase::ONE, - h_tmp, - nbands - ); - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(ps_pointer); - base_device::memory::delete_memory_op, base_device::DEVICE_CPU>()(becp_cpu); -#endif - - } - else if (PARAM.inp.device == "cpu") - { - ModuleBase::gemm_op, base_device::DEVICE_CPU>()( - transa, - transb, - nbands, - nbands, - npm, - &ModuleBase::ONE, - becp_k, - npm, - ps_pointer, - npm, - &ModuleBase::ONE, - h_tmp, - nbands - ); - } - ModuleBase::timer::end("spinconstrain::SpinConstrain", "calculate_delta_hcc"); -} - -/** - * @brief CPU implementation of PW wavefunction and charge density update. - * - * @par Two-stage process: - * Stage 1 - Subspace diagonalization: - * For each k-point, apply DeltaSpin correction to the saved subspace H, - * then diagonalize to rotate the wavefunctions. This is a cheap operation - * in the reduced subspace (nbands x nbands). - * - * Stage 2 - Full-space update: - * Option A (pw_solve=true): Run HSolverPW for iterative refinement in the - * full plane-wave space. This is more accurate but expensive. - * Option B (pw_solve=false): Update weights from new eigenvalues and call - * psiToRho() to build the charge density from current psi. Faster but - * may be less accurate if the subspace rotation was not sufficient. - * - * @par Memory management - * Frees sub_h_save, sub_s_save, becp_save after use. These are allocated - * on the first cal_mw_from_lambda() call and should only be freed here. - * - * @param delta_lambda Lambda change for incremental H correction - * @param pw_solve If true, run full PW solver; if false, just update weights - * @param full_update If true, apply full lambda (not delta) to H correction - */ -template <> -void spinconstrain::SpinConstrain>::update_psi_charge_pw_cpu(const ModuleBase::Vector3* delta_lambda, bool pw_solve, bool full_update) -{ - ModuleBase::TITLE("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); - ModuleBase::timer::start("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); - - psi::Psi>* psi_t = static_cast>*>(this->psi); - hamilt::Hamilt, base_device::DEVICE_CPU>* hamilt_t = static_cast, base_device::DEVICE_CPU>*>(this->p_hamilt); - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - - int nbands = psi_t->get_nbands(); - int npol = psi_t->get_npol(); - int nkb = onsite_p->get_tot_nproj(); - int nk = psi_t->get_nk(); - int size_becp = nbands * nkb * npol; - const int* nh_iat = &onsite_p->get_nh(0); - - std::vector> h_tmp(nbands * nbands), s_tmp(nbands * nbands); - - // CRITICAL: subspace data must have been saved by cal_mw_from_lambda() - assert(this->sub_h_save != nullptr); - assert(this->sub_s_save != nullptr); - assert(this->becp_save != nullptr); - - // Determine which lambda to use for H correction - const ModuleBase::Vector3* lambda_for_hcc = delta_lambda; - std::vector> computed_delta; - if (full_update) - { - lambda_for_hcc = this->lambda_.data(); - } - - // ============================================================= - // STAGE 1: Subspace diagonalization for each k-point - // ============================================================= - for (int ik = 0; ik < nk; ++ik) - { - std::complex* h_k = this->sub_h_save + ik * nbands * nbands; - std::complex* s_k = this->sub_s_save + ik * nbands * nbands; - std::complex* becp_k = this->becp_save + ik * size_becp; - - psi_t->fix_k(ik); - - // Copy saved subspace matrices to temp - memcpy(h_tmp.data(), h_k, sizeof(std::complex) * nbands * nbands); - memcpy(s_tmp.data(), s_k, sizeof(std::complex) * nbands * nbands); - - // Apply DeltaSpin correction: H += becp^† * lambda * becp - this->calculate_delta_hcc(h_tmp.data(), becp_k, lambda_for_hcc, nbands, nkb, nh_iat, ik, full_update); - - // Diagonalize in subspace to update wavefunction coefficients and eigenvalues - hsolver::DiagoIterAssist>::diag_subspace_psi(h_tmp.data(), - s_tmp.data(), - nbands, - psi_t[0], - &this->pelec->ekb(ik, 0)); - } - - // Free saved subspace data (allocated in cal_mw_from_lambda) - delete[] this->sub_h_save; - delete[] this->sub_s_save; - delete[] this->becp_save; - this->sub_h_save = nullptr; - this->sub_s_save = nullptr; - this->becp_save = nullptr; - - // ============================================================= - // STAGE 2: Full-space update - // ============================================================= - if (pw_solve) - { - // Full PW diagonalization: subspace rotation provides a good initial guess, - // then HSolverPW iteratively refines psi in the full plane-wave space and calls psiToRho. - hsolver::HSolverPW, base_device::DEVICE_CPU> hsolver_pw_obj( - this->pw_wfc_, - PARAM.inp.calculation, - PARAM.inp.basis_type, - PARAM.inp.ks_solver, - PARAM.globalv.use_uspp, - PARAM.inp.nspin, - hsolver::DiagoIterAssist>::SCF_ITER, - hsolver::DiagoIterAssist>::PW_DIAG_NMAX, - hsolver::DiagoIterAssist>::PW_DIAG_THR, - hsolver::DiagoIterAssist>::need_subspace, - PARAM.inp.nbands, - PARAM.inp.diago_smooth_ethr, - PARAM.inp.pw_diag_ndim, - PARAM.inp.diag_subspace, - PARAM.inp.nb2d, - PARAM.inp.use_k_continuity); - - hsolver_pw_obj.solve(hamilt_t, psi_t[0], this->pelec, this->pelec->ekb.c, - GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, false, this->tpiba, this->get_nat()); - } - else - { - // No full solver: update weights from new eigenvalues, then build rho from current psi - elecstate::calculate_weights(this->pelec->ekb, - this->pelec->wg, - this->pelec->klist, - this->pelec->eferm, - this->pelec->f_en, - this->pelec->nelec_spin, - PARAM.inp.nbands, - this->pelec->skip_weights); - elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en); - reinterpret_cast, base_device::DEVICE_CPU>*>(this->pelec)->psiToRho(*psi_t); - } - ModuleBase::timer::end("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); -} - -#if ((defined __CUDA) || (defined __ROCM)) -/** - * @brief GPU implementation of PW wavefunction and charge density update. - * - * @details Same algorithm as update_psi_charge_pw_cpu(), but with GPU memory - * management (device allocation, host-device synchronization). - */ -template <> -void spinconstrain::SpinConstrain>::update_psi_charge_pw_gpu(const ModuleBase::Vector3* delta_lambda, bool pw_solve, bool full_update) -{ - ModuleBase::TITLE("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); - ModuleBase::timer::start("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); - - psi::Psi, base_device::DEVICE_GPU>* psi_t = static_cast, base_device::DEVICE_GPU>*>(this->psi); - hamilt::Hamilt, base_device::DEVICE_GPU>* hamilt_t = static_cast, base_device::DEVICE_GPU>*>(this->p_hamilt); - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - - int nbands = psi_t->get_nbands(); - int npol = psi_t->get_npol(); - int nkb = onsite_p->get_tot_nproj(); - int nk = psi_t->get_nk(); - int size_becp = nbands * nkb * npol; - const int* nh_iat = &onsite_p->get_nh(0); - - std::complex* h_tmp = nullptr; - std::complex* s_tmp = nullptr; - base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(h_tmp, nbands * nbands); - base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(s_tmp, nbands * nbands); - - assert(this->sub_h_save != nullptr); - assert(this->sub_s_save != nullptr); - assert(this->becp_save != nullptr); - - const ModuleBase::Vector3* lambda_for_hcc = delta_lambda; - std::vector> computed_delta; - if (full_update) - { - lambda_for_hcc = this->lambda_.data(); - } - - // STAGE 1: Subspace diagonalization for each k-point (GPU) - for (int ik = 0; ik < nk; ++ik) - { - std::complex* h_k = this->sub_h_save + ik * nbands * nbands; - std::complex* s_k = this->sub_s_save + ik * nbands * nbands; - std::complex* becp_k = this->becp_save + ik * size_becp; - - psi_t->fix_k(ik); - - base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(h_tmp, h_k, nbands * nbands); - base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(s_tmp, s_k, nbands * nbands); - - this->calculate_delta_hcc(h_tmp, becp_k, lambda_for_hcc, nbands, nkb, nh_iat, ik, full_update); - - hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::diag_subspace_psi(h_tmp, - s_tmp, - nbands, - psi_t[0], - &this->pelec->ekb(ik, 0)); - } - - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(h_tmp); - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(s_tmp); +// calculate_delta_hcc() has been moved to +// source/source_pw/module_pwdft/deltaspin_pw_impl.cpp +// because it is PW-specific (operates on PW projector subspace Hamiltonian +// and is only called by update_psi_charge_pw_cpu/gpu). - // Free GPU memory for saved subspace data - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(sub_h_save); - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(sub_s_save); - base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(becp_save); - this->sub_h_save = nullptr; - this->sub_s_save = nullptr; - this->becp_save = nullptr; - - // STAGE 2: Full-space update (GPU) - if (pw_solve) - { - hsolver::HSolverPW, base_device::DEVICE_GPU> hsolver_pw_obj( - this->pw_wfc_, - PARAM.inp.calculation, - PARAM.inp.basis_type, - PARAM.inp.ks_solver, - PARAM.globalv.use_uspp, - PARAM.inp.nspin, - hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::SCF_ITER, - hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::PW_DIAG_NMAX, - hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::PW_DIAG_THR, - hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::need_subspace, - PARAM.inp.nbands, - PARAM.inp.diago_smooth_ethr, - PARAM.inp.pw_diag_ndim, - PARAM.inp.diag_subspace, - PARAM.inp.nb2d, - PARAM.inp.use_k_continuity); - - hsolver_pw_obj.solve(hamilt_t, psi_t[0], this->pelec, this->pelec->ekb.c, - GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, false, this->tpiba, this->get_nat()); - } - else - { - elecstate::calculate_weights(this->pelec->ekb, - this->pelec->wg, - this->pelec->klist, - this->pelec->eferm, - this->pelec->f_en, - this->pelec->nelec_spin, - PARAM.inp.nbands, - this->pelec->skip_weights); - elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en); - reinterpret_cast, base_device::DEVICE_GPU>*>(this->pelec)->psiToRho(*psi_t); - } - ModuleBase::timer::end("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); -} -#endif +// update_psi_charge_pw_cpu/gpu() have been moved to +// source/source_pw/module_pwdft/deltaspin_pw_impl.cpp +// because they depend on PW-specific HSolverPW, ElecStatePW, and +// OnsiteProjector. /** * @brief Core workflow: apply lambda -> solve Hamiltonian -> compute magnetic moments. @@ -705,8 +280,9 @@ void spinconstrain::SpinConstrain>::cal_mw_from_lambda( for (int ik = 0; ik < nk; ik++) { const std::complex* becp = &becp_tmp[ik * size_becp]; - this->accumulate_Mi_from_becp(becp, nkb, nbands, this->npol_, ik, - &this->pelec->wg(ik, 0), nh_iat); + const int spin_sign = (this->npol_ == 2) ? 1 : this->get_spin_sign(ik); + accumulate_Mi_from_becp(becp, nkb, nbands, this->npol_, spin_sign, + &this->pelec->wg(ik, 0), nh_iat, this->Mi_); } // MPI reduction: sum Mi across all k-pool ranks Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp index 811d6fef193..42eac35ee41 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.cpp @@ -125,7 +125,8 @@ 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) + const Input_para& inp, + std::ostream& ofs_running) { bool skip_solve = false; @@ -136,14 +137,14 @@ bool run_deltaspin_lambda_loop_lcao(const int iter, if (!sc.mag_converged() && drho > 0 && drho < inp.sc_scf_thr) { /// Charge density is stable enough: optimize lambda for the first time - sc.run_lambda_loop(iter); + sc.run_lambda_loop(iter, true, ofs_running); sc.set_mag_converged(true); skip_solve = true; } else if (sc.mag_converged()) { /// Already converged: refine lambda for the current charge density - sc.run_lambda_loop(iter); + sc.run_lambda_loop(iter, true, ofs_running); skip_solve = true; } } @@ -174,9 +175,11 @@ template void cal_mi_lcao_wrapper>(const int iter, const In template bool run_deltaspin_lambda_loop_lcao(const int iter, const double drho, - const Input_para& inp); + const Input_para& inp, + std::ostream& ofs_running); template bool run_deltaspin_lambda_loop_lcao>(const int iter, const double drho, - const Input_para& inp); + const Input_para& inp, + std::ostream& ofs_running); } // namespace ModuleESolver diff --git a/source/source_lcao/module_deltaspin/deltaspin_lcao.h b/source/source_lcao/module_deltaspin/deltaspin_lcao.h index 959109ece76..dc859590d48 100644 --- a/source/source_lcao/module_deltaspin/deltaspin_lcao.h +++ b/source/source_lcao/module_deltaspin/deltaspin_lcao.h @@ -1,6 +1,8 @@ #ifndef DELTASPIN_LCAO_H #define DELTASPIN_LCAO_H +#include + #include "source_cell/unitcell.h" #include "source_cell/klist.h" #include "source_io/module_parameter/input_parameter.h" @@ -60,7 +62,8 @@ 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); + const Input_para& inp, + std::ostream& ofs_running); } // namespace ModuleESolver diff --git a/source/source_lcao/module_deltaspin/lambda_loop.cpp b/source/source_lcao/module_deltaspin/lambda_loop.cpp index 425e5d74b3d..27a87ec6492 100644 --- a/source/source_lcao/module_deltaspin/lambda_loop.cpp +++ b/source/source_lcao/module_deltaspin/lambda_loop.cpp @@ -1,5 +1,4 @@ #include "spin_constrain.h" - #include #include #include @@ -7,8 +6,9 @@ #include #include "basic_funcs.h" -#include "source_io/module_parameter/parameter.h" +#include "lambda_loop_helper.h" #include "source_base/constants.h" +#include "source_io/module_parameter/parameter.h" /** * @file lambda_loop.cpp @@ -71,7 +71,7 @@ * - decay_grad thresholds are not too aggressive */ template <> -void spinconstrain::SpinConstrain>::run_lambda_loop(int outer_step, bool rerun) +void spinconstrain::SpinConstrain>::run_lambda_loop(int outer_step, bool rerun, std::ostream& ofs_running) { int nat = this->get_nat(); int ntype = this->get_ntype(); @@ -106,7 +106,7 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out double inner_loop_duration = 0.0; - this->print_header(); + print_header(*this, ofs_running); // ============================================================= // MAIN OPTIMIZATION LOOP @@ -128,9 +128,9 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out // Save initial lambda: for unconstrained components (constrain==0), set to 0 where_fill_scalar_else_2d(this->constrain_, 0, zero, this->lambda_, initial_lambda); - print_2d("initial lambda (eV/uB): ", initial_lambda, this->nspin_, ModuleBase::Ry_to_eV); - print_2d("initial spin (uB): ", spin, this->nspin_); - print_2d("target spin (uB): ", this->target_mag_, this->nspin_); + print_2d(" initial lambda (eV/uB): ", initial_lambda, this->nspin_, ModuleBase::Ry_to_eV, ofs_running); + print_2d(" initial spin (uB): ", spin, this->nspin_, 1.0, ofs_running); + print_2d(" target spin (uB): ", this->target_mag_, this->nspin_, 1.0, ofs_running); i_step++; } else @@ -170,7 +170,7 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out new_spin = this->Mi_; // Check if gradient dM/dlambda has decayed below threshold - bool GradLessThanBound = this->check_gradient_decay(new_spin, spin, delta_lambda, dnu_last_step); + bool GradLessThanBound = check_gradient_decay(*this, new_spin, spin, delta_lambda, dnu_last_step, false, ofs_running); if (i_step >= this->nsc_min_ && GradLessThanBound) { // Gradient has decayed: further optimization yields diminishing returns @@ -185,8 +185,8 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out - iterstart)).count() / static_cast(1e6); #endif inner_loop_duration += duration; - std::cout << "Total TIME(s) = " << inner_loop_duration << std::endl; - this->print_termination(); + ofs_running << " Total TIME(s) = " << inner_loop_duration << std::endl; + print_termination(*this, ofs_running); break; } spin = new_spin; @@ -259,8 +259,10 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out - iterstart)).count() / static_cast(1e6); #endif inner_loop_duration += duration; - if (this->check_rms_stop(outer_step, i_step, rms_error, duration, inner_loop_duration)) + if (check_rms_stop(*this, outer_step, i_step, rms_error, duration, inner_loop_duration, ofs_running)) { + // Save RMS for ESolver to display in the SCF iteration table. + this->last_rms_error_ = rms_error; // Converged or max steps reached: final update this->update_psi_charge(dnu_last_step.data(), rerun, true); @@ -280,14 +282,15 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out } mean_error = sum_2d(temp_1) / nat; rms_error = std::sqrt(mean_error); - std::cout<<"Current RMS: "< this->current_sc_thr_ * 10 && rerun == true && this->higher_mag_prec == true) { - std::cout<<"Error: RMS error is too large, rerun the loop"<run_lambda_loop(outer_step, false); + std::cout<<" DeltaSpin: RMS error too large ("<run_lambda_loop(outer_step, false, ofs_running); } } break; @@ -313,7 +316,7 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out } // Cap step size to prevent overshooting - this->check_restriction(search, alpha_trial); + check_restriction(*this, search, alpha_trial, ofs_running); // ============================================================= // CUMULATIVE STEP UPDATE @@ -351,8 +354,8 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out spin_plus = this->Mi_; // Find optimal step size via linear interpolation - alpha_opt = this->cal_alpha_opt(spin, spin_plus, alpha_trial); - this->check_restriction(search, alpha_opt); + alpha_opt = cal_alpha_opt(*this, spin, spin_plus, alpha_trial); + check_restriction(*this, search, alpha_opt, ofs_running); // Correct dnu: dnu += (alpha_opt - alpha_trial) * search alpha_plus = alpha_opt - alpha_trial; @@ -419,7 +422,7 @@ void spinconstrain::SpinConstrain>::run_lambda_loop(int out * step, lambda_eV_uB, Mi_x_0, Mi_y_0, Mi_z_0, Mi_x_1, ... */ template <> -void spinconstrain::SpinConstrain>::run_lambda_linear_scan(int outer_step) +void spinconstrain::SpinConstrain>::run_lambda_linear_scan(int outer_step, std::ostream& ofs_running) { int nat = this->get_nat(); int ntype = this->get_ntype(); @@ -429,7 +432,7 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( int nsteps = PARAM.inp.sc_scan_steps; if (nsteps <= 0) { - std::cout << "[DS-DIAG] linear_scan: sc_scan_steps <= 0, skipping" << std::endl; + ofs_running << " [DS-DIAG] linear_scan: sc_scan_steps <= 0, skipping" << std::endl; return; } @@ -438,15 +441,15 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( double lambda_end_ry = lambda_end / ModuleBase::Ry_to_eV; double lambda_step = (lambda_end_ry - lambda_start_ry) / (nsteps - 1); - std::cout << "\n" << std::string(80, '=') << std::endl; - std::cout << "[DS-DIAG] === LINEAR LAMBDA SCAN START ===" << std::endl; - std::cout << "[DS-DIAG] Scan range: " << lambda_start << " -> " << lambda_end << " eV/uB" << std::endl; - std::cout << "[DS-DIAG] Number of steps: " << nsteps << std::endl; - std::cout << "[DS-DIAG] Lambda step size: " << lambda_step * ModuleBase::Ry_to_eV << " eV/uB" << std::endl; - std::cout << "[DS-DIAG] nat = " << nat << ", ntype = " << ntype << std::endl; - std::cout << "[DS-DIAG] nspin_ = " << this->nspin_ << ", npol_ = " << this->npol_ << std::endl; - std::cout << "[DS-DIAG] p_operator = " << (this->p_operator ? "valid" : "NULL") << std::endl; - std::cout << "[DS-DIAG] constrain_ size = " << this->constrain_.size() << std::endl; + ofs_running << "\n" << std::string(80, '=') << std::endl; + ofs_running << " [DS-DIAG] === LINEAR LAMBDA SCAN START ===" << std::endl; + ofs_running << " [DS-DIAG] Scan range: " << lambda_start << " -> " << lambda_end << " eV/uB" << std::endl; + ofs_running << " [DS-DIAG] Number of steps: " << nsteps << std::endl; + ofs_running << " [DS-DIAG] Lambda step size: " << lambda_step * ModuleBase::Ry_to_eV << " eV/uB" << std::endl; + ofs_running << " [DS-DIAG] nat = " << nat << ", ntype = " << ntype << std::endl; + ofs_running << " [DS-DIAG] nspin_ = " << this->nspin_ << ", npol_ = " << this->npol_ << std::endl; + ofs_running << " [DS-DIAG] p_operator = " << (this->p_operator ? "valid" : "NULL") << std::endl; + ofs_running << " [DS-DIAG] constrain_ size = " << this->constrain_.size() << std::endl; // Check if any constraints are defined; if not, set all atoms as constrained bool has_constraints = false; @@ -458,7 +461,7 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( } if (!has_constraints) { - std::cout << "[DS-DIAG] No constraints found in STRU, setting all atoms as constrained" << std::endl; + ofs_running << " [DS-DIAG] No constraints found in STRU, setting all atoms as constrained" << std::endl; for (int ia = 0; ia < nat; ia++) { if (this->nspin_ == 4) { this->constrain_[ia] = ModuleBase::Vector3(1, 1, 1); @@ -470,11 +473,11 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( } for (int ia = 0; ia < nat; ia++) { - std::cout << "[DS-DIAG] Atom " << ia << " constrain = (" - << this->constrain_[ia].x << ", " << this->constrain_[ia].y << ", " << this->constrain_[ia].z << ")" - << " target_mag = (" << this->target_mag_[ia].x << ", " << this->target_mag_[ia].y << ", " << this->target_mag_[ia].z << ")" << std::endl; + ofs_running << " [DS-DIAG] Atom " << ia << " constrain = (" + << this->constrain_[ia].x << ", " << this->constrain_[ia].y << ", " << this->constrain_[ia].z << ")" + << " target_mag = (" << this->target_mag_[ia].x << ", " << this->target_mag_[ia].y << ", " << this->target_mag_[ia].z << ")" << std::endl; } - std::cout << std::string(80, '=') << "\n" << std::endl; + ofs_running << std::string(80, '=') << "\n" << std::endl; // Save initial lambda to restore after scan std::vector> initial_lambda(nat, 0.0); @@ -526,8 +529,8 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( } } - std::cout << "[DS-DIAG] === Scan step " << istep << "/" << nsteps - << " lambda = " << lambda_val_ev << " eV/uB ===" << std::endl; + ofs_running << " [DS-DIAG] === Scan step " << istep << "/" << nsteps + << " lambda = " << lambda_val_ev << " eV/uB ===" << std::endl; // Compute magnetic moments at current lambda this->cal_mw_from_lambda(istep); @@ -547,14 +550,14 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( } ofs_scan << std::endl; - std::cout << "[DS-DIAG] lambda = " << lambda_val_ev << " eV/uB" << std::endl; + ofs_running << " [DS-DIAG] lambda = " << lambda_val_ev << " eV/uB" << std::endl; for (int ia = 0; ia < nat; ia++) { - std::cout << "[DS-DIAG] Atom " << ia << " Mi = (" - << this->Mi_[ia].x << ", " - << this->Mi_[ia].y << ", " - << this->Mi_[ia].z << ") uB" << std::endl; + ofs_running << " [DS-DIAG] Atom " << ia << " Mi = (" + << this->Mi_[ia].x << ", " + << this->Mi_[ia].y << ", " + << this->Mi_[ia].z << ") uB" << std::endl; } - std::cout << std::endl; + ofs_running << std::endl; } // ============================================================= @@ -562,7 +565,7 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( // to verify that the lambda->Mi mapping is numerically stable // after multiple lambda updates in the scan loop // ============================================================= - std::cout << "[DS-DIAG] === Consistency check: restoring initial lambda ===" << std::endl; + ofs_running << " [DS-DIAG] === Consistency check: restoring initial lambda ===" << std::endl; this->lambda_ = initial_lambda; this->cal_mw_from_lambda(nsteps); @@ -576,12 +579,12 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( } ofs_scan << std::endl; - std::cout << "[DS-DIAG] lambda = " << lambda_start << " eV/uB (restored)" << std::endl; + ofs_running << " [DS-DIAG] lambda = " << lambda_start << " eV/uB (restored)" << std::endl; for (int ia = 0; ia < nat; ia++) { - std::cout << "[DS-DIAG] Atom " << ia << " Mi = (" - << this->Mi_[ia].x << ", " - << this->Mi_[ia].y << ", " - << this->Mi_[ia].z << ") uB" << std::endl; + ofs_running << " [DS-DIAG] Atom " << ia << " Mi = (" + << this->Mi_[ia].x << ", " + << this->Mi_[ia].y << ", " + << this->Mi_[ia].z << ") uB" << std::endl; } // Compare restored Mi with step 0 Mi to check consistency @@ -595,11 +598,11 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( if (diff > max_mi_diff) max_mi_diff = diff; ofs_scan << "# Atom " << ia << " dM = (" << dx << ", " << dy << ", " << dz << ") uB" << std::endl; } - std::cout << "[DS-DIAG] Max Mi difference between step 0 and init_recheck: " << max_mi_diff << " uB" << std::endl; + ofs_running << " [DS-DIAG] Max Mi difference between step 0 and init_recheck: " << max_mi_diff << " uB" << std::endl; if (max_mi_diff > 1e-8) { - std::cout << "[DS-DIAG] WARNING: Mi mapping may be inconsistent after multiple lambda updates!" << std::endl; + ofs_running << " [DS-DIAG] WARNING: Mi mapping may be inconsistent after multiple lambda updates!" << std::endl; } else { - std::cout << "[DS-DIAG] OK: Mi mapping is consistent." << std::endl; + ofs_running << " [DS-DIAG] OK: Mi mapping is consistent." << std::endl; } ofs_scan << "# Max Mi difference: " << max_mi_diff << " uB" << std::endl; @@ -608,10 +611,10 @@ void spinconstrain::SpinConstrain>::run_lambda_linear_scan( // Restore original lambda values (already restored above, but explicit for clarity) this->lambda_ = initial_lambda; - std::cout << std::string(80, '=') << std::endl; - std::cout << "[DS-DIAG] === LINEAR LAMBDA SCAN COMPLETE ===" << std::endl; - std::cout << "[DS-DIAG] Results written to: lambda_scan_results.dat" << std::endl; - std::cout << std::string(80, '=') << "\n" << std::endl; + ofs_running << std::string(80, '=') << std::endl; + ofs_running << " [DS-DIAG] === LINEAR LAMBDA SCAN COMPLETE ===" << std::endl; + ofs_running << " [DS-DIAG] Results written to: lambda_scan_results.dat" << std::endl; + ofs_running << std::string(80, '=') << "\n" << std::endl; return; } diff --git a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp index 9f101bd13f7..b24b29cc9b6 100644 --- a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp +++ b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp @@ -1,40 +1,53 @@ +#include "lambda_loop_helper.h" + +#include +#include +#include +#include +#include +#include +#include + #include "basic_funcs.h" -#include "spin_constrain.h" +#include "source_base/formatter.h" /** * @file lambda_loop_helper.cpp - * @brief Helper/auxiliary methods for the lambda optimization loop. + * @brief Free-function implementations of lambda loop helpers. * - * @par Functions overview - * - print_termination(): Print final spin and lambda values when loop exits - * - check_rms_stop(): Check convergence and print step info - * - print_header(): Print header at loop start - * - check_restriction(): Cap step size to prevent overshooting - * - cal_alpha_opt(): Compute optimal step size via linear interpolation - * - check_gradient_decay(): Check if dM/dlambda has decayed below threshold + * @par History + * Originally these were explicit specializations of SpinConstrain member + * functions: + * template <> void SpinConstrain>::print_termination(...) + * They have been lifted to free function templates so that the SpinConstrain + * class no longer has to carry the lambda-loop workflow as part of its + * interface. The complex-double specialization is provided by instantiating + * the template with TK = std::complex; a double stub is provided + * elsewhere to keep the link surface stable. */ +namespace spinconstrain +{ + /** - * @brief Print final spin and lambda values when lambda loop terminates. + * @brief Print final spin and lambda values when the lambda loop terminates. * * @par Output - * - "after-optimization spin (uB)": Final magnetic moments Mi for each atom - * - "after-optimization lambda (eV/uB)": Final Lagrange multipliers for each atom + * - "after-optimization spin (uB)": Final magnetic moments Mi per atom + * - "after-optimization lambda (eV/uB)": Final Lagrange multipliers per atom * - "Inner optimization for lambda ends.": Termination marker - * - * @par Interpretation - * - Mi close to target_mag: constraint successfully satisfied - * - Mi far from target_mag: constraint not converged (check RMS error in log) - * - lambda ≈ 0: system naturally has the target moment - * - lambda large: system resists the constraint (may indicate unrealistic target) */ -template <> -void spinconstrain::SpinConstrain>::print_termination() +template +void print_termination(const SpinConstrain& sc, std::ostream& ofs_running) { - print_2d("after-optimization spin (uB): (print in the inner loop): ", this->Mi_, this->nspin_); - print_2d("after-optimization lambda (eV/uB): (print in the inner loop): ", this->lambda_, this->nspin_, ModuleBase::Ry_to_eV); - std::cout << "Inner optimization for lambda ends." << std::endl; - std::cout << "===============================================================================" << std::endl; + print_2d(" after-optimization spin (uB): (print in the inner loop): ", sc.get_Mi(), sc.get_nspin(), 1.0, ofs_running); + print_2d(" after-optimization lambda (eV/uB): (print in the inner loop): ", + sc.get_sc_lambda(), + sc.get_nspin(), + ModuleBase::Ry_to_eV, + ofs_running); + ofs_running << " Inner optimization for lambda ends." << std::endl; + ofs_running << " ================================================================================" << std::endl; } /** @@ -50,47 +63,49 @@ void spinconstrain::SpinConstrain>::print_termination() * @par Return value * - true: loop should terminate (either converged or max steps) * - false: continue optimization - * - * @param outer_step Current SCF outer iteration - * @param i_step Current inner lambda step - * @param rms_error Current RMS error of Mi - M_target - * @param duration Time for this step - * @param total_duration Cumulative time for inner loop */ -template <> -bool spinconstrain::SpinConstrain>::check_rms_stop(int outer_step, - int i_step, - double rms_error, - double duration, - double total_duration) +template +bool check_rms_stop(const SpinConstrain& sc, + int outer_step, + int i_step, + double rms_error, + double duration, + double total_duration, + std::ostream& ofs_running) { - std::cout << "Step (Outer -- Inner) = " << outer_step << " -- " << std::left << std::setw(5) << i_step + 1 - << " RMS = " << rms_error << " TIME(s) = " << std::setw(11) << duration << std::endl; - if (rms_error < this->current_sc_thr_ || i_step == this->nsc_ - 1) + ofs_running << " Step (Outer -- Inner) = " << outer_step << " -- " << std::left << std::setw(5) << i_step + 1 + << " RMS = " << rms_error << " TIME(s) = " << std::setw(11) << duration << std::endl; + const double current_sc_thr = sc.get_current_sc_thr(); + const int nsc = sc.get_nsc(); + if (rms_error < current_sc_thr || i_step == nsc - 1) { - if (rms_error < this->current_sc_thr_) + if (rms_error < current_sc_thr) { - std::cout << "Meet convergence criterion ( < " << this->current_sc_thr_ << " ), exit."; - std::cout << " Total TIME(s) = " << total_duration << std::endl; + ofs_running << " DeltaSpin: lambda loop converged ( RMS < " << current_sc_thr + << " ), inner steps = " << (i_step + 1) + << ", Total TIME(s) = " << total_duration << std::endl; + ofs_running << std::endl; } - else if (i_step == this->nsc_ - 1) + else if (i_step == nsc - 1) { - std::cout << "Reach maximum number of steps ( " << this->nsc_ << " ), exit."; - std::cout << " Total TIME(s) = " << total_duration << std::endl; + std::cout << " DeltaSpin: lambda loop reached max steps ( " << nsc + << " ), RMS = " << rms_error + << ", Total TIME(s) = " << total_duration << std::endl; + std::cout << std::endl; } - this->print_termination(); + print_termination(sc, ofs_running); return true; } return false; } /// @brief Print header at start of lambda optimization loop -template <> -void spinconstrain::SpinConstrain>::print_header() +template +void print_header(const SpinConstrain& sc, std::ostream& ofs_running) { - std::cout << "===============================================================================" << std::endl; - std::cout << "Inner optimization for lambda begins ..." << std::endl; - std::cout << "Covergence criterion for the iteration: " << this->sc_thr_ << std::endl; + ofs_running << " ================================================================================" << std::endl; + ofs_running << " Inner optimization for lambda begins ..." << std::endl; + ofs_running << " Covergence criterion for the iteration: " << sc.get_sc_thr() << std::endl; } /** @@ -105,23 +120,22 @@ void spinconstrain::SpinConstrain>::print_header() * @par Output (when restriction is applied) * - "alpha after restrict = X eV/uB^2": The capped step size * - "boundary after = X eV/uB": The actual maximum lambda change - * - * @param search Current search direction - * @param alpha_trial Trial step size (modified in place if capped) */ -template <> -void spinconstrain::SpinConstrain>::check_restriction( - const std::vector>& search, - double& alpha_trial) +template +void check_restriction(const SpinConstrain& sc, + const std::vector>& search, + double& alpha_trial, + std::ostream& ofs_running) { + const double restrict_current = sc.get_sccut(); double boundary = std::abs(alpha_trial * maxval_abs_2d(search)); - if (this->restrict_current_ > 0 && boundary > this->restrict_current_) + if (restrict_current > 0 && boundary > restrict_current) { - alpha_trial = copysign(1.0, alpha_trial) * this->restrict_current_ / maxval_abs_2d(search); + alpha_trial = copysign(1.0, alpha_trial) * restrict_current / maxval_abs_2d(search); boundary = std::abs(alpha_trial * maxval_abs_2d(search)); - std::cout << "alpha after restrict = " << alpha_trial * ModuleBase::Ry_to_eV << std::endl; - std::cout << "boundary after = " << boundary * ModuleBase::Ry_to_eV << std::endl; + ofs_running << " alpha after restrict = " << alpha_trial * ModuleBase::Ry_to_eV << std::endl; + ofs_running << " boundary after = " << boundary * ModuleBase::Ry_to_eV << std::endl; } } @@ -138,26 +152,17 @@ void spinconstrain::SpinConstrain>::check_restriction( * sum_k = sum((target - spin) . (spin_plus - spin)) over constrained components * sum_k2 = sum(|spin - spin_plus|^2) over constrained components * - * This is equivalent to finding the minimum of a quadratic approximation - * to E(lambda) along the search direction. - * * @par Edge case handling * - If |sum_k2| < 1e-30: spin and spin_plus are nearly identical, meaning * the lambda change has no effect on Mi. Return alpha_trial as fallback. - * This can happen if the system is already saturated or if lambda is too small. - * - * @param spin Mi at current lambda - * @param spin_plus Mi at trial lambda (current + alpha_trial * search) - * @param alpha_trial Current trial step size - * @return Optimal step size alpha_opt */ -template <> -double spinconstrain::SpinConstrain>::cal_alpha_opt( - std::vector> spin, - std::vector> spin_plus, - const double alpha_trial) +template +double cal_alpha_opt(const SpinConstrain& sc, + std::vector> spin, + std::vector> spin_plus, + const double alpha_trial) { - int nat = this->get_nat(); + int nat = sc.get_nat(); const bool print = false; const double zero = 0.0; @@ -167,9 +172,9 @@ double spinconstrain::SpinConstrain>::cal_alpha_opt( std::vector> spin_plus_mask(nat, 0.0); std::vector> temp_1(nat, 0.0); std::vector> temp_2(nat, 0.0); - where_fill_scalar_else_2d(this->constrain_, 0, zero, this->target_mag_, target_spin_mask); - where_fill_scalar_else_2d(this->constrain_, 0, zero, spin, spin_mask); - where_fill_scalar_else_2d(this->constrain_, 0, zero, spin_plus, spin_plus_mask); + where_fill_scalar_else_2d(sc.get_constrain(), 0, zero, sc.get_target_mag(), target_spin_mask); + where_fill_scalar_else_2d(sc.get_constrain(), 0, zero, spin, spin_mask); + where_fill_scalar_else_2d(sc.get_constrain(), 0, zero, spin_plus, spin_plus_mask); // Compute dot products for linear interpolation for (int ia = 0; ia < nat; ia++) @@ -227,33 +232,20 @@ double spinconstrain::SpinConstrain>::cal_alpha_opt( * * @par Output (when triggered) * "Reach limitation of current step ( maximum gradient < X uB^2/eV in atom type Y ), exit." - * - * @par Debug output [GRAD-DECAY] - * - WARNING: nu_change too small: indicates delta_lambda and dnu_last_step are - * nearly identical, meaning the optimizer is not making progress. This can happen - * if alpha_trial has become very small or if the search direction is nearly zero. - * Solution: check that alpha_trial is not vanishing; increase sc_thr if target - * is physically unreachable. - * - * @param new_spin Mi at current lambda - * @param spin Mi at previous lambda - * @param delta_lambda Current lambda change - * @param dnu_last_step Previous cumulative step - * @param print Whether to print detailed gradient info - * @return true if gradient decayed below threshold (should terminate), false otherwise */ -template <> -bool spinconstrain::SpinConstrain>::check_gradient_decay( - std::vector> new_spin, - std::vector> spin, - std::vector> delta_lambda, - std::vector> dnu_last_step, - bool print) +template +bool check_gradient_decay(const SpinConstrain& sc, + std::vector> new_spin, + std::vector> spin, + std::vector> delta_lambda, + std::vector> dnu_last_step, + bool print, + std::ostream& ofs_running) { const double one = 1.0; const double zero = 0.0; - int nat = this->get_nat(); - int ntype = this->get_ntype(); + int nat = sc.get_nat(); + int ntype = sc.get_ntype(); // Change in magnetic moments and lambda std::vector> spin_change(nat, 0.0); @@ -273,9 +265,10 @@ bool spinconstrain::SpinConstrain>::check_gradient_decay( subtract_2d(new_spin, spin, spin_change); subtract_2d(delta_lambda, dnu_last_step, nu_change); + const auto& constrain = sc.get_constrain(); // Mask unconstrained components - where_fill_scalar_2d(this->constrain_, 0, zero, spin_change); - where_fill_scalar_2d(this->constrain_, 0, one, nu_change); + where_fill_scalar_2d(constrain, 0, zero, spin_change); + where_fill_scalar_2d(constrain, 0, one, nu_change); // Calculate full gradient matrix for (int ia = 0; ia < nat; ia++) @@ -299,8 +292,12 @@ bool spinconstrain::SpinConstrain>::check_gradient_decay( } } + const auto& atom_counts = sc.get_atomCounts(); + const auto& decay_grad = sc.get_decay_grad(); + const int nspin = sc.get_nspin(); + // Extract diagonal gradient and find max per atom type - for (const auto& sc_elem: this->get_atomCounts()) + for (const auto& sc_elem: atom_counts) { int it = sc_elem.first; int nat_it = sc_elem.second; @@ -322,29 +319,174 @@ bool spinconstrain::SpinConstrain>::check_gradient_decay( if (print) { - print_2d("diagonal gradient: ", spin_nu_gradient_diag, this->nspin_); - std::cout << "maximum gradient appears at: " << std::endl; + print_2d(" diagonal gradient: ", spin_nu_gradient_diag, nspin, 1.0, ofs_running); + ofs_running << " maximum gradient appears at: " << std::endl; for (int it = 0; it < ntype; it++) { - std::cout << "( " << max_gradient_index[it].first << ", " << max_gradient_index[it].second << " )" - << std::endl; + ofs_running << " ( " << max_gradient_index[it].first << ", " << max_gradient_index[it].second << " )" + << std::endl; } - std::cout << "maximum gradient: " << std::endl; + ofs_running << " maximum gradient: " << std::endl; for (int it = 0; it < ntype; it++) { - std::cout << max_gradient[it]/ModuleBase::Ry_to_eV << std::endl; + ofs_running << " " << max_gradient[it]/ModuleBase::Ry_to_eV << std::endl; } } // Check if any atom type's gradient has decayed below threshold for (int it = 0; it < ntype; it++) { - if (this->decay_grad_[it] > 0 && std::abs(max_gradient[it]) < this->decay_grad_[it]) + if (decay_grad[it] > 0 && std::abs(max_gradient[it]) < decay_grad[it]) { - std::cout << "Reach limitation of current step ( maximum gradient < " << this->decay_grad_[it]/ModuleBase::Ry_to_eV // uB^2/Ry to uB^2/eV - << " in atom type " << it << " ), exit." << std::endl; + std::cout << " DeltaSpin: lambda loop early-terminated ( maximum gradient < " + << decay_grad[it]/ModuleBase::Ry_to_eV // uB^2/Ry to uB^2/eV + << " in atom type " << it << " )" << std::endl; + std::cout << std::endl; return true; } } return false; } + +/** + * @brief Print atomic magnetic moments Mi in a formatted table. + * + * @par Output format + * - nspin=2: "Total Magnetism (uB)" with single z-component column + * - nspin=4: three columns (Mx, My, Mz) + * + * Lifted from SpinConstrain::print_Mi; accesses state via getters. + */ +template +void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) +{ + sc.check_atomCounts(); + const int nat = sc.get_nat(); + const int nspin = sc.get_nspin(); + const auto& Mi = sc.get_Mi(); + const auto& atomLabel = sc.get_atomLabels(); + std::vector mag_x(nat, 0.0); + std::vector mag_y(nat, 0.0); + std::vector mag_z(nat, 0.0); + if (nspin == 2) + { + const std::vector title = {"Total Magnetism (uB)", ""}; + const std::vector fmts = {"%-26s", "%20.10f"}; + FmtTable table(/*titles=*/title, + /*nrows=*/nat, + /*formats=*/fmts, + /*indent=*/0, + /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); + for (int iat = 0; iat < nat; ++iat) + { + mag_z[iat] = Mi[iat].z; + } + table << atomLabel << mag_z; + ofs_running << table.str() << std::endl; + } + else if (nspin == 4) + { + const std::vector title = {"Total Magnetism (uB)", "", "", ""}; + const std::vector fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"}; + FmtTable table(/*titles=*/title, + /*nrows=*/nat, + /*formats=*/fmts, + /*indent=*/0, + /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); + for (int iat = 0; iat < nat; ++iat) + { + mag_x[iat] = Mi[iat].x; + mag_y[iat] = Mi[iat].y; + mag_z[iat] = Mi[iat].z; + } + table << atomLabel << mag_x << mag_y << mag_z; + ofs_running << table.str() << std::endl; + } +} + +/** + * @brief Print the magnetic force (-lambda) per atom in eV/uB. + * + * Lifted from SpinConstrain::print_Mag_Force; accesses state via getters. + * lambda is read via get_sc_lambda() and converted from Ry to eV on output. + */ +template +void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) +{ + sc.check_atomCounts(); + const int nat = sc.get_nat(); + const int nspin = sc.get_nspin(); + const auto& lambda = sc.get_sc_lambda(); + const auto& atomLabel = sc.get_atomLabels(); + std::vector mag_force_x(nat, 0.0); + std::vector mag_force_y(nat, 0.0); + std::vector mag_force_z(nat, 0.0); + if (nspin == 2) + { + const std::vector title = {"Magnetic force (eV/uB)", ""}; + const std::vector fmts = {"%-26s", "%20.10f"}; + FmtTable table(/*titles=*/title, + /*nrows=*/nat, + /*formats=*/fmts, + /*indent=*/0, + /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); + for (int iat = 0; iat < nat; ++iat) + { + mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; + } + table << atomLabel << mag_force_z; + ofs_running << table.str() << std::endl; + } + else if (nspin == 4) + { + const std::vector title = {"Magnetic force (eV/uB)", "", "", ""}; + const std::vector fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"}; + FmtTable table(/*titles=*/title, + /*nrows=*/nat, + /*formats=*/fmts, + /*indent=*/0, + /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); + for (int iat = 0; iat < nat; ++iat) + { + mag_force_x[iat] = lambda[iat].x * ModuleBase::Ry_to_eV; + mag_force_y[iat] = lambda[iat].y * ModuleBase::Ry_to_eV; + mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; + } + table << atomLabel << mag_force_x << mag_force_y << mag_force_z; + ofs_running << table.str() << std::endl; + } +} + + +// Explicit instantiation for the only supported TK = std::complex. +// The double stub is provided by template_helpers.cpp via the existing +// specialization mechanism (kept as a separate file to avoid duplicate symbols). +template void print_termination>(const SpinConstrain>&, std::ostream&); +template bool check_rms_stop>(const SpinConstrain>&, + int, int, double, double, double, std::ostream&); +template void print_header>(const SpinConstrain>&, std::ostream&); +template void check_restriction>(const SpinConstrain>&, + const std::vector>&, + double&, std::ostream&); +template double cal_alpha_opt>(const SpinConstrain>&, + std::vector>, + std::vector>, + const double); +template bool check_gradient_decay>(const SpinConstrain>&, + std::vector>, + std::vector>, + std::vector>, + std::vector>, + bool, std::ostream&); + +// print_Mi / print_Mag_Force are generic (no per-TK stub needed): the +// template body works for both TK = std::complex (real usage) and +// TK = double (nspin=2 stub path instantiated by template_helpers_test). +// We instantiate both so callers in PW/LCAO ESolvers that hold either TK +// can link against a single definition. +template void print_Mi>(const SpinConstrain>&, std::ostream&); +template void print_Mag_Force>(const SpinConstrain>&, std::ostream&); +template void print_Mi(const SpinConstrain&, std::ostream&); +template void print_Mag_Force(const SpinConstrain&, std::ostream&); + +} // namespace spinconstrain diff --git a/source/source_lcao/module_deltaspin/lambda_loop_helper.h b/source/source_lcao/module_deltaspin/lambda_loop_helper.h new file mode 100644 index 00000000000..c980452d333 --- /dev/null +++ b/source/source_lcao/module_deltaspin/lambda_loop_helper.h @@ -0,0 +1,165 @@ +#ifndef LAMBDA_LOOP_HELPER_H +#define LAMBDA_LOOP_HELPER_H + +#include +#include + +#include "source_base/vector3.h" +#include "spin_constrain.h" + +/** + * @file lambda_loop_helper.h + * @brief Free-function helpers for the DeltaSpin lambda optimization loop. + * + * @par Background + * Originally these routines were member functions of SpinConstrain. + * They have been lifted to free functions in the spinconstrain namespace to + * shrink the SpinConstrain god class. Each helper takes the SpinConstrain + * instance as its first parameter and accesses internal state through the + * public getters (get_Mi, get_sc_lambda, get_constrain, ...). + * + * @par Template parameter TK + * - std::complex: full implementation (nspin=2 and nspin=4) + * - double: stub specialization (no-ops) provided elsewhere + */ + +namespace spinconstrain +{ + +/** + * @brief Print final spin and lambda values when the lambda loop terminates. + * + * @param sc SpinConstrain instance (read Mi_ and lambda_) + * @param ofs_running Log output stream + */ +template +void print_termination(const SpinConstrain& sc, std::ostream& ofs_running); + +/** + * @brief Check whether RMS error is below the convergence threshold or the + * maximum number of inner steps has been reached. + * + * @param sc SpinConstrain instance + * @param outer_step Current SCF outer iteration + * @param i_step Current inner lambda step + * @param rms_error Current RMS error of Mi - M_target + * @param duration Wall time for this step (s) + * @param total_duration Cumulative wall time for the inner loop (s) + * @param ofs_running Log output stream + * @return true if the inner loop should terminate, false otherwise + */ +template +bool check_rms_stop(const SpinConstrain& sc, + int outer_step, + int i_step, + double rms_error, + double duration, + double total_duration, + std::ostream& ofs_running); + +/** + * @brief Print header at the start of the lambda optimization loop. + * + * @param sc SpinConstrain instance + * @param ofs_running Log output stream + */ +template +void print_header(const SpinConstrain& sc, std::ostream& ofs_running); + +/** + * @brief Cap the step size to prevent the optimizer from overshooting. + * + * @details If |alpha_trial * max(search)| exceeds restrict_current_, the + * trial step is reduced so that the maximum lambda change per step is + * bounded. alpha_trial is modified in place. + * + * @param sc SpinConstrain instance + * @param search Current search direction (per atom, 3 components) + * @param alpha_trial Trial step size, modified in place if capped + * @param ofs_running Log output stream + */ +template +void check_restriction(const SpinConstrain& sc, + const std::vector>& search, + double& alpha_trial, + std::ostream& ofs_running); + +/** + * @brief Compute the optimal step size via two-point linear interpolation. + * + * @par Algorithm + * alpha_opt = sum_k / sum_k2 * alpha_trial + * where + * sum_k = sum((target - spin) . (spin_plus - spin)) over constrained components + * sum_k2 = sum(|spin - spin_plus|^2) over constrained components + * + * @param sc SpinConstrain instance + * @param spin Mi at current lambda + * @param spin_plus Mi at trial lambda (current + alpha_trial * search) + * @param alpha_trial Current trial step size + * @return Optimal step size; falls back to alpha_trial when sum_k2 ~ 0 + */ +template +double cal_alpha_opt(const SpinConstrain& sc, + std::vector> spin, + std::vector> spin_plus, + const double alpha_trial); + +/** + * @brief Check whether the magnetic susceptibility gradient dM/dlambda has + * decayed below the per-atom-type threshold. + * + * @par Algorithm + * 1. Compute spin_change = new_spin - spin + * 2. Compute nu_change = delta_lambda - dnu_last_step + * 3. Build full gradient matrix dM[ia][ic]/dlambda[ja][jc] + * 4. Extract diagonal; pick max abs per atom type + * 5. Return true if max(|diag|) < decay_grad[itype] for any type + * + * @param sc SpinConstrain instance + * @param new_spin Mi at current lambda + * @param spin Mi at previous lambda + * @param delta_lambda Current lambda change + * @param dnu_last_step Previous cumulative step + * @param print Whether to print detailed gradient info + * @param ofs_running Log output stream + * @return true if gradient decayed below threshold for any atom type + */ +template +bool check_gradient_decay(const SpinConstrain& sc, + std::vector> new_spin, + std::vector> spin, + std::vector> delta_lambda, + std::vector> dnu_last_step, + bool print, + std::ostream& ofs_running); + +/** + * @brief Print atomic magnetic moments Mi in a formatted table. + * + * @par Output format + * - nspin=2: "Total Magnetism (uB)" with single z-component column + * - nspin=4: three columns (Mx, My, Mz) + * + * @param sc SpinConstrain instance + * @param ofs_running Log output stream + */ +template +void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running); + +/** + * @brief Print the magnetic force (-lambda) per atom in eV/uB. + * + * @par Physical meaning + * Magnetic force = dL/dMi = -lambda. Large |lambda| means the system + * strongly resists the target moment constraint. + * + * @param sc SpinConstrain instance + * @param ofs_running Log output stream + */ +template +void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running); + +} // namespace spinconstrain + +#endif // LAMBDA_LOOP_HELPER_H diff --git a/source/source_lcao/module_deltaspin/mi_tools.cpp b/source/source_lcao/module_deltaspin/mi_tools.cpp new file mode 100644 index 00000000000..8e38b05e29b --- /dev/null +++ b/source/source_lcao/module_deltaspin/mi_tools.cpp @@ -0,0 +1,92 @@ +#include "mi_tools.h" + +#include "source_base/tool_quit.h" + +namespace spinconstrain +{ + +void accumulate_Mi_from_becp(const std::complex* becp, + int nkb, + int nbands, + int npol, + int spin_sign, + const double* wg_ik, + const int* nh_iat, + std::vector>& mi) +{ + if (becp == nullptr) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "becp is nullptr"); + } + if (wg_ik == nullptr) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "wg_ik is nullptr"); + } + if (nh_iat == nullptr) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "nh_iat is nullptr"); + } + if (nkb <= 0) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "nkb must be positive"); + } + if (nbands <= 0) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "nbands must be positive"); + } + if (npol != 1 && npol != 2) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "npol must be 1 or 2"); + } + if (spin_sign != -1 && spin_sign != 1) + { + ModuleBase::WARNING_QUIT("accumulate_Mi_from_becp", "spin_sign must be -1 or 1"); + } + + const std::complex zero(0.0, 0.0); + if (npol == 2) + { + for (int ib = 0; ib < nbands; ib++) + { + const double weight = wg_ik[ib]; + int begin_ih = 0; + for (int iat = 0; iat < static_cast(mi.size()); iat++) + { + std::complex occ[4] = {zero, zero, zero, zero}; + const int nh = nh_iat[iat]; + for (int ih = 0; ih < nh; ih++) + { + const int index = ib * 2 * nkb + begin_ih + ih; + occ[0] += conj(becp[index]) * becp[index]; + occ[1] += conj(becp[index]) * becp[index + nkb]; + occ[2] += conj(becp[index + nkb]) * becp[index]; + occ[3] += conj(becp[index + nkb]) * becp[index + nkb]; + } + mi[iat] += pauli_to_moment(occ, weight); + begin_ih += nh; + } + } + } + else // npol == 1 + { + for (int ib = 0; ib < nbands; ib++) + { + const double weight = wg_ik[ib]; + int begin_ih = 0; + for (int iat = 0; iat < static_cast(mi.size()); iat++) + { + double occ = 0.0; + const int nh = nh_iat[iat]; + for (int ih = 0; ih < nh; ih++) + { + const int index = ib * nkb + begin_ih + ih; + occ += (conj(becp[index]) * becp[index]).real(); + } + mi[iat].z += weight * occ * spin_sign; + begin_ih += nh; + } + } + } +} + +} // namespace spinconstrain diff --git a/source/source_lcao/module_deltaspin/mi_tools.h b/source/source_lcao/module_deltaspin/mi_tools.h new file mode 100644 index 00000000000..5130aeb5cbd --- /dev/null +++ b/source/source_lcao/module_deltaspin/mi_tools.h @@ -0,0 +1,81 @@ +#ifndef MI_TOOLS_H +#define MI_TOOLS_H + +#include +#include + +#include "source_base/vector3.h" + +/** + * @file mi_tools.h + * @brief Free-function utilities for computing atomic magnetic moments (Mi) + * from intermediate quantities (e.g. becp projector coefficients) + * in the DeltaSpin module. + * + * @par Rationale + * These functions do not depend on the SpinConstrain template state and + * are factored out as free functions to: + * - Reduce coupling to the singleton class members + * - Allow reuse by both LCAO and PW basis paths + * - Be independently unit-testable + */ + +namespace spinconstrain +{ + +/** + * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices. + * + * @details For a two-component spinor wavefunction, the spin density matrix is: + * rho = |a|^2 a*b | = | (1+Mz)/2 (Mx-iMy)/2 | + * |b*a |b|^2 | | (Mx+iMy)/2 (1-Mz)/2 | + * The magnetic moment components are extracted via Pauli matrix traces: + * Mx = Tr(rho * sigma_x) = occ[1] + occ[2] (real part) + * My = Tr(rho * sigma_y) = -Im(occ[1] - occ[2]) (from sigma_y = [[0,-i],[i,0]]) + * Mz = Tr(rho * sigma_z) = occ[0] - occ[3] (real part) + * where occ = {|a|^2, a*b, b*a, |b|^2} from becp coefficients. + * + * @param occ 4-element array of occupation matrix elements (complex) + * @param weight k-point weight for integration + * @return 3D magnetic moment vector (Mx, My, Mz) in Bohr magnetons + */ +inline ModuleBase::Vector3 pauli_to_moment(const std::complex occ[4], double weight) +{ + return ModuleBase::Vector3( + weight * (occ[1] + occ[2]).real(), + weight * (occ[1] - occ[2]).imag(), + weight * (occ[0] - occ[3]).real() + ); +} + +/** + * @brief Accumulate atomic magnetic moments from becp coefficients for one k-point. + * + * @details For npol=2 (nspin=4), computes full Pauli decomposition: + * occ[0] = sum(becp_up^* * becp_up), occ[1] = sum(becp_up^* * becp_dn), + * occ[2] = sum(becp_dn^* * becp_up), occ[3] = sum(becp_dn^* * becp_dn) + * Mi = pauli_to_moment(occ, weight) + * For npol=1 (nspin=2), only z-component: + * occ = sum(|becp|^2), Mi.z += weight * occ * spin_sign + * + * @param becp Projector coefficients + * @param nkb Total number of projectors + * @param nbands Number of bands + * @param npol Number of spinor components (1 for collinear, 2 for non-collinear) + * @param spin_sign +1 for spin-up, -1 for spin-down (nspin=2 only); unused for npol=2 + * @param wg_ik Band occupation weights for this k-point (from Fermi-Dirac) + * @param nh_iat Array of projector counts per atom: nh_iat[iat] = nproj for atom iat + * @param mi [in,out] Magnetic moments vector to accumulate into (size = nat) + */ +void accumulate_Mi_from_becp(const std::complex* becp, + int nkb, + int nbands, + int npol, + int spin_sign, + const double* wg_ik, + const int* nh_iat, + std::vector>& mi); + +} // namespace spinconstrain + +#endif // MI_TOOLS_H diff --git a/source/source_lcao/module_deltaspin/sc_parse_json.cpp b/source/source_lcao/module_deltaspin/sc_parse_json.cpp deleted file mode 100644 index 37f23fa3973..00000000000 --- a/source/source_lcao/module_deltaspin/sc_parse_json.cpp +++ /dev/null @@ -1,4 +0,0 @@ -#include "spin_constrain.h" - -template class spinconstrain::SpinConstrain>; -template class spinconstrain::SpinConstrain; diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp index 5950cc04d9e..fd9d1d83f08 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.cpp +++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp @@ -130,89 +130,6 @@ int SpinConstrain::get_spin_sign(int ik) const return (this->pelec->klist->isk[ik] == 0) ? 1 : -1; } -/** - * @brief Accumulate magnetic moments from projector coefficients (becp) for one k-point. - * - * @par Algorithm (npol=2, non-collinear): - * For each atom, compute the 2x2 occupation matrix from becp coefficients: - * occ[0] = sum_ih becp_up^*(ih) * becp_up(ih) = - * occ[1] = sum_ih becp_up^*(ih) * becp_dn(ih) = - * occ[2] = sum_ih becp_dn^*(ih) * becp_up(ih) = - * occ[3] = sum_ih becp_dn^*(ih) * becp_dn(ih) = - * where P_at = sum_{l,m} |alpha_{l,m}> - * Mz += weight * occ * spin_sign - * where spin_sign = +1 for spin-up, -1 for spin-down. - * - * @param becp Projector coefficients, layout: [ib * npol * nkb + spin * nkb + ih] - * @param nkb Total number of projectors across all atoms - * @param nbands Number of bands (occupied + unoccupied in the subspace) - * @param npol Number of spinor components (1 for collinear, 2 for non-collinear) - * @param ik K-point index (used for spin_sign lookup in collinear mode) - * @param wg_ik Band occupation weights for this k-point (from Fermi-Dirac) - * @param nh_iat Array of projector counts per atom: nh_iat[iat] = nproj for atom iat - */ -template -void SpinConstrain::accumulate_Mi_from_becp(const std::complex* becp, - int nkb, - int nbands, - int npol, - int ik, - const double* wg_ik, - const int* nh_iat) -{ - if (npol == 2) - { - for (int ib = 0; ib < nbands; ib++) - { - const double weight = wg_ik[ib]; - int begin_ih = 0; - for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) - { - std::complex occ[4] = {ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO}; - const int nh = nh_iat[iat]; - for (int ih = 0; ih < nh; ih++) - { - const int index = ib * 2 * nkb + begin_ih + ih; - occ[0] += conj(becp[index]) * becp[index]; - occ[1] += conj(becp[index]) * becp[index + nkb]; - occ[2] += conj(becp[index + nkb]) * becp[index]; - occ[3] += conj(becp[index + nkb]) * becp[index + nkb]; - } - this->Mi_[iat] += pauli_to_moment(occ, weight); - begin_ih += nh; - } - } - } - else // npol == 1 - { - const int sign = this->get_spin_sign(ik); - for (int ib = 0; ib < nbands; ib++) - { - const double weight = wg_ik[ib]; - int begin_ih = 0; - for (int iat = 0; iat < static_cast(this->Mi_.size()); iat++) - { - double occ = 0.0; - const int nh = nh_iat[iat]; - for (int ih = 0; ih < nh; ih++) - { - const int index = ib * nkb + begin_ih + ih; - occ += (conj(becp[index]) * becp[index]).real(); - } - this->Mi_[iat].z += weight * occ * sign; - begin_ih += nh; - } - } - } -} - template int SpinConstrain::get_nw() const { @@ -257,10 +174,10 @@ int SpinConstrain::get_iwt(int itype, int iat, int orbital_index) const /// @brief Get total number of atoms across all element types template -int SpinConstrain::get_nat() +int SpinConstrain::get_nat() const { int nat = 0; - for (std::map::iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) + for (std::map::const_iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) { nat += it->second; } @@ -269,7 +186,7 @@ int SpinConstrain::get_nat() /// @brief Get number of element types template -int SpinConstrain::get_ntype() +int SpinConstrain::get_ntype() const { return this->atomCounts.size(); } @@ -287,7 +204,7 @@ int SpinConstrain::get_ntype() * - "number of atoms <= 0": some element type has no atoms */ template -void SpinConstrain::check_atomCounts() +void SpinConstrain::check_atomCounts() const { if (!this->atomCounts.size()) { @@ -297,7 +214,7 @@ void SpinConstrain::check_atomCounts() { ModuleBase::WARNING_QUIT("SpinConstrain::check_atomCounts", "nat <= 0"); } - for (std::map::iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) + for (std::map::const_iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) { int itype = it->first; if (itype < 0 || itype >= this->get_ntype()) @@ -342,7 +259,7 @@ int SpinConstrain::get_iat(int itype, int atom_index) ModuleBase::WARNING_QUIT("SpinConstrain::get_iat", "atom index out of range [0, nat)"); } int iat = 0; - for (std::map::iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) + for (std::map::const_iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it) { if (it->first == itype) { @@ -618,9 +535,10 @@ void SpinConstrain::zero_Mi() /// this function can only be called by the root process because only /// root process reads the ScDecayGrad from json file template -double SpinConstrain::get_decay_grad(int itype) +double SpinConstrain::get_decay_grad(int itype) const { - return this->ScDecayGrad[itype]; + std::map::const_iterator it = this->ScDecayGrad.find(itype); + return it != this->ScDecayGrad.end() ? it->second : 0.0; } /// set grad_decy @@ -638,7 +556,7 @@ void SpinConstrain::set_decay_grad() /// get decay_grad template -const std::vector& SpinConstrain::get_decay_grad() +const std::vector& SpinConstrain::get_decay_grad() const { return this->decay_grad_; } @@ -684,6 +602,27 @@ double SpinConstrain::get_sc_thr() const return this->sc_thr_; } +/// get current adaptive sc threshold +template +double SpinConstrain::get_current_sc_thr() const +{ + return this->current_sc_thr_; +} + +/// get computed magnetic moments Mi per atom +template +const std::vector>& SpinConstrain::get_Mi() const +{ + return this->Mi_; +} + +/// get human-readable atom labels for table printing +template +const std::vector& SpinConstrain::get_atomLabels() const +{ + return this->atomLabels_; +} + /// get nsc template int SpinConstrain::get_nsc() const @@ -750,127 +689,6 @@ void SpinConstrain::set_ParaV(Parallel_Orbitals* ParaV_in) } } -/** - * @brief Print magnetic moments per atom in formatted table. - * - * @par Output format - * - nspin=2: "ATOM 1 2.0000000000" (z-component only) - * - nspin=4: "ATOM 1 0.0010000000 0.0020000000 1.9990000000" (x, y, z) - * - * @par Interpretation - * - Positive Mi.z: spin aligned with z-axis (spin-up character) - * - Negative Mi.z: spin anti-aligned with z-axis (spin-down character) - * - Non-zero Mi.x/Mi.y: non-collinear spin components - * - Mi close to target_mag: constraint is well-satisfied - * - Mi far from target_mag: constraint is not yet converged - */ -template -void SpinConstrain::print_Mi(std::ofstream& ofs_running) -{ - this->check_atomCounts(); - int nat = this->get_nat(); - std::vector mag_x(nat, 0.0); - std::vector mag_y(nat, 0.0); - std::vector mag_z(nat, 0.0); - if (this->nspin_ == 2) - { - const std::vector title = {"Total Magnetism (uB)", ""}; - const std::vector fmts = {"%-26s", "%20.10f"}; - FmtTable table(/*titles=*/title, - /*nrows=*/nat, - /*formats=*/fmts, - /*indent=*/0, - /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); - for (int iat = 0; iat < nat; ++iat) - { - mag_z[iat] = Mi_[iat].z; - } - table << this->atomLabels_ << mag_z; - ofs_running << table.str() << std::endl; - } - else if (this->nspin_ == 4) - { - const std::vector title = {"Total Magnetism (uB)", "", "", ""}; - const std::vector fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"}; - FmtTable table(/*titles=*/title, - /*nrows=*/nat, - /*formats=*/fmts, - /*indent=*/0, - /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); - for (int iat = 0; iat < nat; ++iat) - { - mag_x[iat] = Mi_[iat].x; - mag_y[iat] = Mi_[iat].y; - mag_z[iat] = Mi_[iat].z; - } - table << this->atomLabels_ << mag_x << mag_y << mag_z; - ofs_running << table.str() << std::endl; - } -} - -/** - * @brief Print the magnetic force (-lambda) per atom in eV/uB. - * - * @par Physical meaning - * The "magnetic force" is the derivative of the constrained Lagrangian - * with respect to the magnetic moment: dL/dMi = -lambda_i. - * It represents how much energy would change if the constraint were relaxed. - * - * @par Interpretation - * - Large |lambda|: The system strongly resists the target moment constraint - * - lambda ≈ 0: The system naturally has the target moment (no constraint needed) - * - Positive lambda.z: The constraint pushes the moment in the +z direction - * - Negative lambda.z: The constraint pushes the moment in the -z direction - * - * @par Typical values - * - Well-converged SCF: lambda ~ 0.01-1 eV/uB - * - Strongly constrained: lambda ~ 1-10 eV/uB - * - Diverging SCF: lambda growing without bound (check target_mag validity) - */ -template -void SpinConstrain::print_Mag_Force(std::ofstream& ofs_running) -{ - this->check_atomCounts(); - int nat = this->get_nat(); - std::vector mag_force_x(nat, 0.0); - std::vector mag_force_y(nat, 0.0); - std::vector mag_force_z(nat, 0.0); - if (this->nspin_ == 2) - { - const std::vector title = {"Magnetic force (eV/uB)", ""}; - const std::vector fmts = {"%-26s", "%20.10f"}; - FmtTable table(/*titles=*/title, - /*nrows=*/nat, - /*formats=*/fmts, - /*indent=*/0, - /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); - for (int iat = 0; iat < nat; ++iat) - { - mag_force_z[iat] = lambda_[iat].z * ModuleBase::Ry_to_eV; - } - table << this->atomLabels_ << mag_force_z; - ofs_running << table.str() << std::endl; - } - else if (this->nspin_ == 4) - { - const std::vector title = {"Magnetic force (eV/uB)", "", "", ""}; - const std::vector fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"}; - FmtTable table(/*titles=*/title, - /*nrows=*/nat, - /*formats=*/fmts, - /*indent=*/0, - /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); - for (int iat = 0; iat < nat; ++iat) - { - mag_force_x[iat] = lambda_[iat].x * ModuleBase::Ry_to_eV; - mag_force_y[iat] = lambda_[iat].y * ModuleBase::Ry_to_eV; - mag_force_z[iat] = lambda_[iat].z * ModuleBase::Ry_to_eV; - } - table << this->atomLabels_ << mag_force_x << mag_force_y << mag_force_z; - ofs_running << table.str() << std::endl; - } -} - /** * @brief Reset DeltaSpin operator initialization state. * diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index 474cef67165..cdac9a615ff 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -37,6 +37,7 @@ #define SPIN_CONSTRAIN_H #include +#include #include #include @@ -59,31 +60,6 @@ namespace spinconstrain { -/** - * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices. - * - * @details For a two-component spinor wavefunction, the spin density matrix is: - * rho = |a|^2 a*b | = | (1+Mz)/2 (Mx-iMy)/2 | - * |b*a |b|^2 | | (Mx+iMy)/2 (1-Mz)/2 | - * The magnetic moment components are extracted via Pauli matrix traces: - * Mx = Tr(rho * sigma_x) = occ[1] + occ[2] (real part) - * My = Tr(rho * sigma_y) = -Im(occ[1] - occ[2]) (from sigma_y = [[0,-i],[i,0]]) - * Mz = Tr(rho * sigma_z) = occ[0] - occ[3] (real part) - * where occ = {|a|^2, a*b, b*a, |b|^2} from becp coefficients. - * - * @param occ 4-element array of occupation matrix elements (complex) - * @param weight k-point weight for integration - * @return 3D magnetic moment vector (Mx, My, Mz) in Bohr magnetons - */ -inline ModuleBase::Vector3 pauli_to_moment(const std::complex occ[4], double weight) -{ - return ModuleBase::Vector3( - weight * (occ[1] + occ[2]).real(), - weight * (occ[1] - occ[2]).imag(), - weight * (occ[0] - occ[3]).real() - ); -} - struct ScAtomData; /** @@ -253,7 +229,11 @@ class SpinConstrain * @param rerun If true, use full PW solver for final charge update */ void run_lambda_loop(int outer_step, - bool rerun = true); + bool rerun, + std::ostream& ofs_running); + + /// @brief RMS error of the most recent lambda optimization loop (-1.0 if none has run). + double get_last_rms_error() const { return last_rms_error_; } /** * @brief Alternative mode: sweep lambda values linearly for energy landscape mapping. @@ -264,7 +244,7 @@ class SpinConstrain * * @param outer_step Current SCF outer iteration number */ - void run_lambda_linear_scan(int outer_step); + void run_lambda_linear_scan(int outer_step, std::ostream& ofs_running); /// @brief Reset DeltaSpin operator initialization state when constraints change void reset_dspin_operator(); @@ -334,40 +314,11 @@ class SpinConstrain int isk); #endif - /// Lambda loop helper: check if RMS error below threshold or max steps reached - bool check_rms_stop(int outer_step, int i_step, double rms_error, double duration, double total_duration); - - /// Lambda loop helper: cap step size via restrict_current_ to prevent overshooting - void check_restriction(const std::vector>& search, double& alpha_trial); - - /** - * @brief Lambda loop helper: check if dM/dlambda gradient has decayed below threshold. - * - * @details Computes the diagonal of the susceptibility matrix dM/dlambda for each - * atom type. If max gradient < decay_grad[itype], the lambda optimization has - * reached diminishing returns and should stop. - * - * @return true if gradient decayed below threshold, false otherwise - */ - bool check_gradient_decay(std::vector> new_spin, - std::vector> old_spin, - std::vector> new_delta_lambda, - std::vector> old_delta_lambda, - bool print = false); - /// @brief Lambda loop helper: calculate optimal step size via linear interpolation - double cal_alpha_opt(std::vector> spin, - std::vector> spin_plus, - const double alpha_trial); - /// Print header at start of lambda loop - void print_header(); - /// Print termination message with final spin and lambda values - void print_termination(); - - /// Print magnetic moments to output stream - void print_Mi(std::ofstream& ofs_running); - - /// Print magnetic force (defined as dL/dMi = -lambda[iat]) in eV/uB - void print_Mag_Force(std::ofstream& ofs_running); + /// Lambda loop helpers (print_rms_stop, check_restriction, check_gradient_decay, + /// cal_alpha_opt, print_header, print_termination) have been lifted to free + /// functions in lambda_loop_helper.h. The class now only carries state and + /// the core lambda-loop driver (run_lambda_loop). + /// print_Mi and print_Mag_Force have also been lifted to lambda_loop_helper.h. /// @brief Use full PW solver (rerun) for higher precision in lambda loop bool higher_mag_prec = false; @@ -444,11 +395,11 @@ class SpinConstrain /// get constrain const std::vector>& get_constrain() const; /// get nat - int get_nat(); + int get_nat() const; /// get ntype - int get_ntype(); + int get_ntype() const; /// check atomCounts - void check_atomCounts(); + void check_atomCounts() const; /// get iat int get_iat(int itype, int atom_index); /// set nspin @@ -458,11 +409,11 @@ class SpinConstrain /// zero atomic magnetic moment void zero_Mi(); /// get decay_grad - double get_decay_grad(int itype); + double get_decay_grad(int itype) const; /// set decay_grad void set_decay_grad(); /// get decay_grad - const std::vector& get_decay_grad(); + const std::vector& get_decay_grad() const; /// set decay_grad from variable void set_decay_grad(const double* decay_grad_in, int ntype_in); /// set decay grad switch @@ -476,6 +427,8 @@ class SpinConstrain double sc_drop_thr_in); /// get sc_thr double get_sc_thr() const; + /// get current adaptive sc threshold (max(initial_rms * sc_drop_thr_, sc_thr_)) + double get_current_sc_thr() const; /// get nsc int get_nsc() const; /// get nsc_min @@ -486,6 +439,10 @@ class SpinConstrain double get_sccut() const; /// get sc_drop_thr double get_sc_drop_thr() const; + /// get computed magnetic moments Mi per atom + const std::vector>& get_Mi() const; + /// get human-readable atom labels ("Fe_0", "Fe_1", ...) for table printing + const std::vector& get_atomLabels() const; /// @brief set orbital parallel info void set_ParaV(Parallel_Orbitals* ParaV_in); /// @brief set parameters for solver @@ -567,31 +524,6 @@ class SpinConstrain int get_iwt(int itype, int iat, int orbital_index) const; ///< Convert (itype, iat, iw) to global orbital index /// @brief Get spin sign for k-point ik: +1 for spin-up, -1 for spin-down (nspin=2 only) int get_spin_sign(int ik) const; - /** - * @brief Accumulate magnetic moments from becp coefficients for a single k-point. - * - * @details For npol=2 (nspin=4), computes full Pauli decomposition: - * occ[0] = sum(becp_up^* * becp_up), occ[1] = sum(becp_up^* * becp_dn), - * occ[2] = sum(becp_dn^* * becp_up), occ[3] = sum(becp_dn^* * becp_dn) - * Mi = pauli_to_moment(occ, weight) - * For npol=1 (nspin=2), only z-component: - * occ = sum(|becp|^2), Mi.z += weight * occ * spin_sign - * - * @param becp Projector coefficients - * @param nkb Total number of projectors - * @param nbands Number of bands - * @param npol Number of spinor components - * @param ik K-point index (for spin_sign lookup in nspin=2) - * @param wg_ik Band weights for this k-point - * @param nh_iat Number of projectors per atom - */ - void accumulate_Mi_from_becp(const std::complex* becp, - int nkb, - int nbands, - int npol, - int ik, - const double* wg_ik, - const int* nh_iat); private: /// DeltaSpin operator pointer for LCAO magnetic moment calculation hamilt::Operator* p_operator = nullptr; @@ -624,6 +556,10 @@ class SpinConstrain TK* sub_s_save = nullptr; ///< Cached subspace overlap matrix for all k-points TK* becp_save = nullptr; ///< Cached becp coefficients for all k-points std::vector> lambda_in_sub_; ///< Lambda values when subspace was saved + + /// RMS error of the most recent lambda optimization loop; -1.0 if no loop has run. + /// Used by ESolver to pass the current DeltaSpin RMS into the SCF iteration table. + double last_rms_error_ = -1.0; }; diff --git a/source/source_lcao/module_deltaspin/template_helpers.cpp b/source/source_lcao/module_deltaspin/template_helpers.cpp index add9ee6398c..4c25ccd4cee 100644 --- a/source/source_lcao/module_deltaspin/template_helpers.cpp +++ b/source/source_lcao/module_deltaspin/template_helpers.cpp @@ -1,3 +1,4 @@ +#include "lambda_loop_helper.h" #include "spin_constrain.h" /** @@ -16,24 +17,17 @@ * - Simple getters/setters (nspin, npol, atomCounts, etc.) still work correctly * because they are not template-specialized (they use the base template) * - * @par Methods stubbed - * - cal_mw_from_lambda: no-op (computed by complex specialization) - * - cal_mi_lcao: no-op (computed by complex specialization) - * - run_lambda_loop: no-op (computed by complex specialization) - * - check_rms_stop: returns false (continue loop) - * - check_restriction: no-op - * - cal_alpha_opt: returns 0.0 - * - print_termination: no-op - * - print_header: no-op - * - check_gradient_decay: returns false (no early termination) - * - run_lambda_linear_scan: no-op - * - reset_dspin_operator: no-op + * @par History + * The following routines were originally member-function stubs of + * SpinConstrain. They have been lifted to free-function stubs in + * the spinconstrain namespace alongside their complex counterparts + * (see lambda_loop_helper.cpp). Each returns the safe no-op value. */ /// @brief cal_mw_from_lambda stub (TK=double): no-op template <> void spinconstrain::SpinConstrain::cal_mw_from_lambda(int i_step, - const ModuleBase::Vector3* delta_lambda) + const ModuleBase::Vector3* delta_lambda) { } @@ -46,65 +40,73 @@ void spinconstrain::SpinConstrain::cal_mi_lcao(const int& step, bool pri /// @brief run_lambda_loop stub (TK=double): no-op template <> void spinconstrain::SpinConstrain::run_lambda_loop(int outer_step, - bool rerun) + bool rerun, + std::ostream& ofs_running) { } /// @brief check_rms_stop stub (TK=double): always return false (continue) template <> -bool spinconstrain::SpinConstrain::check_rms_stop(int outer_step, - int i_step, - double rms_error, - double duration, - double total_duration) +bool spinconstrain::check_rms_stop(const spinconstrain::SpinConstrain& sc, + int outer_step, + int i_step, + double rms_error, + double duration, + double total_duration, + std::ostream& ofs_running) { return false; } /// @brief check_restriction stub (TK=double): no-op template <> -void spinconstrain::SpinConstrain::check_restriction( - const std::vector>& search, - double& alpha_trial) +void spinconstrain::check_restriction(const spinconstrain::SpinConstrain& sc, + const std::vector>& search, + double& alpha_trial, + std::ostream& ofs_running) { } /// @brief cal_alpha_opt stub (TK=double): return 0.0 template <> -double spinconstrain::SpinConstrain::cal_alpha_opt(std::vector> spin, - std::vector> spin_plus, - const double alpha_trial) +double spinconstrain::cal_alpha_opt(const spinconstrain::SpinConstrain& sc, + std::vector> spin, + std::vector> spin_plus, + const double alpha_trial) { return 0.0; } /// @brief print_termination stub (TK=double): no-op template <> -void spinconstrain::SpinConstrain::print_termination() +void spinconstrain::print_termination(const spinconstrain::SpinConstrain& sc, + std::ostream& ofs_running) { } /// @brief print_header stub (TK=double): no-op template <> -void spinconstrain::SpinConstrain::print_header() +void spinconstrain::print_header(const spinconstrain::SpinConstrain& sc, + std::ostream& ofs_running) { } /// @brief check_gradient_decay stub (TK=double): always return false (no early termination) template <> -bool spinconstrain::SpinConstrain::check_gradient_decay( - std::vector> new_spin, - std::vector> old_spin, - std::vector> new_delta_lambda, - std::vector> old_delta_lambda, - bool print) +bool spinconstrain::check_gradient_decay(const spinconstrain::SpinConstrain& sc, + std::vector> new_spin, + std::vector> old_spin, + std::vector> new_delta_lambda, + std::vector> old_delta_lambda, + bool print, + std::ostream& ofs_running) { return false; } /// @brief run_lambda_linear_scan stub (TK=double): no-op template <> -void spinconstrain::SpinConstrain::run_lambda_linear_scan(int outer_step) +void spinconstrain::SpinConstrain::run_lambda_linear_scan(int outer_step, std::ostream& ofs_running) { } diff --git a/source/source_lcao/module_deltaspin/test/CMakeLists.txt b/source/source_lcao/module_deltaspin/test/CMakeLists.txt index c47869123a5..65504369d13 100644 --- a/source/source_lcao/module_deltaspin/test/CMakeLists.txt +++ b/source/source_lcao/module_deltaspin/test/CMakeLists.txt @@ -23,6 +23,8 @@ AddTest( SOURCES template_helpers_test.cpp ../spin_constrain.cpp ../template_helpers.cpp + ../lambda_loop_helper.cpp + ../basic_funcs.cpp ) AddTest( diff --git a/source/source_lcao/module_deltaspin/test/basic_test.cpp b/source/source_lcao/module_deltaspin/test/basic_test.cpp index 0c1c3d0c9b8..3c30dcdddd3 100644 --- a/source/source_lcao/module_deltaspin/test/basic_test.cpp +++ b/source/source_lcao/module_deltaspin/test/basic_test.cpp @@ -181,14 +181,14 @@ TEST_F(BasicFuncsTest, Prin2d) print_2d(info, array, 4); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output,testing::HasSubstr("initial spin")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 1 1.0000000000 2.0000000000 3.0000000000")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 2 4.0000000000 5.0000000000 6.0000000000")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 3 7.0000000000 8.0000000000 9.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 1 1.0000000000 2.0000000000 3.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 2 4.0000000000 5.0000000000 6.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 3 7.0000000000 8.0000000000 9.0000000000")); testing::internal::CaptureStdout(); print_2d(info, array, 2); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output,testing::HasSubstr("initial spin")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 1 3.0000000000")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 2 6.0000000000")); - EXPECT_THAT(output,testing::HasSubstr("ATOM 3 9.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 1 3.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 2 6.0000000000")); + EXPECT_THAT(output,testing::HasSubstr(" ATOM 3 9.0000000000")); } \ No newline at end of file diff --git a/source/source_lcao/module_deltaspin/test/template_helpers_test.cpp b/source/source_lcao/module_deltaspin/test/template_helpers_test.cpp index 900aa56f62b..bc5b71ef95e 100644 --- a/source/source_lcao/module_deltaspin/test/template_helpers_test.cpp +++ b/source/source_lcao/module_deltaspin/test/template_helpers_test.cpp @@ -1,6 +1,7 @@ #include #include +#include "source_lcao/module_deltaspin/lambda_loop_helper.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -30,13 +31,13 @@ TEST_F(SpinConstrainTest, TemplatHelpers) std::vector> Sloc2; EXPECT_NO_THROW(sc.cal_mw_from_lambda(0)); EXPECT_NO_THROW(sc.cal_mi_lcao(0,false)); - EXPECT_NO_THROW(sc.run_lambda_loop(0)); - EXPECT_FALSE(sc.check_rms_stop(0, 0, 0.0, 0.0, 0.0)); - EXPECT_NO_THROW(sc.print_termination()); - EXPECT_NO_THROW(sc.print_header()); + EXPECT_NO_THROW(sc.run_lambda_loop(0, true, std::cout)); + EXPECT_FALSE(spinconstrain::check_rms_stop(sc, 0, 0, 0.0, 0.0, 0.0, std::cout)); + EXPECT_NO_THROW(spinconstrain::print_termination(sc, std::cout)); + EXPECT_NO_THROW(spinconstrain::print_header(sc, std::cout)); std::vector> new_spin, old_spin, new_delta_lambda, old_delta_lambda; - EXPECT_FALSE(sc.check_gradient_decay(new_spin, old_spin, new_delta_lambda, old_delta_lambda, true)); + EXPECT_FALSE(spinconstrain::check_gradient_decay(sc, new_spin, old_spin, new_delta_lambda, old_delta_lambda, true, std::cout)); double alpha = 0.0; - EXPECT_NO_THROW(sc.check_restriction(new_spin, alpha)); - EXPECT_EQ(sc.cal_alpha_opt(new_spin, old_spin, alpha), 0.0); + EXPECT_NO_THROW(spinconstrain::check_restriction(sc, new_spin, alpha, std::cout)); + EXPECT_EQ(spinconstrain::cal_alpha_opt(sc, new_spin, old_spin, alpha), 0.0); } \ No newline at end of file diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index 83884054ccd..4a4b6d71f69 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -21,6 +21,7 @@ list(APPEND objects update_cell_pw.cpp setup_dftu_pw.cpp deltaspin_pw.cpp + deltaspin_pw_impl.cpp force_pw_nl.cpp force_pw_cc.cpp force_pw_scc.cpp diff --git a/source/source_pw/module_pwdft/deltaspin_pw.cpp b/source/source_pw/module_pwdft/deltaspin_pw.cpp index 34bd5156148..53966859998 100644 --- a/source/source_pw/module_pwdft/deltaspin_pw.cpp +++ b/source/source_pw/module_pwdft/deltaspin_pw.cpp @@ -8,7 +8,8 @@ namespace pw bool run_deltaspin_lambda_loop(const int iter, const double drho, - const Input_para& inp) + const Input_para& inp, + std::ostream& ofs_running) { /// Return false if DeltaSpin is not enabled if (!inp.sc_mag_switch) @@ -25,7 +26,7 @@ bool run_deltaspin_lambda_loop(const int iter, /// only records Mi vs lambda to lambda_scan_results.dat. if (inp.sc_lambda_strategy == "linear_scan") { - sc.run_lambda_linear_scan(iter); + sc.run_lambda_linear_scan(iter, ofs_running); return true; } @@ -34,7 +35,7 @@ bool run_deltaspin_lambda_loop(const int iter, if (!sc.mag_converged() && drho > 0 && drho < inp.sc_scf_thr) { /// Optimize lambda to get target magnetic moments - sc.run_lambda_loop(iter - 1); + sc.run_lambda_loop(iter - 1, true, ofs_running); sc.set_mag_converged(true); return true; } @@ -42,7 +43,7 @@ bool run_deltaspin_lambda_loop(const int iter, /// Re-run the lambda loop to update psi and charge density with current lambda. else if (sc.mag_converged()) { - sc.run_lambda_loop(iter - 1); + sc.run_lambda_loop(iter - 1, true, ofs_running); return true; } diff --git a/source/source_pw/module_pwdft/deltaspin_pw.h b/source/source_pw/module_pwdft/deltaspin_pw.h index b965dda53e4..4fecd43f81b 100644 --- a/source/source_pw/module_pwdft/deltaspin_pw.h +++ b/source/source_pw/module_pwdft/deltaspin_pw.h @@ -1,6 +1,7 @@ #ifndef DELTASPIN_PW_H #define DELTASPIN_PW_H +#include class Charge_Mixing; struct Input_para; @@ -17,12 +18,14 @@ namespace pw * @param iter The current iteration number (0-indexed). * @param drho The current charge density difference. * @param inp The input parameters. + * @param ofs_running The output stream for running log (e.g. running_scf.log). * @return true if the solver should be skipped (lambda loop was executed), * false otherwise. */ bool run_deltaspin_lambda_loop(const int iter, const double drho, - const Input_para& inp); + const Input_para& inp, + std::ostream& ofs_running); /** * @brief Check if SCF oscillation occurs for DeltaSpin method. diff --git a/source/source_pw/module_pwdft/deltaspin_pw_impl.cpp b/source/source_pw/module_pwdft/deltaspin_pw_impl.cpp new file mode 100644 index 00000000000..8b8fa5cc09e --- /dev/null +++ b/source/source_pw/module_pwdft/deltaspin_pw_impl.cpp @@ -0,0 +1,522 @@ +#include "source_base/matrix.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_title.h" +#include "source_base/timer.h" +#include "source_base/kernels/math_kernel_op.h" +#include "source_pw/module_pwdft/onsite_proj.h" +#include "source_lcao/module_deltaspin/spin_constrain.h" +#include "source_lcao/module_deltaspin/mi_tools.h" +#include "source_io/module_parameter/parameter.h" +#include "source_hsolver/diago_iter_assist.h" +#include "source_hsolver/hsolver_pw.h" +#include "source_estate/elecstate_pw.h" +#include "source_estate/elecstate_tools.h" + +namespace spinconstrain { + +/** + * @brief Calculate atomic magnetic moments using projector overlap (PW basis). + * + * @details For each k-point: + * 1. Tabulate atomic projectors: set up |alpha_{l,m}> for each atom + * 2. Compute becp = via overlap_proj_psi + * 3. Decompose becp into magnetic moments via accumulate_Mi_from_becp + * + * The magnetic moment is computed as: + * Mi = sum_{k,i} w_{k,i} * + * where P_at is the atomic projector and sigma are the Pauli matrices. + * + * Finally, Mi is summed across all MPI k-pool ranks since each pool only + * has a subset of k-points. + */ +template <> +void SpinConstrain>::cal_mi_pw() +{ + ModuleBase::TITLE("module_deltaspin", "cal_mi_pw"); + ModuleBase::timer::start("spinconstrain::SpinConstrain", "cal_mi_pw"); + + this->zero_Mi(); + if(PARAM.inp.device == "cpu") + { + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + // Loop over k-points to calculate Mi of sum_{k,i,l,m} + std::complex* psi_pointer = nullptr; + psi::Psi, base_device::DEVICE_CPU>* psi_t = static_cast, base_device::DEVICE_CPU>*>(this->psi); + const int nbands = psi_t->get_nbands(); + const int nks = psi_t->get_nk(); + const int npol = psi_t->get_npol(); + for(int ik = 0; ik < nks; ik++) + { + psi_t->fix_k(ik); + psi_pointer = psi_t->get_pointer(); + onsite_p->tabulate_atomic(ik); // Set up atomic projectors for this k-point + onsite_p->overlap_proj_psi(nbands * npol, psi_pointer); // Compute becp = + const std::complex* becp = onsite_p->get_h_becp(); + int nkb = onsite_p->get_tot_nproj(); + const int spin_sign = (npol == 2) ? 1 : this->get_spin_sign(ik); + accumulate_Mi_from_becp(becp, nkb, nbands, npol, spin_sign, + &this->pelec->wg(ik, 0), &onsite_p->get_nh(0), this->Mi_); + } + } +#if ((defined __CUDA) || (defined __ROCM)) + else + { + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + std::complex* psi_pointer = nullptr; + psi::Psi, base_device::DEVICE_GPU>* psi_t = static_cast, base_device::DEVICE_GPU>*>(this->psi); + const int nbands = psi_t->get_nbands(); + const int nks = psi_t->get_nk(); + const int npol = psi_t->get_npol(); + for(int ik = 0; ik < nks; ik++) + { + psi_t->fix_k(ik); + psi_pointer = psi_t->get_pointer(); + onsite_p->tabulate_atomic(ik); + onsite_p->overlap_proj_psi(nbands * npol, psi_pointer); + const std::complex* becp = onsite_p->get_h_becp(); + int nkb = onsite_p->get_size_becp() / nbands / npol; + const int spin_sign = (npol == 2) ? 1 : this->get_spin_sign(ik); + accumulate_Mi_from_becp(becp, nkb, nbands, npol, spin_sign, + &this->pelec->wg(ik, 0), &onsite_p->get_nh(0), this->Mi_); + } + } +#endif + // MPI reduction: sum Mi across all k-pool ranks + Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, GlobalV::NPROC_IN_POOL, &(this->Mi_[0][0]), 3 * this->Mi_.size()); + + ModuleBase::timer::end("spinconstrain::SpinConstrain", "cal_mi_pw"); +} + +/** + * @brief Compute DeltaSpin correction to the subspace Hamiltonian. + * + * @details Adds the constraint term to H in the projector subspace: + * H += becp^† * ps, where ps = delta_lambda * becp + * + * For non-collinear (npol=2), this implements the full 2x2 Pauli matrix: + * H_delta = | lambda_z lambda_x + i*lambda_y | + * | lambda_x - i*lambda_y -lambda_z | + * + * For collinear (npol=1), only the diagonal z-component with spin_sign: + * H_delta = lambda_z * spin_sign + * + * @param h_tmp Subspace Hamiltonian (nbands x nbands, modified in place) + * @param becp_k Projector coefficients for k-point ik + * @param delta_lambda Lambda change per atom (or full lambda if full_update) + * @param nbands Number of bands + * @param nkb Total number of projectors + * @param nh_iat Number of projectors per atom + * @param ik K-point index (for spin_sign lookup in collinear mode) + * @param full_update If true, compute delta = lambda_current - lambda_at_save + */ +template <> +void SpinConstrain>::calculate_delta_hcc(std::complex* h_tmp, const std::complex* becp_k, const ModuleBase::Vector3* delta_lambda, const int nbands, const int nkb, const int* nh_iat, const int ik, bool full_update) +{ + ModuleBase::TITLE("spinconstrain::SpinConstrain", "calculate_delta_hcc"); + ModuleBase::timer::start("spinconstrain::SpinConstrain", "calculate_delta_hcc"); + + // If full_update, compute actual delta = lambda_current - lambda_at_save + // This applies only the CHANGE in lambda, not the full lambda value + std::vector> actual_delta; + const ModuleBase::Vector3* effective_lambda = delta_lambda; + if (full_update) + { + int nat = this->get_nat(); + actual_delta.resize(nat); + for (int iat = 0; iat < nat; iat++) + { + actual_delta[iat] = delta_lambda[iat] - this->lambda_in_sub_[iat]; + } + effective_lambda = actual_delta.data(); + } + + int sum = 0; // Running sum of projectors across atoms + int size_ps = nkb * this->npol_ * nbands; // Total size of ps array + std::complex* becp_cpu = nullptr; + + // Handle GPU/CPU memory for becp + if(PARAM.inp.device == "gpu") + { +#if ((defined __CUDA) || (defined __ROCM)) + base_device::memory::resize_memory_op, base_device::DEVICE_CPU>()(becp_cpu, size_ps); + base_device::memory::synchronize_memory_op, base_device::DEVICE_CPU, base_device::DEVICE_GPU>()(becp_cpu, becp_k, size_ps); +#endif + } + else if (PARAM.inp.device == "cpu") + { + becp_cpu = const_cast*>(becp_k); + } + + // Compute modified projector coefficients: ps = delta_lambda * becp + std::vector> ps(size_ps, 0.0); + if(this->npol_ == 2) + { + // ============================================================= + // nspin=4 (non-collinear): full Pauli matrix treatment + // ============================================================= + // For each atom, construct 2x2 coefficients: + // | lambda_z lambda_x + i*lambda_y | + // | lambda_x - i*lambda_y -lambda_z | + // Then: ps_up = coeff0 * becp_up + coeff2 * becp_dn + // ps_dn = coeff1 * becp_up + coeff3 * becp_dn + for (int iat = 0; iat < this->Mi_.size(); iat++) + { + const int nproj = nh_iat[iat]; + const std::complex coefficients0(effective_lambda[iat][2], 0.0); + const std::complex coefficients1(effective_lambda[iat][0] , effective_lambda[iat][1]); + const std::complex coefficients2(effective_lambda[iat][0] , -1 * effective_lambda[iat][1]); + const std::complex coefficients3(-1 * effective_lambda[iat][2], 0.0); + for (int ib = 0; ib < nbands * this->npol_; ib += this->npol_) + { + for (int ip = 0; ip < nproj; ip++) + { + const int becpind = ib * nkb + sum + ip; + const std::complex becp1 = becp_cpu[becpind]; + const std::complex becp2 = becp_cpu[becpind + nkb]; + ps[becpind] += coefficients0 * becp1 + + coefficients2 * becp2; + ps[becpind + nkb] += coefficients1 * becp1 + + coefficients3 * becp2; + } + } + sum += nproj; + } + } + else if(this->npol_ == 1) + { + // ============================================================= + // nspin=2 (collinear): only z-component with spin_sign + // ============================================================= + // ps = lambda_z * spin_sign * becp + // spin_sign = +1 for spin-up k-points, -1 for spin-down + for (int iat = 0; iat < this->Mi_.size(); iat++) + { + const int nproj = nh_iat[iat]; + double coefficients0 = effective_lambda[iat][2] * this->get_spin_sign(ik); + for (int ib = 0; ib < nbands; ib++) + { + for (int ip = 0; ip < nproj; ip++) + { + const int becpind = ib * nkb + sum + ip; + const std::complex becp1 = becp_cpu[becpind]; + ps[becpind] += coefficients0 * becp1; + } + } + sum += nproj; + } + } + + // Copy ps to GPU if needed + std::complex* ps_pointer = nullptr; + if(PARAM.inp.device == "gpu") + { +#if ((defined __CUDA) || (defined __ROCM)) + base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(ps_pointer, size_ps); + base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_CPU>()(ps_pointer, ps.data(), size_ps); +#endif + } + else if (PARAM.inp.device == "cpu") + { + ps_pointer = ps.data(); + } + + // ============================================================= + // H += becp^† * ps (GEMM: C = alpha * A^† * B + beta * C) + // A = becp_k (npm x nbands), B = ps (npm x nbands), C = h_tmp (nbands x nbands) + // ============================================================= + char transa = 'C'; // Conjugate transpose of becp + char transb = 'N'; // Normal ps + const int npm = nkb * this->npol_; + if (PARAM.inp.device == "gpu") + { +#if ((defined __CUDA) || (defined __ROCM)) + ModuleBase::gemm_op, base_device::DEVICE_GPU>()( + transa, + transb, + nbands, + nbands, + npm, + &ModuleBase::ONE, + becp_k, + npm, + ps_pointer, + npm, + &ModuleBase::ONE, + h_tmp, + nbands + ); + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(ps_pointer); + base_device::memory::delete_memory_op, base_device::DEVICE_CPU>()(becp_cpu); +#endif + + } + else if (PARAM.inp.device == "cpu") + { + ModuleBase::gemm_op, base_device::DEVICE_CPU>()( + transa, + transb, + nbands, + nbands, + npm, + &ModuleBase::ONE, + becp_k, + npm, + ps_pointer, + npm, + &ModuleBase::ONE, + h_tmp, + nbands + ); + } + ModuleBase::timer::end("spinconstrain::SpinConstrain", "calculate_delta_hcc"); +} + +/** + * @brief CPU implementation of PW wavefunction and charge density update. + * + * @par Two-stage process: + * Stage 1 - Subspace diagonalization: + * For each k-point, apply DeltaSpin correction to the saved subspace H, + * then diagonalize to rotate the wavefunctions. This is a cheap operation + * in the reduced subspace (nbands x nbands). + * + * Stage 2 - Full-space update: + * Option A (pw_solve=true): Run HSolverPW for iterative refinement in the + * full plane-wave space. This is more accurate but expensive. + * Option B (pw_solve=false): Update weights from new eigenvalues and call + * psiToRho() to build the charge density from current psi. Faster but + * may be less accurate if the subspace rotation was not sufficient. + * + * @par Memory management + * Frees sub_h_save, sub_s_save, becp_save after use. These are allocated + * on the first cal_mw_from_lambda() call and should only be freed here. + * + * @param delta_lambda Lambda change for incremental H correction + * @param pw_solve If true, run full PW solver; if false, just update weights + * @param full_update If true, apply full lambda (not delta) to H correction + */ +template <> +void SpinConstrain>::update_psi_charge_pw_cpu(const ModuleBase::Vector3* delta_lambda, bool pw_solve, bool full_update) +{ + ModuleBase::TITLE("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); + ModuleBase::timer::start("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); + + psi::Psi>* psi_t = static_cast>*>(this->psi); + hamilt::Hamilt, base_device::DEVICE_CPU>* hamilt_t = static_cast, base_device::DEVICE_CPU>*>(this->p_hamilt); + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + + int nbands = psi_t->get_nbands(); + int npol = psi_t->get_npol(); + int nkb = onsite_p->get_tot_nproj(); + int nk = psi_t->get_nk(); + int size_becp = nbands * nkb * npol; + const int* nh_iat = &onsite_p->get_nh(0); + + std::vector> h_tmp(nbands * nbands), s_tmp(nbands * nbands); + + // CRITICAL: subspace data must have been saved by cal_mw_from_lambda() + assert(this->sub_h_save != nullptr); + assert(this->sub_s_save != nullptr); + assert(this->becp_save != nullptr); + + // Determine which lambda to use for H correction + const ModuleBase::Vector3* lambda_for_hcc = delta_lambda; + std::vector> computed_delta; + if (full_update) + { + lambda_for_hcc = this->lambda_.data(); + } + + // ============================================================= + // STAGE 1: Subspace diagonalization for each k-point + // ============================================================= + for (int ik = 0; ik < nk; ++ik) + { + std::complex* h_k = this->sub_h_save + ik * nbands * nbands; + std::complex* s_k = this->sub_s_save + ik * nbands * nbands; + std::complex* becp_k = this->becp_save + ik * size_becp; + + psi_t->fix_k(ik); + + // Copy saved subspace matrices to temp + memcpy(h_tmp.data(), h_k, sizeof(std::complex) * nbands * nbands); + memcpy(s_tmp.data(), s_k, sizeof(std::complex) * nbands * nbands); + + // Apply DeltaSpin correction: H += becp^† * lambda * becp + this->calculate_delta_hcc(h_tmp.data(), becp_k, lambda_for_hcc, nbands, nkb, nh_iat, ik, full_update); + + // Diagonalize in subspace to update wavefunction coefficients and eigenvalues + hsolver::DiagoIterAssist>::diag_subspace_psi(h_tmp.data(), + s_tmp.data(), + nbands, + psi_t[0], + &this->pelec->ekb(ik, 0)); + } + + // Free saved subspace data (allocated in cal_mw_from_lambda) + delete[] this->sub_h_save; + delete[] this->sub_s_save; + delete[] this->becp_save; + this->sub_h_save = nullptr; + this->sub_s_save = nullptr; + this->becp_save = nullptr; + + // ============================================================= + // STAGE 2: Full-space update + // ============================================================= + if (pw_solve) + { + // Full PW diagonalization: subspace rotation provides a good initial guess, + // then HSolverPW iteratively refines psi in the full plane-wave space and calls psiToRho. + hsolver::HSolverPW, base_device::DEVICE_CPU> hsolver_pw_obj( + this->pw_wfc_, + PARAM.inp.calculation, + PARAM.inp.basis_type, + PARAM.inp.ks_solver, + PARAM.globalv.use_uspp, + PARAM.inp.nspin, + hsolver::DiagoIterAssist>::SCF_ITER, + hsolver::DiagoIterAssist>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist>::PW_DIAG_THR, + hsolver::DiagoIterAssist>::need_subspace, + PARAM.inp.nbands, + PARAM.inp.diago_smooth_ethr, + PARAM.inp.pw_diag_ndim, + PARAM.inp.diag_subspace, + PARAM.inp.nb2d, + PARAM.inp.use_k_continuity); + + hsolver_pw_obj.solve(hamilt_t, psi_t[0], this->pelec, this->pelec->ekb.c, + GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, false, this->tpiba, this->get_nat()); + } + else + { + // No full solver: update weights from new eigenvalues, then build rho from current psi + elecstate::calculate_weights(this->pelec->ekb, + this->pelec->wg, + this->pelec->klist, + this->pelec->eferm, + this->pelec->f_en, + this->pelec->nelec_spin, + PARAM.inp.nbands, + this->pelec->skip_weights); + elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en); + reinterpret_cast, base_device::DEVICE_CPU>*>(this->pelec)->psiToRho(*psi_t); + } + ModuleBase::timer::end("spinconstrain::SpinConstrain", "update_psi_charge_pw_cpu"); +} + +#if ((defined __CUDA) || (defined __ROCM)) +/** + * @brief GPU implementation of PW wavefunction and charge density update. + * + * @details Same algorithm as update_psi_charge_pw_cpu(), but with GPU memory + * management (device allocation, host-device synchronization). + */ +template <> +void SpinConstrain>::update_psi_charge_pw_gpu(const ModuleBase::Vector3* delta_lambda, bool pw_solve, bool full_update) +{ + ModuleBase::TITLE("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); + ModuleBase::timer::start("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); + + psi::Psi, base_device::DEVICE_GPU>* psi_t = static_cast, base_device::DEVICE_GPU>*>(this->psi); + hamilt::Hamilt, base_device::DEVICE_GPU>* hamilt_t = static_cast, base_device::DEVICE_GPU>*>(this->p_hamilt); + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + + int nbands = psi_t->get_nbands(); + int npol = psi_t->get_npol(); + int nkb = onsite_p->get_tot_nproj(); + int nk = psi_t->get_nk(); + int size_becp = nbands * nkb * npol; + const int* nh_iat = &onsite_p->get_nh(0); + + std::complex* h_tmp = nullptr; + std::complex* s_tmp = nullptr; + base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(h_tmp, nbands * nbands); + base_device::memory::resize_memory_op, base_device::DEVICE_GPU>()(s_tmp, nbands * nbands); + + assert(this->sub_h_save != nullptr); + assert(this->sub_s_save != nullptr); + assert(this->becp_save != nullptr); + + const ModuleBase::Vector3* lambda_for_hcc = delta_lambda; + std::vector> computed_delta; + if (full_update) + { + lambda_for_hcc = this->lambda_.data(); + } + + // STAGE 1: Subspace diagonalization for each k-point (GPU) + for (int ik = 0; ik < nk; ++ik) + { + std::complex* h_k = this->sub_h_save + ik * nbands * nbands; + std::complex* s_k = this->sub_s_save + ik * nbands * nbands; + std::complex* becp_k = this->becp_save + ik * size_becp; + + psi_t->fix_k(ik); + + base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(h_tmp, h_k, nbands * nbands); + base_device::memory::synchronize_memory_op, base_device::DEVICE_GPU, base_device::DEVICE_GPU>()(s_tmp, s_k, nbands * nbands); + + this->calculate_delta_hcc(h_tmp, becp_k, lambda_for_hcc, nbands, nkb, nh_iat, ik, full_update); + + hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::diag_subspace_psi(h_tmp, + s_tmp, + nbands, + psi_t[0], + &this->pelec->ekb(ik, 0)); + } + + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(h_tmp); + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(s_tmp); + + // Free GPU memory for saved subspace data + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(sub_h_save); + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(sub_s_save); + base_device::memory::delete_memory_op, base_device::DEVICE_GPU>()(becp_save); + this->sub_h_save = nullptr; + this->sub_s_save = nullptr; + this->becp_save = nullptr; + + // STAGE 2: Full-space update (GPU) + if (pw_solve) + { + hsolver::HSolverPW, base_device::DEVICE_GPU> hsolver_pw_obj( + this->pw_wfc_, + PARAM.inp.calculation, + PARAM.inp.basis_type, + PARAM.inp.ks_solver, + PARAM.globalv.use_uspp, + PARAM.inp.nspin, + hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::SCF_ITER, + hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::PW_DIAG_THR, + hsolver::DiagoIterAssist, base_device::DEVICE_GPU>::need_subspace, + PARAM.inp.nbands, + PARAM.inp.diago_smooth_ethr, + PARAM.inp.pw_diag_ndim, + PARAM.inp.diag_subspace, + PARAM.inp.nb2d, + PARAM.inp.use_k_continuity); + + hsolver_pw_obj.solve(hamilt_t, psi_t[0], this->pelec, this->pelec->ekb.c, + GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, false, this->tpiba, this->get_nat()); + } + else + { + elecstate::calculate_weights(this->pelec->ekb, + this->pelec->wg, + this->pelec->klist, + this->pelec->eferm, + this->pelec->f_en, + this->pelec->nelec_spin, + PARAM.inp.nbands, + this->pelec->skip_weights); + elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en); + reinterpret_cast, base_device::DEVICE_GPU>*>(this->pelec)->psiToRho(*psi_t); + } + ModuleBase::timer::end("spinconstrain::SpinConstrain", "update_psi_charge_pw_gpu"); +} +#endif + +} // namespace spinconstrain