diff --git a/AGENTS.md b/AGENTS.md index e868eceb8ba..57fea4ed84a 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -103,6 +103,9 @@ rules. Read the complete governance document before making or reviewing changes: - Member -> free function: inventory `this->` reads; pass as params (const for config, ref for mutable state); move only when body is `this`-free; keep thin wrapper; compile each step. +- Extract a base-class nested-vector member in three steps (hold + forward, + switch writers, delete legacy) so no commit mixes old-storage writes with + new-storage reads. ## Local Commands diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 0515b73494b..861e716ff16 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -399,13 +399,16 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\ veff_lcao.o\ veff_dh.o\ meta_lcao.o\ - dftu_lcao_op.o\ + dftu_nao_op.o\ deepks_lcao.o\ op_exx_lcao.o\ dspin_lcao.o\ dspin_fs.o\ setup_dftu_lcao.o\ - dftu_fs.o\ + dftu_nao_fs_k.o\ + dftu_nao_for_r.o\ + dftu_nao_fs_r.o\ + dftu_nao_str_r.o\ operator_fs_utils.o\ OBJS_HCONTAINER=base_matrix.o\ @@ -808,9 +811,9 @@ OBJS_SRCPW=h_ewald_pw.o\ setup_pwwfc.o\ update_cell_pw.o\ dftu_base.o\ - dftu_output.o\ - dftu_tools_pw.o\ - dftu_cal_occ_pw.o\ + dftu_base_io.o\ + dftu_base_occ.o\ + dftu_base_tools.o\ setup_dftu_pw.o\ deltaspin_pw.o\ deltaspin_pw_impl.o\ @@ -885,14 +888,16 @@ OBJS_VDW=vdw.o\ vdwd3_autoset_xcname.o\ vdwd3_auto_xcpar.o -OBJS_DFTU=dftu_lcao.o\ - dftu_force.o\ - dftu_yukawa.o\ - dftu_folding.o\ - dftu_lcao_pots.o\ - dftu_lcao_energy.o\ - dftu_lcao_op_legacy.o\ - dftu_lcao_occ.o\ +OBJS_DFTU=dftu_nao.o\ + dftu_nao_fs_k.o\ + dftu_nao_for_r.o\ + dftu_nao_fs_r.o\ + dftu_nao_str_r.o\ + dftu_nao_folding.o\ + dftu_nao_pots.o\ + dftu_nao_energy.o\ + dftu_nao_op_legacy.o\ + dftu_nao_occ.o\ dftu_hamilt.o OBJS_DELTASPIN=basic_funcs.o\ diff --git a/source/source_base/module_container/ATen/kernels/memory.h b/source/source_base/module_container/ATen/kernels/memory.h index da079d7a8c4..ba8ad531a63 100644 --- a/source/source_base/module_container/ATen/kernels/memory.h +++ b/source/source_base/module_container/ATen/kernels/memory.h @@ -81,7 +81,7 @@ struct synchronize_memory_stride { const std::vector& out_size, const std::vector& in_size) { - REQUIRES_OK(in_size.size() == out_size.size() && in_size.size() <= 2); + REQUIRES_OK(in_size.size() == out_size.size() && in_size.size() <= 2, "rank mismatch: in_size and out_size must have the same rank <= 2"); if (in_size.size() == 1) { synchronize_memory()(arr_out, arr_in, in_size[0]); } diff --git a/source/source_base/module_container/base/macros/macros.h b/source/source_base/module_container/base/macros/macros.h index fbcae6fe8a2..597920d5d5a 100644 --- a/source/source_base/module_container/base/macros/macros.h +++ b/source/source_base/module_container/base/macros/macros.h @@ -70,7 +70,7 @@ __func__, \ __FILE__, \ static_cast(__LINE__), \ - CHECK_MSG(expr, ##__VA_ARGS__)); \ + CHECK_MSG(expr, __VA_ARGS__)); \ } // The macro TEMPLATE_1() expands to a switch statement conditioned on diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index eaec4ad8f76..b847dc4024f 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -13,7 +13,7 @@ #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_output/output_log.h" // use write_head #include "source_estate/elecstate_print.h" // print_etot -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 2025-11-07 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 2025-11-07 #include "source_hamilt/module_xc/general_exx_info.h" // for init_general_exx_info namespace ModuleESolver diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index 8efad1be555..b10e4ef7307 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -8,7 +8,7 @@ #include "source_hamilt/hamilt.h" // use Hamiltonian #include "source_hamilt/hamilt_base.h" // use Hamiltonian base class #include "source_hamilt/module_xc/general_exx_info.h" // ESolver owns General_Exx_Info value -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 #include "source_pw/module_pwdft/vnl_pw.h" namespace ModuleESolver diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index f9c7e776fea..293d72ad916 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -202,7 +202,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const // update local occupations for DFT+U // should before lambda loop in DeltaSpin - pw::iter_init_dftu_pw(iter, + DFTU_BASE::iter_init_dftu_pw(iter, istep, this->dftu, this->stp.template get_psi_t(), diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 7eaea3cdc2b..943b6632f42 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -15,7 +15,7 @@ #include "source_lcao/hamilt_lcao.h" #include "source_lcao/lcao_domain.h" #include "source_lcao/module_deltaspin/spin_constrain.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_lcao/module_operator_lcao/op_exx_lcao.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index 6277aea17d5..b747b07d85c 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -42,6 +42,7 @@ list(APPEND objects module_charge/symm_rhog.cpp fp_energy.cpp occupy.cpp + occ_matrix.cpp param_update.cpp setup_estate_pw.cpp update_pot.cpp diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 82ceb15d4b0..5e8b5a9c9cd 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -42,26 +42,57 @@ class ElecState // calculate electronic charge density on grid points or density matrix in real space // the consequence charge density rho saved into rho_out, preparing for charge mixing. + // NOTE: all overloads are intentionally provided (with empty bodies) so that + // template derived classes (e.g. ElecStatePW) can safely mark their + // psiToRho/cal_tau as 'override' regardless of which (T, Device) combo is instantiated. virtual void psiToRho(const psi::Psi>& psi) { return; } + virtual void psiToRho(const psi::Psi, base_device::DEVICE_GPU>& psi) + { + return; + } virtual void psiToRho(const psi::Psi& psi) { return; } + virtual void psiToRho(const psi::Psi& psi) + { + return; + } + virtual void psiToRho(const psi::Psi>& psi) + { + return; + } + virtual void psiToRho(const psi::Psi, base_device::DEVICE_GPU>& psi) + { + return; + } virtual void cal_tau(const psi::Psi>& psi) { return; } + virtual void cal_tau(const psi::Psi, base_device::DEVICE_GPU>& psi) + { + return; + } virtual void cal_tau(const psi::Psi& psi) { return; } + virtual void cal_tau(const psi::Psi& psi) + { + return; + } virtual void cal_tau(const psi::Psi>& psi) { return; } + virtual void cal_tau(const psi::Psi, base_device::DEVICE_GPU>& psi) + { + return; + } // update charge density for next scf step // in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge diff --git a/source/source_estate/elecstate_pw.h b/source/source_estate/elecstate_pw.h index 53e39917a5d..203838878ab 100644 --- a/source/source_estate/elecstate_pw.h +++ b/source/source_estate/elecstate_pw.h @@ -32,9 +32,9 @@ class ElecStatePW : public ElecState ~ElecStatePW(); //! interface for HSolver to calculate rho from Psi - virtual void psiToRho(const psi::Psi& psi); + void psiToRho(const psi::Psi& psi) override; - virtual void cal_tau(const psi::Psi& psi); + void cal_tau(const psi::Psi& psi) override; double get_spin_constrain_energy() override; diff --git a/source/source_estate/module_charge/chgmixing.cpp b/source/source_estate/module_charge/chgmixing.cpp index 617579c4c14..4886e9f5e21 100644 --- a/source/source_estate/module_charge/chgmixing.cpp +++ b/source/source_estate/module_charge/chgmixing.cpp @@ -1,6 +1,6 @@ #include "source_estate/module_charge/chgmixing.h" #include "source_estate/update_pot.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_lcao/module_deltaspin/spin_constrain.h" void module_charge::chgmixing_ks(const int iter, // scf iteration number diff --git a/source/source_estate/occ_matrix.cpp b/source/source_estate/occ_matrix.cpp new file mode 100644 index 00000000000..89283ef5f17 --- /dev/null +++ b/source/source_estate/occ_matrix.cpp @@ -0,0 +1,350 @@ +#include "source_estate/occ_matrix.h" + +#include "source_base/timer.h" +#include "source_cell/unitcell.h" + +void OccupationMatrix::init(const UnitCell& cell, + const std::vector& orbital_corr, + const int nspin, + const int npol) +{ + this->nspin_ = nspin; + this->npol_ = npol; + + this->occ_.resize(cell.nat); + this->occ_save_.resize(cell.nat); + this->iatlnmipol2iwt_.resize(cell.nat); + + for (int it = 0; it < cell.ntype; ++it) + { + for (int ia = 0; ia < cell.atoms[it].na; ia++) + { + const int iat = cell.itia2iat(it, ia); + + occ_[iat].resize(cell.atoms[it].nwl + 1); + occ_save_[iat].resize(cell.atoms[it].nwl + 1); + iatlnmipol2iwt_[iat].resize(cell.atoms[it].nwl + 1); + + if (orbital_corr[it] == -1) + { + continue; + } + + for (int l = 0; l <= cell.atoms[it].nwl; l++) + { + const int N = cell.atoms[it].l_nchi[l]; + + occ_[iat][l].resize(N); + occ_save_[iat][l].resize(N); + + for (int n = 0; n < N; n++) + { + if (nspin == 1 || nspin == 2) + { + occ_[iat][l][n].resize(2); + occ_save_[iat][l][n].resize(2); + + occ_[iat][l][n][0].create(2 * l + 1, 2 * l + 1); + occ_[iat][l][n][1].create(2 * l + 1, 2 * l + 1); + + occ_save_[iat][l][n][0].create(2 * l + 1, 2 * l + 1); + occ_save_[iat][l][n][1].create(2 * l + 1, 2 * l + 1); + } + else if (nspin == 4) + { + occ_[iat][l][n].resize(1); + occ_save_[iat][l][n].resize(1); + + occ_[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol); + occ_save_[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol); + } + } + } + + for (int L = 0; L <= cell.atoms[it].nwl; L++) + { + iatlnmipol2iwt_[iat][L].resize(cell.atoms[it].l_nchi[L]); + + for (int n = 0; n < cell.atoms[it].l_nchi[L]; n++) + { + iatlnmipol2iwt_[iat][L][n].resize(2 * L + 1); + + for (int m = 0; m < 2 * L + 1; m++) + { + iatlnmipol2iwt_[iat][L][n][m].resize(npol); + } + } + } + + for (int iw = 0; iw < cell.atoms[it].nw * npol; iw++) + { + const int iw0 = iw / npol; + const int ipol = iw % npol; + const int iwt = cell.itiaiw2iwt(it, ia, iw); + const int l = cell.atoms[it].iw2l[iw0]; + const int n = cell.atoms[it].iw2n[iw0]; + const int m = cell.atoms[it].iw2m[iw0]; + + iatlnmipol2iwt_[iat][l][n][m][ipol] = iwt; + } + } + } +} + +void OccupationMatrix::get_flat(const int iat, const int l, std::vector& occ) const +{ + const int tlp1 = 2 * l + 1; + const int size = tlp1 * tlp1; + if (nspin_ == 2) + { + for (int is = 0; is < 2; is++) + { + for (int i = 0; i < size; i++) + { + occ[is * size + i] = occ_[iat][l][0][is].c[i]; + } + } + } + else + { + for (int i = 0; i < static_cast(occ.size()); i++) + { + occ[i] = occ_[iat][l][0][0].c[i]; + } + } +} + +void OccupationMatrix::set_flat(const int iat, const int l, const int spin, + const std::vector& occ) +{ + for (int i = 0; i < static_cast(occ.size()); i++) + { + occ_[iat][l][0][spin].c[i] = occ[i]; + } +} + +void OccupationMatrix::zero(const UnitCell& cell, const std::vector& orbital_corr) +{ + for (int T = 0; T < cell.ntype; T++) + { + if (orbital_corr[T] == -1) + { + continue; + } + + for (int I = 0; I < cell.atoms[T].na; I++) + { + const int iat = cell.itia2iat(T, I); + + for (int l = 0; l < cell.atoms[T].nwl + 1; l++) + { + const int N = cell.atoms[T].l_nchi[l]; + + for (int n = 0; n < N; n++) + { + if (nspin_ == 4) + { + occ_[iat][l][n][0].zero_out(); + } + else if (nspin_ == 1 || nspin_ == 2) + { + occ_[iat][l][n][0].zero_out(); + occ_[iat][l][n][1].zero_out(); + } + } + } + } + } +} + +void OccupationMatrix::copy_to_save(const UnitCell& cell, const std::vector& orbital_corr) +{ + ModuleBase::TITLE("OccupationMatrix", "copy_to_save"); + ModuleBase::timer::start("OccupationMatrix", "copy_to_save"); + + for (int T = 0; T < cell.ntype; T++) + { + const int target_l = orbital_corr[T]; + if (target_l == -1) + { + continue; + } + + for (int I = 0; I < cell.atoms[T].na; I++) + { + const int iat = cell.itia2iat(T, I); + + if (nspin_ == 4) + { + occ_save_[iat][target_l][0][0] = occ_[iat][target_l][0][0]; + } + else if (nspin_ == 1 || nspin_ == 2) + { + occ_save_[iat][target_l][0][0] = occ_[iat][target_l][0][0]; + occ_save_[iat][target_l][0][1] = occ_[iat][target_l][0][1]; + } + } + } + ModuleBase::timer::end("OccupationMatrix", "copy_to_save"); +} + +void OccupationMatrix::write_to_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + std::vector& uom) const +{ + if (uom.size() == 0) + { + return; + } + for (int iat = 0; iat < cell.nat; iat++) + { + const int it = cell.iat2it[iat]; + const int target_l = orbital_corr[it]; + if (target_l == -1) + { + continue; + } + const int size = (2 * target_l + 1) * (2 * target_l + 1); + + for (int mm = 0; mm < size; mm++) + { + uom[index[iat] + mm] = occ_[iat][target_l][0][0].c[mm]; + } + if (nspin_ == 2) + { + const int half_size = uom.size() / 2; + for (int mm = 0; mm < size; mm++) + { + uom[half_size + index[iat] + mm] = occ_[iat][target_l][0][1].c[mm]; + } + } + } +} + +void OccupationMatrix::read_from_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + const std::vector& uom) +{ + for (int T = 0; T < cell.ntype; T++) + { + const int l = orbital_corr[T]; + if (l == -1) + { + continue; + } + for (int I = 0; I < cell.atoms[T].na; I++) + { + const int iat = cell.itia2iat(T, I); + if (nspin_ == 4) + { + for (int mm = 0; mm < occ_[iat][l][0][0].nr * occ_[iat][l][0][0].nc; mm++) + { + occ_[iat][l][0][0].c[mm] = uom[index[iat] + mm]; + } + } + else if (nspin_ == 1 || nspin_ == 2) + { + const int half_size = uom.size() / 2; + for (int mm = 0; mm < occ_[iat][l][0][0].nr * occ_[iat][l][0][0].nc; mm++) + { + occ_[iat][l][0][0].c[mm] = uom[index[iat] + mm]; + if (nspin_ == 2) + { + occ_[iat][l][0][1].c[mm] = uom[half_size + index[iat] + mm]; + } + } + } + } + } +} + +void OccupationMatrix::write_save_to_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + std::vector& uom_save) const +{ + if (uom_save.size() == 0) + { + return; + } + for (int T = 0; T < cell.ntype; T++) + { + const int target_l = orbital_corr[T]; + if (target_l == -1) + { + continue; + } + + for (int I = 0; I < cell.atoms[T].na; I++) + { + const int iat = cell.itia2iat(T, I); + const int size = occ_[iat][target_l][0][0].nr * occ_[iat][target_l][0][0].nc; + + if (nspin_ == 4) + { + for (int mm = 0; mm < size; mm++) + { + uom_save[index[iat] + mm] = occ_[iat][target_l][0][0].c[mm]; + } + } + else if (nspin_ == 1 || nspin_ == 2) + { + for (int mm = 0; mm < size; mm++) + { + uom_save[index[iat] + mm] = occ_[iat][target_l][0][0].c[mm]; + } + if (nspin_ == 2) + { + const int half_size = uom_save.size() / 2; + for (int mm = 0; mm < size; mm++) + { + uom_save[half_size + index[iat] + mm] = occ_[iat][target_l][0][1].c[mm]; + } + } + } + } + } +} + +namespace elecstate +{ + +/// occ = beta * occ + (1-beta) * occ_save, applied to the correlated orbital +/// of every atom. nspin-aware: nspin=4 mixes the single Pauli block, +/// nspin=1/2 mixes both spin channels. Replaces the duplicated LCAO +/// k/gamma mixing loops. +void mix_occ_with_save(std::vector>>>& occ_mat, + const std::vector>>>& occ_mat_save, + const UnitCell& cell, + const std::vector& orbital_corr, + const int nspin, + const double beta) +{ + for (int T = 0; T < cell.ntype; T++) + { + const int target_l = orbital_corr[T]; + if (target_l == -1) + { + continue; + } + for (int I = 0; I < cell.atoms[T].na; I++) + { + const int iat = cell.itia2iat(T, I); + const int nchan = (nspin == 4) ? 1 : 2; + for (int is = 0; is < nchan; is++) + { + ModuleBase::matrix& occ = occ_mat[iat][target_l][0][is]; + const ModuleBase::matrix& occ_save = occ_mat_save[iat][target_l][0][is]; + const int size = occ.nr * occ.nc; + for (int mm = 0; mm < size; mm++) + { + occ.c[mm] = occ.c[mm] * beta + occ_save.c[mm] * (1.0 - beta); + } + } + } + } +} + +} // namespace elecstate diff --git a/source/source_estate/occ_matrix.h b/source/source_estate/occ_matrix.h new file mode 100644 index 00000000000..3ac30729fd2 --- /dev/null +++ b/source/source_estate/occ_matrix.h @@ -0,0 +1,129 @@ +#ifndef OCC_MATRIX_H +#define OCC_MATRIX_H + +#include "source_base/matrix.h" + +#include + +class UnitCell; + +/** + * @brief On-site occupation matrices for DFT+U. + * + * Owns the nested occ[iat][l][n][spin] matrices together with their saved + * copy (used by mixing) and the iat->(l,n,m,ipol)->iwt lookup table. + * Layout: + * nspin=1/2: occ[iat][l][n] has 2 spin channels of (2l+1)x(2l+1) + * nspin=4: occ[iat][l][n] has 1 channel of (2l+1)*npol x (2l+1)*npol + * (all Pauli blocks packed together) + */ +class OccupationMatrix +{ + public: + /// allocate occ/occ_save/iatlnmipol2iwt according to the cell + void init(const UnitCell& cell, + const std::vector& orbital_corr, + int nspin, + int npol); + + // --- element access --- + double get(int iat, int l, int n, int spin, int m1, int m2) const + { + return occ_[iat][l][n][spin](m1, m2); + } + double get_save(int iat, int l, int n, int spin, int m1, int m2) const + { + return occ_save_[iat][l][n][spin](m1, m2); + } + void set(int iat, int l, int n, int spin, int m1, int m2, double val) + { + occ_[iat][l][n][spin](m1, m2) = val; + } + + /// direct matrix access for kernels that operate on whole blocks + ModuleBase::matrix& mat(int iat, int l, int n, int spin) + { + return occ_[iat][l][n][spin]; + } + const ModuleBase::matrix& mat(int iat, int l, int n, int spin) const + { + return occ_[iat][l][n][spin]; + } + ModuleBase::matrix& mat_save(int iat, int l, int n, int spin) + { + return occ_save_[iat][l][n][spin]; + } + const ModuleBase::matrix& mat_save(int iat, int l, int n, int spin) const + { + return occ_save_[iat][l][n][spin]; + } + + // --- bulk data access (used by IO and legacy call sites) --- + std::vector>>>& data() { return occ_; } + const std::vector>>>& data() const { return occ_; } + std::vector>>>& data_save() { return occ_save_; } + const std::vector>>>& data_save() const { return occ_save_; } + + // --- lookup table --- + int iwt(int iat, int l, int n, int m, int ipol) const + { + return iatlnmipol2iwt_[iat][l][n][m][ipol]; + } + const std::vector>>>>& iatlnmipol2iwt() const + { + return iatlnmipol2iwt_; + } + + // --- flat (de)serialization of one atom's correlated orbital --- + /// nspin=1: fills occ with occ[iat][l][0][0] data + /// nspin=2: fills occ with interleaved spin-up then spin-down data + /// nspin=4: fills occ with occ[iat][l][0][0] data (all Pauli blocks) + void get_flat(int iat, int l, std::vector& occ) const; + void set_flat(int iat, int l, int spin, const std::vector& occ); + + // --- whole-array operations --- + void zero(const UnitCell& cell, const std::vector& orbital_corr); + void copy_to_save(const UnitCell& cell, const std::vector& orbital_corr); + + // --- flat mixing buffer (de)serialization over all atoms --- + /// write occ into uom at offsets given by index (split spin layout) + void write_to_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + std::vector& uom) const; + /// read occ from uom at offsets given by index (split spin layout) + void read_from_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + const std::vector& uom); + /// write occ_save into uom_save (skips when uom_save is empty) + void write_save_to_flat(const UnitCell& cell, + const std::vector& orbital_corr, + const std::vector& index, + std::vector& uom_save) const; + + int nspin() const { return nspin_; } + int npol() const { return npol_; } + + private: + std::vector>>> occ_; + std::vector>>> occ_save_; + std::vector>>>> iatlnmipol2iwt_; + int nspin_ = 0; + int npol_ = 0; +}; + +namespace elecstate +{ +/// occ = beta * occ + (1-beta) * occ_save on every atom's correlated orbital. +/// nspin-aware: nspin=4 mixes the single Pauli block, nspin=1/2 mixes both +/// spin channels. Replaces the duplicated LCAO k/gamma mixing loops. +void mix_occ_with_save(std::vector>>>& occ_mat, + const std::vector>>>& occ_mat_save, + const UnitCell& cell, + const std::vector& orbital_corr, + const int nspin, + const double beta); +} // namespace elecstate + +#endif diff --git a/source/source_hsolver/hsolver_pw_sdft.cpp b/source/source_hsolver/hsolver_pw_sdft.cpp index 339ac576a1c..eee59358463 100644 --- a/source/source_hsolver/hsolver_pw_sdft.cpp +++ b/source/source_hsolver/hsolver_pw_sdft.cpp @@ -1,6 +1,7 @@ #include "hsolver_pw_sdft.h" #include "source_base/global_function.h" +#include "source_base/parallel_comm.h" #include "source_base/parallel_device.h" #include "source_base/timer.h" #include "source_base/tool_title.h" diff --git a/source/source_hsolver/test/diago_bpcg_test.cpp b/source/source_hsolver/test/diago_bpcg_test.cpp index 37529ec60e7..dbd309041e3 100644 --- a/source/source_hsolver/test/diago_bpcg_test.cpp +++ b/source/source_hsolver/test/diago_bpcg_test.cpp @@ -7,6 +7,7 @@ #include "../diago_bpcg.h" #include "diago_mock.h" #include "mpi.h" +#include "source_base/parallel_comm.h" #include "source_basis/module_pw/test/test_tool.h" #include diff --git a/source/source_hsolver/test/test_hsolver_sdft.cpp b/source/source_hsolver/test/test_hsolver_sdft.cpp index 0820922019d..570f3ff9227 100644 --- a/source/source_hsolver/test/test_hsolver_sdft.cpp +++ b/source/source_hsolver/test/test_hsolver_sdft.cpp @@ -388,6 +388,7 @@ class TestHSolverPW_SDFT : public ::testing::Test // } #ifdef __MPI +#include "source_base/parallel_comm.h" #include "source_base/timer.h" #include "mpi.h" int main(int argc, char** argv) diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.h b/source/source_io/module_ctrl/ctrl_scf_lcao.h index b658253f12e..63688c08452 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.h +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.h @@ -9,7 +9,7 @@ #include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04 #include "source_hamilt/module_surchem/surchem.h" // use surchem (for dH veff pots) #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 #include "source_lcao/module_rdmft/rdmft.h" // use RDMFT codes #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008 diff --git a/source/source_io/module_hs/output_mat_sparse.h b/source/source_io/module_hs/output_mat_sparse.h index 0ce5e9ccf2c..7ee4ecd9068 100644 --- a/source/source_io/module_hs/output_mat_sparse.h +++ b/source/source_io/module_hs/output_mat_sparse.h @@ -6,7 +6,7 @@ #include "source_cell/klist.h" #include "source_hamilt/hamilt.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 namespace ModuleIO { diff --git a/source/source_io/module_hs/write_hs_r.h b/source/source_io/module_hs/write_hs_r.h index e0a0eab995a..0efac9b0f82 100644 --- a/source/source_io/module_hs/write_hs_r.h +++ b/source/source_io/module_hs/write_hs_r.h @@ -2,11 +2,13 @@ #define WRITE_HS_R_H #include "source_base/matrix.h" +#include "source_basis/module_ao/parallel_orbitals.h" #include "source_basis/module_nao/two_center_bundle.h" #include "source_cell/klist.h" +#include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_hamilt/hamilt.h" #include "source_lcao/lcao_hs_arrays.hpp" -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 #ifdef __EXX #include "RI/global/Tensor.h" // for RI::Tensor diff --git a/source/source_io/module_hs/write_vxc.hpp b/source/source_io/module_hs/write_vxc.hpp index 313d211cc05..ff794ea1488 100644 --- a/source/source_io/module_hs/write_vxc.hpp +++ b/source/source_io/module_hs/write_vxc.hpp @@ -4,7 +4,7 @@ #include "source_base/parallel_reduce.h" #include "source_base/module_container/base/third_party/blas.h" #include "source_base/module_external/scalapack_connector.h" -#include "source_lcao/module_dftu/dftu_lcao_op_legacy.h" +#include "source_lcao/module_dftu/dftu_nao_op_legacy.h" #include "source_lcao/module_operator_lcao/veff_lcao.h" #include "source_hamilt/module_xc/exx_info.h" #ifdef __EXX diff --git a/source/source_io/module_hs/write_vxc_r.hpp b/source/source_io/module_hs/write_vxc_r.hpp index 1c31008a2ea..3f132b16de5 100644 --- a/source/source_io/module_hs/write_vxc_r.hpp +++ b/source/source_io/module_hs/write_vxc_r.hpp @@ -2,7 +2,7 @@ #define __WRITE_VXC_R_H_ #include "source_io/module_parameter/parameter.h" #include "source_io/module_hs/write_hs_sparse.h" -#include "source_lcao/module_dftu/dftu_lcao_op_legacy.h" +#include "source_lcao/module_dftu/dftu_nao_op_legacy.h" #include "source_lcao/module_operator_lcao/veff_lcao.h" #include "source_lcao/spar_hsr.h" #ifdef __EXX diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h index f7f1c208444..233d07ea21c 100644 --- a/source/source_io/test/for_testing_input_conv.h +++ b/source/source_io/test/for_testing_input_conv.h @@ -13,7 +13,7 @@ #include "source_io/module_restart/restart.h" #include "source_io/module_unk/berryphase.h" #include "source_lcao/force_stress_lcao.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_md/md_func.h" #include "source_pw/module_pwdft/stru_fac.h" #include "source_pw/module_pwdft/vnl_pw.h" diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 4bdd62a3605..895e03c1378 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -1,8 +1,8 @@ #include "force_stress_lcao.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_dftu/dftu_lcao.h" //Quxin add for DFT+U on 20201029 -#include "source_lcao/module_dftu/dftu_force.h" +#include "source_lcao/module_dftu/dftu_nao.h" //Quxin add for DFT+U on 20201029 +#include "source_lcao/module_dftu/dftu_nao_fs_k.h" #include "source_io/module_output/output_log.h" #include "source_io/module_parameter/parameter.h" // new @@ -21,7 +21,7 @@ #include "source_lcao/module_deepks/lcao_deepks_io.h" // mohan add 2024-07-22 #include "source_lcao/module_deepks/deepks_force.h" #endif -#include "source_lcao/module_dftu/dftu_lcao_op.h" +#include "source_lcao/module_dftu/dftu_nao_op.h" #include "source_lcao/module_operator_lcao/dspin_lcao.h" #include "source_lcao/module_operator_lcao/nonlocal.h" #include "source_lcao/module_operator_lcao/ekinetic.h" diff --git a/source/source_lcao/force_stress_lcao.h b/source/source_lcao/force_stress_lcao.h index 02f62090547..0df810ad5d3 100644 --- a/source/source_lcao/force_stress_lcao.h +++ b/source/source_lcao/force_stress_lcao.h @@ -16,7 +16,8 @@ #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251010 #include "source_lcao/setup_dm.h" // mohan add 2025-11-03 -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 2025-11-07 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 2025-11-07 +#include "source_hamilt/hamilt.h" #include "source_hamilt/module_xc/exx_info.h" namespace vdw diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index c56c6e45ed1..8073516745b 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -3,7 +3,7 @@ #include "source_base/global_variable.h" #include "source_base/memory_recorder.h" #include "source_base/timer.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_lcao/setup_exx.h" #include "source_lcao/setup_deepks.h" #include "source_estate/module_dm/density_matrix.h" @@ -31,12 +31,12 @@ #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include "source_hsolver/hsolver_lcao.h" -#include "module_dftu/dftu_lcao_op.h" +#include "module_dftu/dftu_nao_op.h" #include "module_operator_lcao/dspin_lcao.h" #include "module_operator_lcao/ekinetic.h" #include "module_operator_lcao/meta_lcao.h" #include "module_operator_lcao/nonlocal.h" -#include "module_dftu/dftu_lcao_op_legacy.h" +#include "module_dftu/dftu_nao_op_legacy.h" #include "module_operator_lcao/op_exx_lcao.h" #include "module_operator_lcao/overlap.h" #include "module_operator_lcao/td_ekinetic_lcao.h" diff --git a/source/source_lcao/hamilt_lcao.h b/source/source_lcao/hamilt_lcao.h index 529b998c2e0..62ed73df3ce 100644 --- a/source/source_lcao/hamilt_lcao.h +++ b/source/source_lcao/hamilt_lcao.h @@ -22,7 +22,7 @@ namespace elecstate { template class DensityMatrix; } // Setup_DeePKS forward declaration, full definition in setup_deepks.h (moved to .cpp) // mohan add 20260605 template class Setup_DeePKS; -// Plus_U forward declaration, full definition in module_dftu/dftu_lcao.h (moved to .cpp) +// Plus_U forward declaration, full definition in module_dftu/dftu_nao.h (moved to .cpp) // mohan add 20260605 class Plus_U; diff --git a/source/source_lcao/lcao_set.h b/source/source_lcao/lcao_set.h index 6b9c5d81da6..7ed7a4cc9a0 100644 --- a/source/source_lcao/lcao_set.h +++ b/source/source_lcao/lcao_set.h @@ -12,7 +12,7 @@ #include "source_basis/module_pw/pw_basis.h" #include "source_hamilt/module_surchem/surchem.h" #include "source_pw/module_pwdft/vl_pw.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_lcao/setup_exx.h" #include "source_lcao/setup_deepks.h" diff --git a/source/source_lcao/module_dftu/CMakeLists.txt b/source/source_lcao/module_dftu/CMakeLists.txt index 370aa6b6262..479cda17c84 100644 --- a/source/source_lcao/module_dftu/CMakeLists.txt +++ b/source/source_lcao/module_dftu/CMakeLists.txt @@ -1,19 +1,20 @@ list(APPEND objects - dftu_lcao.cpp - dftu_force.cpp - dftu_yukawa.cpp - dftu_folding.cpp - dftu_lcao_pots.cpp - dftu_lcao_occ.cpp - dftu_lcao_energy.cpp + dftu_nao.cpp + dftu_nao_fs_k.cpp + dftu_nao_folding.cpp + dftu_nao_pots.cpp + dftu_nao_occ.cpp + dftu_nao_energy.cpp dftu_hamilt.cpp ) if(ENABLE_LCAO) list(APPEND objects - dftu_lcao_op.cpp - dftu_fs.cpp - dftu_lcao_op_legacy.cpp + dftu_nao_op.cpp + dftu_nao_fs_r.cpp + dftu_nao_for_r.cpp + dftu_nao_str_r.cpp + dftu_nao_op_legacy.cpp ) endif() diff --git a/source/source_lcao/module_dftu/dftu_fs.cpp b/source/source_lcao/module_dftu/dftu_fs.cpp deleted file mode 100644 index 477d8e25e49..00000000000 --- a/source/source_lcao/module_dftu/dftu_fs.cpp +++ /dev/null @@ -1,501 +0,0 @@ -#include "dftu_lcao_op.h" -#include "source_base/parallel_reduce.h" -#include "source_base/timer.h" - -namespace hamilt -{ - -template -void DFTU>::cal_force_stress(const bool cal_force, - const bool cal_stress, - ModuleBase::matrix& force, - ModuleBase::matrix& stress) -{ - ModuleBase::TITLE("DFTU", "cal_force_stress"); - if (this->dftu->get_dmr(0) == nullptr) - { - ModuleBase::WARNING_QUIT("DFTU", "dmr is not set"); - } - - // try to get the density matrix, if the density matrix is empty, skip the calculation and return - const hamilt::HContainer* dmR_tmp[this->nspin]; - dmR_tmp[0] = this->dftu->get_dmr(0); - - if (this->nspin == 2) - { - dmR_tmp[1] = this->dftu->get_dmr(1); - } - if (dmR_tmp[0]->size_atom_pairs() == 0) - { - return; - } - - // begin the calculation of force and stress - ModuleBase::timer::start("DFTU", "cal_force_stress"); - - const Parallel_Orbitals* pv = dmR_tmp[0]->get_paraV(); - const int npol = this->ucell->get_npol(); - std::vector stress_tmp(6, 0); - if (cal_force) - { - force.zero_out(); - } - // calculate atom_index for adjs_all, induced by omp parallel - int atom_index = 0; - std::vector atom_index_all(this->ucell->nat, -1); - for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) - { - int T0=0; - int I0=0; - ucell->iat2iait(iat0, &I0, &T0); - if(!this->dftu->has_correlated_orbital(T0)) - { - continue; - } - atom_index_all[iat0] = atom_index; - atom_index++; - } - - // 1. calculate for each pair of atoms - // loop over all on-site atoms - #pragma omp parallel - { - std::vector stress_local(6, 0); - ModuleBase::matrix force_local(force.nr, force.nc); - #pragma omp for schedule(dynamic) - for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) - { - // skip the atoms without plus-U - auto tau0 = ucell->get_tau(iat0); - int T0=0; - int I0=0; - ucell->iat2iait(iat0, &I0, &T0); - if (!this->dftu->has_correlated_orbital(T0)) - { - continue; - } - const int target_L = this->dftu->get_orbital_corr(T0); - const int tlp1 = 2 * target_L + 1; - AdjacentAtomInfo& adjs = this->adjs_all[atom_index_all[iat0]]; - - std::vector>> nlm_tot; - nlm_tot.resize(adjs.adj_num + 1); - - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int T1 = adjs.ntype[ad]; - const int I1 = adjs.natom[ad]; - const int iat1 = ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; - const Atom* atom1 = &ucell->atoms[T1]; - - auto all_indexes = pv->get_indexes_row(iat1); - auto col_indexes = pv->get_indexes_col(iat1); - // insert col_indexes into all_indexes to get universal set with no repeat elements - all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end()); - std::sort(all_indexes.begin(), all_indexes.end()); - all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end()); - for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += npol) - { - const int iw1 = all_indexes[iw1l] / npol; - std::vector> nlm; - // nlm is a vector of vectors, but size of outer vector is only 1 here - // If we are calculating force, we need also to store the gradient - // and size of outer vector is then 4 - // inner loop : all projectors (L0,M0) - int L1 = atom1->iw2l[iw1]; - int N1 = atom1->iw2n[iw1]; - int m1 = atom1->iw2m[iw1]; - - // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) - int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2; - - ModuleBase::Vector3 dtau = tau0 - tau1; - intor_->snap(T1, L1, N1, M1, T0, dtau * this->ucell->lat0, 1 /*cal_deri*/, nlm); - - // select the elements of nlm with target_L - std::vector nlm_target(tlp1 * 4); - for (int iw = 0; iw < this->ucell->atoms[T0].nw; iw++) - { - const int L0 = this->ucell->atoms[T0].iw2l[iw]; - if (L0 == target_L) - { - for (int m = 0; m < tlp1; m++) //-l, -l+1, ..., l-1, l - { - for (int n = 0; n < 4; n++) // value, deri_x, deri_y, deri_z - { - nlm_target[m + n * tlp1] = nlm[n][iw + m]; - // if(dtau.norm2 == 0.0) std::cout<<__FILE__<<__LINE__<<" "< occ(tlp1 * tlp1 * this->nspin, 0); - this->dftu->get_occ_mat_flat(iat0, target_L, occ); - - // calculate pot_onsite - const double u_value = this->dftu->get_u_current(T0); - std::vector pot_onsite(occ.size()); - double eu_tmp = 0; - this->cal_pot_onsite(occ, tlp1, u_value, &pot_onsite[0], eu_tmp); - - // second iteration to calculate force and stress - // calculate Force for atom J - // DMR_{I,J,R'-R} * U*(1/2*delta(m, m')-occ(m, m')) - // \frac{\partial }{\partial \tau_J} for each pair of atoms - // calculate Stress for strain tensor \varepsilon_{\alpha\beta} - // -1/Omega * DMR_{I,J,R'-R} * [ \frac{\partial }{\partial \tau_{J,\alpha}}\tau_{J,\beta} - // U*(1/2*delta(m, m')-occ(m, m')) - // + U*(1/2*delta(m, m')-occ(m, m')) - // \frac{\partial }{\partial \tau_{J,\alpha}}\tau_{J,\beta}] for each pair of atoms - for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) - { - const int T1 = adjs.ntype[ad1]; - const int I1 = adjs.natom[ad1]; - const int iat1 = ucell->itia2iat(T1, I1); - double* force_tmp1 = (cal_force) ? &force_local(iat1, 0) : nullptr; - double* force_tmp2 = (cal_force) ? &force_local(iat0, 0) : nullptr; - ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - ModuleBase::Vector3 dis1 = adjs.adjacent_tau[ad1] - tau0; - for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) - { - const int T2 = adjs.ntype[ad2]; - const int I2 = adjs.natom[ad2]; - const int iat2 = ucell->itia2iat(T2, I2); - ModuleBase::Vector3& R_index2 = adjs.box[ad2]; - ModuleBase::Vector3 dis2 = adjs.adjacent_tau[ad2] - tau0; - ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], - R_index2[1] - R_index1[1], - R_index2[2] - R_index1[2]); - const hamilt::BaseMatrix* tmp[this->nspin]; - tmp[0] = dmR_tmp[0]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); - if (this->nspin == 2) - { - tmp[1] = dmR_tmp[1]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); - } - // if not found , skip this pair of atoms - if (tmp[0] != nullptr) - { - // calculate force - if (cal_force) { - this->cal_force_IJR(iat1, - iat2, - pv, - nlm_tot[ad1], - nlm_tot[ad2], - pot_onsite, - tmp, - this->nspin, - force_tmp1, - force_tmp2); - } - - // calculate stress - if (cal_stress) { - this->cal_stress_IJR(iat1, - iat2, - pv, - nlm_tot[ad1], - nlm_tot[ad2], - pot_onsite, - tmp, - this->nspin, - dis1, - dis2, - stress_local.data()); - } - } - } - } - } - #pragma omp critical - { - if(cal_force) - { - force += force_local; - } - if(cal_stress) - { - for(int i = 0; i < 6; i++) - { - stress_tmp[i] += stress_local[i]; - } - } - } - } - - if (cal_force) - { -#ifdef __MPI - Parallel_Reduce::reduce_all(force.c, force.nr * force.nc); -#endif - if (this->nspin != 4) - { - for (int i = 0; i < force.nr * force.nc; i++) - { - force.c[i] *= 2.0; - } - } - } - - // stress renormalization - if (cal_stress) - { -#ifdef __MPI - // sum up the occupation matrix - Parallel_Reduce::reduce_all(stress_tmp.data(), 6); -#endif - const double weight = this->ucell->lat0 / this->ucell->omega; - for (int i = 0; i < 6; i++) - { - stress.c[i] = stress_tmp[i] * weight; - } - stress.c[8] = stress.c[5]; // stress(2,2) - stress.c[7] = stress.c[4]; // stress(2,1) - stress.c[6] = stress.c[2]; // stress(2,0) - stress.c[5] = stress.c[4]; // stress(1,2) - stress.c[4] = stress.c[3]; // stress(1,1) - stress.c[3] = stress.c[1]; // stress(1,0) - } - - ModuleBase::timer::end("DFTU", "cal_force_stress"); -} - - -template -void DFTU>::cal_force_IJR(const int& iat1, - const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - double* force1, - double* force2) -{ - // npol is the number of polarizations, - // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), - // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) - const int npol = this->ucell->get_npol(); - // --------------------------------------------- - // calculate the Nonlocal matrix for each pair of orbitals - // --------------------------------------------- - auto row_indexes = pv->get_indexes_row(iat1); - auto col_indexes = pv->get_indexes_col(iat2); - const int m_size = int(sqrt(pot_onsite_in.size() / nspin)); - const int m_size2 = m_size * m_size; - - // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 - std::vector step_trace(npol * npol, 0); - - if (npol == 2) - { - step_trace[1] = 1; - step_trace[2] = col_indexes.size(); - step_trace[3] = col_indexes.size() + 1; - } - - double tmp[3] = {0.0}; - // calculate the local matrix - for (int is = 0; is < nspin; is++) - { - const int is0 = nspin==2 ? is : 0; - const int step_is = nspin==4 ? is : 0; - const double* dm_pointer = dmR_pointer[is0]->get_pointer(); - for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) - { - const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; - for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) - { - const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; -#ifdef __DEBUG - assert(nlm1.size() == nlm2.size()); -#endif - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - tmp[0] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size] - * nlm2[m2] * dm_pointer[step_trace[step_is]]; - tmp[1] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 2] - * nlm2[m2] * dm_pointer[step_trace[step_is]]; - tmp[2] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 3] - * nlm2[m2] * dm_pointer[step_trace[step_is]]; - // force1 = - pot_onsite * * - // force2 = - pot_onsite * * } - force1[0] += tmp[0]; - force1[1] += tmp[1]; - force1[2] += tmp[2]; - force2[0] -= tmp[0]; - force2[1] -= tmp[1]; - force2[2] -= tmp[2]; - } - } - dm_pointer += npol; - } - dm_pointer += (npol - 1) * col_indexes.size(); - } - } -} - -template -void DFTU>::cal_stress_IJR(const int& iat1, - const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress) -{ - // npol is the number of polarizations, - // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), - // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) - const int npol = this->ucell->get_npol(); - // --------------------------------------------- - // calculate the Nonlocal matrix for each pair of orbitals - // --------------------------------------------- - auto row_indexes = pv->get_indexes_row(iat1); - auto col_indexes = pv->get_indexes_col(iat2); - const int m_size = int(sqrt(pot_onsite_in.size() / nspin)); - const int m_size2 = m_size * m_size; - - // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 - std::vector step_trace(npol * npol, 0); - - if (npol == 2) - { - step_trace[1] = 1; - step_trace[2] = col_indexes.size(); - step_trace[3] = col_indexes.size() + 1; - } - - // calculate the local matrix - for (int is = 0; is < nspin; is++) - { - const int is0 = nspin==2 ? is : 0; - const int step_is = nspin==4 ? is : 0; - const double* dm_pointer = dmR_pointer[is0]->get_pointer(); - for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) - { - const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; - for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) - { - const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; -#ifdef __DEBUG - assert(nlm1.size() == nlm2.size()); -#endif - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - double tmp = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * dm_pointer[step_trace[step_is]]; - // std::cout<<__FILE__<<__LINE__<<" "<>::cal_force_stress( - const bool cal_force, const bool cal_stress, - ModuleBase::matrix& force, ModuleBase::matrix& stress); -template void DFTU, double>>::cal_force_stress( - const bool cal_force, const bool cal_stress, - ModuleBase::matrix& force, ModuleBase::matrix& stress); -template void DFTU, std::complex>>::cal_force_stress( - const bool cal_force, const bool cal_stress, - ModuleBase::matrix& force, ModuleBase::matrix& stress); - -template void DFTU>::cal_force_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - double* force1, double* force2); -template void DFTU, double>>::cal_force_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - double* force1, double* force2); -template void DFTU, std::complex>>::cal_force_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - double* force1, double* force2); - -template void DFTU>::cal_stress_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress); -template void DFTU, double>>::cal_stress_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress); -template void DFTU, std::complex>>::cal_stress_IJR( - const int& iat1, const int& iat2, - const Parallel_Orbitals* pv, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& pot_onsite_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress); - -} // namespace hamilt diff --git a/source/source_lcao/module_dftu/dftu_hamilt.cpp b/source/source_lcao/module_dftu/dftu_hamilt.cpp index 787ef9acc57..11251f85bf3 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/source_lcao/module_dftu/dftu_hamilt.cpp @@ -1,8 +1,11 @@ -#include "dftu_lcao.h" +#include "dftu_nao.h" #include "dftu_hamilt.h" -#include "dftu_lcao_pots.h" +#include "dftu_nao_pots.h" +#include "source_base/global_function.h" #include "source_base/module_external/scalapack_connector.h" #include "source_base/timer.h" +#include "source_base/tool_title.h" +#include "source_basis/module_ao/parallel_orbitals.h" #ifdef __LCAO @@ -125,9 +128,9 @@ void pot_uterm_real(Plus_U& dftu, return; } -} // namespace DFTU_LCAO - -void Plus_U::cal_eff_pot_mat_R_double(const UnitCell& ucell, const Parallel_Orbitals* pv, const int ispin, double* SR, double* HR, const int npol) +/// @brief Accumulate the DFT+U term into the real-space HR (double). +/// Wraps pot_onsite_real plus the (pot_onsite*SR + SR*pot_onsite)/2 GEMM pair. +void pot_uterm_HR_real(const Plus_U& dftu, const UnitCell& ucell, const Parallel_Orbitals* pv, const int ispin, double* SR, double* HR, const int npol) { const char transN = 'N', transT = 'T'; const int one_int = 1; @@ -135,7 +138,7 @@ void Plus_U::cal_eff_pot_mat_R_double(const UnitCell& ucell, const Parallel_Orbi const int nlocal = pv->get_global_row_size(); std::vector pot_onsite(pv->nloc); - DFTU_LCAO::pot_onsite_real(*this, ucell, pv, ispin, true, &pot_onsite[0], npol); + pot_onsite_real(dftu, ucell, pv, ispin, true, &pot_onsite[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, @@ -158,7 +161,9 @@ void Plus_U::cal_eff_pot_mat_R_double(const UnitCell& ucell, const Parallel_Orbi return; } -void Plus_U::cal_eff_pot_mat_R_complex_double(const UnitCell& ucell, const Parallel_Orbitals* pv, const int ispin, std::complex* SR, std::complex* HR, const int npol) +/// @brief Accumulate the DFT+U term into the real-space HR (complex). +/// Wraps pot_onsite_complex plus the (pot_onsite*SR + SR*pot_onsite)/2 GEMM pair. +void pot_uterm_HR_complex(const Plus_U& dftu, const UnitCell& ucell, const Parallel_Orbitals* pv, const int ispin, std::complex* SR, std::complex* HR, const int npol) { const char transN = 'N', transT = 'T'; const int one_int = 1; @@ -166,7 +171,7 @@ void Plus_U::cal_eff_pot_mat_R_complex_double(const UnitCell& ucell, const Paral const int nlocal = pv->get_global_row_size(); std::vector> pot_onsite(pv->nloc); - DFTU_LCAO::pot_onsite_complex(*this, ucell, pv, ispin, true, &pot_onsite[0], npol); + pot_onsite_complex(dftu, ucell, pv, ispin, true, &pot_onsite[0], npol); #ifdef __MPI ScalapackConnector::gemm(transN, transN, @@ -189,4 +194,6 @@ void Plus_U::cal_eff_pot_mat_R_complex_double(const UnitCell& ucell, const Paral return; } +} // namespace DFTU_LCAO + #endif diff --git a/source/source_lcao/module_dftu/dftu_hamilt.h b/source/source_lcao/module_dftu/dftu_hamilt.h index 70e3c0cf34e..7b207c08f47 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.h +++ b/source/source_lcao/module_dftu/dftu_hamilt.h @@ -33,6 +33,26 @@ void pot_uterm_real(Plus_U& dftu, const double* sk, const int npol); +/// @brief Accumulate the DFT+U term into the real-space HR (double). +/// Wraps pot_onsite_real plus the (pot_onsite*SR + SR*pot_onsite)/2 GEMM pair. +void pot_uterm_HR_real(const Plus_U& dftu, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + const int ispin, + double* SR, + double* HR, + const int npol); + +/// @brief Accumulate the DFT+U term into the real-space HR (complex). +/// Wraps pot_onsite_complex plus the (pot_onsite*SR + SR*pot_onsite)/2 GEMM pair. +void pot_uterm_HR_complex(const Plus_U& dftu, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + const int ispin, + std::complex* SR, + std::complex* HR, + const int npol); + } // namespace DFTU_LCAO #endif diff --git a/source/source_lcao/module_dftu/dftu_lcao.cpp b/source/source_lcao/module_dftu/dftu_nao.cpp similarity index 92% rename from source/source_lcao/module_dftu/dftu_lcao.cpp rename to source/source_lcao/module_dftu/dftu_nao.cpp index a7833e001a8..f4a3bd6eed7 100644 --- a/source/source_lcao/module_dftu/dftu_lcao.cpp +++ b/source/source_lcao/module_dftu/dftu_nao.cpp @@ -1,8 +1,13 @@ -#include "dftu_lcao.h" +#include "dftu_nao.h" #include "source_base/tool_quit.h" #include "source_base/tool_title.h" #include "source_base/timer.h" +#include "source_basis/module_ao/parallel_orbitals.h" +#ifdef __LCAO +#include "source_basis/module_ao/orb_read.h" +#include "source_estate/module_dm/density_matrix.h" +#endif #include #include @@ -38,8 +43,6 @@ void Plus_U::init(UnitCell& cell, { ModuleBase::TITLE("Plus_U", "init"); - this->yukawa_lambda = yukawa_lambda; - #ifdef __LCAO ptr_orb_ = orb; if(ptr_orb_ != nullptr) @@ -67,6 +70,7 @@ void Plus_U::init(UnitCell& cell, nspin, orbital_corr, yukawa_potential, + yukawa_lambda, global_readin_dir, global_out_dir, init_chg, diff --git a/source/source_lcao/module_dftu/dftu_lcao.h b/source/source_lcao/module_dftu/dftu_nao.h similarity index 71% rename from source/source_lcao/module_dftu/dftu_lcao.h rename to source/source_lcao/module_dftu/dftu_nao.h index 780c7bfc9a7..fb0bef416c1 100644 --- a/source/source_lcao/module_dftu/dftu_lcao.h +++ b/source/source_lcao/module_dftu/dftu_nao.h @@ -1,21 +1,31 @@ #ifndef DFTU_LCAO_H #define DFTU_LCAO_H -#include "source_cell/klist.h" -#include "source_cell/unitcell.h" -#include "source_basis/module_ao/parallel_orbitals.h" #include "source_pw/module_pwdft/dftu_base.h" -#ifdef __LCAO -#include "source_basis/module_ao/orb_read.h" -#include "source_hamilt/hamilt.h" -#include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_estate/module_dm/density_matrix.h" -#endif +#include #include #include +class UnitCell; +class Parallel_Orbitals; + +#ifdef __LCAO +class LCAO_Orbitals; +namespace hamilt +{ +template +class HContainer; +} // namespace hamilt +namespace elecstate +{ +template +class DensityMatrix; +} // namespace elecstate +#endif + + class Plus_U : public Plus_U_Base { @@ -49,32 +59,15 @@ class Plus_U : public Plus_U_Base private: - double yukawa_lambda = 0.0; - #ifdef __LCAO const LCAO_Orbitals* ptr_orb_ = nullptr; std::vector orb_cutoff_; //============================================================= - // In dftu_hamilt.cpp + // In dftu_hamilt.cpp (DFTU_LCAO free functions) // For calculating contribution to Hamiltonian matrices //============================================================= public: - void cal_eff_pot_mat_R_double(const UnitCell& ucell, - const Parallel_Orbitals* pv, - const int ispin, - double* SR, - double* HR, - const int npol); - - void cal_eff_pot_mat_R_complex_double(const UnitCell& ucell, - const Parallel_Orbitals* pv, - const int ispin, - std::complex* SR, - std::complex* HR, - const int npol); - - /** * @brief get the density matrix of target spin * nspin = 1 and 4 : ispin should be 0 @@ -90,7 +83,6 @@ class Plus_U : public Plus_U_Base /// read-only accessors for state needed by DFTU_LCAO free functions const std::vector& get_orb_cutoff() const { return orb_cutoff_; } - double get_yukawa_lambda() const { return yukawa_lambda; } const LCAO_Orbitals* get_ptr_orb() const { return ptr_orb_; } private: diff --git a/source/source_lcao/module_dftu/dftu_lcao_energy.cpp b/source/source_lcao/module_dftu/dftu_nao_energy.cpp similarity index 85% rename from source/source_lcao/module_dftu/dftu_lcao_energy.cpp rename to source/source_lcao/module_dftu/dftu_nao_energy.cpp index 9abd260f36e..dcb1f5949d6 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_energy.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_energy.cpp @@ -1,8 +1,9 @@ -#include "dftu_lcao.h" -#include "dftu_lcao_energy.h" -#include "dftu_lcao_pots.h" +#include "dftu_nao.h" +#include "dftu_nao_energy.h" +#include "dftu_nao_pots.h" #include "source_base/timer.h" #include "source_base/tool_title.h" +#include "source_cell/unitcell.h" #include "source_io/module_parameter/parameter.h" #ifdef __LCAO @@ -68,16 +69,16 @@ void DFTU_LCAO::cal_energy_correction(Plus_U& dftu, const UnitCell& ucell) for (int m0 = 0; m0 < 2 * l + 1; m0++) { - nm_trace += dftu.get_occ_mat(iat, l, n, spin, m0, m0); + nm_trace += dftu.occmat().get(iat, l, n, spin, m0, m0); for (int m1 = 0; m1 < 2 * l + 1; m1++) { - nm2_trace += dftu.get_occ_mat(iat, l, n, spin, m0, m1) - * dftu.get_occ_mat(iat, l, n, spin, m1, m0); + nm2_trace += dftu.occmat().get(iat, l, n, spin, m0, m1) + * dftu.occmat().get(iat, l, n, spin, m1, m0); } } if (dftu.use_yukawa()) { - energy_u += 0.5 * (dftu.get_U_Yukawa(T, l, n) - dftu.get_J_Yukawa(T, l, n)) + energy_u += 0.5 * (dftu.yukawa().get_U(T, l, n) - dftu.yukawa().get_J(T, l, n)) * (nm_trace - nm2_trace); } else @@ -96,7 +97,7 @@ void DFTU_LCAO::cal_energy_correction(Plus_U& dftu, const UnitCell& ucell) for (int ipol0 = 0; ipol0 < npol; ipol0++) { const int m0_all = m0 + (2 * l + 1) * ipol0; - nm_trace += dftu.get_occ_mat(iat, l, n, 0, m0_all, m0_all); + nm_trace += dftu.occmat().get(iat, l, n, 0, m0_all, m0_all); for (int m1 = 0; m1 < 2 * l + 1; m1++) { @@ -104,15 +105,15 @@ void DFTU_LCAO::cal_energy_correction(Plus_U& dftu, const UnitCell& ucell) { int m1_all = m1 + (2 * l + 1) * ipol1; - nm2_trace += dftu.get_occ_mat(iat, l, n, 0, m0_all, m1_all) - * dftu.get_occ_mat(iat, l, n, 0, m1_all, m0_all); + nm2_trace += dftu.occmat().get(iat, l, n, 0, m0_all, m1_all) + * dftu.occmat().get(iat, l, n, 0, m1_all, m0_all); } } } } if (dftu.use_yukawa()) { - energy_u += 0.5 * (dftu.get_U_Yukawa(T, l, n) - dftu.get_J_Yukawa(T, l, n)) + energy_u += 0.5 * (dftu.yukawa().get_U(T, l, n) - dftu.yukawa().get_J(T, l, n)) * (nm_trace - nm2_trace); } else @@ -138,14 +139,14 @@ void DFTU_LCAO::cal_energy_correction(Plus_U& dftu, const UnitCell& ucell) { double pot_onsite = 0.0; pot_onsite = get_onsite_pot(dftu, T, iat, l, n, is, m1_all, m2_all, false); - energy_dc += pot_onsite * dftu.get_occ_mat(iat, l, n, is, m1_all, m2_all); + energy_dc += pot_onsite * dftu.occmat().get(iat, l, n, is, m1_all, m2_all); } } else if (nspin == 4) { double pot_onsite = 0.0; pot_onsite = get_onsite_pot(dftu, T, iat, l, n, 0, m1_all, m2_all, false); - energy_dc += pot_onsite * dftu.get_occ_mat(iat, l, n, 0, m1_all, m2_all); + energy_dc += pot_onsite * dftu.occmat().get(iat, l, n, 0, m1_all, m2_all); } } } diff --git a/source/source_lcao/module_dftu/dftu_lcao_energy.h b/source/source_lcao/module_dftu/dftu_nao_energy.h similarity index 100% rename from source/source_lcao/module_dftu/dftu_lcao_energy.h rename to source/source_lcao/module_dftu/dftu_nao_energy.h diff --git a/source/source_lcao/module_dftu/dftu_folding.cpp b/source/source_lcao/module_dftu/dftu_nao_folding.cpp similarity index 99% rename from source/source_lcao/module_dftu/dftu_folding.cpp rename to source/source_lcao/module_dftu/dftu_nao_folding.cpp index fc09e23ae84..71a50fab8ba 100644 --- a/source/source_lcao/module_dftu/dftu_folding.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_folding.cpp @@ -1,6 +1,6 @@ #ifdef __LCAO -#include "dftu_folding.h" -#include "dftu_lcao.h" +#include "dftu_nao_folding.h" +#include "dftu_nao.h" #include "source_base/timer.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_lcao/hamilt_lcao.h" diff --git a/source/source_lcao/module_dftu/dftu_folding.h b/source/source_lcao/module_dftu/dftu_nao_folding.h similarity index 100% rename from source/source_lcao/module_dftu/dftu_folding.h rename to source/source_lcao/module_dftu/dftu_nao_folding.h diff --git a/source/source_lcao/module_dftu/dftu_nao_for_r.cpp b/source/source_lcao/module_dftu/dftu_nao_for_r.cpp new file mode 100644 index 00000000000..452ddf663a0 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_for_r.cpp @@ -0,0 +1,126 @@ +/// @file dftu_nao_for_r.cpp +/// @brief DFT+U force calculation in real space (r-space) - implementation +/// +/// See dftu_nao_for_r.h for the mathematical formula and detailed documentation. + +#include "dftu_nao_for_r.h" +#include "dftu_nao_op.h" +#include "source_base/timer.h" + +namespace hamilt +{ + +template +void cal_for_IJR_nao_r(const DFTU>* dftu_op, + const int& iat1, + const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, + double* force2) +{ + // npol is the number of polarizations, + // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), + // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) + const int npol = dftu_op->get_ucell()->get_npol(); + + // --------------------------------------------- + // calculate the Nonlocal matrix for each pair of orbitals + // --------------------------------------------- + auto row_indexes = pv->get_indexes_row(iat1); + auto col_indexes = pv->get_indexes_col(iat2); + const int m_size = int(sqrt(pot_onsite_in.size() / nspin)); + const int m_size2 = m_size * m_size; + + // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 + std::vector step_trace(npol * npol, 0); + + if (npol == 2) + { + step_trace[1] = 1; + step_trace[2] = col_indexes.size(); + step_trace[3] = col_indexes.size() + 1; + } + + double tmp[3] = {0.0}; + // calculate the local matrix + for (int is = 0; is < nspin; is++) + { + const int is0 = nspin == 2 ? is : 0; + const int step_is = nspin == 4 ? is : 0; + const double* dm_pointer = dmR_pointer[is0]->get_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; +#ifdef __DEBUG + assert(nlm1.size() == nlm2.size()); +#endif + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + tmp[0] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size] + * nlm2[m2] * dm_pointer[step_trace[step_is]]; + tmp[1] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 2] + * nlm2[m2] * dm_pointer[step_trace[step_is]]; + tmp[2] = pot_onsite_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 3] + * nlm2[m2] * dm_pointer[step_trace[step_is]]; + // force1 = - pot_onsite * * + // force2 = - pot_onsite * * + force1[0] += tmp[0]; + force1[1] += tmp[1]; + force1[2] += tmp[2]; + force2[0] -= tmp[0]; + force2[1] -= tmp[1]; + force2[2] -= tmp[2]; + } + } + dm_pointer += npol; + } + dm_pointer += (npol - 1) * col_indexes.size(); + } + } +} + +// explicit template instantiation +template void cal_for_IJR_nao_r( + const DFTU>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); + +template void cal_for_IJR_nao_r, double>( + const DFTU, double>>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); + +template void cal_for_IJR_nao_r, std::complex>( + const DFTU, std::complex>>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); + +} // namespace hamilt diff --git a/source/source_lcao/module_dftu/dftu_nao_for_r.h b/source/source_lcao/module_dftu/dftu_nao_for_r.h new file mode 100644 index 00000000000..8c3f46896f6 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_for_r.h @@ -0,0 +1,90 @@ +/// @file dftu_nao_for_r.h +/// @brief DFT+U force calculation in real space (r-space) +/// +/// This file provides the real-space implementation of DFT+U force contribution +/// from a single atom pair (I,J,R). It is independent of k-point sampling because +/// the real-space density matrix (DMR) already contains the Brillouin-zone integration. +/// +/// Naming convention: _r suffix denotes real-space implementation, +/// corresponding to _k suffix for k-space (legacy) implementation. +/// +/// The force formula for atom pair (I,J,R) is: +/// +/// F_{J1} += sum_{m,m'} V_U_{mm'}(I) * +/// * d/d tau_{J1} * DMR_{mu,nu}(J1,J2,R) +/// +/// F_{J2} -= sum_{m,m'} V_U_{mm'}(I) * d/d tau_{J2} +/// * * DMR_{mu,nu}(J1,J2,R) +/// +/// where V_U_{mm'}(I) = U_eff * (delta_{mm'}/2 - n_{m'm}(I)) is the on-site +/// Hubbard potential, and the two-center integrals are pre-computed +/// by TwoCenterIntegrator::snap(). + +#ifndef DFTU_NAO_FOR_R_H +#define DFTU_NAO_FOR_R_H + +#include "source_basis/module_ao/parallel_orbitals.h" +#include "source_hamilt/module_hcontainer/hcontainer.h" + +#include +#include + +namespace hamilt +{ + +// Forward declarations to avoid circular dependency with dftu_lcao_op.h +template +class OperatorLCAO; + +template +class DFTU; + +/** + * @brief Compute DFT+U force contribution from a single atom pair (I,J,R) in real space + * + * For a given on-site atom I with correlated orbital l, and a pair of basis atoms (J1, J2) + * separated by lattice vector R, the force contribution is: + * + * F_{J1} += sum_{m,m'} V_U_{mm'}(I) * + * * d/d tau_{J1} * DMR_{mu,nu}(J1,J2,R) + * + * F_{J2} -= sum_{m,m'} V_U_{mm'}(I) * d/d tau_{J2} + * * * DMR_{mu,nu}(J1,J2,R) + * + * where mu runs over orbitals on J1, nu over orbitals on J2, and m,m' over the 2l+1 + * magnetic quantum numbers of the correlated orbital. + * + * The two-center integrals and their derivatives are pre-computed by + * TwoCenterIntegrator::snap() and stored in nlm_tot. + * + * @param iat1 [in] global atom index of J1 + * @param iat2 [in] global atom index of J2 + * @param pv [in] Parallel_Orbitals for basis index mapping + * @param nlm1_all [in] pre-computed and derivatives for atom J1 + * @param nlm2_all [in] pre-computed and derivatives for atom J2 + * @param pot_onsite [in] flattened V_U matrix: [m1*m_size + m2 + is*m_size2] + * @param dmR_pointer [in] pointer to DMR matrix blocks for each spin + * @param nspin [in] number of spin channels (1, 2, or 4) + * @param force1 [out] force accumulator for atom J1 (3 components) + * @param force2 [out] force accumulator for atom J2 (3 components) + * + * @note For nspin=1, the force is scaled by 2.0 at the caller level to account + * for spin degeneracy. For nspin=2, spin-up and spin-down are summed explicitly. + * For nspin=4 (non-collinear), the spinor structure is handled via npol=2 indexing. + */ +template +void cal_for_IJR_nao_r(const DFTU>* dftu_op, + const int& iat1, + const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, + double* force2); + +} // namespace hamilt + +#endif // DFTU_NAO_FOR_R_H diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_nao_fs_k.cpp similarity index 98% rename from source/source_lcao/module_dftu/dftu_force.cpp rename to source/source_lcao/module_dftu/dftu_nao_fs_k.cpp index 79f10f3a3d6..6fe8c0a68af 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_fs_k.cpp @@ -1,8 +1,8 @@ #ifdef __LCAO -#include "dftu_force.h" -#include "dftu_folding.h" -#include "dftu_lcao.h" -#include "dftu_lcao_pots.h" +#include "dftu_nao_fs_k.h" +#include "dftu_nao_folding.h" +#include "dftu_nao.h" +#include "dftu_nao_pots.h" #include "source_base/global_function.h" #include "source_base/module_external/scalapack_connector.h" #include "source_base/parallel_reduce.h" @@ -135,7 +135,7 @@ void force_stress(Plus_U& dftu, if (cal_force) { cal_force_gamma(nlocal, npol, - dftu.get_orbital_corr_vec(), dftu.get_iatlnmipol2iwt(), + dftu.get_orbital_corr_vec(), dftu.occmat().iatlnmipol2iwt(), ucell, &rho_pot_onsite[0], pv, fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z, force_dftu); } @@ -182,7 +182,7 @@ void force_stress(Plus_U& dftu, { cal_force_k(nlocal, npol, PARAM.inp.ks_solver, dftu.get_orb_cutoff(), - dftu.get_orbital_corr_vec(), dftu.get_iatlnmipol2iwt(), + dftu.get_orbital_corr_vec(), dftu.occmat().iatlnmipol2iwt(), ucell, gd, fsr, pv, ik, &rho_pot_onsite[0], force_dftu, kv.kvec_d[ik]); } if (cal_stress) diff --git a/source/source_lcao/module_dftu/dftu_force.h b/source/source_lcao/module_dftu/dftu_nao_fs_k.h similarity index 100% rename from source/source_lcao/module_dftu/dftu_force.h rename to source/source_lcao/module_dftu/dftu_nao_fs_k.h diff --git a/source/source_lcao/module_dftu/dftu_nao_fs_r.cpp b/source/source_lcao/module_dftu/dftu_nao_fs_r.cpp new file mode 100644 index 00000000000..35dfa562ee3 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_fs_r.cpp @@ -0,0 +1,283 @@ +/// @file dftu_nao_fs_r.cpp +/// @brief DFT+U force and stress unified entry in real space (r-space) - implementation +/// +/// See dftu_nao_fs_r.h for the mathematical formula and detailed documentation. + +#include "dftu_nao_fs_r.h" +#include "dftu_nao_for_r.h" +#include "dftu_nao_str_r.h" +#include "dftu_nao_op.h" +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" + +namespace hamilt +{ + +template +void cal_fs_nao_r(DFTU>* dftu_op, + const bool cal_force, + const bool cal_stress, + ModuleBase::matrix& force, + ModuleBase::matrix& stress) +{ + ModuleBase::TITLE("DFTU", "cal_fs_nao_r"); + if (dftu_op->get_dftu()->get_dmr(0) == nullptr) + { + ModuleBase::WARNING_QUIT("DFTU", "dmr is not set"); + } + + // try to get the density matrix, if the density matrix is empty, skip the calculation and return + const hamilt::HContainer* dmR_tmp[dftu_op->get_nspin()]; + dmR_tmp[0] = dftu_op->get_dftu()->get_dmr(0); + + if (dftu_op->get_nspin() == 2) + { + dmR_tmp[1] = dftu_op->get_dftu()->get_dmr(1); + } + if (dmR_tmp[0]->size_atom_pairs() == 0) + { + return; + } + + // begin the calculation of force and stress + ModuleBase::timer::start("DFTU", "cal_fs_nao_r"); + + const Parallel_Orbitals* pv = dmR_tmp[0]->get_paraV(); + const int npol = dftu_op->get_ucell()->get_npol(); + std::vector stress_tmp(6, 0); + if (cal_force) + { + force.zero_out(); + } + // calculate atom_index for adjs_all, induced by omp parallel + int atom_index = 0; + std::vector atom_index_all(dftu_op->get_ucell()->nat, -1); + for (int iat0 = 0; iat0 < dftu_op->get_ucell()->nat; iat0++) + { + int T0 = 0; + int I0 = 0; + dftu_op->get_ucell()->iat2iait(iat0, &I0, &T0); + if (!dftu_op->get_dftu()->has_correlated_orbital(T0)) + { + continue; + } + atom_index_all[iat0] = atom_index; + atom_index++; + } + + // 1. calculate for each pair of atoms + // loop over all on-site atoms +#pragma omp parallel + { + std::vector stress_local(6, 0); + ModuleBase::matrix force_local(force.nr, force.nc); +#pragma omp for schedule(dynamic) + for (int iat0 = 0; iat0 < dftu_op->get_ucell()->nat; iat0++) + { + // skip the atoms without plus-U + auto tau0 = dftu_op->get_ucell()->get_tau(iat0); + int T0 = 0; + int I0 = 0; + dftu_op->get_ucell()->iat2iait(iat0, &I0, &T0); + if (!dftu_op->get_dftu()->has_correlated_orbital(T0)) + { + continue; + } + const int target_L = dftu_op->get_dftu()->get_orbital_corr(T0); + const int tlp1 = 2 * target_L + 1; + AdjacentAtomInfo& adjs = dftu_op->get_adjs_all()[atom_index_all[iat0]]; + + std::vector>> nlm_tot; + nlm_tot.resize(adjs.adj_num + 1); + + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T1 = adjs.ntype[ad]; + const int I1 = adjs.natom[ad]; + const int iat1 = dftu_op->get_ucell()->itia2iat(T1, I1); + const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; + const Atom* atom1 = &dftu_op->get_ucell()->atoms[T1]; + + auto all_indexes = pv->get_indexes_row(iat1); + auto col_indexes = pv->get_indexes_col(iat1); + // insert col_indexes into all_indexes to get universal set with no repeat elements + all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end()); + std::sort(all_indexes.begin(), all_indexes.end()); + all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end()); + for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += npol) + { + const int iw1 = all_indexes[iw1l] / npol; + std::vector> nlm; + // nlm is a vector of vectors, but size of outer vector is only 1 here + // If we are calculating force, we need also to store the gradient + // and size of outer vector is then 4 + // inner loop : all projectors (L0,M0) + int L1 = atom1->iw2l[iw1]; + int N1 = atom1->iw2n[iw1]; + int m1 = atom1->iw2m[iw1]; + + // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) + int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2; + + ModuleBase::Vector3 dtau = tau0 - tau1; + dftu_op->get_intor()->snap(T1, L1, N1, M1, T0, dtau * dftu_op->get_ucell()->lat0, + 1 /*cal_deri*/, nlm); + + // select the elements of nlm with target_L + std::vector nlm_target(tlp1 * 4); + for (int iw = 0; iw < dftu_op->get_ucell()->atoms[T0].nw; iw++) + { + const int L0 = dftu_op->get_ucell()->atoms[T0].iw2l[iw]; + if (L0 == target_L) + { + for (int m = 0; m < tlp1; m++) //-l, -l+1, ..., l-1, l + { + for (int n = 0; n < 4; n++) // value, deri_x, deri_y, deri_z + { + nlm_target[m + n * tlp1] = nlm[n][iw + m]; + } + } + break; + } + } + nlm_tot[ad].insert({all_indexes[iw1l], nlm_target}); + } + } + // first iteration to calculate occupation matrix + std::vector occ(tlp1 * tlp1 * dftu_op->get_nspin(), 0); + dftu_op->get_dftu()->occmat().get_flat(iat0, target_L, occ); + + // calculate pot_onsite + const double u_value = dftu_op->get_dftu()->get_u_current(T0); + std::vector pot_onsite(occ.size()); + double eu_tmp = 0; + dftu_op->cal_pot_onsite(occ, tlp1, u_value, &pot_onsite[0], eu_tmp); + + // second iteration to calculate force and stress + // calculate Force for atom J + // DMR_{I,J,R'-R} * U*(1/2*delta(m, m')-occ(m, m')) + // d/d tau_J for each pair of atoms + // calculate Stress for strain tensor epsilon_{alpha,beta} + // -1/Omega * DMR_{I,J,R'-R} * [ d/d tau_{J,alpha} * tau_{J,beta} + // U*(1/2*delta(m, m')-occ(m, m')) + // + U*(1/2*delta(m, m')-occ(m, m')) + // d/d tau_{J,alpha} * tau_{J,beta} ] for each pair of atoms + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T1 = adjs.ntype[ad1]; + const int I1 = adjs.natom[ad1]; + const int iat1 = dftu_op->get_ucell()->itia2iat(T1, I1); + double* force_tmp1 = (cal_force) ? &force_local(iat1, 0) : nullptr; + double* force_tmp2 = (cal_force) ? &force_local(iat0, 0) : nullptr; + ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + ModuleBase::Vector3 dis1 = adjs.adjacent_tau[ad1] - tau0; + for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) + { + const int T2 = adjs.ntype[ad2]; + const int I2 = adjs.natom[ad2]; + const int iat2 = dftu_op->get_ucell()->itia2iat(T2, I2); + ModuleBase::Vector3& R_index2 = adjs.box[ad2]; + ModuleBase::Vector3 dis2 = adjs.adjacent_tau[ad2] - tau0; + ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], + R_index2[1] - R_index1[1], + R_index2[2] - R_index1[2]); + const hamilt::BaseMatrix* tmp[dftu_op->get_nspin()]; + tmp[0] = dmR_tmp[0]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); + if (dftu_op->get_nspin() == 2) + { + tmp[1] = dmR_tmp[1]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); + } + // if not found , skip this pair of atoms + if (tmp[0] != nullptr) + { + // calculate force + if (cal_force) + { + cal_for_IJR_nao_r(dftu_op, iat1, iat2, pv, + nlm_tot[ad1], nlm_tot[ad2], + pot_onsite, tmp, dftu_op->get_nspin(), + force_tmp1, force_tmp2); + } + + // calculate stress + if (cal_stress) + { + cal_str_IJR_nao_r(dftu_op, iat1, iat2, pv, + nlm_tot[ad1], nlm_tot[ad2], + pot_onsite, tmp, dftu_op->get_nspin(), + dis1, dis2, stress_local.data()); + } + } + } + } + } +#pragma omp critical + { + if (cal_force) + { + force += force_local; + } + if (cal_stress) + { + for (int i = 0; i < 6; i++) + { + stress_tmp[i] += stress_local[i]; + } + } + } + } + + if (cal_force) + { +#ifdef __MPI + Parallel_Reduce::reduce_all(force.c, force.nr * force.nc); +#endif + if (dftu_op->get_nspin() != 4) + { + for (int i = 0; i < force.nr * force.nc; i++) + { + force.c[i] *= 2.0; + } + } + } + + // stress renormalization + if (cal_stress) + { +#ifdef __MPI + // sum up the occupation matrix + Parallel_Reduce::reduce_all(stress_tmp.data(), 6); +#endif + const double weight = dftu_op->get_ucell()->lat0 / dftu_op->get_ucell()->omega; + for (int i = 0; i < 6; i++) + { + stress.c[i] = stress_tmp[i] * weight; + } + stress.c[8] = stress.c[5]; // stress(2,2) + stress.c[7] = stress.c[4]; // stress(2,1) + stress.c[6] = stress.c[2]; // stress(2,0) + stress.c[5] = stress.c[4]; // stress(1,2) + stress.c[4] = stress.c[3]; // stress(1,1) + stress.c[3] = stress.c[1]; // stress(1,0) + } + + ModuleBase::timer::end("DFTU", "cal_fs_nao_r"); +} + +// explicit template instantiation +template void cal_fs_nao_r( + DFTU>* dftu_op, + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void cal_fs_nao_r, double>( + DFTU, double>>* dftu_op, + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void cal_fs_nao_r, std::complex>( + DFTU, std::complex>>* dftu_op, + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +} // namespace hamilt diff --git a/source/source_lcao/module_dftu/dftu_nao_fs_r.h b/source/source_lcao/module_dftu/dftu_nao_fs_r.h new file mode 100644 index 00000000000..899f34d50a0 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_fs_r.h @@ -0,0 +1,79 @@ +/// @file dftu_nao_fs_r.h +/// @brief DFT+U force and stress unified entry in real space (r-space) +/// +/// This file provides the unified entry for DFT+U force/stress using the +/// real-space density matrix (DMR). It is independent of k-point sampling +/// because DMR already contains the Brillouin-zone integration. +/// +/// Naming convention: _r suffix denotes real-space implementation, +/// corresponding to _k suffix for k-space (legacy) implementation. +/// +/// The DFT+U force on atom J is derived from the Hubbard correction energy: +/// +/// E_U = (U_eff/2) * sum_{I,m,m',sigma} [ n^sigma_{mm'}(I) * (delta_{mm'} - n^sigma_{m'm}(I)) ] +/// +/// where n^sigma_{mm'}(I) is the on-site occupation matrix for correlated orbital l on atom I. +/// +/// The force on atom J is: +/// +/// F_J = -dE_U/d tau_J +/// = -sum_{I,R} sum_{m,m'} V_U_{mm'}(I) * [ +/// sum_{mu,nu} DMR_{mu,nu}(I,R) * d/d tau_J * +/// ] +/// +/// For stress, the derivative is with respect to strain tensor epsilon_{alpha,beta}: +/// +/// sigma_{alpha,beta} = -(1/Omega) * dE_U/d epsilon_{alpha,beta} +/// = -(1/Omega) * sum_{I,R} sum_{m,m'} V_U_{mm'}(I) * [ +/// sum_{mu,nu} DMR_{mu,nu}(I,R) * ( +/// d/d epsilon_{alpha,beta} * * R_beta +/// + * d/d epsilon_{alpha,beta} * R_beta +/// ) +/// ] + +#ifndef DFTU_NAO_FS_R_H +#define DFTU_NAO_FS_R_H + +#include "source_base/matrix.h" + +namespace hamilt +{ + +// Forward declarations to avoid circular dependency with dftu_lcao_op.h +template +class OperatorLCAO; + +template +class DFTU; + +/** + * @brief Calculate DFT+U force and stress in real space (unified for gamma-only and multik) + * + * This is the unified entry for DFT+U force/stress calculation. It loops over all + * on-site atoms with correlated orbitals, computes the two-center integrals + * via TwoCenterIntegrator, and accumulates force/stress contributions from all + * atom pairs (I,J,R) using OpenMP parallelization. + * + * @note This implementation uses the real-space density matrix DMR (HContainer) + * and two-center integrals computed by TwoCenterIntegrator. It is + * independent of k-point sampling because DMR already contains the BZ integration. + * + * @param dftu_op [in] pointer to the DFTU operator object (for accessing ucell, dftu, intor_) + * @param cal_force [in] whether to compute force + * @param cal_stress [in] whether to compute stress + * @param force [out] force matrix (nat, 3), accumulated + * @param stress [out] stress matrix (3, 3), accumulated + * + * @warning The density matrix must be set via Plus_U::set_dmr() before calling this. + * If get_dmr(0) returns nullptr, the function aborts with WARNING_QUIT. + */ +template +void cal_fs_nao_r(DFTU>* dftu_op, + const bool cal_force, + const bool cal_stress, + ModuleBase::matrix& force, + ModuleBase::matrix& stress); + +} // namespace hamilt + +#endif // DFTU_NAO_FS_R_H diff --git a/source/source_lcao/module_dftu/dftu_lcao_occ.cpp b/source/source_lcao/module_dftu/dftu_nao_occ.cpp similarity index 58% rename from source/source_lcao/module_dftu/dftu_lcao_occ.cpp rename to source/source_lcao/module_dftu/dftu_nao_occ.cpp index 247c8b1174f..291b21a7c20 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_occ.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_occ.cpp @@ -1,16 +1,17 @@ -#include "dftu_lcao_occ.h" -#include "dftu_lcao.h" -#include "dftu_folding.h" +#include "dftu_nao_occ.h" +#include "dftu_nao.h" +#include "dftu_nao_folding.h" #include "source_base/timer.h" #include "source_base/module_external/scalapack_connector.h" +#include "source_estate/occ_matrix.h" #include "source_io/module_parameter/parameter.h" #ifdef __LCAO #include "source_lcao/hamilt_lcao.h" #endif -// copy_occ_mat(), zero_occ_mat(), mix_occ_mat(), set_occ_mat(ucell), -// get_occ_mat_flat(), set_occ_mat_flat() -// are now implemented in dftu_base.cpp as Plus_U_Base methods (inherited by Plus_U). +// cal_occ_mat_k / cal_occ_mat_gamma take Plus_U& dftu directly and read all +// occupation-matrix state (occ/save arrays, lookup table, nspin/npol, and the +// occ_mat_initialized flag) from dftu.occmat() and the Plus_U_Base accessors. #ifdef __LCAO @@ -22,63 +23,21 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, const double& mixing_beta, hamilt::Hamilt>* p_ham, const bool gamma_only_local, - const int nspin, - const int npol, - const int nlocal, - const std::string& ks_solver, - const std::vector>>>>& iatlnmipol2iwt, - const std::vector& orbital_corr, - std::vector>>>& occ_mat, - std::vector>>>& occ_mat_save, - bool& occ_mat_initialized) + Plus_U& dftu) { ModuleBase::TITLE("DFTU_LCAO", "cal_occ_mat_k"); ModuleBase::timer::start("DFTU_LCAO", "cal_occ_mat_k"); - // copy occ_mat to occ_mat_save - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = orbital_corr[T]; - if (target_l == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - if (nspin == 4) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - } - else if (nspin == 1 || nspin == 2) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - occ_mat_save[iat][target_l][0][1] = occ_mat[iat][target_l][0][1]; - } - } - } - // zero occ_mat - for (int T = 0; T < ucell.ntype; T++) - { - if (orbital_corr[T] == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - for (int l = 0; l < ucell.atoms[T].nwl + 1; l++) - { - const int N = ucell.atoms[T].l_nchi[l]; - for (int n = 0; n < N; n++) - { - if (nspin == 4) - { - occ_mat[iat][l][n][0].zero_out(); - } - else if (nspin == 1 || nspin == 2) - { - occ_mat[iat][l][n][0].zero_out(); - occ_mat[iat][l][n][1].zero_out(); - } - } - } - } - } + const int nspin = dftu.occmat().nspin(); + const int npol = dftu.occmat().npol(); + const int nlocal = pv->get_global_row_size(); + const std::string& ks_solver = PARAM.inp.ks_solver; + const auto& iatlnmipol2iwt = dftu.occmat().iatlnmipol2iwt(); + const std::vector& orbital_corr = dftu.get_orbital_corr_vec(); + + // copy occ_mat to occ_mat_save, then zero occ_mat + dftu.occmat().copy_to_save(ucell, orbital_corr); + dftu.occmat().zero(ucell, orbital_corr); //=================Part 1====================== // call SCALAPACK routine to calculate the product of the S and density matrix @@ -160,6 +119,7 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, } // Calculate the local occupation number matrix + ModuleBase::matrix& occ = dftu.occmat().mat(iat, l, n, spin); for (int m0 = 0; m0 < 2 * l + 1; m0++) { for (int ipol0 = 0; ipol0 < npol; ipol0++) @@ -184,12 +144,12 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, if ((nu >= 0) && (mu >= 0)) { - occ_mat[iat][l][n][spin](m0_all, m1_all) += (srho[irc]).real() / 4.0; + occ(m0_all, m1_all) += (srho[irc]).real() / 4.0; } if ((nu_prime >= 0) && (mu_prime >= 0)) { - occ_mat[iat][l][n][spin](m0_all, m1_all) + occ(m0_all, m1_all) += (std::conj(srho[irc_prime])).real() / 4.0; } } // ipol1 @@ -237,9 +197,10 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, #ifdef __MPI if (nspin == 1 || nspin == 4) { - ModuleBase::matrix temp(occ_mat[iat][l][n][0]); + ModuleBase::matrix& occ0 = dftu.occmat().mat(iat, l, n, 0); + ModuleBase::matrix temp(occ0); MPI_Allreduce(&temp(0, 0), - &occ_mat[iat][l][n][0](0, 0), + &occ0(0, 0), (2 * l + 1) * npol * (2 * l + 1) * npol, MPI_DOUBLE, MPI_SUM, @@ -247,17 +208,19 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, } else if (nspin == 2) { - ModuleBase::matrix temp0(occ_mat[iat][l][n][0]); + ModuleBase::matrix& occ0 = dftu.occmat().mat(iat, l, n, 0); + ModuleBase::matrix temp0(occ0); MPI_Allreduce(&temp0(0, 0), - &occ_mat[iat][l][n][0](0, 0), + &occ0(0, 0), (2 * l + 1) * (2 * l + 1), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - ModuleBase::matrix temp1(occ_mat[iat][l][n][1]); + ModuleBase::matrix& occ1 = dftu.occmat().mat(iat, l, n, 1); + ModuleBase::matrix temp1(occ1); MPI_Allreduce(&temp1(0, 0), - &occ_mat[iat][l][n][1](0, 0), + &occ1(0, 0), (2 * l + 1) * (2 * l + 1), MPI_DOUBLE, MPI_SUM, @@ -268,19 +231,28 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, switch (nspin) { case 1: - occ_mat[iat][l][n][0] += transpose(occ_mat[iat][l][n][0]); - occ_mat[iat][l][n][0] *= 0.5; - occ_mat[iat][l][n][1] += occ_mat[iat][l][n][0]; + { + ModuleBase::matrix& occ0 = dftu.occmat().mat(iat, l, n, 0); + occ0 += transpose(occ0); + occ0 *= 0.5; + dftu.occmat().mat(iat, l, n, 1) += occ0; break; + } case 2: for (int is = 0; is < nspin; is++) - occ_mat[iat][l][n][is] += transpose(occ_mat[iat][l][n][is]); + { + ModuleBase::matrix& occ_is = dftu.occmat().mat(iat, l, n, is); + occ_is += transpose(occ_is); + } break; case 4: - occ_mat[iat][l][n][0] += transpose(occ_mat[iat][l][n][0]); + { + ModuleBase::matrix& occ0 = dftu.occmat().mat(iat, l, n, 0); + occ0 += transpose(occ0); break; + } default: std::cout << "Not supported NSPIN parameter" << std::endl; @@ -291,41 +263,13 @@ void DFTU_LCAO::cal_occ_mat_k(const Parallel_Orbitals* pv, } // end ia } // end it - if(PARAM.inp.mixing_dftu && occ_mat_initialized) + if(PARAM.inp.mixing_dftu && dftu.is_occ_mat_initialized()) { - double beta = mixing_beta; - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = orbital_corr[T]; - if (target_l == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - if (nspin == 4) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta - + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - } - } - else if (nspin == 1 || nspin == 2) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta - + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - occ_mat[iat][target_l][0][1].c[mm] = occ_mat[iat][target_l][0][1].c[mm] * beta - + occ_mat_save[iat][target_l][0][1].c[mm] * (1.0 - beta); - } - } - } - } + elecstate::mix_occ_with_save(dftu.occmat().data(), dftu.occmat().data_save(), + ucell, orbital_corr, nspin, mixing_beta); } - occ_mat_initialized = true; + dftu.mark_occ_mat_initialized(); ModuleBase::timer::end("DFTU_LCAO", "cal_occ_mat_k"); return; } @@ -336,61 +280,20 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, const std::vector> &dm_gamma, const double& mixing_beta, hamilt::Hamilt* p_ham, - const int nspin, - const int npol, - const int nlocal, - const std::vector>>>>& iatlnmipol2iwt, - const std::vector& orbital_corr, - std::vector>>>& occ_mat, - std::vector>>>& occ_mat_save, - bool& occ_mat_initialized) + Plus_U& dftu) { ModuleBase::TITLE("DFTU_LCAO", "cal_occ_mat_gamma"); ModuleBase::timer::start("DFTU_LCAO", "cal_occ_mat_gamma"); - // copy occ_mat to occ_mat_save - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = orbital_corr[T]; - if (target_l == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - if (nspin == 4) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - } - else if (nspin == 1 || nspin == 2) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - occ_mat_save[iat][target_l][0][1] = occ_mat[iat][target_l][0][1]; - } - } - } - // zero occ_mat - for (int T = 0; T < ucell.ntype; T++) - { - if (orbital_corr[T] == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - for (int l = 0; l < ucell.atoms[T].nwl + 1; l++) - { - const int N = ucell.atoms[T].l_nchi[l]; - for (int n = 0; n < N; n++) - { - if (nspin == 4) - { - occ_mat[iat][l][n][0].zero_out(); - } - else if (nspin == 1 || nspin == 2) - { - occ_mat[iat][l][n][0].zero_out(); - occ_mat[iat][l][n][1].zero_out(); - } - } - } - } - } + + const int nspin = dftu.occmat().nspin(); + const int npol = dftu.occmat().npol(); + const int nlocal = pv->get_global_row_size(); + const auto& iatlnmipol2iwt = dftu.occmat().iatlnmipol2iwt(); + const std::vector& orbital_corr = dftu.get_orbital_corr_vec(); + + // copy occ_mat to occ_mat_save, then zero occ_mat + dftu.occmat().copy_to_save(ucell, orbital_corr); + dftu.occmat().zero(ucell, orbital_corr); //=================Part 1====================== // call PBLAS routine to calculate the product of the S and density matrix @@ -456,6 +359,7 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, } // Calculate the local occupation number matrix + ModuleBase::matrix& occ_is = dftu.occmat().mat(iat, l, n, is); for (int m0 = 0; m0 < 2 * l + 1; m0++) { for (int ipol0 = 0; ipol0 < npol; ipol0++) @@ -480,7 +384,7 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, int m0_all = m0 + (2 * l + 1) * ipol0; int m1_all = m0 + (2 * l + 1) * ipol1; - occ_mat[iat][l][n][is](m0, m1) += srho[irc] / 4.0; + occ_is(m0, m1) += srho[irc] / 4.0; } if ((nu_prime >= 0) && (mu_prime >= 0)) @@ -488,18 +392,18 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, int m0_all = m0 + (2 * l + 1) * ipol0; int m1_all = m0 + (2 * l + 1) * ipol1; - occ_mat[iat][l][n][is](m0, m1) += srho[irc_prime] / 4.0; + occ_is(m0, m1) += srho[irc_prime] / 4.0; } } } } } - ModuleBase::matrix temp(occ_mat[iat][l][n][is]); + ModuleBase::matrix temp(occ_is); #ifdef __MPI MPI_Allreduce(&temp(0, 0), - &occ_mat[iat][l][n][is](0, 0), + &occ_is(0, 0), (2 * l + 1) * npol * (2 * l + 1) * npol, MPI_DOUBLE, MPI_SUM, @@ -510,13 +414,16 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, switch (nspin) { case 1: - occ_mat[iat][l][n][0] += transpose(occ_mat[iat][l][n][0]); - occ_mat[iat][l][n][0] *= 0.5; - occ_mat[iat][l][n][1] += occ_mat[iat][l][n][0]; + { + ModuleBase::matrix& occ0 = dftu.occmat().mat(iat, l, n, 0); + occ0 += transpose(occ0); + occ0 *= 0.5; + dftu.occmat().mat(iat, l, n, 1) += occ0; break; + } case 2: - occ_mat[iat][l][n][is] += transpose(occ_mat[iat][l][n][is]); + occ_is += transpose(occ_is); break; default: @@ -530,41 +437,13 @@ void DFTU_LCAO::cal_occ_mat_gamma(const Parallel_Orbitals* pv, } // it } // is - if(PARAM.inp.mixing_dftu && occ_mat_initialized) + if(PARAM.inp.mixing_dftu && dftu.is_occ_mat_initialized()) { - double beta = mixing_beta; - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = orbital_corr[T]; - if (target_l == -1) continue; - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - if (nspin == 4) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta - + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - } - } - else if (nspin == 1 || nspin == 2) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta - + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - occ_mat[iat][target_l][0][1].c[mm] = occ_mat[iat][target_l][0][1].c[mm] * beta - + occ_mat_save[iat][target_l][0][1].c[mm] * (1.0 - beta); - } - } - } - } + elecstate::mix_occ_with_save(dftu.occmat().data(), dftu.occmat().data_save(), + ucell, orbital_corr, nspin, mixing_beta); } - occ_mat_initialized = true; + dftu.mark_occ_mat_initialized(); ModuleBase::timer::end("DFTU_LCAO", "cal_occ_mat_gamma"); return; } @@ -584,13 +463,7 @@ void cal_occ_mat(const Parallel_Orbitals* pv, const bool gamma_only_local, const int nspin) { - bool occ_mat_initialized = dftu.get_occ_mat_initialized(); - DFTU_LCAO::cal_occ_mat_gamma(pv, iter, ucell, dm, mixing_beta, p_ham, nspin, - ucell.get_npol(), pv->get_global_row_size(), dftu.get_iatlnmipol2iwt(), - dftu.get_orbital_corr_vec(), - dftu.get_occ_mat_data(), dftu.get_occ_mat_save_data(), - occ_mat_initialized); - dftu.set_occ_mat_initialized(occ_mat_initialized); + DFTU_LCAO::cal_occ_mat_gamma(pv, iter, ucell, dm, mixing_beta, p_ham, dftu); } //! dftu occupation matrix for multiple k-points using dm(complex) @@ -606,13 +479,7 @@ void cal_occ_mat(const Parallel_Orbitals* pv, const bool gamma_only_local, const int nspin) { - bool occ_mat_initialized = dftu.get_occ_mat_initialized(); - DFTU_LCAO::cal_occ_mat_k(pv, iter, ucell, dm, kv, mixing_beta, p_ham, gamma_only_local, nspin, - ucell.get_npol(), pv->get_global_row_size(), PARAM.inp.ks_solver, dftu.get_iatlnmipol2iwt(), - dftu.get_orbital_corr_vec(), - dftu.get_occ_mat_data(), dftu.get_occ_mat_save_data(), - occ_mat_initialized); - dftu.set_occ_mat_initialized(occ_mat_initialized); + DFTU_LCAO::cal_occ_mat_k(pv, iter, ucell, dm, kv, mixing_beta, p_ham, gamma_only_local, dftu); } } // namespace DFTU_LCAO diff --git a/source/source_lcao/module_dftu/dftu_lcao_occ.h b/source/source_lcao/module_dftu/dftu_nao_occ.h similarity index 65% rename from source/source_lcao/module_dftu/dftu_lcao_occ.h rename to source/source_lcao/module_dftu/dftu_nao_occ.h index afcb4c45b13..eb2619499eb 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_occ.h +++ b/source/source_lcao/module_dftu/dftu_nao_occ.h @@ -43,15 +43,7 @@ void cal_occ_mat_k(const Parallel_Orbitals* pv, const double& mixing_beta, hamilt::Hamilt>* p_ham, const bool gamma_only_local, - const int nspin, - const int npol, - const int nlocal, - const std::string& ks_solver, - const std::vector>>>>& iatlnmipol2iwt, - const std::vector& orbital_corr, - std::vector>>>& occ_mat, - std::vector>>>& occ_mat_save, - bool& occ_mat_initialized); + Plus_U& dftu); // calculate the local occupation number matrix (gamma-point version) void cal_occ_mat_gamma(const Parallel_Orbitals* pv, @@ -60,14 +52,7 @@ void cal_occ_mat_gamma(const Parallel_Orbitals* pv, const std::vector>& dm_gamma, const double& mixing_beta, hamilt::Hamilt* p_ham, - const int nspin, - const int npol, - const int nlocal, - const std::vector>>>>& iatlnmipol2iwt, - const std::vector& orbital_corr, - std::vector>>>& occ_mat, - std::vector>>>& occ_mat_save, - bool& occ_mat_initialized); + Plus_U& dftu); } // namespace DFTU_LCAO #endif diff --git a/source/source_lcao/module_dftu/dftu_lcao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp similarity index 96% rename from source/source_lcao/module_dftu/dftu_lcao_op.cpp rename to source/source_lcao/module_dftu/dftu_nao_op.cpp index 0f4fbd15ffe..c23ed2be5e0 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -1,4 +1,4 @@ -#include "dftu_lcao_op.h" +#include "dftu_nao_op.h" #include "source_base/timer.h" #include "source_base/tool_title.h" @@ -7,6 +7,9 @@ #include "source_io/module_parameter/parameter.h" #include "source_base/parallel_reduce.h" +// Include the free function implementations for force/stress in real space +#include "dftu_nao_fs_r.h" + template hamilt::DFTU>::DFTU(HS_Matrix_K* hsk_in, const std::vector>& kvec_d_in, @@ -324,7 +327,7 @@ void hamilt::DFTU>::contributeHR() { for (auto& v : occ) { v *= 0.5; } } - this->dftu->set_occ_mat_flat(iat0, target_L, this->current_spin, occ); + this->dftu->occmat().set_flat(iat0, target_L, this->current_spin, occ); } // ============================================================ // BRANCH 2: Occ_mat IS initialized (use pre-read data) @@ -343,7 +346,7 @@ void hamilt::DFTU>::contributeHR() // For nspin=4, occ_mat is stored as 4 stacked tlp1^2 blocks // at offsets 0, tlp1^2, 2*tlp1^2, 3*tlp1^2 for the 4 Pauli channels. // Use get_occ_mat_flat to read the stacked blocks directly - this->dftu->get_occ_mat_flat(iat0, target_L, occ); + this->dftu->occmat().get_flat(iat0, target_L, occ); } // nspin=1 or nspin=2: Collinear spin case // Occ_mat stored separately for each spin channel @@ -354,7 +357,7 @@ void hamilt::DFTU>::contributeHR() // TODO: UNSAFE - current_spin must be correct for nspin=2. // If current_spin is not toggled properly, wrong spin channel's occ_mat is read. // This can happen if contributeHR() is called out of expected order. - occ[i] = this->dftu->get_occ_mat(iat0, target_L, 0, this->current_spin, + occ[i] = this->dftu->occmat().get(iat0, target_L, 0, this->current_spin, i / (2 * target_L + 1), i % (2 * target_L + 1)); } } @@ -682,6 +685,17 @@ void hamilt::DFTU>::cal_pot_onsite(const std::vecto } } +// cal_force_stress(): thin wrapper calling the real-space free function implementation +// See dftu_nao_fs_r.cpp for the actual implementation and mathematical formulas +template +void hamilt::DFTU>::cal_force_stress(const bool cal_force, + const bool cal_stress, + ModuleBase::matrix& force, + ModuleBase::matrix& stress) +{ + cal_fs_nao_r(this, cal_force, cal_stress, force, stress); +} + template class hamilt::DFTU>; template class hamilt::DFTU, double>>; template class hamilt::DFTU, std::complex>>; diff --git a/source/source_lcao/module_dftu/dftu_lcao_op.h b/source/source_lcao/module_dftu/dftu_nao_op.h similarity index 93% rename from source/source_lcao/module_dftu/dftu_lcao_op.h rename to source/source_lcao/module_dftu/dftu_nao_op.h index d244ebf37d5..d7f6c51d16f 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_op.h +++ b/source/source_lcao/module_dftu/dftu_nao_op.h @@ -1,11 +1,11 @@ -#ifndef DFTPLUSU_H -#define DFTPLUSU_H +#ifndef DFTU_NAO_OP_H +#define DFTU_NAO_OP_H #include "source_basis/module_ao/parallel_orbitals.h" #include "source_basis/module_nao/two_center_integrator.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_cell/unitcell.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_nao.h" #include "source_hamilt/module_hcontainer/hcontainer.h" #include @@ -60,6 +60,20 @@ class DFTU> : public OperatorLCAO ModuleBase::matrix& force, ModuleBase::matrix& stress); + // Getters for free functions in dftu_nao_fs_r/dftu_nao_for_r/dftu_nao_str_r + const UnitCell* get_ucell() const { return ucell; } + Plus_U* get_dftu() const { return dftu; } + const TwoCenterIntegrator* get_intor() const { return intor_; } + int get_nspin() const { return nspin; } + std::vector& get_adjs_all() { return adjs_all; } + + /// pot_onsite_{m, m'} = sum_{m,m'} (1/2*delta_{m, m'} - occ_{m, m'}) * U + /// EU = sum_{m,m'} 1/2 * U * occ_{m, m'} * occ_{m', m} + void cal_pot_onsite(const std::vector& occ, const int m_size, const double u_value, double* pot_onsite, double& eu); + + /// transfer pot_onsite format from pauli matrix to normal for non-collinear spin case + void transfer_pot_onsite(std::vector& pot_onsite_tmp, std::vector& pot_onsite); + private: const UnitCell* ucell = nullptr; @@ -98,12 +112,6 @@ class DFTU> : public OperatorLCAO const double* data_pointer, std::vector& occupations); - /// transfer pot_onsite format from pauli matrix to normal for non-collinear spin case - void transfer_pot_onsite(std::vector& pot_onsite_tmp, std::vector& pot_onsite); - /// pot_onsite_{m, m'} = sum_{m,m'} (1/2*delta_{m, m'} - occ_{m, m'}) * U - /// EU = sum_{m,m'} 1/2 * U * occ_{m, m'} * occ_{m', m} - void cal_pot_onsite(const std::vector& occ, const int m_size, const double u_value, double* pot_onsite, double& eu); - /** * @brief calculate the HR local matrix of atom pair */ diff --git a/source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp b/source/source_lcao/module_dftu/dftu_nao_op_legacy.cpp similarity index 98% rename from source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp rename to source/source_lcao/module_dftu/dftu_nao_op_legacy.cpp index 392c08cccd4..a1851efc3c4 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op_legacy.cpp @@ -1,4 +1,4 @@ -#include "dftu_lcao_op_legacy.h" +#include "dftu_nao_op_legacy.h" #include "dftu_hamilt.h" #include "source_base/timer.h" #include "source_base/tool_title.h" diff --git a/source/source_lcao/module_dftu/dftu_lcao_op_legacy.h b/source/source_lcao/module_dftu/dftu_nao_op_legacy.h similarity index 94% rename from source/source_lcao/module_dftu/dftu_lcao_op_legacy.h rename to source/source_lcao/module_dftu/dftu_nao_op_legacy.h index 3e78207095b..e6e9168fa66 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_op_legacy.h +++ b/source/source_lcao/module_dftu/dftu_nao_op_legacy.h @@ -2,7 +2,7 @@ #define OPDFTULCAO_H #include "source_lcao/module_operator_lcao/operator_lcao.h" -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 namespace hamilt { diff --git a/source/source_lcao/module_dftu/dftu_lcao_pots.cpp b/source/source_lcao/module_dftu/dftu_nao_pots.cpp similarity index 84% rename from source/source_lcao/module_dftu/dftu_lcao_pots.cpp rename to source/source_lcao/module_dftu/dftu_nao_pots.cpp index c244508cfc2..1f7486593e7 100644 --- a/source/source_lcao/module_dftu/dftu_lcao_pots.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_pots.cpp @@ -1,5 +1,10 @@ -#include "dftu_lcao.h" -#include "dftu_lcao_pots.h" +#include "dftu_nao.h" +#include "dftu_nao_pots.h" + +#include "source_base/global_function.h" +#include "source_base/tool_title.h" +#include "source_basis/module_ao/parallel_orbitals.h" +#include "source_cell/unitcell.h" #ifdef __LCAO void DFTU_LCAO::pot_onsite_complex(const Plus_U& dftu, @@ -13,7 +18,7 @@ void DFTU_LCAO::pot_onsite_complex(const Plus_U& dftu, ModuleBase::TITLE("DFTU_LCAO", "pot_onsite_complex"); ModuleBase::GlobalFunc::ZEROS(pot_onsite, pv->nloc); - const auto& iatlnmipol2iwt = dftu.get_iatlnmipol2iwt(); + const auto& iatlnmipol2iwt = dftu.occmat().iatlnmipol2iwt(); for (int it = 0; it < ucell.ntype; ++it) { @@ -86,7 +91,7 @@ void DFTU_LCAO::pot_onsite_real(const Plus_U& dftu, ModuleBase::TITLE("DFTU_LCAO", "pot_onsite_real"); ModuleBase::GlobalFunc::ZEROS(pot_onsite, pv->nloc); - const auto& iatlnmipol2iwt = dftu.get_iatlnmipol2iwt(); + const auto& iatlnmipol2iwt = dftu.occmat().iatlnmipol2iwt(); for (int it = 0; it < ucell.ntype; ++it) { @@ -178,13 +183,13 @@ double DFTU_LCAO::get_onsite_pot(const Plus_U& dftu, { if (m0 == m1) { - pot_onsite = (dftu.get_U_Yukawa(T, L, N) - dftu.get_J_Yukawa(T, L, N)) - * (0.5 - dftu.get_occ_mat(iat, L, N, spin, m0, m1)); + pot_onsite = (dftu.yukawa().get_U(T, L, N) - dftu.yukawa().get_J(T, L, N)) + * (0.5 - dftu.occmat().get(iat, L, N, spin, m0, m1)); } else { - pot_onsite = -(dftu.get_U_Yukawa(T, L, N) - dftu.get_J_Yukawa(T, L, N)) - * dftu.get_occ_mat(iat, L, N, spin, m0, m1); + pot_onsite = -(dftu.yukawa().get_U(T, L, N) - dftu.yukawa().get_J(T, L, N)) + * dftu.occmat().get(iat, L, N, spin, m0, m1); } } else @@ -192,12 +197,12 @@ double DFTU_LCAO::get_onsite_pot(const Plus_U& dftu, if (m0 == m1) { pot_onsite = dftu.get_u_current(T) - * (0.5 - dftu.get_occ_mat(iat, L, N, spin, m0, m1)); + * (0.5 - dftu.occmat().get(iat, L, N, spin, m0, m1)); } else { pot_onsite = -dftu.get_u_current(T) - * dftu.get_occ_mat(iat, L, N, spin, m0, m1); + * dftu.occmat().get(iat, L, N, spin, m0, m1); } } } @@ -207,13 +212,13 @@ double DFTU_LCAO::get_onsite_pot(const Plus_U& dftu, { if (m0 == m1) { - pot_onsite = (dftu.get_U_Yukawa(T, L, N) - dftu.get_J_Yukawa(T, L, N)) - * (0.5 - dftu.get_occ_mat_save(iat, L, N, spin, m0, m1)); + pot_onsite = (dftu.yukawa().get_U(T, L, N) - dftu.yukawa().get_J(T, L, N)) + * (0.5 - dftu.occmat().get_save(iat, L, N, spin, m0, m1)); } else { - pot_onsite = -(dftu.get_U_Yukawa(T, L, N) - dftu.get_J_Yukawa(T, L, N)) - * dftu.get_occ_mat_save(iat, L, N, spin, m0, m1); + pot_onsite = -(dftu.yukawa().get_U(T, L, N) - dftu.yukawa().get_J(T, L, N)) + * dftu.occmat().get_save(iat, L, N, spin, m0, m1); } } else @@ -221,12 +226,12 @@ double DFTU_LCAO::get_onsite_pot(const Plus_U& dftu, if (m0 == m1) { pot_onsite = dftu.get_u_current(T) - * (0.5 - dftu.get_occ_mat_save(iat, L, N, spin, m0, m1)); + * (0.5 - dftu.occmat().get_save(iat, L, N, spin, m0, m1)); } else { pot_onsite = -dftu.get_u_current(T) - * dftu.get_occ_mat_save(iat, L, N, spin, m0, m1); + * dftu.occmat().get_save(iat, L, N, spin, m0, m1); } } } diff --git a/source/source_lcao/module_dftu/dftu_lcao_pots.h b/source/source_lcao/module_dftu/dftu_nao_pots.h similarity index 100% rename from source/source_lcao/module_dftu/dftu_lcao_pots.h rename to source/source_lcao/module_dftu/dftu_nao_pots.h diff --git a/source/source_lcao/module_dftu/dftu_nao_str_r.cpp b/source/source_lcao/module_dftu/dftu_nao_str_r.cpp new file mode 100644 index 00000000000..5b27943490a --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_str_r.cpp @@ -0,0 +1,135 @@ +/// @file dftu_nao_str_r.cpp +/// @brief DFT+U stress calculation in real space (r-space) - implementation +/// +/// See dftu_nao_str_r.h for the mathematical formula and detailed documentation. + +#include "dftu_nao_str_r.h" +#include "dftu_nao_op.h" +#include "source_base/timer.h" + +namespace hamilt +{ + +template +void cal_str_IJR_nao_r(const DFTU>* dftu_op, + const int& iat1, + const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress) +{ + // npol is the number of polarizations, + // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), + // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) + const int npol = dftu_op->get_ucell()->get_npol(); + + // --------------------------------------------- + // calculate the Nonlocal matrix for each pair of orbitals + // --------------------------------------------- + auto row_indexes = pv->get_indexes_row(iat1); + auto col_indexes = pv->get_indexes_col(iat2); + const int m_size = int(sqrt(pot_onsite_in.size() / nspin)); + const int m_size2 = m_size * m_size; + + // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 + std::vector step_trace(npol * npol, 0); + + if (npol == 2) + { + step_trace[1] = 1; + step_trace[2] = col_indexes.size(); + step_trace[3] = col_indexes.size() + 1; + } + + // calculate the local matrix + for (int is = 0; is < nspin; is++) + { + const int is0 = nspin == 2 ? is : 0; + const int step_is = nspin == 4 ? is : 0; + const double* dm_pointer = dmR_pointer[is0]->get_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; +#ifdef __DEBUG + assert(nlm1.size() == nlm2.size()); +#endif + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + double tmp = pot_onsite_in[m1 * m_size + m2 + is * m_size2] + * dm_pointer[step_trace[step_is]]; + // Voigt notation: stress[0]=xx, stress[1]=xy, stress[2]=xz, + // stress[3]=yy, stress[4]=yz, stress[5]=zz + stress[0] += tmp * (nlm1[m1 + m_size] * dis1.x * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size] * dis2.x); + stress[1] += tmp * (nlm1[m1 + m_size] * dis1.y * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size] * dis2.y); + stress[2] += tmp * (nlm1[m1 + m_size] * dis1.z * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size] * dis2.z); + + stress[3] += tmp * (nlm1[m1 + m_size * 2] * dis1.y * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size * 2] * dis2.y); + stress[4] += tmp * (nlm1[m1 + m_size * 2] * dis1.z * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size * 2] * dis2.z); + stress[5] += tmp * (nlm1[m1 + m_size * 3] * dis1.z * nlm2[m2] + + nlm1[m1] * nlm2[m2 + m_size * 3] * dis2.z); + } + } + dm_pointer += npol; + } + dm_pointer += (npol - 1) * col_indexes.size(); + } + } +} + +// explicit template instantiation +template void cal_str_IJR_nao_r( + const DFTU>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +template void cal_str_IJR_nao_r, double>( + const DFTU, double>>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +template void cal_str_IJR_nao_r, std::complex>( + const DFTU, std::complex>>* dftu_op, + const int& iat1, const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +} // namespace hamilt diff --git a/source/source_lcao/module_dftu/dftu_nao_str_r.h b/source/source_lcao/module_dftu/dftu_nao_str_r.h new file mode 100644 index 00000000000..2c230534065 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_nao_str_r.h @@ -0,0 +1,87 @@ +/// @file dftu_nao_str_r.h +/// @brief DFT+U stress calculation in real space (r-space) +/// +/// This file provides the real-space implementation of DFT+U stress contribution +/// from a single atom pair (I,J,R). It is independent of k-point sampling because +/// the real-space density matrix (DMR) already contains the Brillouin-zone integration. +/// +/// Naming convention: _r suffix denotes real-space implementation, +/// corresponding to _k suffix for k-space (legacy) implementation. +/// +/// The stress formula for atom pair (I,J,R) is: +/// +/// sigma_{alpha,beta} += -(1/Omega) * sum_{m,m'} V_U_{mm'}(I) * DMR_{mu,nu}(J1,J2,R) * [ +/// d/d tau_{J1,alpha} * * R_{J1,beta} +/// + * d/d tau_{J2,alpha} * R_{J2,beta} +/// ] +/// +/// where R_{J1} and R_{J2} are the position vectors of atoms J1 and J2 relative to +/// the on-site atom I, and Omega is the unit cell volume. + +#ifndef DFTU_NAO_STR_R_H +#define DFTU_NAO_STR_R_H + +#include "source_basis/module_ao/parallel_orbitals.h" +#include "source_base/vector3.h" +#include "source_hamilt/module_hcontainer/hcontainer.h" + +#include +#include + +namespace hamilt +{ + +// Forward declarations to avoid circular dependency with dftu_lcao_op.h +template +class OperatorLCAO; + +template +class DFTU; + +/** + * @brief Compute DFT+U stress contribution from a single atom pair (I,J,R) in real space + * + * The stress tensor component (alpha,beta) receives contributions from both the derivative + * of with respect to the strain and the explicit R_beta factor: + * + * sigma_{alpha,beta} += -(1/Omega) * sum_{m,m'} V_U_{mm'}(I) * DMR_{mu,nu}(J1,J2,R) * [ + * d/d tau_{J1,alpha} * * R_{J1,beta} + * + * d/d tau_{J2,alpha} * R_{J2,beta} + * ] + * + * The strain derivative of the two-center integral is approximated by: + * d/d epsilon_{alpha,beta} ~ d/d tau_alpha * R_beta + * + * @note The stress is accumulated in Voigt notation (6 components: xx, xy, xz, yy, yz, zz) + * and symmetrized at the caller level. The final scaling by lat0/Omega is applied + * after MPI reduction. + * + * @param iat1 [in] global atom index of J1 + * @param iat2 [in] global atom index of J2 + * @param pv [in] Parallel_Orbitals for basis index mapping + * @param nlm1_all [in] pre-computed and derivatives for atom J1 + * @param nlm2_all [in] pre-computed and derivatives for atom J2 + * @param pot_onsite [in] flattened V_U matrix + * @param dmR_pointer [in] pointer to DMR matrix blocks for each spin + * @param nspin [in] number of spin channels + * @param dis1 [in] position vector of J1 relative to I + * @param dis2 [in] position vector of J2 relative to I + * @param stress [out] stress accumulator (6 components in Voigt notation) + */ +template +void cal_str_IJR_nao_r(const DFTU>* dftu_op, + const int& iat1, + const int& iat2, + const Parallel_Orbitals* pv, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& pot_onsite_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +} // namespace hamilt + +#endif // DFTU_NAO_STR_R_H diff --git a/source/source_lcao/module_dftu/dftu_yukawa.h b/source/source_lcao/module_dftu/dftu_yukawa.h deleted file mode 100644 index b249effa979..00000000000 --- a/source/source_lcao/module_dftu/dftu_yukawa.h +++ /dev/null @@ -1,84 +0,0 @@ -#ifndef DFTU_YUKAWA_H -#define DFTU_YUKAWA_H - -class Plus_U; -class UnitCell; - -#ifdef __LCAO -namespace DFTU_LCAO { - -/** - * @brief Spherical modified Bessel function of the first kind (Yukawa kernel). - * - * Evaluates i_k(r*lambda) for even orders k = 0, 2, 4, 6 used in the Slater - * integral construction for the Yukawa-screened DFT+U. The implementation uses - * piecewise small-x expansions and large-x closed forms, so the caller must - * pass an order handled by those branches; other orders return 0. - * - * @param k order (supported: 0, 2, 4, 6) - * @param r radial distance - * @param lambda Yukawa screening length - * @return function value - */ -double spherical_Bessel(const int k, const double r, const double lambda); - -/** - * @brief Spherical modified Hankel function of the second kind (Yukawa kernel). - * - * Evaluates k_k(r*lambda) for even orders k = 0, 2, 4, 6 used together with - * spherical_Bessel in the radial Slater integrals. Same piecewise scheme as - * spherical_Bessel; unsupported orders return 0. - * - * @param k order (supported: 0, 2, 4, 6) - * @param r radial distance - * @param lambda Yukawa screening length - * @return function value - */ -double spherical_Hankel(const int k, const double r, const double lambda); - -/** - * @brief Determine the Yukawa screening length lambda and store it on dftu. - * - * If dftu carries a positive yukawa_lambda configuration value, that value is - * used directly; otherwise lambda is estimated from the charge density - * (Thomas-Fermi-like) and rescaled by 1.6. The spin channel count is read - * from the global PARAM.inp.nspin rather than a Plus_U member. - * - * @param dftu Plus_U state (lambda is written back via set_lambda) - * @param rho charge density per spin channel - * @param nrxx number of real-space grid points - */ -void cal_yukawa_lambda(Plus_U& dftu, double** rho, const int& nrxx); - -/** - * @brief Compute the Slater integrals Fk for the correlated orbital of atom type T. - * - * Accumulates into dftu.Fk[T][L][chi][k] using the radial orbital grid and the - * Yukawa-screened spherical Bessel/Hankel kernels. Only acts when Yukawa - * screening is enabled on dftu. - * - * @param dftu Plus_U state (provides ptr_orb, lambda, Fk, use_yukawa) - * @param ucell unit cell - * @param L angular momentum - * @param T atom type - */ -void cal_slater_Fk(Plus_U& dftu, const UnitCell& ucell, const int L, const int T); - -/** - * @brief Compute Yukawa-screened Slater integrals and derive U/J for all atoms. - * - * Drives cal_yukawa_lambda then cal_slater_Fk over correlated orbitals and - * writes the resulting U_Yukawa/J_Yukawa/u_current back onto dftu. No-op when - * Yukawa screening is disabled on dftu. - * - * @param dftu Plus_U state (U_Yukawa, J_Yukawa, u_current, lambda written back) - * @param ucell unit cell - * @param rho charge density per spin channel - * @param nrxx number of real-space grid points - */ -void cal_slater_UJ(Plus_U& dftu, const UnitCell& ucell, double** rho, const int& nrxx); - -} // namespace DFTU_LCAO -#endif - -#endif diff --git a/source/source_lcao/module_dftu/test/CMakeLists.txt b/source/source_lcao/module_dftu/test/CMakeLists.txt index bad024242c6..94738d40d60 100644 --- a/source/source_lcao/module_dftu/test/CMakeLists.txt +++ b/source/source_lcao/module_dftu/test/CMakeLists.txt @@ -4,7 +4,7 @@ AddTest( TARGET dftu_pw_test LIBS base device parameter SOURCES dftu_pw_test.cpp - ../../../source_pw/module_pwdft/dftu_tools_pw.cpp + ../../../source_pw/module_pwdft/dftu_base_tools.cpp ) AddTest( @@ -23,8 +23,11 @@ if(ENABLE_LCAO AND ENABLE_MPI) AddTest( TARGET dftu_lcao_test LIBS parameter psi base device container - SOURCES dftu_lcao_test.cpp ../dftu_lcao_op.cpp ../dftu_fs.cpp + SOURCES dftu_lcao_test.cpp ../dftu_nao_op.cpp ../dftu_nao_fs_r.cpp ../dftu_nao_for_r.cpp ../dftu_nao_str_r.cpp ../../../source_pw/module_pwdft/dftu_base.cpp + ../../../source_pw/module_pwdft/dftu_base_io.cpp + ../../../source_pw/module_pwdft/yukawa_screening.cpp + ../../../source_estate/occ_matrix.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp diff --git a/source/source_lcao/module_dftu/test/dftu_lcao_test.cpp b/source/source_lcao/module_dftu/test/dftu_lcao_test.cpp index f89a4351d72..6c6db0902ce 100644 --- a/source/source_lcao/module_dftu/test/dftu_lcao_test.cpp +++ b/source/source_lcao/module_dftu/test/dftu_lcao_test.cpp @@ -5,8 +5,8 @@ #define private public #include "source_io/module_parameter/parameter.h" #undef private -#include "../dftu_lcao_op.h" -#include "source_lcao/module_dftu/dftu_lcao.h" +#include "../dftu_nao_op.h" +#include "source_lcao/module_dftu/dftu_nao.h" Plus_U::Plus_U(){}; Plus_U::~Plus_U(){}; @@ -85,16 +85,16 @@ class DFTUTest : public ::testing::Test tmp_DMR = DMR; // setting of DFTU - dftu.occ_mat.resize(test_size); + dftu.occmat().data().resize(test_size); for (int iat = 0; iat < test_size; iat++) { - dftu.occ_mat[iat].resize(3); + dftu.occmat().data()[iat].resize(3); for (int l = 0; l < 3; l++) { - dftu.occ_mat[iat][l].resize(1); - dftu.occ_mat[iat][l][0].resize(2); - dftu.occ_mat[iat][l][0][0].create(2 * l + 1, 2 * l + 1); - dftu.occ_mat[iat][l][0][1].create(2 * l + 1, 2 * l + 1); + dftu.occmat().data()[iat][l].resize(1); + dftu.occmat().data()[iat][l][0].resize(2); + dftu.occmat().data()[iat][l][0][0].create(2 * l + 1, 2 * l + 1); + dftu.occmat().data()[iat][l][0][1].create(2 * l + 1, 2 * l + 1); } } dftu.u_current = {U_test}; @@ -112,12 +112,12 @@ class DFTUTest : public ::testing::Test } // Helper for TEST_F bodies: gtest-derived classes do not inherit - // the friend declaration, so direct dftu.occ_mat[...] access from - // TestBody would fail to compile. This wrapper runs inside + // the friend declaration, so direct dftu.occmat().data()[...] access + // from TestBody would fail to compile. This wrapper runs inside // DFTUTest, which is a friend of Plus_U_Base. double occ_mat_c(int iat, int spin, int icc) const { - return dftu.occ_mat[iat][2][0][spin].c[icc]; + return dftu.occmat().data()[iat][2][0][spin].c[icc]; } #ifdef __MPI diff --git a/source/source_lcao/module_dftu/test/dftu_pw_test.cpp b/source/source_lcao/module_dftu/test/dftu_pw_test.cpp index bea26626428..fcca0a5a75c 100644 --- a/source/source_lcao/module_dftu/test/dftu_pw_test.cpp +++ b/source/source_lcao/module_dftu/test/dftu_pw_test.cpp @@ -5,7 +5,7 @@ #include "source_io/module_parameter/parameter.h" #undef private #include "source_base/matrix.h" -#include "source_pw/module_pwdft/dftu_tools_pw.h" +#include "source_pw/module_pwdft/dftu_base_tools.h" /*********************************************************************** * Unit tests for DFT+U PW nspin=1/2/4 support (PR-2) @@ -83,7 +83,7 @@ TEST_F(DftuPwTest, PotOnsitePotNspin1_DiagonalLocale) occ_mat_c[m * m_size + m] = 0.3; // diagonal std::vector> pot_onsite(size, {0.0, 0.0}); - dftu_pw::compute_pot_onsite_scalar(pot_onsite.data(), occ_mat_c.data(), U_val, 0.5, 1.0, m_size); + DFTU_BASE::compute_pot_onsite_scalar(pot_onsite.data(), occ_mat_c.data(), U_val, 0.5, 1.0, m_size); // diagonal: U*(0.5 - 0.3) = 4.0*0.2 = 0.8 for (int m = 0; m < m_size; m++) @@ -107,8 +107,8 @@ TEST_F(DftuPwTest, PotOnsitePotNspin2_TwoSpinChannels) std::vector> pot_onsite_up(size, {0.0, 0.0}); std::vector> pot_onsite_dn(size, {0.0, 0.0}); - dftu_pw::compute_pot_onsite_scalar(pot_onsite_up.data(), occ_mat_up.data(), U_val, 0.5, 0.5, m_size); - dftu_pw::compute_pot_onsite_scalar(pot_onsite_dn.data(), occ_mat_dn.data(), U_val, 0.5, 0.5, m_size); + DFTU_BASE::compute_pot_onsite_scalar(pot_onsite_up.data(), occ_mat_up.data(), U_val, 0.5, 0.5, m_size); + DFTU_BASE::compute_pot_onsite_scalar(pot_onsite_dn.data(), occ_mat_dn.data(), U_val, 0.5, 0.5, m_size); // pot_onsite_up[0,0] = U*(0.5 - 0.4) = 0.5 EXPECT_DOUBLE_EQ(pot_onsite_up[0].real(), 0.5); @@ -133,7 +133,7 @@ TEST_F(DftuPwTest, PotOnsitePotNspin4_PauliTransform) pot_onsite[2] = {0.3, 0.0}; // sigma_y pot_onsite[3] = {0.2, 0.0}; // sigma_z - dftu_pw::pauli_to_spin_basis(pot_onsite, m_size); + DFTU_BASE::pauli_to_spin_basis(pot_onsite, m_size); EXPECT_DOUBLE_EQ(pot_onsite[0].real(), 0.6); // 0.5*(1.0+0.2) EXPECT_DOUBLE_EQ(pot_onsite[0].imag(), 0.0); @@ -164,7 +164,7 @@ TEST_F(DftuPwTest, EnergyNspin12_DiagonalLocale) // nspin=1: E = U * 1.0 * (0.5^2 + 0.3^2 + 0.2^2) = 4 * 0.38 = 1.52 std::vector> pot_onsite_nspin1(size, {0.0, 0.0}); - double energy_u = dftu_pw::compute_pot_onsite_scalar( + double energy_u = DFTU_BASE::compute_pot_onsite_scalar( pot_onsite_nspin1.data(), occ_mat_c.data(), U_val, 0.5, 1.0, m_size); EXPECT_DOUBLE_EQ(energy_u, 1.52); @@ -174,9 +174,9 @@ TEST_F(DftuPwTest, EnergyNspin12_DiagonalLocale) std::vector> pot_onsite_up(size, {0.0, 0.0}); std::vector> pot_onsite_dn(size, {0.0, 0.0}); energy_u = 0.0; - energy_u += dftu_pw::compute_pot_onsite_scalar( + energy_u += DFTU_BASE::compute_pot_onsite_scalar( pot_onsite_up.data(), occ_mat_up.data(), U_val, 0.5, 0.5, m_size); - energy_u += dftu_pw::compute_pot_onsite_scalar( + energy_u += DFTU_BASE::compute_pot_onsite_scalar( pot_onsite_dn.data(), occ_mat_dn.data(), U_val, 0.5, 0.5, m_size); // E = U*0.5*(0.4^2 + 0.6^2) = 4*0.5*(0.16+0.36) = 1.04 EXPECT_DOUBLE_EQ(energy_u, 1.04); @@ -200,7 +200,7 @@ TEST_F(DftuPwTest, EnergyNspin4_WithOffDiagonal) occ_mat_c[size + 2] = 0.0; occ_mat_c[size + 3] = 0.2; std::vector> pot_onsite(size * 4, {0.0, 0.0}); - double energy_u = dftu_pw::compute_pot_onsite_spinor( + double energy_u = DFTU_BASE::compute_pot_onsite_spinor( pot_onsite.data(), occ_mat_c.data(), U_val, 1.0, weight_eu, m_size); // is=0: 2*0.25*(0.5*0.5 + 0.1*0.1 + 0.1*0.1 + 0.5*0.5) = 0.26 @@ -227,7 +227,7 @@ TEST_F(DftuPwTest, LocaleAccumNspin12) wg(0, 1) = 0.5; std::vector occ_mat_c(m_size * m_size, 0.0); - dftu_pw::accumulate_occ_scalar( + DFTU_BASE::accumulate_occ_scalar( occ_mat_c.data(), becp.data(), nbands, nkb, begin_ih, m_begin, m_size, wg, ik); @@ -261,7 +261,7 @@ TEST_F(DftuPwTest, LocaleAccumNspin4_PauliComponents) ModuleBase::matrix wg(1, nbands); wg(0, 0) = 1.0; - dftu_pw::accumulate_occ_spinor( + DFTU_BASE::accumulate_occ_spinor( occ_mat_c.data(), becp.data(), nbands, npol, nkb, 0, 0, m_size, wg, ik); diff --git a/source/source_lcao/module_lr/potentials/xc_kernel.h b/source/source_lcao/module_lr/potentials/xc_kernel.h index a52b0c3c03e..06cb78c9e3d 100644 --- a/source/source_lcao/module_lr/potentials/xc_kernel.h +++ b/source/source_lcao/module_lr/potentials/xc_kernel.h @@ -3,8 +3,8 @@ #include "source_cell/unitcell.h" #include "source_base/parallel_grid.h" #include "source_estate/module_charge/charge.h" -#define CREF(x) const std::vector& x = x##_; -#define CREF3(x) const std::vector>& x = x##_; +#define CREF(x) const std::vector& x = x##_ +#define CREF3(x) const std::vector>& x = x##_ namespace LR { /// @brief Calculate the exchange-correlation (XC) kernel ($f_{xc}=\delta^2E_xc/\delta\rho^2$) and store its components. @@ -21,7 +21,7 @@ namespace LR const std::string& kernel_name, const std::vector& lr_init_xc_kernel, const bool openshell = false); - ~KernelXC() {}; + ~KernelXC() {} // const references CREF(vrho);CREF(vsigma); CREF(v2rho2); CREF(v2rhosigma); CREF(v2sigma2); diff --git a/source/source_lcao/setup_dftu_lcao.cpp b/source/source_lcao/setup_dftu_lcao.cpp index 771797208a6..4d9d62195e7 100644 --- a/source/source_lcao/setup_dftu_lcao.cpp +++ b/source/source_lcao/setup_dftu_lcao.cpp @@ -1,9 +1,9 @@ #include "setup_dftu_lcao.h" -#include "source_lcao/module_dftu/dftu_lcao.h" -#include "source_lcao/module_dftu/dftu_lcao_occ.h" -#include "source_lcao/module_dftu/dftu_lcao_energy.h" -#include "source_lcao/module_dftu/dftu_yukawa.h" -#include "source_pw/module_pwdft/dftu_output.h" // mohan add 2025-11-08 +#include "source_lcao/module_dftu/dftu_nao.h" +#include "source_lcao/module_dftu/dftu_nao_occ.h" +#include "source_lcao/module_dftu/dftu_nao_energy.h" +#include "source_pw/module_pwdft/dftu_base_io.h" // mohan add 2025-11-08 +#include "source_io/module_parameter/parameter.h" #include "source_estate/module_dm/density_matrix.h" #include "source_lcao/hamilt_lcao.h" @@ -34,7 +34,18 @@ void init_dftu_lcao(const int istep, } /// Calculate U and J if Yukawa potential is used - DFTU_LCAO::cal_slater_UJ(*dftu_ptr, ucell, rho, nrxx); + if (dftu_ptr->use_yukawa()) + { + dftu_ptr->yukawa().cal_slater_UJ(ucell, rho, nrxx, PARAM.inp.nspin, dftu_ptr->get_ptr_orb()); + // update current U with calculated U-J from Slater integrals + for (int T = 0; T < ucell.ntype; T++) + { + if (dftu_ptr->has_correlated_orbital(T)) + { + dftu_ptr->set_u_current(T, dftu_ptr->yukawa().get_Ueff(T)); + } + } + } } template @@ -73,7 +84,7 @@ void finish_dftu_lcao(const int iter, } DFTU_LCAO::cal_energy_correction(*dftu_ptr, ucell); } - dftu_io::output(*dftu_ptr, ucell, out_chg, global_out_dir, nspin, npol); + DFTU_BASE::output(*dftu_ptr, ucell, out_chg, global_out_dir, nspin, npol); /// use the converged occupation matrix for next MD/Relax SCF calculation if (conv_esolver) diff --git a/source/source_lcao/spar_u.cpp b/source/source_lcao/spar_u.cpp index 0f8e4b2aac2..4a04fe25c31 100644 --- a/source/source_lcao/spar_u.cpp +++ b/source/source_lcao/spar_u.cpp @@ -1,7 +1,11 @@ #include "spar_u.h" +#include "source_base/global_function.h" #include "source_base/parallel_reduce.h" +#include "source_base/tool_title.h" +#include "source_basis/module_ao/parallel_orbitals.h" #include "source_io/module_parameter/parameter.h" #include "source_base/timer.h" +#include "source_lcao/module_dftu/dftu_hamilt.h" void sparse_format::cal_HR_dftu( Plus_U &dftu, // mohan add 2025-11-07 @@ -74,7 +78,7 @@ void sparse_format::cal_HR_dftu( } } - dftu.cal_eff_pot_mat_R_double(ucell, &pv, current_spin, SR_tmp, HR_tmp, PARAM.globalv.npol); + DFTU_LCAO::pot_uterm_HR_real(dftu, ucell, &pv, current_spin, SR_tmp, HR_tmp, PARAM.globalv.npol); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { @@ -194,7 +198,7 @@ void sparse_format::cal_HR_dftu_soc( } } - dftu.cal_eff_pot_mat_R_complex_double(ucell, &pv, current_spin, SR_soc_tmp, HR_soc_tmp, PARAM.globalv.npol); + DFTU_LCAO::pot_uterm_HR_complex(dftu, ucell, &pv, current_spin, SR_soc_tmp, HR_soc_tmp, PARAM.globalv.npol); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { diff --git a/source/source_lcao/spar_u.h b/source/source_lcao/spar_u.h index 5ad5c661ceb..b111e117cde 100644 --- a/source/source_lcao/spar_u.h +++ b/source/source_lcao/spar_u.h @@ -2,7 +2,7 @@ #define SPARSE_FORMAT_U_H #include "source_lcao/module_ri/abfs_vector3_order.h" -#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_nao.h" // mohan add 20251107 namespace sparse_format { diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index 208e1cb2770..691ec986f98 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -13,9 +13,10 @@ list(APPEND objects op_pw_exx_ace.cpp op_pw_exx_pot.cpp dftu_base.cpp - dftu_output.cpp - dftu_tools_pw.cpp - dftu_cal_occ_pw.cpp + dftu_base_io.cpp + dftu_base_occ.cpp + dftu_base_tools.cpp + yukawa_screening.cpp setup_pot.cpp setup_pwrho.cpp setup_pwwfc.cpp diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index d54b8dc30e6..a2005daf21d 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -1,5 +1,7 @@ #include "source_pw/module_pwdft/dftu_base.h" +#include "source_cell/unitcell.h" +#include "source_pw/module_pwdft/dftu_base_io.h" #include "source_base/global_function.h" #include "source_base/memory_recorder.h" #include "source_base/parallel_global.h" @@ -12,7 +14,7 @@ #include // local inline helpers for eigenvalue calculation (JacobiRotate, CalculateEigenvalues) -// have been migrated to dftu_output.cpp, where they are used by dftu_io::write_occup_m. +// have been migrated to dftu_base_io.cpp, where they are used by DFTU_BASE::write_occup_m. // mohan refactored 2025-11-08 // All members are now non-static; default values are in the header. @@ -32,6 +34,7 @@ void Plus_U_Base::init_base(UnitCell& cell, const int nspin, const std::vector& orbital_corr, const bool yukawa_potential, + const double yukawa_lambda, const std::string& global_readin_dir, const std::string& global_out_dir, const std::string& init_chg, @@ -51,7 +54,6 @@ void Plus_U_Base::init_base(UnitCell& cell, this->nspin = nspin; this->orbital_corr = orbital_corr; - this->use_yukawa_ = yukawa_potential; this->uramping = uramping; this->occ_mat_ctrl = occ_mat_ctrl; this->mixing_dftu = mixing_dftu; @@ -68,13 +70,11 @@ void Plus_U_Base::init_base(UnitCell& cell, this->energy_u = 0.0; - this->occ_mat.resize(cell.nat); - this->occ_mat_save.resize(cell.nat); + this->occmat_.init(cell, orbital_corr, nspin, npol); + this->pot_uterm_pw_index.resize(cell.nat); int pot_index = 0; - this->iatlnmipol2iwt.resize(cell.nat); - int num_locale = 0; for (int it = 0; it < cell.ntype; ++it) { @@ -82,11 +82,6 @@ void Plus_U_Base::init_base(UnitCell& cell, { const int iat = cell.itia2iat(it, ia); - occ_mat[iat].resize(cell.atoms[it].nwl + 1); - occ_mat_save[iat].resize(cell.atoms[it].nwl + 1); - - this->iatlnmipol2iwt[iat].resize(cell.atoms[it].nwl + 1); - if(!has_correlated_orbital(it)) { continue; @@ -110,62 +105,18 @@ void Plus_U_Base::init_base(UnitCell& cell, { const int N = cell.atoms[it].l_nchi[l]; - occ_mat[iat][l].resize(N); - occ_mat_save[iat][l].resize(N); - for (int n = 0; n < N; n++) { if (nspin == 1 || nspin == 2) { - occ_mat[iat][l][n].resize(2); - occ_mat_save[iat][l][n].resize(2); - - occ_mat[iat][l][n][0].create(2 * l + 1, 2 * l + 1); - occ_mat[iat][l][n][1].create(2 * l + 1, 2 * l + 1); - - occ_mat_save[iat][l][n][0].create(2 * l + 1, 2 * l + 1); - occ_mat_save[iat][l][n][1].create(2 * l + 1, 2 * l + 1); num_locale += (2 * l + 1) * (2 * l + 1) * 2; } else if (nspin == 4) { - occ_mat[iat][l][n].resize(1); - occ_mat_save[iat][l][n].resize(1); - - occ_mat[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol); - occ_mat_save[iat][l][n][0].create((2 * l + 1) * npol, (2 * l + 1) * npol); num_locale += (2 * l + 1) * (2 * l + 1) * npol * npol; } } } - - this->iatlnmipol2iwt[iat].resize(cell.atoms[it].nwl + 1); - for (int L = 0; L <= cell.atoms[it].nwl; L++) - { - this->iatlnmipol2iwt[iat][L].resize(cell.atoms[it].l_nchi[L]); - - for (int n = 0; n < cell.atoms[it].l_nchi[L]; n++) - { - this->iatlnmipol2iwt[iat][L][n].resize(2 * L + 1); - - for (int m = 0; m < 2 * L + 1; m++) - { - this->iatlnmipol2iwt[iat][L][n][m].resize(npol); - } - } - } - - for (int iw = 0; iw < cell.atoms[it].nw * npol; iw++) - { - int iw0 = iw / npol; - int ipol = iw % npol; - int iwt = cell.itiaiw2iwt(it, ia, iw); - int l = cell.atoms[it].iw2l[iw0]; - int n = cell.atoms[it].iw2n[iw0]; - int m = cell.atoms[it].iw2m[iw0]; - - this->iatlnmipol2iwt[iat][l][n][m][ipol] = iwt; - } } } @@ -175,48 +126,26 @@ void Plus_U_Base::init_base(UnitCell& cell, this->uom_array.resize(pot_index, 0.0); this->uom_save.resize(pot_index, 0.0); - if (use_yukawa_) + if (yukawa_potential) { - this->Fk.resize(cell.ntype); - - this->U_Yukawa.resize(cell.ntype); - this->J_Yukawa.resize(cell.ntype); - - for (int it = 0; it < cell.ntype; it++) - { - const int NL = cell.atoms[it].nwl + 1; - - this->Fk[it].resize(NL); - this->U_Yukawa[it].resize(NL); - this->J_Yukawa[it].resize(NL); - - for (int l = 0; l < NL; l++) - { - int N = cell.atoms[it].l_nchi[l]; - - this->Fk[it][l].resize(N); - for (int n = 0; n < N; n++) - { - this->Fk[it][l][n].resize(l + 1, 0.0); - } - - this->U_Yukawa[it][l].resize(N, 0.0); - this->J_Yukawa[it][l].resize(N, 0.0); - } - } + this->yukawa_.reset(new YukawaScreening()); + this->yukawa_->init(cell, orbital_corr, yukawa_lambda); } if (occ_mat_ctrl != 0) { std::stringstream sst; sst << global_readin_dir << "dm_onsite_ini.txt"; - this->read_occup_m(cell, sst.str(), init_chg, nspin, npol); + DFTU_BASE::read_occup_m(cell, this->occmat_, this->orbital_corr, this->occ_mat_ctrl, + sst.str(), init_chg, nspin, npol); #ifdef __MPI - this->local_occup_bcast(cell, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->orbital_corr, nspin, npol); #endif mark_occ_mat_initialized(); - this->copy_occ_mat(cell); + this->occmat_.copy_to_save(cell, this->orbital_corr); + this->occmat_.write_save_to_flat(cell, this->orbital_corr, + this->pot_uterm_pw_index, this->uom_save); } else { @@ -224,15 +153,16 @@ void Plus_U_Base::init_base(UnitCell& cell, { std::stringstream sst; sst << global_readin_dir << "dm_onsite.txt"; - this->read_occup_m(cell, sst.str(), init_chg, nspin, npol); + DFTU_BASE::read_occup_m(cell, this->occmat_, this->orbital_corr, this->occ_mat_ctrl, + sst.str(), init_chg, nspin, npol); #ifdef __MPI - this->local_occup_bcast(cell, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->orbital_corr, nspin, npol); #endif mark_occ_mat_initialized(); } else { - this->zero_occ_mat(cell); + this->occmat_.zero(cell, this->orbital_corr); } } @@ -244,7 +174,7 @@ void Plus_U_Base::init_base(UnitCell& cell, void Plus_U_Base::uramping_update() { // Yukawa calculates U directly every iteration, no need for ramping - if (use_yukawa_) { + if (use_yukawa()) { return; } // if uramping < 0.1, use the original U @@ -269,7 +199,7 @@ void Plus_U_Base::uramping_update() bool Plus_U_Base::u_converged() { // Yukawa calculates U directly every iteration, always considered converged - if (use_yukawa_) { + if (use_yukawa()) { return true; } for (int i = 0; i < static_cast(this->u_target.size()); i++) @@ -283,453 +213,6 @@ bool Plus_U_Base::u_converged() } -// copy_locale — save current occ_mat to occ_mat_save and uom_save -void Plus_U_Base::copy_occ_mat(const UnitCell& ucell) -{ - ModuleBase::TITLE("Plus_U_Base", "copy_occ_mat"); - ModuleBase::timer::start("Plus_U_Base", "copy_occ_mat"); - - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = get_orbital_corr(T); - if (target_l == -1) - continue; - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - - if (this->nspin == 4) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - if(this->uom_save.size() != 0) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for(int mm=0; mmuom_save[pot_uterm_pw_index[iat]+mm] = occ_mat[iat][target_l][0][0].c[mm]; - } - } - } - else if (this->nspin == 1 || this->nspin == 2) - { - occ_mat_save[iat][target_l][0][0] = occ_mat[iat][target_l][0][0]; - occ_mat_save[iat][target_l][0][1] = occ_mat[iat][target_l][0][1]; - if(this->uom_save.size() != 0) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - const int half_size = this->uom_save.size() / 2; - for(int mm=0; mmuom_save[pot_uterm_pw_index[iat]+mm] = occ_mat[iat][target_l][0][0].c[mm]; - this->uom_save[half_size + pot_uterm_pw_index[iat]+mm] = occ_mat[iat][target_l][0][1].c[mm]; - } - } - } - } - } - ModuleBase::timer::end("Plus_U_Base", "copy_occ_mat"); -} - - -void Plus_U_Base::zero_occ_mat(const UnitCell& ucell) -{ - ModuleBase::TITLE("Plus_U_Base", "zero_occ_mat"); - ModuleBase::timer::start("Plus_U_Base", "zero_occ_mat"); - - for (int T = 0; T < ucell.ntype; T++) - { - if (!has_correlated_orbital(T)) - { - continue; - } - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - - for (int l = 0; l < ucell.atoms[T].nwl + 1; l++) - { - const int N = ucell.atoms[T].l_nchi[l]; - - for (int n = 0; n < N; n++) - { - if (this->nspin == 4) - { - occ_mat[iat][l][n][0].zero_out(); - } - else if (this->nspin == 1 || this->nspin == 2) - { - occ_mat[iat][l][n][0].zero_out(); - occ_mat[iat][l][n][1].zero_out(); - } - } - } - } - } - ModuleBase::timer::end("Plus_U_Base", "zero_occ_mat"); -} - - -void Plus_U_Base::mix_occ_mat(const UnitCell& ucell, - const double& mixing_beta) -{ - ModuleBase::TITLE("Plus_U_Base", "mix_occ_mat"); - ModuleBase::timer::start("Plus_U_Base", "mix_occ_mat"); - - double beta = mixing_beta; - - for (int T = 0; T < ucell.ntype; T++) - { - int target_l = get_orbital_corr(T); - if (target_l == -1) - continue; - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - - if (this->nspin == 4) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - } - if (this->uom_save.size() != 0) - { - for (int mm = 0; mm < size; mm++) - { - this->uom_save[pot_uterm_pw_index[iat] + mm] = occ_mat[iat][target_l][0][0].c[mm]; - } - } - } - else if (this->nspin == 1 || this->nspin == 2) - { - const int size = occ_mat[iat][target_l][0][0].nr * occ_mat[iat][target_l][0][0].nc; - const int half_size = this->uom_save.size() / 2; - for (int mm = 0; mm < size; mm++) - { - occ_mat[iat][target_l][0][0].c[mm] = occ_mat[iat][target_l][0][0].c[mm] * beta + occ_mat_save[iat][target_l][0][0].c[mm] * (1.0 - beta); - occ_mat[iat][target_l][0][1].c[mm] = occ_mat[iat][target_l][0][1].c[mm] * beta + occ_mat_save[iat][target_l][0][1].c[mm] * (1.0 - beta); - } - if (this->uom_save.size() != 0) - { - for (int mm = 0; mm < size; mm++) - { - this->uom_save[pot_uterm_pw_index[iat] + mm] = occ_mat[iat][target_l][0][0].c[mm]; - this->uom_save[half_size + pot_uterm_pw_index[iat] + mm] = occ_mat[iat][target_l][0][1].c[mm]; - } - } - } - } - } - ModuleBase::timer::end("Plus_U_Base", "mix_occ_mat"); -} - - -void Plus_U_Base::set_occ_mat(const UnitCell& ucell) -{ - ModuleBase::TITLE("Plus_U_Base", "set_occ_mat"); - ModuleBase::timer::start("Plus_U_Base", "set_occ_mat"); - - for (int T = 0; T < ucell.ntype; T++) - { - if (!has_correlated_orbital(T)) continue; - const int l = get_orbital_corr(T); - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - if (this->nspin == 4) - { - for(int mm = 0; mm < occ_mat[iat][l][0][0].nr * occ_mat[iat][l][0][0].nc; mm++) - occ_mat[iat][l][0][0].c[mm] = this->uom_array[pot_uterm_pw_index[iat] + mm]; - } - else if (this->nspin == 1 || this->nspin == 2) - { - const int half_size = this->uom_array.size() / 2; - for(int mm = 0; mm < occ_mat[iat][l][0][0].nr * occ_mat[iat][l][0][0].nc; mm++) - { - occ_mat[iat][l][0][0].c[mm] = this->uom_array[pot_uterm_pw_index[iat] + mm]; - if (this->nspin == 2) - { - occ_mat[iat][l][0][1].c[mm] = this->uom_array[half_size + pot_uterm_pw_index[iat] + mm]; - } - } - } - } - } - - ModuleBase::timer::end("Plus_U_Base", "set_occ_mat"); -} - - -void Plus_U_Base::get_occ_mat_flat(const int iat, const int l, std::vector& occ) const -{ - const int tlp1 = 2 * l + 1; - const int size = tlp1 * tlp1; - if (nspin == 2) - { - for (int is = 0; is < 2; is++) - { - for (int i = 0; i < size; i++) - { - occ[is * size + i] = occ_mat[iat][l][0][is].c[i]; - } - } - } - else - { - for (int i = 0; i < static_cast(occ.size()); i++) - { - occ[i] = occ_mat[iat][l][0][0].c[i]; - } - } -} - - -void Plus_U_Base::set_occ_mat_flat(const int iat, const int l, const int spin, - const std::vector& occ) -{ - for (int i = 0; i < static_cast(occ.size()); i++) - { - occ_mat[iat][l][0][spin].c[i] = occ[i]; - } -} - - -void Plus_U_Base::read_occup_m(const UnitCell& ucell, - const std::string& fn, - const std::string& init_chg, - int nspin, - int npol) -{ - ModuleBase::TITLE("Plus_U_Base", "read_occup_m"); - - if (GlobalV::MY_RANK != 0) - { - return; - } - - std::ifstream ifdftu(fn.c_str(), std::ios::in); - - if (!ifdftu) - { - if (occ_mat_ctrl > 0) - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); - } - else - { - if (init_chg == "file") - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "Can not find the file dm_onsite.txt. Please do scf calculation first"); - } - } - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "Can not open dm_onsite.txt file"); - } - - ifdftu.clear(); - ifdftu.seekg(0); - - char word[20]; - - int T = 0; - int iat = 0; - int spin = 0; - int L = 0; - int zeta = 0; - - ifdftu.rdstate(); - - while (ifdftu.good()) - { - ifdftu >> word; - if (ifdftu.eof()) - { - break; - } - - if (strcmp("Atom=", word) == 0) - { - ifdftu >> iat; - iat -= 1; - ifdftu >> word; - - if (strcmp("L=", word) != 0) - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> L; - ifdftu >> word; - - if (strcmp("ORBITAL=", word) != 0) - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> zeta; - ifdftu.ignore(150, '\n'); - - T = ucell.iat2it[iat]; - const int NL = ucell.atoms[T].nwl + 1; - const int LC = get_orbital_corr(T); - - for (int l = 0; l < NL; l++) - { - if (l != get_orbital_corr(T)) - { - continue; - } - - if (nspin == 1 || nspin == 2) - { - for (int is = 0; is < 2; is++) - { - ifdftu >> word; - if (strcmp("spin=", word) == 0) - { - ifdftu >> spin; - spin -= 1; - ifdftu.ignore(150, '\n'); - - double value = 0.0; - for (int m0 = 0; m0 < 2 * L + 1; m0++) - { - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - ifdftu >> value; - occ_mat[iat][L][zeta][spin](m0, m1) = value; - } - ifdftu.ignore(150, '\n'); - } - } - else - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - } - } - else if (nspin == 4) // SOC - { - double value = 0.0; - for (int m0 = 0; m0 < 2 * L + 1; m0++) - { - for (int ipol0 = 0; ipol0 < npol; ipol0++) - { - const int m0_all = m0 + (2 * L + 1) * ipol0; - - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - for (int ipol1 = 0; ipol1 < npol; ipol1++) - { - int m1_all = m1 + (2 * L + 1) * ipol1; - ifdftu >> value; - occ_mat[iat][L][zeta][0](m0_all, m1_all) = value; - } - } - ifdftu.ignore(150, '\n'); - } - } - } - } - } - else - { - ModuleBase::WARNING_QUIT("Plus_U_Base::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - - ifdftu.rdstate(); - - if (ifdftu.eof() != 0) - { - break; - } - } - - return; -} - - -void Plus_U_Base::local_occup_bcast(const UnitCell& ucell, - int nspin, - int npol) -{ - ModuleBase::TITLE("Plus_U_Base", "local_occup_bcast"); - - for (int T = 0; T < ucell.ntype; T++) - { - if (!has_correlated_orbital(T)) - { - continue; - } - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - const int L = get_orbital_corr(T); - - for (int l = 0; l <= ucell.atoms[T].nwl; l++) - { - if (l != get_orbital_corr(T)) - { - continue; - } - - for (int n = 0; n < ucell.atoms[T].l_nchi[l]; n++) - { - if (n != 0) - { - continue; - } - - if (nspin == 1 || nspin == 2) - { - for (int spin = 0; spin < 2; spin++) - { - for (int m0 = 0; m0 < 2 * l + 1; m0++) - { - for (int m1 = 0; m1 < 2 * l + 1; m1++) - { -#ifdef __MPI - MPI_Bcast(&occ_mat[iat][l][n][spin](m0, m1), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); -#endif - } - } - } - } - else if (nspin == 4) // SOC - { - for (int m0 = 0; m0 < 2 * L + 1; m0++) - { - for (int ipol0 = 0; ipol0 < npol; ipol0++) - { - const int m0_all = m0 + (2 * L + 1) * ipol0; - - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - for (int ipol1 = 0; ipol1 < npol; ipol1++) - { - int m1_all = m1 + (2 * L + 1) * ipol1; -#ifdef __MPI - MPI_Bcast(&occ_mat[iat][l][n][0](m0_all, m1_all), - 1, - MPI_DOUBLE, - 0, - MPI_COMM_WORLD); -#endif - } - } - } - } - } - } - } - } - } - return; -} - - -// cal_occ_pw() is implemented in source_pw/module_pwdft/dftu_cal_occ_pw.cpp -// as a Plus_U_Base method. Pure per-atom kernels live in dftu_tools_pw.{h,cpp} -// as free functions in namespace dftu_pw. +// cal_occ_pw() is implemented in source_pw/module_pwdft/dftu_base_occ.cpp +// as a Plus_U_Base method. Pure per-atom kernels live in dftu_base_tools.{h,cpp} +// as free functions in namespace DFTU_BASE. diff --git a/source/source_pw/module_pwdft/dftu_base.h b/source/source_pw/module_pwdft/dftu_base.h index cb2ae65a186..f91ecf6ad3a 100644 --- a/source/source_pw/module_pwdft/dftu_base.h +++ b/source/source_pw/module_pwdft/dftu_base.h @@ -1,22 +1,25 @@ #ifndef DFTU_BASE_H #define DFTU_BASE_H -#include "source_cell/unitcell.h" -#include "source_estate/module_charge/charge_mixing.h" +#include "source_base/matrix.h" +#include "source_estate/occ_matrix.h" +#include "source_pw/module_pwdft/yukawa_screening.h" +#include +#include #include #include +class UnitCell; +class Charge_Mixing; + class DFTUTest; class Plus_U_Base { friend class DFTUTest; - //============================================================= - // public section - //============================================================= public: Plus_U_Base(); ~Plus_U_Base(); @@ -27,6 +30,7 @@ class Plus_U_Base const int nspin, const std::vector& orbital_corr, const bool yukawa_potential, + const double yukawa_lambda, const std::string& global_readin_dir, const std::string& global_out_dir, const std::string& init_chg, @@ -42,32 +46,23 @@ class Plus_U_Base // --- Accessors for U values and orbital configuration --- double get_u_current(int it) const { return u_current[it]; } - double get_u_target(int it) const { return u_target[it]; } int get_num_u_types() const { return static_cast(u_current.size()); } int get_orbital_corr(int it) const { return orbital_corr[it]; } bool has_correlated_orbital(int it) const { return orbital_corr[it] != -1; } - const int* get_orbital_corr_data() const { return orbital_corr.data(); } /// read-only access to the orbital_corr vector (length ntype) const std::vector& get_orbital_corr_vec() const { return orbital_corr; } - /// read-only access to the iat->(l,n,m,ipol)->iwt lookup table - const std::vector>>>>& - get_iatlnmipol2iwt() const { return iatlnmipol2iwt; } - // --- Accessors for DFT+U configuration --- double get_uramping() const { return uramping; } int get_occ_mat_ctrl() const { return occ_mat_ctrl; } int get_cal_type() const { return cal_type; } - bool use_yukawa() const { return use_yukawa_; } - - double get_U_Yukawa(int it, int l, int n) const { return U_Yukawa[it][l][n]; } - double get_J_Yukawa(int it, int l, int n) const { return J_Yukawa[it][l][n]; } - void set_U_Yukawa(int it, int l, int n, double val) { U_Yukawa[it][l][n] = val; } - void set_J_Yukawa(int it, int l, int n, double val) { J_Yukawa[it][l][n] = val; } - double get_lambda() const { return lambda; } - void set_lambda(double l) { lambda = l; } - std::vector>>>& get_Fk_data() { return Fk; } + bool use_yukawa() const { return yukawa_ != nullptr; } + + /// access the Yukawa screening object (non-null only when use_yukawa()) + YukawaScreening& yukawa() { return *yukawa_; } + const YukawaScreening& yukawa() const { return *yukawa_; } + void set_u_current(int it, double val) { u_current[it] = val; } double get_energy() const { return energy_u; } @@ -108,14 +103,6 @@ class Plus_U_Base : static_cast(pot_uterm_pw.size()); } - /// get effective potential matrix for PW base (per-atom, raw index) - /// @deprecated Use get_pot_uterm_pw_spin() for nspin-aware access. - [[deprecated("Use get_pot_uterm_pw_spin() for nspin-aware access")]] - const std::complex* get_pot_uterm_pw(const int iat) const - { - return &(pot_uterm_pw[pot_uterm_pw_index[iat]]); - } - int get_size_pot_uterm_pw() const { return pot_uterm_pw.size(); @@ -129,45 +116,9 @@ class Plus_U_Base bool is_mixing_enabled() const { return mixing_dftu != 0; } void enable_mixing() { mixing_dftu = 1; } - /// get occupation matrix element occ_mat[iat][l][n][spin](m1,m2) - double get_occ_mat(const int iat, const int l, const int n, const int spin, - const int m1, const int m2) const - { - return occ_mat[iat][l][n][spin](m1, m2); - } - - /// get saved occupation matrix element occ_mat_save[iat][l][n][spin](m1,m2) - double get_occ_mat_save(const int iat, const int l, const int n, const int spin, - const int m1, const int m2) const - { - return occ_mat_save[iat][l][n][spin](m1, m2); - } - - /// set occupation matrix element occ_mat[iat][l][n][spin](m1,m2) - void set_occ_mat(const int iat, const int l, const int n, const int spin, - const int m1, const int m2, const double val) - { - occ_mat[iat][l][n][spin](m1, m2) = val; - } - - /// get reference to occ_mat data - std::vector>>>& get_occ_mat_data() { return occ_mat; } - /// get reference to occ_mat_save data - std::vector>>>& get_occ_mat_save_data() { return occ_mat_save; } - /// get occ_mat_initialized flag - bool get_occ_mat_initialized() const { return occ_mat_initialized; } - /// set occ_mat_initialized flag - void set_occ_mat_initialized(bool val) { occ_mat_initialized = val; } - - /// get flat occupation matrix for an atom's correlated orbital. - /// nspin=1: fills occ with occ_mat[iat][l][0][0] data - /// nspin=2: fills occ with interleaved occ_mat[iat][l][0][0] and [1] data - /// nspin=4: fills occ with occ_mat[iat][l][0][0] data (all 4 Pauli blocks) - void get_occ_mat_flat(const int iat, const int l, std::vector& occ) const; - - /// set flat occupation matrix for an atom's correlated orbital (write-back) - void set_occ_mat_flat(const int iat, const int l, const int spin, - const std::vector& occ); + /// direct access to the occupation matrix object (new write path) + OccupationMatrix& occmat() { return occmat_; } + const OccupationMatrix& occmat() const { return occmat_; } protected: // --- U values and orbital configuration (set in init_base) --- @@ -180,15 +131,13 @@ class Plus_U_Base int occ_mat_ctrl = 0; int mixing_dftu = 0; int nspin = 0; - bool use_yukawa_ = false; // --- State flags --- // dftu can be calculated only after occ_mat has been initialized bool occ_mat_initialized = false; // --- Occupation matrices --- - std::vector>>> occ_mat; - std::vector>>> occ_mat_save; + OccupationMatrix occmat_; // --- Internal state --- double energy_u = 0.0; @@ -197,51 +146,14 @@ class Plus_U_Base std::string device; int kpar = 1; - // transform between iwt index and it, ia, L, N and m index - std::vector>>>> - iatlnmipol2iwt; - - void copy_occ_mat(const UnitCell& ucell); - void zero_occ_mat(const UnitCell& ucell); - void mix_occ_mat(const UnitCell& ucell, const double& mixing_beta); - void set_occ_mat(const UnitCell& ucell); - - /// accumulate occ_mat from psi for all k-points (per-device template) - template - void accumulate_occ_one_k(const void* psi_in, - const ModuleBase::matrix& wg_in, - const UnitCell& cell, - const int* isk); - - /// reduce occ_mat across k-pools (per-atom, nspin-aware) - void reduce_occ_mat(const UnitCell& cell); - - /// copy occ_mat to uom_array for mixing (nspin-aware split layout) - void sync_occ_to_uom(const UnitCell& cell); - - /// compute effective potential pot_onsite and DFT+U energy from occ_mat - /// (assumes occ_mat has already been reduced across k-pools) - void compute_eff_pot_and_energy(const UnitCell& cell); - std::vector> pot_uterm_pw; std::vector pot_uterm_pw_index; std::vector uom_array; std::vector uom_save; - // Yukawa-related members (base part, no LCAO dependency) - double lambda = 0.0; - std::vector>>> Fk; - std::vector>> U_Yukawa; - std::vector>> J_Yukawa; - - void read_occup_m(const UnitCell& ucell, - const std::string& fn, - const std::string& init_chg, - int nspin, - int npol); - void local_occup_bcast(const UnitCell& ucell, - int nspin, - int npol); + // Yukawa screening object; constructed only when use_yukawa() is true. + // Owns the screening length, Slater integrals and derived U/J. + std::unique_ptr yukawa_; }; diff --git a/source/source_pw/module_pwdft/dftu_output.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp similarity index 56% rename from source/source_pw/module_pwdft/dftu_output.cpp rename to source/source_pw/module_pwdft/dftu_base_io.cpp index 0f2aee7615f..04ca5278877 100644 --- a/source/source_pw/module_pwdft/dftu_output.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -1,18 +1,25 @@ -#include "source_pw/module_pwdft/dftu_output.h" +#include "source_pw/module_pwdft/dftu_base_io.h" +#include "source_cell/unitcell.h" #include "source_pw/module_pwdft/dftu_base.h" #include "source_base/constants.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include "source_base/parallel_common.h" +#include "source_base/parallel_global.h" #include "source_base/timer.h" #include +#include #include #include #include -// local inline helpers for eigenvalue calculation +// local helpers for eigenvalue calculation // migrated from dftu_base.cpp, mohan 2025-11-08 +namespace +{ + inline void JacobiRotate(std::vector>& A, int p, int q, int n) { if (std::abs(A[p][q]) > 1e-10) @@ -79,9 +86,226 @@ inline std::vector CalculateEigenvalues(std::vector> return eigenvalues; } +} // namespace + + +namespace DFTU_BASE +{ + +void read_occup_m(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& orbital_corr, + const int occ_mat_ctrl, + const std::string& fn, + const std::string& init_chg, + int nspin, + int npol) +{ + ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); + + if (GlobalV::MY_RANK != 0) + { + return; + } + + std::ifstream ifdftu(fn.c_str(), std::ios::in); + + if (!ifdftu) + { + if (occ_mat_ctrl > 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); + } + else + { + if (init_chg == "file") + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite.txt. Please do scf calculation first"); + } + } + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not open dm_onsite.txt file"); + } + + ifdftu.clear(); + ifdftu.seekg(0); + + char word[20]; + + int T = 0; + int iat = 0; + int spin = 0; + int L = 0; + int zeta = 0; + + ifdftu.rdstate(); + + while (ifdftu.good()) + { + ifdftu >> word; + if (ifdftu.eof()) + { + break; + } + + if (strcmp("Atom=", word) == 0) + { + ifdftu >> iat; + iat -= 1; + ifdftu >> word; + + if (strcmp("L=", word) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + ifdftu >> L; + ifdftu >> word; + + if (strcmp("ORBITAL=", word) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + ifdftu >> zeta; + ifdftu.ignore(150, '\n'); + + T = ucell.iat2it[iat]; + const int NL = ucell.atoms[T].nwl + 1; + + for (int l = 0; l < NL; l++) + { + if (l != orbital_corr[T]) + { + continue; + } + + if (nspin == 1 || nspin == 2) + { + for (int is = 0; is < 2; is++) + { + ifdftu >> word; + if (strcmp("spin=", word) == 0) + { + ifdftu >> spin; + spin -= 1; + ifdftu.ignore(150, '\n'); + + double value = 0.0; + for (int m0 = 0; m0 < 2 * L + 1; m0++) + { + for (int m1 = 0; m1 < 2 * L + 1; m1++) + { + ifdftu >> value; + occ.set(iat, L, zeta, spin, m0, m1, value); + } + ifdftu.ignore(150, '\n'); + } + } + else + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + } + } + else if (nspin == 4) // SOC + { + double value = 0.0; + for (int m0 = 0; m0 < 2 * L + 1; m0++) + { + for (int ipol0 = 0; ipol0 < npol; ipol0++) + { + const int m0_all = m0 + (2 * L + 1) * ipol0; + + for (int m1 = 0; m1 < 2 * L + 1; m1++) + { + for (int ipol1 = 0; ipol1 < npol; ipol1++) + { + int m1_all = m1 + (2 * L + 1) * ipol1; + ifdftu >> value; + occ.set(iat, L, zeta, 0, m0_all, m1_all, value); + } + } + ifdftu.ignore(150, '\n'); + } + } + } + } + } + else + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } -namespace dftu_io + ifdftu.rdstate(); + + if (ifdftu.eof() != 0) + { + break; + } + } + + return; +} + +#ifdef __MPI +/// Broadcast the local occupation number matrices from rank 0 to all ranks. +/// +/// Each occupation matrix is broadcast as one contiguous block +/// (matrix::c stores nr * nc consecutive doubles) instead of element by +/// element. +void local_occup_bcast(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& orbital_corr, + int nspin, + int npol) { + ModuleBase::TITLE("DFTU_BASE", "local_occup_bcast"); + + for (int T = 0; T < ucell.ntype; T++) + { + if (orbital_corr[T] == -1) + { + continue; + } + + for (int I = 0; I < ucell.atoms[T].na; I++) + { + const int iat = ucell.itia2iat(T, I); + const int L = orbital_corr[T]; + + for (int l = 0; l <= ucell.atoms[T].nwl; l++) + { + if (l != orbital_corr[T]) + { + continue; + } + + for (int n = 0; n < ucell.atoms[T].l_nchi[l]; n++) + { + if (n != 0) + { + continue; + } + + if (nspin == 1 || nspin == 2) + { + for (int spin = 0; spin < 2; spin++) + { + Parallel_Common::bcast_double(occ.mat(iat, l, n, spin).c, + occ.mat(iat, l, n, spin).nr * occ.mat(iat, l, n, spin).nc); + } + } + else if (nspin == 4) // SOC + { + Parallel_Common::bcast_double(occ.mat(iat, l, n, 0).c, + occ.mat(iat, l, n, 0).nr * occ.mat(iat, l, n, 0).nc); + } + } + } + } + } + return; +} +#endif + void output(const Plus_U_Base& dftu, const UnitCell& ucell, @@ -90,7 +314,7 @@ void output(const Plus_U_Base& dftu, int nspin, int npol) { - ModuleBase::TITLE("dftu_io", "output"); + ModuleBase::TITLE("DFTU_BASE", "output"); GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; GlobalV::ofs_running << " | #DFT+U INFORMATION# |" << std::endl; @@ -124,10 +348,10 @@ void output(const Plus_U_Base& dftu, { continue; } - double Ueff = (dftu.get_U_Yukawa(T, L, n) - dftu.get_J_Yukawa(T, L, n)) * ModuleBase::Ry_to_eV; + double Ueff = (dftu.yukawa().get_U(T, L, n) - dftu.yukawa().get_J(T, L, n)) * ModuleBase::Ry_to_eV; GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << n - << " U=" << dftu.get_U_Yukawa(T, L, n) * ModuleBase::Ry_to_eV << " eV" - << " J=" << dftu.get_J_Yukawa(T, L, n) * ModuleBase::Ry_to_eV << " eV" + << " U=" << dftu.yukawa().get_U(T, L, n) * ModuleBase::Ry_to_eV << " eV" + << " J=" << dftu.yukawa().get_J(T, L, n) * ModuleBase::Ry_to_eV << " eV" << std::endl; } } @@ -136,7 +360,7 @@ void output(const Plus_U_Base& dftu, } GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; - dftu_io::write_occup_m(dftu, ucell, GlobalV::ofs_running, true, nspin, npol); + write_occup_m(dftu, ucell, GlobalV::ofs_running, true, nspin, npol); // Write dm_onsite.txt if (out_chg && GlobalV::MY_RANK == 0) @@ -145,9 +369,9 @@ void output(const Plus_U_Base& dftu, ofdftu.open(global_out_dir + "dm_onsite.txt"); if (!ofdftu) { - ModuleBase::WARNING_QUIT("dftu_io::output", "Can't create file dm_onsite.txt"); + ModuleBase::WARNING_QUIT("DFTU_BASE::output", "Can't create file dm_onsite.txt"); } - dftu_io::write_occup_m(dftu, ucell, ofdftu, false, nspin, npol); + write_occup_m(dftu, ucell, ofdftu, false, nspin, npol); ofdftu.close(); } @@ -166,7 +390,7 @@ void write_occup_m(const Plus_U_Base& dftu, int nspin, int npol) { - ModuleBase::TITLE("dftu_io", "write_occup_m"); + ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); if (GlobalV::MY_RANK != 0) { @@ -218,7 +442,7 @@ void write_occup_m(const Plus_U_Base& dftu, { for (int m1 = 0; m1 < 2 * l + 1; m1++) { - A[m0][m1] = dftu.get_occ_mat(iat, l, n, is, m0, m1); + A[m0][m1] = dftu.occmat().get(iat, l, n, is, m0, m1); } } std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); @@ -240,7 +464,7 @@ void write_occup_m(const Plus_U_Base& dftu, for (int m1 = 0; m1 < 2 * l + 1; m1++) { ofs << std::setw(12) - << dftu.get_occ_mat(iat, l, n, is, m0, m1); + << dftu.occmat().get(iat, l, n, is, m0, m1); } ofs << std::endl; } @@ -265,7 +489,7 @@ void write_occup_m(const Plus_U_Base& dftu, { for (int m1 = 0; m1 < 2 * l + 1; m1++) { - A[m0][m1] = dftu.get_occ_mat(iat, l, n, 0, m0, m1); + A[m0][m1] = dftu.occmat().get(iat, l, n, 0, m0, m1); index++; } } @@ -299,7 +523,7 @@ void write_occup_m(const Plus_U_Base& dftu, { int m1_all = m1 + (2 * l + 1) * ipol1; ofs << std::setw(12) << std::setprecision(8) << std::fixed - << dftu.get_occ_mat(iat, l, n, 0, m0_all, m1_all); + << dftu.occmat().get(iat, l, n, 0, m0_all, m1_all); } } ofs << std::endl; @@ -316,4 +540,4 @@ void write_occup_m(const Plus_U_Base& dftu, } -} // namespace dftu_io +} // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h new file mode 100644 index 00000000000..ed48701f2e5 --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -0,0 +1,70 @@ +#ifndef DFTU_BASE_IO_H +#define DFTU_BASE_IO_H + +#include "source_base/matrix.h" +#include "source_estate/occ_matrix.h" + +#include +#include +#include + +class Plus_U_Base; +class UnitCell; + +namespace DFTU_BASE +{ + +/// nested occupation-matrix type used by DFT+U: occ_mat[iat][l][n][spin](m0, m1) +using OccMatData = std::vector>>>; + +/// Read the local occupation number matrix from file (rank 0 only). +/// +/// The file format matches the output of write_occup_m(). When the file can +/// not be opened, the run quits with an error message that depends on +/// occ_mat_ctrl and init_chg. +void read_occup_m(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& orbital_corr, + const int occ_mat_ctrl, + const std::string& fn, + const std::string& init_chg, + int nspin, + int npol); + +/// Broadcast the local occupation number matrices from rank 0 to all ranks. +/// +/// Implemented in dftu_base_io.cpp (only available in MPI builds). +void local_occup_bcast(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& orbital_corr, + int nspin, + int npol); + +/// Output DFT+U information (Hubbard U/J, local occupation matrices) to the +/// running log and, when out_chg is set, to the dm_onsite.txt file. +/// +/// Extracted from Plus_U_Base::output as a free function so that IO logic is +/// decoupled from the Plus_U_Base class. The function only reads the +/// Plus_U_Base state via public accessors; no friend declaration needed. +void output(const Plus_U_Base& dftu, + const UnitCell& ucell, + bool out_chg, + const std::string& global_out_dir, + int nspin, + int npol); + +/// Write local occupation matrices to the given stream. +/// +/// Extracted from Plus_U_Base::write_occup_m. When diag is true, eigenvalues +/// and magnetism are also printed; otherwise only raw matrix elements. +/// Caller is responsible for opening/closing the stream. +void write_occup_m(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ofstream& ofs, + bool diag, + int nspin, + int npol); + +} // namespace DFTU_BASE + +#endif diff --git a/source/source_pw/module_pwdft/dftu_base_occ.cpp b/source/source_pw/module_pwdft/dftu_base_occ.cpp new file mode 100644 index 00000000000..42b40e96a43 --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_occ.cpp @@ -0,0 +1,247 @@ +#include "source_pw/module_pwdft/dftu_base.h" +#include "source_pw/module_pwdft/dftu_base_io.h" +#include "source_pw/module_pwdft/dftu_base_tools.h" +#include "source_pw/module_pwdft/onsite_proj.h" +#include "source_cell/unitcell.h" +#include "source_estate/module_charge/charge_mixing.h" +#include "source_base/parallel_reduce.h" +#include "source_base/global_variable.h" +#include "source_base/timer.h" +#include "source_base/parallel_global.h" + + + + +/// calculate occupation matrix for DFT+U (PW basis) +/// +/// nspin=1 (npol=1): single spin channel; occ_mat[iat][l][n][0] only; +/// pot_uterm_pw has one block of tlp1^2 per atom. +/// +/// nspin=2 (npol=1): two spin channels stored separately: +/// occ_mat[iat][l][n][0] = spin-up, occ_mat[iat][l][n][1] = spin-down; +/// becp indices: ib*nkb + begin_ih + m (same formula for both spins); +/// spin channel selected by `isk[ik]` (not ik >= nk/2, which fails for kpar>1); +/// +/// nspin=4 (npol=2): spinor calculation; +/// occ_mat has a single matrix of size (2*tlp1) x (2*tlp1) per atom +/// storing all 4 Pauli blocks contiguously. +void Plus_U_Base::cal_occ_pw(const void* psi_in, + const ModuleBase::matrix& wg_in, + const UnitCell& cell, + Charge_Mixing* p_chgmix, + const int* isk) +{ + ModuleBase::timer::start("Plus_U_Base", "cal_occ_pw"); + this->occmat_.copy_to_save(cell, this->orbital_corr); + this->occmat_.write_save_to_flat(cell, this->orbital_corr, + this->pot_uterm_pw_index, this->uom_save); + this->occmat_.zero(cell, this->orbital_corr); + + if(this->device == "cpu") + { + DFTU_BASE::accumulate_occ_one_k( + psi_in, wg_in, cell, isk, this->nspin, this->orbital_corr, this->occmat_); + } +#if defined(__CUDA) || defined(__ROCM) + else + { + DFTU_BASE::accumulate_occ_one_k( + psi_in, wg_in, cell, isk, this->nspin, this->orbital_corr, this->occmat_); + } +#endif + + // reduce occ_mat across k-pools, then copy to uom_array for mixing + DFTU_BASE::reduce_occ_mat(cell, this->nspin, this->kpar, + this->orbital_corr, this->occmat_); + this->occmat_.write_to_flat(cell, this->orbital_corr, + this->pot_uterm_pw_index, this->uom_array); + + // mixing + if(is_mixing_enabled() && p_chgmix != nullptr) + { + p_chgmix->mix_uom(this->uom_array, this->uom_save); + this->occmat_.read_from_flat(cell, this->orbital_corr, + this->pot_uterm_pw_index, this->uom_array); + } + + DFTU_BASE::compute_pot_uterm_and_energy(cell, this->nspin, + this->u_current, this->orbital_corr, this->pot_uterm_pw_index, + this->occmat_, this->pot_uterm_pw, this->energy_u); + + ModuleBase::timer::end("Plus_U_Base", "cal_occ_pw"); +} + +namespace DFTU_BASE { + +void reduce_occ_mat(const UnitCell& cell, + const int nspin, + const int kpar, + const std::vector& orbital_corr, + OccupationMatrix& occmat) +{ + for(int iat = 0; iat < cell.nat; iat++) + { + const int it = cell.iat2it[iat]; + const int target_l = orbital_corr[it]; + if(target_l == -1) + { + continue; + } + const int size = (2 * target_l + 1) * (2 * target_l + 1); + + if(nspin != 4) + { + Parallel_Reduce::reduce_double_allpool(kpar, + GlobalV::NPROC_IN_POOL, + occmat.mat(iat, target_l, 0, 0).c, + size); + if(nspin == 2) + { + Parallel_Reduce::reduce_double_allpool(kpar, + GlobalV::NPROC_IN_POOL, + occmat.mat(iat, target_l, 0, 1).c, + size); + } + } + else + { + Parallel_Reduce::reduce_double_allpool(kpar, + GlobalV::NPROC_IN_POOL, + occmat.mat(iat, target_l, 0, 0).c, + size * 4); + } + } +} + +void compute_pot_uterm_and_energy(const UnitCell& cell, + const int nspin, + const std::vector& u_current, + const std::vector& orbital_corr, + const std::vector& pot_uterm_pw_index, + const OccupationMatrix& occmat, + std::vector>& pot_uterm_pw, + double& energy_u) +{ + energy_u = 0.0; + const double weight_eu = (nspin == 1) ? 1.0 : (nspin == 2) ? 0.5 : 0.25; + const double diag_coeff = (nspin == 4) ? 1.0 : 0.5; + // calculate pot_onsite and energy (occ_mat already reduced above) + for(int iat = 0; iat < cell.nat; iat++) + { + const int it = cell.iat2it[iat]; + const int target_l = orbital_corr[it]; + if(target_l == -1) + { + continue; + } + const int size = (2 * target_l + 1) * (2 * target_l + 1); + + //update effective potential + const double u_value = u_current[it]; + std::complex* pot_onsite_iat = &(pot_uterm_pw[pot_uterm_pw_index[iat]]); + const int m_size = 2 * target_l + 1; + + if(nspin == 4) + { + // pot_onsite is stored as 4 contiguous Pauli blocks per atom: + // is=0: charge channel (identity), Hubbard U contributes the + // diagonal term diag_coeff*delta(m1,m2) + // is=1,2,3: spin channels (sigma_x/y/z), no U diagonal term + // The occupation matrix occ_mat[...][0][0].c packs all 4 blocks + // contiguously, each of size m_size*m_size. + energy_u += compute_pot_onsite_spinor( + pot_onsite_iat, + occmat.mat(iat, target_l, 0, 0).c, + u_value, diag_coeff, weight_eu, m_size); + } + else // nspin=1 or nspin=2 + { + // spin-up channel + energy_u += compute_pot_onsite_scalar( + pot_onsite_iat, + occmat.mat(iat, target_l, 0, 0).c, + u_value, diag_coeff, weight_eu, m_size); + // spin-down channel for nspin=2 + if(nspin == 2) + { + std::complex* pot_onsite_iat1 = &(pot_uterm_pw[pot_uterm_pw.size()/2 + pot_uterm_pw_index[iat]]); + energy_u += compute_pot_onsite_scalar( + pot_onsite_iat1, + occmat.mat(iat, target_l, 0, 1).c, + u_value, diag_coeff, weight_eu, m_size); + } + } + } +} + +} // namespace DFTU_BASE + +namespace DFTU_BASE { + +template +void accumulate_occ_one_k(const void* psi_in, + const ModuleBase::matrix& wg_in, + const UnitCell& cell, + const int* isk, + const int nspin, + const std::vector& orbital_corr, + OccupationMatrix& occmat) +{ + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + const psi::Psi, Device>* psi_p = + (const psi::Psi, Device>*)psi_in; + const int nbands = psi_p->get_nbands(); + const int npol = psi_p->get_npol(); + for(int ik = 0; ik < psi_p->get_nk(); ik++) + { + int is = (nspin == 2) ? isk[ik] : 0; + psi_p->fix_k(ik); + onsite_p->tabulate_atomic(ik); + + onsite_p->overlap_proj_psi(nbands*npol, psi_p->get_pointer()); + const std::complex* becp = onsite_p->get_h_becp(); + int nkb = onsite_p->get_size_becp() / nbands / npol; + + int begin_ih = 0; + for(int iat = 0; iat < cell.nat; iat++) + { + const int it = cell.iat2it[iat]; + const int nh = onsite_p->get_nh(iat); + const int target_l = orbital_corr[it]; + if(target_l == -1) + { + begin_ih += nh; + continue; + } + const int m_begin = target_l * target_l; + const int tlp1 = 2 * target_l + 1; + if(nspin == 4) + { + accumulate_occ_spinor( + occmat.mat(iat, target_l, 0, 0).c, + becp, nbands, npol, nkb, begin_ih, m_begin, tlp1, + wg_in, ik); + } + else // nspin=1 or nspin=2 + { + accumulate_occ_scalar( + occmat.mat(iat, target_l, 0, is).c, + becp, nbands, nkb, begin_ih, m_begin, tlp1, + wg_in, ik); + } + begin_ih += nh; + } + } +} + +} // namespace DFTU_BASE + +// explicit instantiations +template void DFTU_BASE::accumulate_occ_one_k( + const void*, const ModuleBase::matrix&, const UnitCell&, const int*, + const int, const std::vector&, OccupationMatrix&); +#if defined(__CUDA) || defined(__ROCM) +template void DFTU_BASE::accumulate_occ_one_k( + const void*, const ModuleBase::matrix&, const UnitCell&, const int*, + const int, const std::vector&, OccupationMatrix&); +#endif diff --git a/source/source_pw/module_pwdft/dftu_tools_pw.cpp b/source/source_pw/module_pwdft/dftu_base_tools.cpp similarity index 97% rename from source/source_pw/module_pwdft/dftu_tools_pw.cpp rename to source/source_pw/module_pwdft/dftu_base_tools.cpp index 036d2d8625e..d253eeed62e 100644 --- a/source/source_pw/module_pwdft/dftu_tools_pw.cpp +++ b/source/source_pw/module_pwdft/dftu_base_tools.cpp @@ -1,6 +1,6 @@ -#include "source_pw/module_pwdft/dftu_tools_pw.h" +#include "source_pw/module_pwdft/dftu_base_tools.h" -namespace dftu_pw { +namespace DFTU_BASE { void pauli_to_spin_basis(std::complex* pot_onsite, int m_size) { @@ -146,4 +146,4 @@ void accumulate_occ_scalar( } } -} // namespace dftu_pw +} // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_tools_pw.h b/source/source_pw/module_pwdft/dftu_base_tools.h similarity index 58% rename from source/source_pw/module_pwdft/dftu_tools_pw.h rename to source/source_pw/module_pwdft/dftu_base_tools.h index 188a5127ce9..57aca9226c4 100644 --- a/source/source_pw/module_pwdft/dftu_tools_pw.h +++ b/source/source_pw/module_pwdft/dftu_base_tools.h @@ -1,16 +1,20 @@ -#ifndef DFTU_TOOLS_PW_H -#define DFTU_TOOLS_PW_H +#ifndef DFTU_BASE_TOOLS_H +#define DFTU_BASE_TOOLS_H #include +#include #include "source_base/matrix.h" +class UnitCell; +class OccupationMatrix; + /// Free functions for DFT+U PW basis calculations. /// /// These functions are pure (no access to Plus_U_Base members) so they can be /// unit-tested directly by including this header. The member functions in -/// dftu_pw.cpp call them after computing per-atom offsets and fetching the -/// relevant member state (occ_mat, pot_uterm_pw, u_current, etc.). -namespace dftu_pw { +/// dftu_base_occ.cpp call them after computing per-atom offsets and fetching +/// the relevant member state (occ_mat, pot_uterm_pw, u_current, etc.). +namespace DFTU_BASE { /// transform pot_onsite from Pauli basis to spin basis (in-place, nspin==4 only). /// @@ -96,6 +100,53 @@ void accumulate_occ_scalar( const ModuleBase::matrix& wg, int ik); -} // namespace dftu_pw +/// reduce occ_mat across all k-pools (per-atom, nspin-aware). +/// +/// Each k-pool only accumulates occ_mat contributions from the k-points it +/// owns; this sums them across pools so occmat holds the full result. +/// nspin=1: single channel, size elements +/// nspin=2: two channels (spin-up/down) reduced separately +/// nspin=4: 4 Pauli blocks packed contiguously, reduced in one shot +void reduce_occ_mat(const UnitCell& cell, + const int nspin, + const int kpar, + const std::vector& orbital_corr, + OccupationMatrix& occmat); + +/// compute effective potential pot_onsite and DFT+U energy from occ_mat. +/// +/// Preconditions: +/// - occmat has been accumulated from psi and reduced across k-pools. +/// +/// Outputs: +/// - pot_uterm_pw: pot_onsite = U * (diag*delta - occ) written per atom +/// nspin=4: 4 Pauli blocks per atom, then transformed to spin basis +/// nspin=1: single channel +/// nspin=2: two channels in split layout [all_up | all_dn] +/// - energy_u (out): E_U = sum U * weight_eu * occ(m2,m1) * occ(m1,m2), +/// overwritten with the total energy of this call +void compute_pot_uterm_and_energy(const UnitCell& cell, + const int nspin, + const std::vector& u_current, + const std::vector& orbital_corr, + const std::vector& pot_uterm_pw_index, + const OccupationMatrix& occmat, + std::vector>& pot_uterm_pw, + double& energy_u); + +/// accumulate occ_mat from psi for all k-points (per-device template). +/// +/// Explicitly instantiated for DEVICE_CPU (and DEVICE_GPU when available) +/// in dftu_base_occ.cpp. +template +void accumulate_occ_one_k(const void* psi_in, + const ModuleBase::matrix& wg_in, + const UnitCell& cell, + const int* isk, + const int nspin, + const std::vector& orbital_corr, + OccupationMatrix& occmat); + +} // namespace DFTU_BASE #endif diff --git a/source/source_pw/module_pwdft/dftu_cal_occ_pw.cpp b/source/source_pw/module_pwdft/dftu_cal_occ_pw.cpp deleted file mode 100644 index 54848f456ca..00000000000 --- a/source/source_pw/module_pwdft/dftu_cal_occ_pw.cpp +++ /dev/null @@ -1,269 +0,0 @@ -#include "source_pw/module_pwdft/dftu_base.h" -#include "source_pw/module_pwdft/dftu_tools_pw.h" -#include "source_pw/module_pwdft/onsite_proj.h" -#include "source_base/parallel_reduce.h" -#include "source_base/global_variable.h" -#include "source_base/timer.h" -#include "source_base/parallel_global.h" - -/// calculate occupation matrix for DFT+U (PW basis) -/// -/// nspin=1 (npol=1): single spin channel; occ_mat[iat][l][n][0] only; -/// pot_uterm_pw has one block of tlp1^2 per atom. -/// -/// nspin=2 (npol=1): two spin channels stored separately: -/// occ_mat[iat][l][n][0] = spin-up, occ_mat[iat][l][n][1] = spin-down; -/// becp indices: ib*nkb + begin_ih + m (same formula for both spins); -/// spin channel selected by `isk[ik]` (not ik >= nk/2, which fails for kpar>1); -/// -/// nspin=4 (npol=2): spinor calculation; -/// occ_mat has a single matrix of size (2*tlp1) x (2*tlp1) per atom -/// storing all 4 Pauli blocks contiguously. -void Plus_U_Base::cal_occ_pw(const void* psi_in, - const ModuleBase::matrix& wg_in, - const UnitCell& cell, - Charge_Mixing* p_chgmix, - const int* isk) -{ - ModuleBase::timer::start("Plus_U_Base", "cal_occ_pw"); - this->copy_occ_mat(cell); - this->zero_occ_mat(cell); - - if(this->device == "cpu") - { - this->accumulate_occ_one_k(psi_in, wg_in, cell, isk); - } -#if defined(__CUDA) || defined(__ROCM) - else - { - this->accumulate_occ_one_k(psi_in, wg_in, cell, isk); - } -#endif - - // reduce occ_mat across k-pools, then copy to uom_array for mixing - this->reduce_occ_mat(cell); - this->sync_occ_to_uom(cell); - - // mixing - if(is_mixing_enabled() && p_chgmix != nullptr) - { - p_chgmix->mix_uom(this->uom_array, this->uom_save); - this->set_occ_mat(cell); - } - - this->compute_eff_pot_and_energy(cell); - - ModuleBase::timer::end("Plus_U_Base", "cal_occ_pw"); -} - -/// reduce occ_mat across all k-pools. -/// -/// Each k-pool only accumulates occ_mat contributions from the k-points it -/// owns; this sums them across pools so occ_mat holds the full result. -/// nspin=1: single channel, size elements -/// nspin=2: two channels (spin-up/down) reduced separately -/// nspin=4: 4 Pauli blocks packed contiguously, reduced in one shot -void Plus_U_Base::reduce_occ_mat(const UnitCell& cell) -{ - for(int iat = 0; iat < cell.nat; iat++) - { - const int it = cell.iat2it[iat]; - const int target_l = get_orbital_corr(it); - if(!has_correlated_orbital(it)) - { - continue; - } - const int size = (2 * target_l + 1) * (2 * target_l + 1); - - if(this->nspin != 4) - { - Parallel_Reduce::reduce_double_allpool(this->kpar, - GlobalV::NPROC_IN_POOL, - this->occ_mat[iat][target_l][0][0].c, - size); - if(this->nspin == 2) - { - Parallel_Reduce::reduce_double_allpool(this->kpar, - GlobalV::NPROC_IN_POOL, - this->occ_mat[iat][target_l][0][1].c, - size); - } - } - else - { - Parallel_Reduce::reduce_double_allpool(this->kpar, - GlobalV::NPROC_IN_POOL, - this->occ_mat[iat][target_l][0][0].c, - size * 4); - } - } -} - -/// copy occ_mat to uom_array for mixing. -/// -/// Layout: -/// nspin=1: uom_array[pot_uterm_pw_index[iat] + mm] = occ_mat[...][0][0] -/// nspin=2: split layout [all_up | all_dn], each atom's spin-up in the -/// first half and spin-down in the second half, both indexed by -/// pot_uterm_pw_index[iat] -/// nspin=4: not used here (uom_array mixing only covers nspin=1/2 in the -/// current code path; the nspin=4 branch is a no-op) -void Plus_U_Base::sync_occ_to_uom(const UnitCell& cell) -{ - if(this->uom_array.size() == 0) - { - return; - } - for(int iat = 0; iat < cell.nat; iat++) - { - const int it = cell.iat2it[iat]; - const int target_l = get_orbital_corr(it); - if(!has_correlated_orbital(it)) - { - continue; - } - const int size = (2 * target_l + 1) * (2 * target_l + 1); - - for(int mm = 0; mm < size; mm++) - { - this->uom_array[pot_uterm_pw_index[iat] + mm] = - this->occ_mat[iat][target_l][0][0].c[mm]; - } - if(this->nspin == 2) - { - const int half_size = this->uom_array.size() / 2; - for(int mm = 0; mm < size; mm++) - { - this->uom_array[half_size + pot_uterm_pw_index[iat] + mm] = - this->occ_mat[iat][target_l][0][1].c[mm]; - } - } - } -} - -/// compute effective potential pot_onsite and DFT+U energy from occ_mat. -/// -/// Preconditions: -/// - occ_mat has been accumulated from psi and reduced across k-pools -/// (cal_occ_pw calls this after the reduce + mixing steps). -/// -/// Outputs: -/// - pot_uterm_pw: pot_onsite = U * (diag*delta - occ) written per atom -/// nspin=4: 4 Pauli blocks per atom, then transformed to spin basis -/// nspin=1: single channel -/// nspin=2: two channels in split layout [all_up | all_dn] -/// - energy_u: E_U = sum U * weight_eu * occ(m2,m1) * occ(m1,m2) -void Plus_U_Base::compute_eff_pot_and_energy(const UnitCell& cell) -{ - this->energy_u = 0.0; - const double weight_eu = (this->nspin == 1) ? 1.0 : (this->nspin == 2) ? 0.5 : 0.25; - const double diag_coeff = (this->nspin == 4) ? 1.0 : 0.5; - // calculate pot_onsite and energy (occ_mat already reduced above) - for(int iat = 0; iat < cell.nat; iat++) - { - const int it = cell.iat2it[iat]; - const int target_l = get_orbital_corr(it); - if(!has_correlated_orbital(it)) - { - continue; - } - const int size = (2 * target_l + 1) * (2 * target_l + 1); - - //update effective potential - const double u_value = this->u_current[it]; - std::complex* pot_onsite_iat = &(this->pot_uterm_pw[this->pot_uterm_pw_index[iat]]); - const int m_size = 2 * target_l + 1; - - if(this->nspin == 4) - { - // pot_onsite is stored as 4 contiguous Pauli blocks per atom: - // is=0: charge channel (identity), Hubbard U contributes the - // diagonal term diag_coeff*delta(m1,m2) - // is=1,2,3: spin channels (sigma_x/y/z), no U diagonal term - // The occupation matrix occ_mat[...][0][0].c packs all 4 blocks - // contiguously, each of size m_size*m_size. - this->energy_u += dftu_pw::compute_pot_onsite_spinor( - pot_onsite_iat, - this->occ_mat[iat][target_l][0][0].c, - u_value, diag_coeff, weight_eu, m_size); - } - else // nspin=1 or nspin=2 - { - // spin-up channel - this->energy_u += dftu_pw::compute_pot_onsite_scalar( - pot_onsite_iat, - this->occ_mat[iat][target_l][0][0].c, - u_value, diag_coeff, weight_eu, m_size); - // spin-down channel for nspin=2 - if(this->nspin == 2) - { - std::complex* pot_onsite_iat1 = &(this->pot_uterm_pw[this->pot_uterm_pw.size()/2 + this->pot_uterm_pw_index[iat]]); - this->energy_u += dftu_pw::compute_pot_onsite_scalar( - pot_onsite_iat1, - this->occ_mat[iat][target_l][0][1].c, - u_value, diag_coeff, weight_eu, m_size); - } - } - } -} - -template -void Plus_U_Base::accumulate_occ_one_k(const void* psi_in, - const ModuleBase::matrix& wg_in, - const UnitCell& cell, - const int* isk) -{ - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - const psi::Psi, Device>* psi_p = - (const psi::Psi, Device>*)psi_in; - const int nbands = psi_p->get_nbands(); - const int npol = psi_p->get_npol(); - for(int ik = 0; ik < psi_p->get_nk(); ik++) - { - int is = (this->nspin == 2) ? isk[ik] : 0; - psi_p->fix_k(ik); - onsite_p->tabulate_atomic(ik); - - onsite_p->overlap_proj_psi(nbands*npol, psi_p->get_pointer()); - const std::complex* becp = onsite_p->get_h_becp(); - int nkb = onsite_p->get_size_becp() / nbands / npol; - - int begin_ih = 0; - for(int iat = 0; iat < cell.nat; iat++) - { - const int it = cell.iat2it[iat]; - const int nh = onsite_p->get_nh(iat); - const int target_l = get_orbital_corr(it); - if(!has_correlated_orbital(it)) - { - begin_ih += nh; - continue; - } - const int m_begin = target_l * target_l; - const int tlp1 = 2 * target_l + 1; - if(this->nspin == 4) - { - dftu_pw::accumulate_occ_spinor( - this->occ_mat[iat][target_l][0][0].c, - becp, nbands, npol, nkb, begin_ih, m_begin, tlp1, - wg_in, ik); - } - else // nspin=1 or nspin=2 - { - dftu_pw::accumulate_occ_scalar( - this->occ_mat[iat][target_l][0][is].c, - becp, nbands, nkb, begin_ih, m_begin, tlp1, - wg_in, ik); - } - begin_ih += nh; - } - } -} - -// explicit instantiations -template void Plus_U_Base::accumulate_occ_one_k( - const void*, const ModuleBase::matrix&, const UnitCell&, const int*); -#if defined(__CUDA) || defined(__ROCM) -template void Plus_U_Base::accumulate_occ_one_k( - const void*, const ModuleBase::matrix&, const UnitCell&, const int*); -#endif diff --git a/source/source_pw/module_pwdft/dftu_output.h b/source/source_pw/module_pwdft/dftu_output.h deleted file mode 100644 index 98275d63a0a..00000000000 --- a/source/source_pw/module_pwdft/dftu_output.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef DFTU_OUTPUT_H -#define DFTU_OUTPUT_H - -#include -#include - -class Plus_U_Base; -class UnitCell; - -namespace dftu_io -{ - -/// Output DFT+U information (Hubbard U/J, local occupation matrices) to the -/// running log and, when out_chg is set, to the dm_onsite.txt file. -/// -/// Extracted from Plus_U_Base::output as a free function so that IO logic is -/// decoupled from the Plus_U_Base class. The function only reads the -/// Plus_U_Base state via public accessors; no friend declaration needed. -void output(const Plus_U_Base& dftu, - const UnitCell& ucell, - bool out_chg, - const std::string& global_out_dir, - int nspin, - int npol); - -/// Write local occupation matrices to the given stream. -/// -/// Extracted from Plus_U_Base::write_occup_m. When diag is true, eigenvalues -/// and magnetism are also printed; otherwise only raw matrix elements. -/// Caller is responsible for opening/closing the stream. -void write_occup_m(const Plus_U_Base& dftu, - const UnitCell& ucell, - std::ofstream& ofs, - bool diag, - int nspin, - int npol); - -} // namespace dftu_io - -#endif diff --git a/source/source_pw/module_pwdft/onsite_proj_force_stress.cpp b/source/source_pw/module_pwdft/onsite_proj_force_stress.cpp index 47ef14a8067..e794839ffe8 100644 --- a/source/source_pw/module_pwdft/onsite_proj_force_stress.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_force_stress.cpp @@ -12,7 +12,7 @@ void projectors::OnsiteProjector::cal_force_onsite_dftu(int ik, int n const std::complex* pot_onsite_ptr = dftu.get_pot_uterm_pw_spin(isk_val); const int pot_onsite_size = dftu.get_size_pot_uterm_pw_spin(); this->fs_tools->cal_force_dftu(ik, npm, force, - dftu.get_orbital_corr_data(), pot_onsite_ptr, pot_onsite_size, wg_ik); + dftu.get_orbital_corr_vec().data(), pot_onsite_ptr, pot_onsite_size, wg_ik); } template @@ -24,7 +24,7 @@ double projectors::OnsiteProjector::cal_stress_onsite_dftu(int ik, in const std::complex* pot_onsite_ptr = dftu.get_pot_uterm_pw_spin(isk_val); const int pot_onsite_size = dftu.get_size_pot_uterm_pw_spin(); return this->fs_tools->cal_stress_dftu(ik, npm, - dftu.get_orbital_corr_data(), pot_onsite_ptr, pot_onsite_size, wg_ik); + dftu.get_orbital_corr_vec().data(), pot_onsite_ptr, pot_onsite_size, wg_ik); } template diff --git a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp index 672b7d48b78..216617eaa9a 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp @@ -4,6 +4,13 @@ namespace hamilt { +extern template class OperatorEXXPW, base_device::DEVICE_CPU>; +extern template class OperatorEXXPW, base_device::DEVICE_CPU>; +#if ((defined __CUDA) || (defined __ROCM)) +extern template class OperatorEXXPW, base_device::DEVICE_GPU>; +extern template class OperatorEXXPW, base_device::DEVICE_GPU>; +#endif + template void get_exx_potential(const K_Vectors* kv, const ModulePW::PW_Basis_K* wfcpw, @@ -531,8 +538,6 @@ double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, return exx_div; } -template class OperatorEXXPW, base_device::DEVICE_CPU>; -template class OperatorEXXPW, base_device::DEVICE_CPU>; template void get_exx_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, @@ -576,8 +581,6 @@ template void get_exx_stress_potential(const K_ int, const CoulombParam&); #if ((defined __CUDA) || (defined __ROCM)) -template class OperatorEXXPW, base_device::DEVICE_GPU>; -template class OperatorEXXPW, base_device::DEVICE_GPU>; template void get_exx_potential(const K_Vectors*, const ModulePW::PW_Basis_K*, ModulePW::PW_Basis*, diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index a49a99f4d6c..f330a6c6982 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -1,9 +1,9 @@ #include "source_pw/module_pwdft/setup_dftu_pw.h" #include "source_pw/module_pwdft/dftu_base.h" // mohan add 2025-11-06 -#include "source_pw/module_pwdft/dftu_output.h" // mohan add 2025-11-08 +#include "source_pw/module_pwdft/dftu_base_io.h" // mohan add 2025-11-08 #include "source_io/module_parameter/parameter.h" -namespace pw +namespace DFTU_BASE { void iter_init_dftu_pw(const int iter, @@ -29,7 +29,7 @@ void iter_init_dftu_pw(const int iter, { dftu.cal_occ_pw(psi, wg, ucell, p_chgmix, isk); } - dftu_io::output(dftu, ucell, PARAM.inp.out_chg[0], PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol); + DFTU_BASE::output(dftu, ucell, PARAM.inp.out_chg[0], PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol); } } diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.h b/source/source_pw/module_pwdft/setup_dftu_pw.h index 9eabec7bd9c..c9c4731cc49 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.h +++ b/source/source_pw/module_pwdft/setup_dftu_pw.h @@ -8,7 +8,7 @@ struct Input_para; class Plus_U_Base; // mohan add 2025-11-06 -namespace pw +namespace DFTU_BASE { void iter_init_dftu_pw(const int iter, diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 738664d700f..98e46340f1d 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -124,6 +124,7 @@ void pw::setup_pot(const int istep, dftu.init_base(ucell, PARAM.globalv.npol, inp.nspin, inp.orbital_corr, inp.yukawa_potential, + inp.yukawa_lambda, PARAM.globalv.global_readin_dir, PARAM.globalv.global_out_dir, inp.init_chg, diff --git a/source/source_lcao/module_dftu/dftu_yukawa.cpp b/source/source_pw/module_pwdft/yukawa_screening.cpp similarity index 54% rename from source/source_lcao/module_dftu/dftu_yukawa.cpp rename to source/source_pw/module_pwdft/yukawa_screening.cpp index 1d0f2a2ceb3..9f5892d28fc 100644 --- a/source/source_lcao/module_dftu/dftu_yukawa.cpp +++ b/source/source_pw/module_pwdft/yukawa_screening.cpp @@ -1,26 +1,62 @@ -#ifdef __LCAO +#include "yukawa_screening.h" + #include "source_base/constants.h" -#include "source_base/global_function.h" -#include "dftu_lcao.h" -#include "dftu_yukawa.h" -#include "source_io/module_parameter/parameter.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_quit.h" +#include "source_base/tool_title.h" +#include "source_cell/unitcell.h" +#ifdef __LCAO +#include "source_basis/module_ao/orb_read.h" +#endif +#include #include #include #include - -void DFTU_LCAO::cal_yukawa_lambda(Plus_U& dftu, double** rho, const int& nrxx) +void YukawaScreening::init(const UnitCell& cell, + const std::vector& orbital_corr, + double yukawa_lambda_cfg) { - ModuleBase::TITLE("DFTU_LCAO", "cal_yukawa_lambda"); + this->yukawa_lambda_cfg_ = yukawa_lambda_cfg; + this->lambda_ = 0.0; + this->orbital_corr_ = orbital_corr; + + this->Fk_.resize(cell.ntype); + this->U_Yukawa_.resize(cell.ntype); + this->J_Yukawa_.resize(cell.ntype); - // read from the global PARAM.inp.nspin instead of a Plus_U member; - // the member indirection is being removed during the refactor - const int nspin = PARAM.inp.nspin; + for (int it = 0; it < cell.ntype; it++) + { + const int NL = cell.atoms[it].nwl + 1; + + this->Fk_[it].resize(NL); + this->U_Yukawa_[it].resize(NL); + this->J_Yukawa_[it].resize(NL); + + for (int l = 0; l < NL; l++) + { + const int N = cell.atoms[it].l_nchi[l]; + + this->Fk_[it][l].resize(N); + for (int n = 0; n < N; n++) + { + this->Fk_[it][l][n].resize(l + 1, 0.0); + } + + this->U_Yukawa_[it][l].resize(N, 0.0); + this->J_Yukawa_[it][l].resize(N, 0.0); + } + } +} + +void YukawaScreening::cal_lambda(double** rho, int nrxx, int nspin) +{ + ModuleBase::TITLE("YukawaScreening", "cal_lambda"); - if (dftu.get_yukawa_lambda() > 0) + if (this->yukawa_lambda_cfg_ > 0) { - dftu.set_lambda(dftu.get_yukawa_lambda()); + this->lambda_ = this->yukawa_lambda_cfg_; return; } @@ -36,9 +72,9 @@ void DFTU_LCAO::cal_yukawa_lambda(Plus_U& dftu, double** rho, const int& nrxx) double min_rho = std::numeric_limits::max(); for (int is = 0; is < nspin; is++) { - if(nspin == 4 && is > 0) + if (nspin == 4 && is > 0) { - continue;// for non-collinear spin case, first spin contains the charge density + continue; // for non-collinear spin case, first spin contains the charge density } for (int ir = 0; ir < nrxx; ir++) { @@ -86,80 +122,77 @@ void DFTU_LCAO::cal_yukawa_lambda(Plus_U& dftu, double** rho, const int& nrxx) << " min_rho=" << min_rho_global << " (need finite sum_rho > 0); nspin=" << nspin << " nrxx=" << nrxx; - ModuleBase::WARNING_QUIT("DFTU_LCAO::cal_yukawa_lambda", oss.str()); + ModuleBase::WARNING_QUIT("YukawaScreening::cal_lambda", oss.str()); } - dftu.set_lambda(val2 / val1); + this->lambda_ = val2 / val1; // rescaling - dftu.set_lambda(dftu.get_lambda() / 1.6); - - return; + this->lambda_ /= 1.6; } -void DFTU_LCAO::cal_slater_Fk(Plus_U& dftu, const UnitCell& ucell, const int L, const int T) +void YukawaScreening::cal_slater_Fk(const UnitCell& ucell, int L, int T, const LCAO_Orbitals* orb) { - ModuleBase::TITLE("DFTU_LCAO", "cal_slater_Fk"); + ModuleBase::TITLE("YukawaScreening", "cal_slater_Fk"); + +#ifdef __LCAO + const double lambda_val = this->lambda_; - if (dftu.use_yukawa()) + for (int chi = 0; chi < ucell.atoms[T].l_nchi[L]; chi++) { - const LCAO_Orbitals* orb = dftu.get_ptr_orb(); - const double lambda_val = dftu.get_lambda(); - auto& Fk = dftu.get_Fk_data(); + const int mesh = orb->Phi[T].PhiLN(L, chi).getNr(); - for (int chi = 0; chi < ucell.atoms[T].l_nchi[L]; chi++) + for (int k = 0; k <= L; k++) { - // if(chi!=0) continue; - const int mesh = orb->Phi[T].PhiLN(L, chi).getNr(); - - for (int k = 0; k <= L; k++) + for (int ir0 = 1; ir0 < mesh; ir0++) { - for (int ir0 = 1; ir0 < mesh; ir0++) + double r0 = orb->Phi[T].PhiLN(L, chi).getRadial(ir0); + const double rab0 = orb->Phi[T].PhiLN(L, chi).getRab(ir0); + const double R_L0 = orb->Phi[T].PhiLN(L, chi).getPsi(ir0); + + for (int ir1 = 1; ir1 < mesh; ir1++) { - double r0 = orb->Phi[T].PhiLN(L, chi).getRadial(ir0); - const double rab0 = orb->Phi[T].PhiLN(L, chi).getRab(ir0); - const double R_L0 = orb->Phi[T].PhiLN(L, chi).getPsi(ir0); + double bslval, hnkval; + double r1 = orb->Phi[T].PhiLN(L, chi).getRadial(ir1); + const double rab1 = orb->Phi[T].PhiLN(L, chi).getRab(ir1); + const double R_L1 = orb->Phi[T].PhiLN(L, chi).getPsi(ir1); - for (int ir1 = 1; ir1 < mesh; ir1++) + int l = 2 * k; + if (ir0 < ir1) // less than { - double bslval, hnkval; - double r1 = orb->Phi[T].PhiLN(L, chi).getRadial(ir1); - const double rab1 = orb->Phi[T].PhiLN(L, chi).getRab(ir1); - const double R_L1 = orb->Phi[T].PhiLN(L, chi).getPsi(ir1); - - int l = 2 * k; - if (ir0 < ir1) // less than - { - bslval = DFTU_LCAO::spherical_Bessel(l, r0, lambda_val); - hnkval = DFTU_LCAO::spherical_Hankel(l, r1, lambda_val); - } - else // greater than - { - bslval = DFTU_LCAO::spherical_Bessel(l, r1, lambda_val); - hnkval = DFTU_LCAO::spherical_Hankel(l, r0, lambda_val); - } - Fk[T][L][chi][k] -= (4 * k + 1) * lambda_val * pow(R_L0, 2) * bslval * hnkval * pow(R_L1, 2) - * pow(r0, 2) * pow(r1, 2) * rab0 * rab1; + bslval = spherical_Bessel(l, r0, lambda_val); + hnkval = spherical_Hankel(l, r1, lambda_val); } + else // greater than + { + bslval = spherical_Bessel(l, r1, lambda_val); + hnkval = spherical_Hankel(l, r0, lambda_val); + } + this->Fk_[T][L][chi][k] -= (4 * k + 1) * lambda_val * pow(R_L0, 2) * bslval * hnkval + * pow(R_L1, 2) * pow(r0, 2) * pow(r1, 2) * rab0 * rab1; } } } } - - return; +#else + (void)ucell; + (void)L; + (void)T; + (void)orb; + ModuleBase::WARNING_QUIT("YukawaScreening::cal_slater_Fk", + "Slater integrals require numerical orbitals; compile with __LCAO"); +#endif } -void DFTU_LCAO::cal_slater_UJ(Plus_U& dftu, const UnitCell& ucell, double** rho, const int& nrxx) +void YukawaScreening::cal_slater_UJ(const UnitCell& ucell, + double** rho, + int nrxx, + int nspin, + const LCAO_Orbitals* orb) { - ModuleBase::TITLE("DFTU_LCAO", "cal_slater_UJ"); - if (!dftu.use_yukawa()) - { - return; - } - - cal_yukawa_lambda(dftu, rho, nrxx); + ModuleBase::TITLE("YukawaScreening", "cal_slater_UJ"); - auto& Fk = dftu.get_Fk_data(); + this->cal_lambda(rho, nrxx, nspin); for (int it = 0; it < ucell.ntype; it++) { @@ -170,7 +203,7 @@ void DFTU_LCAO::cal_slater_UJ(Plus_U& dftu, const UnitCell& ucell, double** rho, int N = ucell.atoms[it].l_nchi[l]; for (int n = 0; n < N; n++) { - ModuleBase::GlobalFunc::ZEROS(ModuleBase::GlobalFunc::VECTOR_TO_PTR(Fk[it][l][n]), l + 1); + std::fill(this->Fk_[it][l][n].begin(), this->Fk_[it][l][n].end(), 0.0); } } } @@ -181,52 +214,43 @@ void DFTU_LCAO::cal_slater_UJ(Plus_U& dftu, const UnitCell& ucell, double** rho, for (int L = 0; L < NL; L++) { - const int N = ucell.atoms[T].l_nchi[L]; - - if (L >= dftu.get_orbital_corr(T) && dftu.get_orbital_corr(T) != -1) + if (L >= this->orbital_corr_[T] && this->orbital_corr_[T] != -1) { - if (L != dftu.get_orbital_corr(T)) + if (L != this->orbital_corr_[T]) { continue; } - cal_slater_Fk(dftu, ucell, L, T); - + this->cal_slater_Fk(ucell, L, T, orb); - if( L == 1) + if (L == 1) { - dftu.set_U_Yukawa(T, L, 0, Fk[T][L][0][0]); - dftu.set_J_Yukawa(T, L, 0, Fk[T][L][0][1] / 5.0); + this->U_Yukawa_[T][L][0] = this->Fk_[T][L][0][0]; + this->J_Yukawa_[T][L][0] = this->Fk_[T][L][0][1] / 5.0; } - else if( L == 2) + else if (L == 2) { - dftu.set_U_Yukawa(T, L, 0, Fk[T][L][0][0]); - dftu.set_J_Yukawa(T, L, 0, (Fk[T][L][0][1] + Fk[T][L][0][2]) / 14.0); + this->U_Yukawa_[T][L][0] = this->Fk_[T][L][0][0]; + this->J_Yukawa_[T][L][0] = (this->Fk_[T][L][0][1] + this->Fk_[T][L][0][2]) / 14.0; } - else if( L == 3) + else if (L == 3) { - dftu.set_U_Yukawa(T, L, 0, Fk[T][L][0][0]); - dftu.set_J_Yukawa(T, L, 0, (286.0 * Fk[T][L][0][1] + 195.0 * Fk[T][L][0][2] - + 250.0 * Fk[T][L][0][3]) - / 6435.0); + this->U_Yukawa_[T][L][0] = this->Fk_[T][L][0][0]; + this->J_Yukawa_[T][L][0] = (286.0 * this->Fk_[T][L][0][1] + + 195.0 * this->Fk_[T][L][0][2] + + 250.0 * this->Fk_[T][L][0][3]) / 6435.0; } // Hartree to Rydeberg - dftu.set_U_Yukawa(T, L, 0, dftu.get_U_Yukawa(T, L, 0) * 2.0); - dftu.set_J_Yukawa(T, L, 0, dftu.get_J_Yukawa(T, L, 0) * 2.0); - // update current U with calculated U-J from Slater integrals - dftu.set_u_current(T, dftu.get_U_Yukawa(T, L, 0) - dftu.get_J_Yukawa(T, L, 0)); + this->U_Yukawa_[T][L][0] *= 2.0; + this->J_Yukawa_[T][L][0] *= 2.0; } // end if } // end L } // end T - - return; } -double DFTU_LCAO::spherical_Bessel(const int k, const double r, const double lambda) +double YukawaScreening::spherical_Bessel(const int k, const double r, const double lambda) { - ModuleBase::TITLE("DFTU_LCAO", "spherical_Bessel"); - - double val=0.0; + double val = 0.0; double x = r * lambda; if (k == 0) { @@ -277,11 +301,9 @@ double DFTU_LCAO::spherical_Bessel(const int k, const double r, const double lam return val; } -double DFTU_LCAO::spherical_Hankel(const int k, const double r, const double lambda) +double YukawaScreening::spherical_Hankel(const int k, const double r, const double lambda) { - ModuleBase::TITLE("DFTU_LCAO", "spherical_Hankel"); - - double val=0.0; + double val = 0.0; double x = r * lambda; if (k == 0) { @@ -334,5 +356,3 @@ double DFTU_LCAO::spherical_Hankel(const int k, const double r, const double lam } return val; } - -#endif diff --git a/source/source_pw/module_pwdft/yukawa_screening.h b/source/source_pw/module_pwdft/yukawa_screening.h new file mode 100644 index 00000000000..1cf8e32cfdc --- /dev/null +++ b/source/source_pw/module_pwdft/yukawa_screening.h @@ -0,0 +1,72 @@ +#ifndef YUKAWA_SCREENING_H +#define YUKAWA_SCREENING_H + +#include + +class UnitCell; +class LCAO_Orbitals; + +/** + * @brief Yukawa-screened DFT+U: self-consistent U/J from Slater integrals. + * + * Encapsulates the Yukawa screening length (lambda), the Slater integrals Fk, + * and the derived U_Yukawa / J_Yukawa values. All Yukawa-related state that + * used to live on Plus_U_Base / Plus_U is owned here, so the DFT+U classes + * only hold an instance of this class when Yukawa screening is enabled. + * + * Currently only the LCAO path drives the calculation (it needs the radial + * orbitals from LCAO_Orbitals), but the class itself has no LCAO-only data. + */ +class YukawaScreening +{ + public: + YukawaScreening() = default; + ~YukawaScreening() = default; + + /// allocate Fk / U_Yukawa / J_Yukawa according to the cell and record the + /// user-provided screening length (yukawa_lambda_cfg > 0 means fixed). + void init(const UnitCell& cell, + const std::vector& orbital_corr, + double yukawa_lambda_cfg); + + /// determine lambda: use the fixed config value when positive, otherwise + /// estimate from the charge density (Thomas-Fermi-like) and rescale by 1.6. + void cal_lambda(double** rho, int nrxx, int nspin); + + /// compute Slater integrals Fk for the correlated orbital of atom type T. + void cal_slater_Fk(const UnitCell& ucell, int L, int T, const LCAO_Orbitals* orb); + + /// drive cal_lambda + cal_slater_Fk over all correlated orbitals and derive + /// U_Yukawa / J_Yukawa. Returns via get_U/get_J; u_current of the owning + /// DFT+U object is updated by the caller. + void cal_slater_UJ(const UnitCell& ucell, + double** rho, + int nrxx, + int nspin, + const LCAO_Orbitals* orb); + + double get_lambda() const { return lambda_; } + double get_U(int it, int l, int n) const { return U_Yukawa_[it][l][n]; } + double get_J(int it, int l, int n) const { return J_Yukawa_[it][l][n]; } + /// effective U-J of the correlated orbital (n = 0) for atom type it + double get_Ueff(int it) const + { + const int l = orbital_corr_[it]; + return U_Yukawa_[it][l][0] - J_Yukawa_[it][l][0]; + } + + private: + /// spherical modified Bessel function of the first kind, orders 0/2/4/6 + static double spherical_Bessel(int k, double r, double lambda); + /// spherical modified Hankel function of the second kind, orders 0/2/4/6 + static double spherical_Hankel(int k, double r, double lambda); + + double lambda_ = 0.0; + double yukawa_lambda_cfg_ = 0.0; + std::vector orbital_corr_; + std::vector>>> Fk_; + std::vector>> U_Yukawa_; + std::vector>> J_Yukawa_; +}; + +#endif diff --git a/source/source_pw/module_stodft/sto_dos.cpp b/source/source_pw/module_stodft/sto_dos.cpp index 245a23f8fd4..884679a8b83 100644 --- a/source/source_pw/module_stodft/sto_dos.cpp +++ b/source/source_pw/module_stodft/sto_dos.cpp @@ -1,5 +1,6 @@ #include "sto_dos.h" +#include "source_base/parallel_comm.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_title.h" diff --git a/source/source_pw/module_stodft/sto_elecond.cpp b/source/source_pw/module_stodft/sto_elecond.cpp index 0e9e2aa97bd..ba60a6114a1 100644 --- a/source/source_pw/module_stodft/sto_elecond.cpp +++ b/source/source_pw/module_stodft/sto_elecond.cpp @@ -5,6 +5,7 @@ #include "source_base/memory_recorder.h" #include "source_base/module_container/ATen/tensor.h" #include "source_base/parallel_device.h" +#include "source_base/parallel_comm.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/vector3.h" diff --git a/source/source_pw/module_stodft/sto_iter.cpp b/source/source_pw/module_stodft/sto_iter.cpp index dcc2e187a78..fc1f28d6e4b 100644 --- a/source/source_pw/module_stodft/sto_iter.cpp +++ b/source/source_pw/module_stodft/sto_iter.cpp @@ -2,6 +2,7 @@ #include "source_base/kernels/math_kernel_op.h" #include "source_base/para_gemm.h" +#include "source_base/parallel_comm.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" diff --git a/source/source_pw/module_stodft/sto_tool.cpp b/source/source_pw/module_stodft/sto_tool.cpp index d21851fa9a3..4c432a687c6 100644 --- a/source/source_pw/module_stodft/sto_tool.cpp +++ b/source/source_pw/module_stodft/sto_tool.cpp @@ -1,6 +1,7 @@ #include "sto_tool.h" #include "source_base/math_chebyshev.h" +#include "source_base/parallel_comm.h" #include "source_base/parallel_device.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h"