From eb03771cb3e37188e49287316bb6d9617a860a81 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 10:25:28 +0800 Subject: [PATCH 1/3] Refactor DFT+U onsite Hamiltonian assembly --- .../source_lcao/module_dftu/dftu_nao_op.cpp | 277 +++++++++++------- source/source_lcao/module_dftu/dftu_nao_op.h | 13 + 2 files changed, 180 insertions(+), 110 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index 2dc5c9cb6b5..da980da15ad 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -12,7 +12,8 @@ #include "source_cell/klist.h" #include "source_cell/module_symmetry/symmetry.h" -#include +#include +#include #include "dftu_nao_adj.h" #include "dftu_nao_fs_r.h" @@ -49,13 +50,170 @@ hamilt::DFTU_onsite>::DFTU_onsite(HS_Matrix_K* // and are computed once here; the operator is rebuilt every ionic step. // Kept in the constructor body (not the initializer list) so failures // inside these calls are easy to debug. - this->adjs_all = DFTU_LCAO::build_adjacent_atoms(this->ucell, this->dftu, GridD_in, this->orb_cutoff_, onsite_radius); + this->adjs_all = DFTU_LCAO::build_adjacent_atoms( + this->ucell, this->dftu, GridD_in, this->orb_cutoff_, onsite_radius); const Parallel_Orbitals* pv = this->hR->get_atom_pair(0).get_paraV(); this->nlm_tot = DFTU_LCAO::cal_nlm_all(*this->ucell, *this->dftu, *this->intor_, this->adjs_all, *pv); ModuleBase::timer::end("DFTU_onsite", "DFTU_onsite"); } +template +bool hamilt::DFTU_onsite>::should_restore_symmetric_dmr() const +{ + return !this->dftu->is_occmat_ready() && this->kv_ != nullptr + && ModuleSymmetry::Symmetry::symm_flag == 1 && !this->kv_->kstars.empty(); +} + +template +void hamilt::DFTU_onsite>::ensure_symmetry_rotation( + const Parallel_Orbitals& pv) +{ + if (this->symrot_built_) + { + return; + } + + const std::array period{this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2]}; + const std::vector bvk_cells + = ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period); + this->symrot_.find_irred_sector(this->ucell->symm, + this->ucell->atoms, + this->ucell->nat, + this->ucell->ntype, + this->ucell->iat2it, + this->ucell->iat2ia, + bvk_cells, + period, + this->ucell->lat); + this->symrot_.cal_Ms(*this->kv_, + this->ucell->symm, + this->ucell->atoms, + this->ucell->latvec, + this->ucell->lmax, + this->ucell->nat, + this->ucell->iat2it, + this->ucell->iat2ia, + pv, + this->nspin); + this->symrot_built_ = true; +} + +template +std::vector> +hamilt::DFTU_onsite>::build_full_bz_kvectors(const int nspin0) const +{ + // restore_dm() returns the stars owned by this k-point pool. Build the + // matching k-vector sequence using the same local-to-global IBZ mapping. + const int nk_local = this->kv_->get_nks() / nspin0; + const int nks_ibz_global = static_cast(this->kv_->kstars.size()); + std::vector> kvec_d_full; + for (int ik_local = 0; ik_local < nk_local; ++ik_local) + { + const int ik_ibz = this->kv_->ik2iktot[ik_local] % nks_ibz_global; + for (const std::pair>& isym_kvd + : this->kv_->kstars[ik_ibz]) + { + kvec_d_full.push_back(isym_kvd.second); + } + } + return kvec_d_full; +} + +template +void hamilt::DFTU_onsite>::contribute_from_symmetric_dmr( + const Parallel_Orbitals& pv) +{ + this->ensure_symmetry_rotation(pv); + + const int nspin0 = (this->nspin == 2) ? 2 : 1; + const std::vector> kvec_d_full + = this->build_full_bz_kvectors(nspin0); + const std::vector> dmk_full + = this->symrot_.restore_dm(*this->kv_, this->dm_->get_dmk_vec(), pv); + const int nks_full = static_cast(kvec_d_full.size()); + module_dm::DensityMatrix dmr_sym(&pv, nspin0, kvec_d_full, nks_full); + dmr_sym.init_dmr(*this->dm_->get_dmr_ptr(1)); + dmr_sym.get_dmk_vec() = dmk_full; + dmr_sym.cal_dmr(-1); + + const hamilt::HContainer* dmr = dmr_sym.get_dmr_ptr(this->current_spin + 1); + this->contribute_hubbard_atoms(pv, dmr); +} + +template +void hamilt::DFTU_onsite>::contribute_hubbard_atoms( + const Parallel_Orbitals& pv, + const hamilt::HContainer* dmr) +{ + int atom_index = 0; + for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + { + int T0 = 0; + int I0 = 0; + this->ucell->iat2iait(iat0, &I0, &T0); + if (!this->dftu->has_l_channel(T0)) + { + continue; + } + + const int target_L = this->dftu->get_l_channel(T0); + const int tlp1 = 2 * target_L + 1; + AdjacentAtomInfo& adjs = this->adjs_all[atom_index++]; + const int spin_fold = (this->nspin == 4) ? 4 : 1; + std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); + + if (!this->dftu->is_occmat_ready()) + { + assert(dmr != nullptr); + DFTU_LCAO::compute_occ_from_dmr(*this->ucell, + *this->dftu, + iat0, + target_L, + this->current_spin, + this->nspin, + adjs, + pv, + this->nlm_tot, + *dmr, + occ); + } + else + { + DFTU_LCAO::load_occ_from_file( + *this->dftu, iat0, target_L, this->nspin, this->current_spin, occ); + } + + const double u_value = this->dftu->get_u_current(T0); + std::vector pot_onsite_tmp(occ.size()); + double u_energy = this->dftu->get_energy(); + DFTU_LCAO::cal_pot_onsite(occ, tlp1, u_value, pot_onsite_tmp.data(), u_energy); + this->dftu->set_energy(u_energy); + + std::vector pot_onsite(occ.size()); + DFTU_LCAO::transfer_pot_onsite(pot_onsite_tmp, pot_onsite); + DFTU_LCAO::accumulate_hr_for_iat0( + *this->ucell, this->hR, this->nlm_tot, iat0, adjs, pv, pot_onsite); + } +} + +template +void hamilt::DFTU_onsite>::finish_spin_step() +{ + if (this->nspin == 1) + { + this->dftu->set_double_energy(); + } + if (this->current_spin == this->nspin - 1 || this->nspin == 4) + { + this->dftu->set_occmat_stale(); + } + if (this->nspin == 2) + { + this->current_spin = 1 - this->current_spin; + } +} + // contributeHR() /** * @brief Contribute DFT+U Hamiltonian to real-space HR matrix @@ -128,122 +286,21 @@ void hamilt::DFTU_onsite>::contributeHR() const Parallel_Orbitals* pv = this->hR->get_atom_pair(0).get_paraV(); // nlm_tot is precomputed in the constructor (structure snapshot) - // (symmetry) when crystal symmetry reduces the k-mesh, dm_->get_dmr_ptr() - // was Fourier-transformed from the irreducible k-points only and is not - // actually symmetric; reconstruct the full-BZ DMR once here (reused by every - // atom below) via the same D(k) restoration EXX already uses for its own - // real-space density matrix (ModuleSymmetry::Symmetry_rotation::restore_dm). - std::unique_ptr> dmr_sym; - if (!this->dftu->is_occmat_ready() && this->kv_ != nullptr && ModuleSymmetry::Symmetry::symm_flag == 1 - && !this->kv_->kstars.empty()) + if (this->should_restore_symmetric_dmr()) { - if (!this->symrot_built_) - { - const std::array period{ this->kv_->nmp[0], this->kv_->nmp[1], this->kv_->nmp[2] }; - // for return_lattice to calculate Ms - this->symrot_.find_irred_sector(this->ucell->symm, this->ucell->atoms, - this->ucell->nat, this->ucell->ntype, this->ucell->iat2it, this->ucell->iat2ia, - ModuleSymmetry::Symmetry_rotation_k::get_bvk_cells(period), period, this->ucell->lat); - this->symrot_.cal_Ms(*this->kv_, this->ucell->symm, this->ucell->atoms, - this->ucell->latvec, this->ucell->lmax, - this->ucell->nat, this->ucell->iat2it, this->ucell->iat2ia, - *pv, this->nspin); - this->symrot_built_ = true; - } - const int nspin0 = (this->nspin == 2) ? 2 : 1; - // (k-point pools, KPAR>1) restore_dm() now returns only the stars of THIS pool's own - // local irreducible k-points (see its definition for why that's enough); kvec_d_full - // must be built the same way -- one entry per star member of each local ibz-k, - // enumerated in the same order restore_dm uses for its spin-0 block (kv.ik2iktot maps - // the spin-0 and spin-1 blocks to the same sequence of global ibz indices, so a single - // list built from the spin-0 mapping is valid for the whole nspin0-block DensityMatrix). - const int nk_local = this->kv_->get_nks() / nspin0; - const int nks_ibz_global = static_cast(this->kv_->kstars.size()); - std::vector> kvec_d_full; - for (int ik_local = 0; ik_local < nk_local; ++ik_local) - { - const int ik_ibz = this->kv_->ik2iktot[ik_local] % nks_ibz_global; - for (const std::pair>& isym_kvd : this->kv_->kstars[ik_ibz]) { kvec_d_full.push_back(isym_kvd.second); } - } - const std::vector> dmk_full = this->symrot_.restore_dm(*this->kv_, this->dm_->get_dmk_vec(), *pv); - dmr_sym.reset(new module_dm::DensityMatrix(pv, nspin0, kvec_d_full, static_cast(kvec_d_full.size()))); - dmr_sym->init_dmr(*this->dm_->get_dmr_ptr(1)); - dmr_sym->get_dmk_vec() = dmk_full; - dmr_sym->cal_dmr(-1); + this->contribute_from_symmetric_dmr(*pv); } - - // loop over all Hubbard-projector center atoms (iat0) - int atom_index = 0; - for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + else { - int T0 = 0; - int I0 = 0; - ucell->iat2iait(iat0, &I0, &T0); - if (!this->dftu->has_l_channel(T0)) - { - continue; - } - const int target_L = this->dftu->get_l_channel(T0); - const int tlp1 = 2 * target_L + 1; - AdjacentAtomInfo& adjs = this->adjs_all[atom_index++]; - - const int spin_fold = (this->nspin == 4) ? 4 : 1; - std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); - - // compute or load occupation matrix + const hamilt::HContainer* dmr = nullptr; if (!this->dftu->is_occmat_ready()) { - // DMR is guaranteed ready here: otherwise the early exit above - // would have returned. DMR index is 1-based, hence +1. - const hamilt::HContainer* dmr = this->dm_->get_dmr_ptr(this->current_spin + 1); - if (dmr_sym) { dmr = dmr_sym->get_dmr_ptr(this->current_spin + 1); } - DFTU_LCAO::compute_occ_from_dmr(*this->ucell, - *this->dftu, - iat0, - target_L, - this->current_spin, - this->nspin, - adjs, - *pv, - this->nlm_tot, - *dmr, - occ); + dmr = this->dm_->get_dmr_ptr(this->current_spin + 1); } - else - { - DFTU_LCAO::load_occ_from_file(*this->dftu, iat0, target_L, this->nspin, this->current_spin, occ); - } - - // compute Hubbard potential and energy - const double u_value = this->dftu->get_u_current(T0); - std::vector pot_onsite_tmp(occ.size()); - double u_energy = this->dftu->get_energy(); - DFTU_LCAO::cal_pot_onsite(occ, tlp1, u_value, pot_onsite_tmp.data(), u_energy); - this->dftu->set_energy(u_energy); - - std::vector pot_onsite(occ.size()); - DFTU_LCAO::transfer_pot_onsite(pot_onsite_tmp, pot_onsite); - - // accumulate HR contributions from neighbor pairs - DFTU_LCAO::accumulate_hr_for_iat0(*this->ucell, this->hR, this->nlm_tot, iat0, adjs, *pv, pot_onsite); - } - - // post-processing: energy doubling for nspin=1 - if (this->nspin == 1) - { - this->dftu->set_double_energy(); - } - // mark occ_mat stale for next iteration - if (this->current_spin == this->nspin - 1 || this->nspin == 4) - { - this->dftu->set_occmat_stale(); - } - // toggle spin channel for nspin=2 - if (this->nspin == 2) - { - this->current_spin = 1 - this->current_spin; + this->contribute_hubbard_atoms(*pv, dmr); } + this->finish_spin_step(); ModuleBase::timer::end("DFTU_onsite", "contributeHR"); } diff --git a/source/source_lcao/module_dftu/dftu_nao_op.h b/source/source_lcao/module_dftu/dftu_nao_op.h index 51fbb1dd9bd..54483d904b6 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.h +++ b/source/source_lcao/module_dftu/dftu_nao_op.h @@ -71,6 +71,19 @@ class DFTU_onsite> : public OperatorLCAO void set_kv(const K_Vectors* kv_in) { this->kv_ = kv_in; } private: + bool should_restore_symmetric_dmr() const; + + void ensure_symmetry_rotation(const Parallel_Orbitals& pv); + + std::vector> build_full_bz_kvectors(const int nspin0) const; + + void contribute_from_symmetric_dmr(const Parallel_Orbitals& pv); + + void contribute_hubbard_atoms(const Parallel_Orbitals& pv, + const hamilt::HContainer* dmr); + + void finish_spin_step(); + const UnitCell* ucell = nullptr; Plus_U_Base* dftu = nullptr; From 2d30e9e60c4ef154392409a23fc493c623788625 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 19:52:36 +0800 Subject: [PATCH 2/3] Avoid copying restored DFT+U density matrix --- source/source_lcao/module_dftu/dftu_nao_op.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index da980da15ad..5dedc370d33 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -129,12 +129,11 @@ void hamilt::DFTU_onsite>::contribute_from_symmetri const int nspin0 = (this->nspin == 2) ? 2 : 1; const std::vector> kvec_d_full = this->build_full_bz_kvectors(nspin0); - const std::vector> dmk_full - = this->symrot_.restore_dm(*this->kv_, this->dm_->get_dmk_vec(), pv); const int nks_full = static_cast(kvec_d_full.size()); module_dm::DensityMatrix dmr_sym(&pv, nspin0, kvec_d_full, nks_full); dmr_sym.init_dmr(*this->dm_->get_dmr_ptr(1)); - dmr_sym.get_dmk_vec() = dmk_full; + dmr_sym.get_dmk_vec() + = this->symrot_.restore_dm(*this->kv_, this->dm_->get_dmk_vec(), pv); dmr_sym.cal_dmr(-1); const hamilt::HContainer* dmr = dmr_sym.get_dmr_ptr(this->current_spin + 1); From 943e514c20c5713ddf99c620142126cabed45045 Mon Sep 17 00:00:00 2001 From: lanshuyue <140165754+lanshuyue@users.noreply.github.com> Date: Wed, 21 Oct 2026 20:04:34 +0800 Subject: [PATCH 3/3] ci: retrigger cancelled documentation build