diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index d5dc5f18ddf..3e764e161f7 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.h" // mohan add 2025-11-07 +#include "source_lcao/module_dftu/dftu_lcao.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 c0b69a87984..a86b205155c 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.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 #include "source_pw/module_pwdft/vnl_pw.h" namespace ModuleESolver diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 9f958d68a79..eaa03faf030 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -506,6 +506,14 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // 1) calculate the local occupation number matrix and energy correction in DFT+U finish_dftu_lcao(iter, conv_esolver, this->inp_->dft_plus_u, this->inp_->out_chg[0], &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, this->inp_->nspin, PARAM.globalv.npol); + // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState. + // Covers both dft_plus_u==1 (new method, energy accumulated by DFTU::contributeHR + // via cal_v_of_u) and dft_plus_u==2 (old method, energy from cal_energy_correction). + if (this->inp_->dft_plus_u) + { + this->pelec->set_dftu_energy(this->dftu.get_energy()); + } + // 2) for deepks, calculate delta_e, output labels during electronic steps this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, *this->inp_); diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index f2fccf7e3c3..cf131ef96c1 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -209,6 +209,12 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const ucell, this->p_chgmix, this->kv.isk.data()); + + // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState + if (this->inp_->dft_plus_u) + { + this->pelec->set_dftu_energy(this->dftu.get_energy()); + } } // Temporary, it should be replaced by hsolver later. diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 961c0302e5c..e4aa10e4da7 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -2,7 +2,7 @@ #include "source_cell/cal_ux.h" #include "source_estate/module_charge/symm_rho.h" #include "source_lcao/hamilt_lcao.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/module_dftu/dftu_lcao.h" #include "source_hamilt/module_gint/gint.h" #include "source_base/formatter.h" #include "source_base/timer.h" diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index 228d5248368..82ceb15d4b0 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -131,6 +131,7 @@ class ElecState } double get_dftu_energy(); + void set_dftu_energy(double e) { dftu_energy_ = e; } double get_local_pp_energy(); fenergy f_en; ///< energies contribute to the total free energy @@ -148,6 +149,9 @@ class ElecState public: bool skip_weights = false; + + private: + double dftu_energy_ = 0.0; ///< DFT+U energy, set by ESolver via set_dftu_energy() }; /** diff --git a/source/source_estate/estate_e_terms.cpp b/source/source_estate/estate_e_terms.cpp index f5514a6b8ed..997c8fa2c98 100644 --- a/source/source_estate/estate_e_terms.cpp +++ b/source/source_estate/estate_e_terms.cpp @@ -3,7 +3,6 @@ #include "source_estate/module_pot/h_hartree_pw.h" #include "source_estate/module_pot/efield.h" #include "source_estate/module_pot/gatefield.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 namespace elecstate { @@ -35,7 +34,7 @@ double ElecState::get_solvent_model_Acav() double ElecState::get_dftu_energy() { - return Plus_U::get_energy(); + return dftu_energy_; } double ElecState::get_local_pp_energy() diff --git a/source/source_estate/module_charge/chgmixing.cpp b/source/source_estate/module_charge/chgmixing.cpp index 1fd48fac5d3..b15f28a8d7d 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.h" +#include "source_lcao/module_dftu/dftu_lcao.h" #include "source_lcao/module_deltaspin/spin_constrain.h" void module_charge::chgmixing_ks(const int iter, // scf iteration number @@ -119,7 +119,7 @@ void module_charge::chgmixing_ks(const int iter, // scf iteration number void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number Charge_Mixing* p_chgmix, // charge mixing class - Plus_U &dftu, // mohan add 2025-11-06 + Plus_U_Base& dftu, const Input_para& inp) // input parameters { ModuleBase::TITLE("module_charge", "chgmixing_ks_pw"); @@ -167,7 +167,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number std::cout << " U-Ramping! Current U = "; for (int i = 0; i < dftu.get_num_u_types(); i++) { - std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } @@ -180,7 +180,7 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number Charge_Mixing* p_chgmix, // charge mixing class - Plus_U &dftu, // mohan add 2025-11-06 + Plus_U_Base& dftu, const int nnr, // dimension of density matrix const Input_para& inp) // input parameters { @@ -202,7 +202,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number std::cout << " U-Ramping! Current U = "; for (int i = 0; i < dftu.get_num_u_types(); i++) { - std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } @@ -221,7 +221,7 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number std::cout << " U-Ramping! Current U = "; for (int i = 0; i < dftu.get_num_u_types(); i++) { - std::cout << dftu.get_hubbard_u(i) * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.get_u_current(i) * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } diff --git a/source/source_estate/module_charge/chgmixing.h b/source/source_estate/module_charge/chgmixing.h index 4a04a0880d0..6b751128e36 100644 --- a/source/source_estate/module_charge/chgmixing.h +++ b/source/source_estate/module_charge/chgmixing.h @@ -6,7 +6,7 @@ #include "source_estate/module_charge/charge_mixing.h" // use p_chgmix #include "source_io/module_parameter/input_parameter.h" // use Input_para #include "source_cell/unitcell.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 +#include "source_pw/module_pwdft/dftu_base.h" namespace module_charge { @@ -28,12 +28,12 @@ void chgmixing_ks(const int iter, // scf iteration number void chgmixing_ks_pw(const int iter, Charge_Mixing* p_chgmix, - Plus_U &dftu, // mohan add 2025-11-06 + Plus_U_Base& dftu, const Input_para& inp); // input parameters void chgmixing_ks_lcao(const int iter, // scf iteration number Charge_Mixing* p_chgmix, // charge mixing class - Plus_U &dftu, // mohan add 2025-11-06 + Plus_U_Base& dftu, const int nnr, // dimension of density matrix const Input_para& inp); // input parameters diff --git a/source/source_hsolver/test/diago_mock.h b/source/source_hsolver/test/diago_mock.h index cbab29db8d6..279b5dc566e 100644 --- a/source/source_hsolver/test/diago_mock.h +++ b/source/source_hsolver/test/diago_mock.h @@ -577,7 +577,7 @@ template<> hamilt::HamiltPW::HamiltPW( ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl* ppcell, - Plus_U* p_dftu, // mohan add 20251108 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info) { @@ -599,7 +599,7 @@ template<> hamilt::HamiltPW>::HamiltPW( ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl* ppcell, - Plus_U* p_dftu, // mohan add 20251108 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info) { @@ -621,7 +621,7 @@ template<> hamilt::HamiltPW>::HamiltPW( ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl* ppcell, - Plus_U* p_dftu, // mohan add 20251108 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info) { diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.h b/source/source_io/module_ctrl/ctrl_scf_lcao.h index 53749f3c97b..b658253f12e 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.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.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 cfea51ca73e..0ce5e9ccf2c 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.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.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 00a4466b56a..e0a0eab995a 100644 --- a/source/source_io/module_hs/write_hs_r.h +++ b/source/source_io/module_hs/write_hs_r.h @@ -6,7 +6,7 @@ #include "source_cell/klist.h" #include "source_hamilt/hamilt.h" #include "source_lcao/lcao_hs_arrays.hpp" -#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.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 b8368903f09..37b76864a55 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/op_dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_lcao_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 b2e91a7437b..1c31008a2ea 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/op_dftu_lcao.h" +#include "source_lcao/module_dftu/dftu_lcao_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/module_parameter/input_conv.cpp b/source/source_io/module_parameter/input_conv.cpp index a0eecbe45a1..adeac3353c1 100644 --- a/source/source_io/module_parameter/input_conv.cpp +++ b/source/source_io/module_parameter/input_conv.cpp @@ -20,7 +20,6 @@ #include "source_lcao/module_ri/exx_abfs_jle.h" #endif -#include "source_lcao/module_dftu/dftu.h" #ifdef __LCAO #include "source_basis/module_ao/orb_read.h" #include "source_lcao/force_stress_lcao.h" @@ -93,19 +92,9 @@ void Input_Conv::Convert() // iteration (1/3) //---------------------------------------------------------- - if (PARAM.inp.dft_plus_u) - { - Plus_U::omc = PARAM.inp.omc; - Plus_U::orbital_corr = PARAM.inp.orbital_corr; - Plus_U::uramping = PARAM.globalv.uramping; - Plus_U::mixing_dftu = PARAM.inp.mixing_dftu; - Plus_U::U = PARAM.globalv.hubbard_u; - Plus_U::U0 = PARAM.globalv.hubbard_u; - if (PARAM.globalv.uramping > 0.01) - { - ModuleBase::GlobalFunc::ZEROS(Plus_U::U.data(), PARAM.inp.ntype); - } - } + // Note: DFT+U static members (u_current, u_target, occ_mat_ctrl, etc.) + // are now initialized inside Plus_U_Base::init_base() which is called + // from setup_pot.cpp (PW) and lcao_set.cpp (LCAO). Mohan refactor 2025-11. //---------------------------------------------------------- // Yu Liu add 2022-05-18 diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 196b11fad9f..c5890b4dc3b 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -626,7 +626,7 @@ struct Input_para bool yukawa_potential = false; ///< default: false double yukawa_lambda = -1.0; ///< default: -1.0, which means we calculate lambda double uramping_eV = -1.0; ///< U-Ramping method (eV) - int omc = 0; ///< the mode of occupation matrix control + int occ_mat_ctrl = 0; ///< the mode of occupation matrix control double onsite_radius = 0.0; ///< radius of the sphere for onsite projection (Bohr) std::vector hubbard_u_eV = {}; ///< Hubbard Coulomb interaction parameter U(ev) std::vector orbital_corr = {}; ///< which correlated orbitals need corrected ; d:2 ,f:3, do not diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index ad43f1b966c..13515d403ee 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -818,7 +818,7 @@ void ReadInput::item_dftu() [NOTE] The easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. The format of the file is rather straight-forward.)"; item.default_value = "0"; item.unit = ""; - read_sync_int(input.omc); + read_sync_int(input.occ_mat_ctrl); this->add_item(item); } { diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h index 95e9b80a1df..f7f1c208444 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.h" +#include "source_lcao/module_dftu/dftu_lcao.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_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index a8678774d89..239f1650df6 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -358,7 +358,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_FALSE(param.inp.yukawa_potential); EXPECT_DOUBLE_EQ(param.inp.yukawa_lambda, -1.0); EXPECT_EQ(param.inp.onsite_radius, 0.0); - EXPECT_EQ(param.inp.omc, 0); + EXPECT_EQ(param.inp.occ_mat_ctrl, 0); EXPECT_FALSE(param.inp.dft_plus_dmft); EXPECT_FALSE(param.inp.rpa); EXPECT_FALSE(param.inp.rpa_out_vel); diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index f2b4a04435b..ae5635f5ce8 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -1,7 +1,7 @@ #include "force_stress_lcao.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_dftu/dftu.h" //Quxin add for DFT+U on 20201029 +#include "source_lcao/module_dftu/dftu_lcao.h" //Quxin add for DFT+U on 20201029 #include "source_io/module_output/output_log.h" #include "source_io/module_parameter/parameter.h" // new @@ -20,7 +20,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.h" +#include "source_lcao/module_dftu/dftu_lcao_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" @@ -426,7 +426,32 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, } if (PARAM.inp.dft_plus_u == 2) { - // Old DFT+U implementation (dft_plus_u==2) still needs ForceStressArrays + // The legacy dft_plus_u==2 force/stress path is currently broken. + // + // Background: Plus_U::force_stress relies on ForceStressArrays + // members DSloc_x/y/z (gamma_only) or DSloc_Rx/Ry/Rz (multik) + // and DH_r being pre-allocated and filled with dS/dR data by the + // main force flow (formerly ForceLcaoGamma::ftable). The DFT+U + // step 2 refactor (commit 70c54c9d5a, 2026-01-23) removed the + // main-flow ForceStressArrays because the operator-based force + // calculation no longer needs it, but the legacy dft_plus_u==2 + // path still depends on it. The local fsr_dftu below is declared + // without allocating those arrays, so any call into + // cal_force_gamma / cal_stress_gamma / folding_matrix_k would + // pass nullptr to pdgemm_ and crash with SIGSEGV. + // + // Until the legacy path is restored or re-implemented, we + // explicitly reject dft_plus_u==2 with cal_force or cal_stress + // enabled. SCF-only runs (no force/stress) are unaffected + // because the energy is computed in cal_energy_correction, + // which does not touch DSloc arrays. Use dft_plus_u=1 for + // force/stress calculations. + if (isforce || isstress) + { + ModuleBase::WARNING_QUIT("Force_Stress_LCAO::getForceStress", + "dft_plus_u==2 with cal_force or cal_stress is currently broken; " + "please use dft_plus_u=1 instead. See notes in source/source_lcao/force_stress_lcao.cpp."); + } ForceStressArrays fsr_dftu; std::vector>* dmk_d = nullptr; std::vector>>* dmk_c = nullptr; diff --git a/source/source_lcao/force_stress_lcao.h b/source/source_lcao/force_stress_lcao.h index af89b500bad..02f62090547 100644 --- a/source/source_lcao/force_stress_lcao.h +++ b/source/source_lcao/force_stress_lcao.h @@ -16,7 +16,7 @@ #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.h" // mohan add 2025-11-07 +#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 2025-11-07 #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 8d9fcbbb62c..16ef0547aaa 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.h" +#include "source_lcao/module_dftu/dftu_lcao.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.h" +#include "module_dftu/dftu_lcao_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/op_dftu_lcao.h" +#include "module_dftu/dftu_lcao_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 375fa864acd..529b998c2e0 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.h (moved to .cpp) +// Plus_U forward declaration, full definition in module_dftu/dftu_lcao.h (moved to .cpp) // mohan add 20260605 class Plus_U; diff --git a/source/source_lcao/lcao_set.cpp b/source/source_lcao/lcao_set.cpp index dda20f80580..e2996b070b1 100644 --- a/source/source_lcao/lcao_set.cpp +++ b/source/source_lcao/lcao_set.cpp @@ -94,6 +94,10 @@ void LCAO_domain::set_pot( inp.cal_stress, inp.device, inp.kpar, + PARAM.globalv.hubbard_u, + PARAM.globalv.uramping, + inp.occ_mat_ctrl, + inp.mixing_dftu, &orb); } diff --git a/source/source_lcao/lcao_set.h b/source/source_lcao/lcao_set.h index f3defdb9823..6b9c5d81da6 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.h" +#include "source_lcao/module_dftu/dftu_lcao.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 bc5af4c68e8..343bc2077b2 100644 --- a/source/source_lcao/module_dftu/CMakeLists.txt +++ b/source/source_lcao/module_dftu/CMakeLists.txt @@ -1,5 +1,5 @@ list(APPEND objects - dftu.cpp + dftu_lcao.cpp dftu_force.cpp dftu_yukawa.cpp dftu_folding.cpp @@ -10,9 +10,9 @@ list(APPEND objects if(ENABLE_LCAO) list(APPEND objects - dftu_lcao.cpp + dftu_lcao_op.cpp dftu_fs.cpp - op_dftu_lcao.cpp + dftu_lcao_op_legacy.cpp ) endif() diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp deleted file mode 100644 index f387f09409c..00000000000 --- a/source/source_lcao/module_dftu/dftu.cpp +++ /dev/null @@ -1,317 +0,0 @@ -#include "dftu.h" - -#include "source_io/module_parameter/parameter.h" -#include "source_base/constants.h" -#include "source_base/global_function.h" -#include "source_base/inverse_matrix.h" -#include "source_base/memory_recorder.h" -#include "source_base/timer.h" -#include "source_cell/magnetism.h" -#include "source_estate/module_charge/charge.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - // mohan add 2025-11-06 -// Static member definitions moved to dftu_base.cpp (Plus_U_Base::) -// Plus_U inherits these from Plus_U_Base. - -Plus_U::Plus_U() -{} - -Plus_U::~Plus_U() -{} - -void Plus_U::init(UnitCell& cell, - const Parallel_Orbitals* pv, - const int npol, - 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, - const int nlocal, - const bool gamma_only_local, - const std::string& ks_solver, - const bool cal_force, - const bool cal_stress, - const std::string& device, - const int kpar -#ifdef __LCAO - , const LCAO_Orbitals* orb -#endif - ) -{ - ModuleBase::TITLE("Plus_U", "init"); - - this->paraV = pv; - -#ifdef __LCAO - ptr_orb_ = orb; - if(ptr_orb_ != nullptr) - { - orb_cutoff_ = orb->cutoffs(); - } - ucell = &cell; -#endif - - if (pv != nullptr) - { - const int global_rows = pv->get_global_row_size(); - const int global_cols = pv->get_global_col_size(); - if (global_rows != global_cols) - { - ModuleBase::WARNING_QUIT("Plus_U::init", "Global row and column dimensions do not match"); - } - if (nlocal != global_rows) - { - ModuleBase::WARNING_QUIT("Plus_U::init", "nlocal does not match global matrix dimension"); - } - } - - this->init_base(cell, - npol, - nspin, - orbital_corr, - yukawa_potential, - yukawa_lambda, - global_readin_dir, - global_out_dir, - init_chg, - nlocal, - gamma_only_local, - ks_solver, - cal_force, - cal_stress, - device, - kpar); - return; -} - -#ifdef __LCAO - -void Plus_U::cal_energy_correction(const UnitCell& ucell, - const int istep) -{ - ModuleBase::TITLE("Plus_U", "cal_energy_correction"); - ModuleBase::timer::start("Plus_U", "cal_energy_correction"); - if (!is_occ_mat_initialized()) - { - ModuleBase::timer::end("Plus_U", "cal_energy_correction"); - return; - } - - // mohan update 20251106 - Plus_U::energy_u = 0.0; - - double energy_dc = 0.0; - - for (int T = 0; T < ucell.ntype; T++) - { - const int NL = ucell.atoms[T].nwl + 1; - const int LC = get_orbital_corr(T); - for (int I = 0; I < ucell.atoms[T].na; I++) - { - if (LC == -1) - { - continue; - } - - const int iat = ucell.itia2iat(T, I); - const int L = get_orbital_corr(T); - - for (int l = 0; l < NL; l++) - { - if (l != get_orbital_corr(T)) - { - continue; - } - - const int N = ucell.atoms[T].l_nchi[l]; - - const int m_tot = 2 * l + 1; - - // part 1: calculate the DFT+U energy correction - for (int n = 0; n < N; n++) - { - if (n != 0) - { - continue; - } - - if (Plus_U::nspin == 1 || Plus_U::nspin == 2) - { - for (int spin = 0; spin < 2; spin++) - { - double nm_trace = 0.0; - double nm2_trace = 0.0; - - for (int m0 = 0; m0 < 2 * l + 1; m0++) - { - nm_trace += this->occ_mat[iat][l][n][spin](m0, m0); - for (int m1 = 0; m1 < 2 * l + 1; m1++) - { - nm2_trace += this->occ_mat[iat][l][n][spin](m0, m1) - * this->occ_mat[iat][l][n][spin](m1, m0); - } - } - if (Yukawa) - { - Plus_U::energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) - * (nm_trace - nm2_trace); - } - else - { - Plus_U::energy_u += 0.5 * this->U[T] * (nm_trace - nm2_trace); - } - } - } - else if (Plus_U::nspin == 4) - { - double nm_trace = 0.0; - double nm2_trace = 0.0; - - for (int m0 = 0; m0 < 2 * l + 1; m0++) - { - for (int ipol0 = 0; ipol0 < this->npol; ipol0++) - { - const int m0_all = m0 + (2 * l + 1) * ipol0; - nm_trace += this->occ_mat[iat][l][n][0](m0_all, m0_all); - - for (int m1 = 0; m1 < 2 * l + 1; m1++) - { - for (int ipol1 = 0; ipol1 < this->npol; ipol1++) - { - int m1_all = m1 + (2 * l + 1) * ipol1; - - nm2_trace += this->occ_mat[iat][l][n][0](m0_all, m1_all) - * this->occ_mat[iat][l][n][0](m1_all, m0_all); - } - } - } - } - if (Yukawa) - { - Plus_U::energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) - * (nm_trace - nm2_trace); - } - else - { - Plus_U::energy_u += 0.5 * this->U[T] * (nm_trace - nm2_trace); - } - } - - for (int m1 = 0; m1 < 2 * l + 1; m1++) - { - for (int ipol1 = 0; ipol1 < this->npol; ipol1++) - { - const int m1_all = m1 + ipol1 * (2 * l + 1); - for (int m2 = 0; m2 < 2 * l + 1; m2++) - { - for (int ipol2 = 0; ipol2 < this->npol; ipol2++) - { - const int m2_all = m2 + ipol2 * (2 * l + 1); - - if (Plus_U::nspin == 1 || Plus_U::nspin == 2) - { - for (int is = 0; is < 2; is++) - { - double VU = 0.0; - VU = get_onebody_eff_pot(T, iat, l, n, is, m1_all, m2_all, false); - energy_dc += VU * this->occ_mat[iat][l][n][is](m1_all, m2_all); - } - } - else if (Plus_U::nspin == 4) - { - double VU = 0.0; - VU = get_onebody_eff_pot(T, iat, l, n, 0, m1_all, m2_all, false); - energy_dc += VU * this->occ_mat[iat][l][n][0](m1_all, m2_all); - } - } - } - } - } - } // end n - } // end L - } // end I - } // end T - - // substract the double counting energy_dc included in band energy eband - Plus_U::energy_u -= energy_dc; - - ModuleBase::timer::end("Plus_U", "cal_energy_correction"); - return; -} - -#endif - -// uramping_update() and u_converged() are now implemented in -// dftu_base.cpp as Plus_U_Base methods (inherited by Plus_U). - -#ifdef __LCAO - -void Plus_U::set_dmr(const elecstate::DensityMatrix, double>* dmr) -{ - this->dm_in_dftu_cd = dmr; - return; -} - -void Plus_U::set_dmr(const elecstate::DensityMatrix* dmr) -{ - this->dm_in_dftu_d = dmr; - return; -} - -const hamilt::HContainer* Plus_U::get_dmr(int ispin) const -{ - if (this->dm_in_dftu_d != nullptr) - { - return this->dm_in_dftu_d->get_DMR_pointer(ispin + 1); - } - else if (this->dm_in_dftu_cd != nullptr) - { - return this->dm_in_dftu_cd->get_DMR_pointer(ispin + 1); - } - else - { - return nullptr; - } -} - -//! dftu occupation matrix for gamma only using dm(double) -template <> -void dftu_cal_occup_m(const int iter, - const UnitCell& ucell, - const std::vector>& dm, - const K_Vectors& kv, - const double& mixing_beta, - hamilt::Hamilt* p_ham, - Plus_U &dftu) -{ - dftu.cal_occup_m_gamma(iter, ucell ,dm, mixing_beta, p_ham); -} - -//! dftu occupation matrix for multiple k-points using dm(complex) -template <> -void dftu_cal_occup_m(const int iter, - const UnitCell& ucell, - const std::vector>>& dm, - const K_Vectors& kv, - const double& mixing_beta, - hamilt::Hamilt>* p_ham, - Plus_U &dftu) -{ - dftu.cal_occup_m_k(iter,ucell, dm, kv, mixing_beta, p_ham); -} - -#endif diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h deleted file mode 100644 index a4c37798ef5..00000000000 --- a/source/source_lcao/module_dftu/dftu.h +++ /dev/null @@ -1,268 +0,0 @@ -#ifndef DFTU_H -#define DFTU_H - -#include "source_cell/klist.h" -#include "source_cell/unitcell.h" -#include "source_basis/module_ao/parallel_orbitals.h" -#include "source_estate/module_charge/charge_mixing.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" -#include "source_lcao/force_stress_arrays.h" // mohan add 2024-06-15 -#endif - -#include -#include - - -class Plus_U : public Plus_U_Base -{ - - public: - Plus_U(); - ~Plus_U(); - - public: - // allocate relevant data strcutures - void init(UnitCell& cell, - const Parallel_Orbitals* pv, - const int npol, - 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, - const int nlocal, - const bool gamma_only_local, - const std::string& ks_solver, - const bool cal_force, - const bool cal_stress, - const std::string& device, - const int kpar -#ifdef __LCAO - , const LCAO_Orbitals* orb = nullptr -#endif - ); - - // calculate the energy correction - void cal_energy_correction(const UnitCell& ucell, const int istep); - - private: - - const Parallel_Orbitals* paraV = nullptr; - -#ifdef __LCAO - const LCAO_Orbitals* ptr_orb_ = nullptr; - std::vector orb_cutoff_; -#endif - -#ifdef __LCAO - //============================================================= - // In dftu_hamilt.cpp - // For calculating contribution to Hamiltonian matrices - //============================================================= - public: - void cal_eff_pot_mat_complex(const int ik, - std::complex* eff_pot, - const std::vector& isk, - const std::complex* sk, - const int npol); - - void cal_eff_pot_mat_real(const int ik, - double* eff_pot, - const std::vector& isk, - const double* sk, - const int npol); - - void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, const int npol); - - void cal_eff_pot_mat_R_complex_double(const int ispin, - std::complex* SR, - std::complex* HR, - const int npol); -#endif - -#ifdef __LCAO - // calculate the local occupation number matrix - void cal_occup_m_k(const int iter, - const UnitCell& ucell, - const std::vector>>& dm_k, - const K_Vectors& kv, - const double& mixing_beta, - hamilt::Hamilt>* p_ham); - - void cal_occup_m_gamma(const int iter, - const UnitCell& ucell, - const std::vector>& dm_gamma, - const double& mixing_beta, - hamilt::Hamilt* p_ham); -#endif - -#ifdef __LCAO -private: - //============================================================= - // In dftu_tools.cpp - // For calculating onsite potential, which is used - // for both Hamiltonian and force/stress - //============================================================= - - void cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU, const int npol); - void cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU, const int npol); - - double get_onebody_eff_pot(const int T, - const int iat, - const int L, - const int N, - const int spin, - const int m0, - const int m1, - const bool newlocale); - - //============================================================= - // In dftu_folding.cpp - // Subroutines for folding S and dS matrix - //============================================================= - - void fold_dSR_gamma(const UnitCell& ucell, - const Parallel_Orbitals& pv, - const Grid_Driver* gd, - double* dsloc_x, - double* dsloc_y, - double* dsloc_z, - double* dh_r, - const int dim1, - const int dim2, - double* dSR_gamma); - - // dim = 0 : S, for Hamiltonian - // dim = 1-3 : dS, for force - // dim = 4-6 : dS * dR, for stress - - void folding_matrix_k(const UnitCell& ucell, - const Grid_Driver& gd, - ForceStressArrays& fsr, - const Parallel_Orbitals& pv, - const int ik, - const int dim1, - const int dim2, - std::complex* mat_k, - const ModuleBase::Vector3& kvec_d); - - /** - * @brief new function of folding_S_matrix - * only for Hamiltonian now, for force and stress will be developed later - * use HContainer as input and output in mat_k - */ - void folding_matrix_k_new(const int ik, - hamilt::Hamilt>* p_ham); - - //============================================================= - // In dftu_force.cpp - // For calculating force and stress fomr DFT+U - //============================================================= - public: - void force_stress(const UnitCell& ucell, - const Grid_Driver& gd, - std::vector>* dmk_d, - std::vector>>* dmk_c, - const Parallel_Orbitals& pv, - ForceStressArrays& fsr, - ModuleBase::matrix& force_dftu, - ModuleBase::matrix& stress_dftu, - const K_Vectors& kv, - const int npol); - - private: - void cal_force_k(const UnitCell& ucell, - const Grid_Driver& gd, - ForceStressArrays& fsr, - const Parallel_Orbitals& pv, - const int ik, - const std::complex* rho_VU, - ModuleBase::matrix& force_dftu, - const ModuleBase::Vector3& kvec_d); - - void cal_stress_k(const UnitCell& ucell, - const Grid_Driver& gd, - ForceStressArrays& fsr, - const Parallel_Orbitals& pv, - const int ik, - const std::complex* rho_VU, - ModuleBase::matrix& stress_dftu, - const ModuleBase::Vector3& kvec_d); - - void cal_force_gamma(const UnitCell& ucell, - const double* rho_VU, - const Parallel_Orbitals& pv, - double* dsloc_x, - double* dsloc_y, - double* dsloc_z, - ModuleBase::matrix& force_dftu); - - void cal_stress_gamma(const UnitCell& ucell, - const Parallel_Orbitals& pv, - const Grid_Driver* gd, - double* dsloc_x, - double* dsloc_y, - double* dsloc_z, - double* dh_r, - const double* rho_VU, - ModuleBase::matrix& stress_dftu); -#endif - - //============================================================= - // In dftu_yukawa.cpp - // Relevant for calculating U using Yukawa potential - //============================================================= - - public: - void cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx); - - private: - void cal_slater_Fk(const UnitCell& ucell,const int L, const int T); // L:angular momnet, T:atom type - void cal_yukawa_lambda(double** rho, const int& nrxx); - - double spherical_Bessel(const int k, const double r, const double lambda); - double spherical_Hankel(const int k, const double r, const double lambda); - -#ifdef __LCAO - public: - /** - * @brief get the density matrix of target spin - * nspin = 1 and 4 : ispin should be 0 - * nspin = 2 : ispin should be 0/1 - */ - const hamilt::HContainer* get_dmr(int ispin) const; - /** - * @brief set the density matrix for DFT+U calculation - * if the density matrix is not set or set to nullptr, the DFT+U calculation will not be performed - */ - void set_dmr(const elecstate::DensityMatrix* dm_in_dftu_d); - void set_dmr(const elecstate::DensityMatrix, double>* dm_in_dftu_cd); - - private: - const UnitCell* ucell = nullptr; - const elecstate::DensityMatrix* dm_in_dftu_d = nullptr; - const elecstate::DensityMatrix, double>* dm_in_dftu_cd = nullptr; -#endif -}; - - -#ifdef __LCAO -template -void dftu_cal_occup_m(const int iter, - const UnitCell& ucell, - const std::vector>& dm, - const K_Vectors& kv, - const double& mixing_beta, - hamilt::Hamilt* p_ham, - Plus_U &dftu); -#endif - - -#endif diff --git a/source/source_lcao/module_dftu/dftu.hpp b/source/source_lcao/module_dftu/dftu.hpp deleted file mode 100644 index 3f8478f7abc..00000000000 --- a/source/source_lcao/module_dftu/dftu.hpp +++ /dev/null @@ -1,18 +0,0 @@ -namespace hamilt -{ - -#ifndef __DFTUTEMPLATE -#define __DFTUTEMPLATE - -/// The DFTU class template inherits from class T -/// it is used to calculate the non-local pseudopotential of wavefunction basis -/// Template parameters: -/// - T: base class, it would be OperatorLCAO or OperatorPW -template -class DFTU : public T -{ -}; - -#endif - -} \ No newline at end of file diff --git a/source/source_lcao/module_dftu/dftu_folding.cpp b/source/source_lcao/module_dftu/dftu_folding.cpp index 07f0b0d0456..4958f9ad594 100644 --- a/source/source_lcao/module_dftu/dftu_folding.cpp +++ b/source/source_lcao/module_dftu/dftu_folding.cpp @@ -1,5 +1,5 @@ #ifdef __LCAO -#include "dftu.h" +#include "dftu_lcao.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" @@ -298,7 +298,7 @@ void Plus_U::folding_matrix_k_new(const int ik, } else { - if(Plus_U::nspin != 4) + if(this->nspin != 4) { dynamic_cast, double>*>(p_ham) ->updateSk(ik, hk_type); diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_force.cpp index 1f58d094d9f..3f119721bb2 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_force.cpp @@ -1,7 +1,7 @@ #include "source_io/module_parameter/parameter.h" #ifdef __LCAO -#include "dftu.h" +#include "dftu_lcao.h" #include "source_base/constants.h" #include "source_base/global_function.h" #include "source_base/inverse_matrix.h" @@ -36,6 +36,52 @@ void Plus_U::force_stress(const UnitCell& ucell, ModuleBase::TITLE("Plus_U", "force_stress"); ModuleBase::timer::start("Plus_U", "force_stress"); + // Defensive null check: the legacy dft_plus_u==2 force/stress path + // requires fsr.DSloc_x/y/z (gamma_only) or fsr.DSloc_Rx/Ry/Rz (multik) + // and fsr.DH_r to be allocated and filled by the caller. If the caller + // forgot to allocate them (as in force_stress_lcao.cpp where the local + // fsr_dftu is created without allocation), we fail early with a clear + // message instead of letting pdgemm_ dereference nullptr and crash. + // See force_stress_lcao.cpp for the historical background. + if (this->gamma_only_local) + { + if (this->cal_force + && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr)) + { + ModuleBase::WARNING_QUIT("Plus_U::force_stress", + "fsr.DSloc_x/y/z are nullptr in gamma_only path; the caller must allocate and fill them. " + "See notes in source/source_lcao/force_stress_lcao.cpp."); + } + if (this->cal_stress + && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr + || fsr.DH_r == nullptr)) + { + ModuleBase::WARNING_QUIT("Plus_U::force_stress", + "fsr.DSloc_x/y/z or fsr.DH_r is nullptr in gamma_only path; " + "the caller must allocate and fill them. " + "See notes in source/source_lcao/force_stress_lcao.cpp."); + } + } + else + { + if (this->cal_force + && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr)) + { + ModuleBase::WARNING_QUIT("Plus_U::force_stress", + "fsr.DSloc_Rx/Ry/Rz are nullptr in multik path; the caller must allocate and fill them. " + "See notes in source/source_lcao/force_stress_lcao.cpp."); + } + if (this->cal_stress + && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr + || fsr.DH_r == nullptr)) + { + ModuleBase::WARNING_QUIT("Plus_U::force_stress", + "fsr.DSloc_Rx/Ry/Rz or fsr.DH_r is nullptr in multik path; " + "the caller must allocate and fill them. " + "See notes in source/source_lcao/force_stress_lcao.cpp."); + } + } + const int nlocal = this->nlocal; if (this->cal_force) diff --git a/source/source_lcao/module_dftu/dftu_fs.cpp b/source/source_lcao/module_dftu/dftu_fs.cpp index 7db4896ed06..72f2ffcce0f 100644 --- a/source/source_lcao/module_dftu/dftu_fs.cpp +++ b/source/source_lcao/module_dftu/dftu_fs.cpp @@ -1,4 +1,4 @@ -#include "dftu_lcao.h" +#include "dftu_lcao_op.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" @@ -141,7 +141,7 @@ void DFTU>::cal_force_stress(const bool cal_force, this->dftu->get_occ_mat_flat(iat0, target_L, occ); // calculate VU - const double u_value = this->dftu->U[T0]; + const double u_value = this->dftu->u_current[T0]; std::vector VU(occ.size()); double eu_tmp = 0; this->cal_v_of_u(occ, tlp1, u_value, &VU[0], eu_tmp); diff --git a/source/source_lcao/module_dftu/dftu_hamilt.cpp b/source/source_lcao/module_dftu/dftu_hamilt.cpp index fc1afdbc8e5..c8362bc749b 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/source_lcao/module_dftu/dftu_hamilt.cpp @@ -1,4 +1,4 @@ -#include "dftu.h" +#include "dftu_lcao.h" #include "source_base/module_external/scalapack_connector.h" #include "source_io/module_parameter/parameter.h" #include "source_base/timer.h" diff --git a/source/source_lcao/module_dftu/dftu_lcao.cpp b/source/source_lcao/module_dftu/dftu_lcao.cpp index 3d02c7b3e5e..67e915e0713 100644 --- a/source/source_lcao/module_dftu/dftu_lcao.cpp +++ b/source/source_lcao/module_dftu/dftu_lcao.cpp @@ -1,692 +1,327 @@ #include "dftu_lcao.h" -#include "source_base/timer.h" -#include "source_base/tool_title.h" -#include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_lcao/module_operator_lcao/operator_lcao.h" -#include "source_hamilt/module_hcontainer/hcontainer_funcs.h" #include "source_io/module_parameter/parameter.h" -#ifdef _OPENMP -#include +#include "source_base/constants.h" +#include "source_base/global_function.h" +#include "source_base/inverse_matrix.h" +#include "source_base/memory_recorder.h" +#include "source_base/timer.h" +#include "source_cell/magnetism.h" +#include "source_estate/module_charge/charge.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // mohan add 2025-11-06 +// Static member definitions moved to dftu_base.cpp (Plus_U_Base::) +// Plus_U inherits these from Plus_U_Base. + +Plus_U::Plus_U() +{} + +Plus_U::~Plus_U() +{} + +void Plus_U::init(UnitCell& cell, + const Parallel_Orbitals* pv, + const int npol, + 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, + const int nlocal, + const bool gamma_only_local, + const std::string& ks_solver, + const bool cal_force, + const bool cal_stress, + const std::string& device, + const int kpar, + const std::vector& hubbard_u, + const double uramping, + const int occ_mat_ctrl, + const int mixing_dftu +#ifdef __LCAO + , const LCAO_Orbitals* orb #endif -#include "source_base/parallel_reduce.h" - -template -hamilt::DFTU>::DFTU(HS_Matrix_K* hsk_in, - const std::vector>& kvec_d_in, - hamilt::HContainer* hR_in, - const UnitCell& ucell_in, - const Grid_Driver* GridD_in, - const TwoCenterIntegrator* intor, - const std::vector& orb_cutoff, - Plus_U* p_dftu) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor), orb_cutoff_(orb_cutoff) + ) { - this->cal_type = calculation_type::lcao_dftu; - this->ucell = &ucell_in; - this->dftu = p_dftu; -#ifdef __DEBUG - assert(this->ucell != nullptr); -#endif - // initialize HR to allocate sparse Nonlocal matrix memory - this->initialize_HR(GridD_in); - // set nspin - this->nspin = PARAM.inp.nspin; -} + ModuleBase::TITLE("Plus_U", "init"); -// destructor -template -hamilt::DFTU>::~DFTU() -{ -} + this->paraV = pv; -// initialize_HR() -template -void hamilt::DFTU>::initialize_HR(const Grid_Driver* GridD) -{ - ModuleBase::TITLE("DFTU", "initialize_HR"); - ModuleBase::timer::start("DFTU", "initialize_HR"); + this->yukawa_lambda = yukawa_lambda; + this->npol = npol; + this->nlocal = nlocal; + this->gamma_only_local = gamma_only_local; + this->ks_solver = ks_solver; + this->cal_force = cal_force; + this->cal_stress = cal_stress; - this->adjs_all.clear(); - this->adjs_all.reserve(this->ucell->nat); - for (int iat0 = 0; iat0 < ucell->nat; iat0++) +#ifdef __LCAO + ptr_orb_ = orb; + if(ptr_orb_ != nullptr) { - 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); - - AdjacentAtomInfo adjs; - GridD->Find_atom(*ucell, tau0, T0, I0, &adjs); - std::vector is_adj(adjs.adj_num + 1, false); - 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); - const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; - const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - // choose the real adjacent atoms - // Note: the distance of atoms should less than the cutoff radius, - // When equal, the theoretical value of matrix element is zero, - // but the calculated value is not zero due to the numerical error, which would lead to result changes. - if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb_cutoff_[T1] + PARAM.inp.onsite_radius) - { - is_adj[ad1] = true; - } - } - filter_adjs(is_adj, adjs); - this->adjs_all.push_back(adjs); + orb_cutoff_ = orb->cutoffs(); } + ucell = &cell; +#endif - ModuleBase::timer::end("DFTU", "initialize_HR"); -} - -template -void hamilt::DFTU>::cal_nlm_all(const Parallel_Orbitals* paraV) -{ - ModuleBase::TITLE("DFTU", "cal_nlm_all"); - if (this->precal_nlm_done) - { - return; - } - - ModuleBase::timer::start("DFTU", "cal_nlm_all"); - nlm_tot.resize(this->ucell->nat); - const int npol = this->ucell->get_npol(); - int atom_index = 0; - for (int iat0 = 0; iat0 < ucell->nat; iat0++) + if (pv != nullptr) { - auto tau0 = ucell->get_tau(iat0); - int T0=0; - int I0=0; - ucell->iat2iait(iat0, &I0, &T0); - if (!this->dftu->has_correlated_orbital(T0)) + const int global_rows = pv->get_global_row_size(); + const int global_cols = pv->get_global_col_size(); + if (global_rows != global_cols) { - continue; + ModuleBase::WARNING_QUIT("Plus_U::init", "Global row and column dimensions do not match"); } - const int target_L = this->dftu->get_orbital_corr(T0); - const int tlp1 = 2 * target_L + 1; - AdjacentAtomInfo& adjs = this->adjs_all[atom_index++]; - - // calculate and save the table of two-center integrals - nlm_tot[iat0].resize(adjs.adj_num + 1); - - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + if (nlocal != global_rows) { - 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 = paraV->get_indexes_row(iat1); - auto col_indexes = paraV->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; - // only first zeta orbitals in target L of atom iat0 are needed - std::vector nlm_target(tlp1); - const int L1 = atom1->iw2l[iw1]; - const int N1 = atom1->iw2n[iw1]; - const int m1 = atom1->iw2m[iw1]; - 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) - - // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) - const 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, false /*cal_deri*/, nlm); - // select the elements of nlm with target_L - 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 < 2 * L0 + 1; m++) - { - nlm_target[m] = nlm[0][iw + m]; - } - break; - } - } - nlm_tot[iat0][ad].insert({all_indexes[iw1l], nlm_target}); - } + ModuleBase::WARNING_QUIT("Plus_U::init", "nlocal does not match global matrix dimension"); } } - this->precal_nlm_done = true; - ModuleBase::timer::end("DFTU", "cal_nlm_all"); + + this->init_base(cell, + npol, + nspin, + orbital_corr, + yukawa_potential, + global_readin_dir, + global_out_dir, + init_chg, + device, + kpar, + hubbard_u, + uramping, + occ_mat_ctrl, + mixing_dftu); + return; } -// contributeHR() -/** - * @brief Contribute DFT+U Hamiltonian to real-space HR matrix - * - * @details This function handles different scenarios based on: - * 1. Whether occ_mat (occupation matrix) is read from file (is_occ_mat_initialized) - * 2. Spin configuration (nspin=1, 2, or 4) - * 3. SCF iteration stage (first vs subsequent iterations) - * - * Case 1: Occ_mat NOT initialized (!is_occ_mat_initialized) - * - First electronic iteration: calculates occupation matrix from density matrix (DMR) - * * Uses get_dmr(current_spin) to get real-space density matrix - * * Accumulates contributions from all atom pairs via cal_occ() - * * Performs MPI reduction to sum occ across processes - * * Stores result via set_occ_mat_flat() for use in VU calculation - * * For nspin=1: occ is scaled by 0.5 (since only one spin channel computed) - * - Subsequent iterations: occ_mat is computed fresh each iteration from updated DMR - * - * Case 2: Occ_mat IS initialized (is_occ_mat_initialized, i.e., read from dm_onsite.txt file) - * - First electronic iteration: uses pre-read occ_mat directly without DMR calculation - * * Skips DMR-based occ calculation entirely - * * Reads locale from stored data via get_occ_mat() - * * Different indexing for nspin=4 vs nspin=1/2 (see below) - * - After first iteration: mark_occ_mat_dirty() is called to force recomputation - * - * Spin configurations: - * nspin=1 (non-spin-polarized): - * - Single spin channel, occ computed once - * - Energy correction doubled at end (set_double_energy) - * - current_spin always 0 - * - * nspin=2 (collinear spin-polarized): - * - Two separate spin channels (spin-up: 0, spin-down: 1) - * - current_spin toggles between 0 and 1 across iterations - * - mark_occ_mat_dirty() called when current_spin == 1 (last spin) - * - HR accumulated separately for each spin - * - * nspin=4 (non-collinear/SOC): - * - Single 4x4 Pauli matrix representation per atom - * - occ has 4*(2l+1)^2 elements (spin_fold=4) - * - get_occ_mat uses spin=0, ipol indices for Pauli blocks - * - mark_occ_mat_dirty() always called (current_spin check always true) - * - No current_spin toggling (all spins handled simultaneously) - * - * @warning THREAD SAFETY: cal_HR_IJR() updates shared HR matrix entries. - * Different iat0 may contribute to same HR(iat1, iat2, R), requiring - * critical section protection for multithreaded correctness. - * TODO: Consider refactoring to atom_row_list pattern (see nonlocal.cpp) - * for better parallel performance instead of critical section. - */ -template -void hamilt::DFTU>::contributeHR() +#ifdef __LCAO + +void Plus_U::cal_energy_correction(const UnitCell& ucell, + const int istep) { - ModuleBase::TITLE("DFTU", "contributeHR"); - // Early exit conditions: - // - get_dmr(0) == nullptr: DMR not available (typical in first iteration without file input) - // - !is_occ_mat_initialized(): occ_mat not read from file AND not yet computed from DMR - // When both true, skip DFT+U contribution entirely (first iteration, no file input) - const bool dmr_null = (this->dftu->get_dmr(0) == nullptr); - const bool occ_mat_not_init = !this->dftu->is_occ_mat_initialized(); - - if (dmr_null && occ_mat_not_init) + ModuleBase::TITLE("Plus_U", "cal_energy_correction"); + ModuleBase::timer::start("Plus_U", "cal_energy_correction"); + if (!is_occ_mat_initialized()) { + ModuleBase::timer::end("Plus_U", "cal_energy_correction"); return; } - else - { - // Reset DFT+U energy at start of each spin cycle - // For nspin=2: reset when current_spin==0 (start of spin-up calculation) - // For nspin=1/4: reset once (current_spin always 0) - if (this->current_spin == 0) - { - this->dftu->set_energy(0.0); - } - } - ModuleBase::timer::start("DFTU", "contributeHR"); - - const Parallel_Orbitals* paraV = this->hR->get_atom_pair(0).get_paraV(); - const int npol = this->ucell->get_npol(); - // 1. Calculate two-center integrals for all atom pairs - // This is reused in both occ and HR calculations - this->cal_nlm_all(paraV); - - // 2. Loop over all Hubbard-projector center atoms (iat0) - int atom_index = 0; - for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + + // mohan update 20251106 + this->energy_u = 0.0; + + double energy_dc = 0.0; + + for (int T = 0; T < ucell.ntype; T++) { - auto tau0 = ucell->get_tau(iat0); - int T0, I0; - 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++]; - - ModuleBase::timer::start("DFTU", "cal_occ"); - // spin_fold: number of spin components in occ array - // nspin=4: 4 (Pauli matrix blocks), nspin=1/2: 1 (single spin channel) - const int spin_fold = (this->nspin == 4) ? 4 : 1; - std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); - - // ============================================================ - // BRANCH 1: Occ_mat NOT initialized (compute from DMR) - // ============================================================ - // This branch is taken when: - // - is_occ_mat_initialized() == false (no file read or omc != 0) - // - DMR is available (get_dmr() != nullptr) - // Typical scenario: normal SCF iterations after first step - if (!this->dftu->is_occ_mat_initialized()) + const int NL = ucell.atoms[T].nwl + 1; + const int LC = get_orbital_corr(T); + for (int I = 0; I < ucell.atoms[T].na; I++) { - // TODO: UNSAFE - get_dmr(current_spin) assumes DMR has correct spin indexing. - // For nspin=2, current_spin must be correctly toggled (0 then 1). - // If current_spin is wrong, wrong spin channel's DMR is used. - const hamilt::HContainer* dmR_current = this->dftu->get_dmr(this->current_spin); - for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + if (LC == -1) { - const int T1 = adjs.ntype[ad1]; - const int I1 = adjs.natom[ad1]; - const int iat1 = ucell->itia2iat(T1, I1); - ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - const std::unordered_map>& nlm1 = nlm_tot[iat0][ad1]; - 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); - const std::unordered_map>& nlm2 = nlm_tot[iat0][ad2]; - ModuleBase::Vector3& R_index2 = adjs.box[ad2]; - 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 - = dmR_current->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); - if (tmp != nullptr) - { - this->cal_occ(iat1, iat2, paraV, nlm1, nlm2, tmp->get_pointer(), occ); - } - } - } -#ifdef __MPI - // CRITICAL: MPI reduction required for distributed DMR calculations. - // Each process computes partial occ from its local DMR blocks. - // Without this, occ would be incomplete and DFT+U potential wrong. - // TODO: Verify that occ size is consistent across processes before reduction. - // TODO: Consider using MPI_IN_PLACE to avoid extra buffer allocation. - Parallel_Reduce::reduce_all(occ.data(), occ.size()); -#endif - // For nspin=1: occ computed from single spin channel, but should represent - // total occupation (both spins). Scale by 0.5 to account for this. - if (this->nspin == 1) - { - for (auto& v : occ) { v *= 0.5; } + continue; } - this->dftu->set_occ_mat_flat(iat0, target_L, this->current_spin, occ); - } - // ============================================================ - // BRANCH 2: Occ_mat IS initialized (use pre-read data) - // ============================================================ - // This branch is taken when: - // - is_occ_mat_initialized() == true (occ_mat read from dm_onsite.txt file) - // - OR omc != 0 (occupation matrix control with dm_onsite_ini.txt) - // Typical scenario: first SCF iteration with file input, or restart calculation - else - { - // nspin=4: Non-collinear case with Pauli matrix representation - // Occ_mat stored as single 4x4 block per atom, with spin indices embedded - // in the matrix indices (ipol0, ipol1 for Pauli block indices) - if (this->nspin == 4) - { - // 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); - } - // nspin=1 or nspin=2: Collinear spin case - // Occ_mat stored separately for each spin channel - else + + const int iat = ucell.itia2iat(T, I); + const int L = get_orbital_corr(T); + + for (int l = 0; l < NL; l++) { - for (int i = 0; i < static_cast(occ.size()); i++) + if (l != get_orbital_corr(T)) { - // 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, - i / (2 * target_L + 1), i % (2 * target_L + 1)); + continue; } - } - } - ModuleBase::timer::end("DFTU", "cal_occ"); - - // 3. Calculate Hubbard potential VU from occupation matrix - // VU = U * (1/2 * delta(m,m') - occ(m,m')) for each spin channel - // Energy: EU = U * 1/2 * occ(m,m') * occ(m',m) - ModuleBase::timer::start("DFTU", "cal_vu"); - const double u_value = this->dftu->U[T0]; - std::vector VU_tmp(occ.size()); - - // TODO: GLOBAL STATE - Plus_U::get_energy()/set_energy() uses static member variable. - // This is NOT thread-safe for parallel SCF calculations. - // TODO: Refactor to use instance member or pass energy by reference. - double u_energy = Plus_U::get_energy(); - this->cal_v_of_u(occ, tlp1, u_value, VU_tmp.data(), u_energy); - Plus_U::set_energy(u_energy); - - // 4. Convert VU to appropriate data type (real or complex) - // For nspin=4 with complex Hamiltonian, VU needs Pauli matrix transformation - std::vector VU(occ.size()); - this->transfer_vu(VU_tmp, VU); - - // 5. Second iteration: Calculate Hamiltonian matrix contribution - // HR += * VU(m,m') * - // for all atom pairs within cutoff - // Note: different iat0 may contribute to the same HR(iat1, iat2, R), so we need to protect the update - // to avoid race conditions in multithreading. Reference: nonlocal.cpp for the atom_row_list pattern. - // TODO: CRITICAL SECTION PERFORMANCE - This critical section serializes HR updates. - // For systems with many Hubbard atoms, this becomes a bottleneck. - // Consider refactoring to atom_row_list pattern (see nonlocal.cpp lines 127-220): - // 1. Use #pragma omp for to distribute iat0 across threads - // 2. Each thread records its assigned iat0 in thread-local atom_row_list - // 3. When updating HR(iat1, iat2, R), skip if iat1 not in thread's atom_row_list - // 4. This eliminates race conditions without critical section - 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); - ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - const std::unordered_map>& nlm1 = nlm_tot[iat0][ad1]; - 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); - const std::unordered_map>& nlm2 = nlm_tot[iat0][ad2]; - ModuleBase::Vector3& R_index2 = adjs.box[ad2]; - ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], - R_index2[1] - R_index1[1], - R_index2[2] - R_index1[2]); - hamilt::BaseMatrix* tmp = this->hR->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); - if (tmp != nullptr) + + const int N = ucell.atoms[T].l_nchi[l]; + + const int m_tot = 2 * l + 1; + + // part 1: calculate the DFT+U energy correction + for (int n = 0; n < N; n++) { -#ifdef _OPENMP -#pragma omp critical(dftu_hr_update) -#endif + if (n != 0) { - this->cal_HR_IJR(iat1, iat2, paraV, nlm1, nlm2, VU, tmp->get_pointer()); + continue; } - } - } - } - ModuleBase::timer::end("DFTU", "cal_vu"); - } - // 6. Post-processing: Energy correction and occ_mat state management - // For nspin=1: DFT+U energy computed for single spin channel, but should count both spins - // set_double_energy() doubles the energy to account for degenerate spin-up/down - if (this->nspin == 1) - { - this->dftu->set_double_energy(); - } - - // 7. Mark occ_mat as dirty to force recomputation in next iteration - // This is called when: - // - nspin=4: Always (all spins handled simultaneously, current_spin==0==nspin-1) - // - nspin=2: When current_spin==1 (after spin-down calculation, last spin channel) - // - nspin=1: When current_spin==0==nspin-1 (always called) - // - // Purpose: Ensure occ_mat is recomputed from updated DMR in next SCF iteration, - // rather than using stale pre-read data from file. - // TODO: This logic is confusing. Consider explicit variable like `is_last_spin_channel`. - if (this->current_spin == this->nspin - 1 || this->nspin == 4) - { - this->dftu->mark_occ_mat_dirty(); - } - - // 8. Spin channel toggling for nspin=2 - // nspin=2 requires separate HR updates for spin-up (current_spin=0) and spin-down (current_spin=1) - // The HR matrix is updated twice per SCF iteration, once for each spin channel - // current_spin toggles: 0 -> 1 -> 0 -> 1 ... - // For nspin=1: current_spin always 0 (no toggling needed) - // For nspin=4: current_spin always 0 (all spins handled simultaneously via Pauli matrices) - // TODO: UNSAFE - This assumes contributeHR() is called in strict alternating order. - // If called out of order (e.g., due to parallel k-point distribution), current_spin may be wrong. - // TODO: Consider deriving current_spin from ik or explicit parameter instead of toggling. - if (this->nspin == 2) - { - this->current_spin = 1 - this->current_spin; - } - - ModuleBase::timer::end("DFTU", "contributeHR"); -} + if (this->nspin == 1 || this->nspin == 2) + { + for (int spin = 0; spin < 2; spin++) + { + double nm_trace = 0.0; + double nm2_trace = 0.0; + + for (int m0 = 0; m0 < 2 * l + 1; m0++) + { + nm_trace += this->occ_mat[iat][l][n][spin](m0, m0); + for (int m1 = 0; m1 < 2 * l + 1; m1++) + { + nm2_trace += this->occ_mat[iat][l][n][spin](m0, m1) + * this->occ_mat[iat][l][n][spin](m1, m0); + } + } + if (use_yukawa) + { + this->energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) + * (nm_trace - nm2_trace); + } + else + { + this->energy_u += 0.5 * this->u_current[T] * (nm_trace - nm2_trace); + } + } + } + else if (this->nspin == 4) + { + double nm_trace = 0.0; + double nm2_trace = 0.0; -// cal_HR_IJR() -template -void hamilt::DFTU>::cal_HR_IJR( - const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& VU, - TR* data_pointer) -{ + for (int m0 = 0; m0 < 2 * l + 1; m0++) + { + for (int ipol0 = 0; ipol0 < this->npol; ipol0++) + { + const int m0_all = m0 + (2 * l + 1) * ipol0; + nm_trace += this->occ_mat[iat][l][n][0](m0_all, m0_all); + + for (int m1 = 0; m1 < 2 * l + 1; m1++) + { + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) + { + int m1_all = m1 + (2 * l + 1) * ipol1; + + nm2_trace += this->occ_mat[iat][l][n][0](m0_all, m1_all) + * this->occ_mat[iat][l][n][0](m1_all, m0_all); + } + } + } + } + if (use_yukawa) + { + this->energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) + * (nm_trace - nm2_trace); + } + else + { + this->energy_u += 0.5 * this->u_current[T] * (nm_trace - nm2_trace); + } + } - // 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 = paraV->get_indexes_row(iat1); - auto col_indexes = paraV->get_indexes_col(iat2); - const int m_size = int(sqrt(VU.size()) / npol); - // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 - std::vector step_trace(npol * npol, 0); - for (int is = 0; is < npol; is++) - { - for (int is2 = 0; is2 < npol; is2++) - { - step_trace[is * npol + is2] = paraV->get_ncol_atom(iat2) * is + is2; - } - } - // calculate the local matrix - const TR* tmp_d = nullptr; - 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 is = 0; is < npol * npol; ++is) - { - int start = is * m_size * m_size; - TR nlm_tmp = TR(0); - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) + for (int m1 = 0; m1 < 2 * l + 1; m1++) { - nlm_tmp += nlm1[m1] * nlm2[m2] * VU[m1 * m_size + m2 + start]; + for (int ipol1 = 0; ipol1 < this->npol; ipol1++) + { + const int m1_all = m1 + ipol1 * (2 * l + 1); + for (int m2 = 0; m2 < 2 * l + 1; m2++) + { + for (int ipol2 = 0; ipol2 < this->npol; ipol2++) + { + const int m2_all = m2 + ipol2 * (2 * l + 1); + + if (this->nspin == 1 || this->nspin == 2) + { + for (int is = 0; is < 2; is++) + { + double VU = 0.0; + VU = get_onebody_eff_pot(T, iat, l, n, is, m1_all, m2_all, false); + energy_dc += VU * this->occ_mat[iat][l][n][is](m1_all, m2_all); + } + } + else if (this->nspin == 4) + { + double VU = 0.0; + VU = get_onebody_eff_pot(T, iat, l, n, 0, m1_all, m2_all, false); + energy_dc += VU * this->occ_mat[iat][l][n][0](m1_all, m2_all); + } + } + } + } } - } - data_pointer[step_trace[is]] += nlm_tmp; - } - data_pointer += npol; - } - data_pointer += (npol - 1) * col_indexes.size(); - } + } // end n + } // end L + } // end I + } // end T + + // substract the double counting energy_dc included in band energy eband + this->energy_u -= energy_dc; + + ModuleBase::timer::end("Plus_U", "cal_energy_correction"); + return; } -template -void hamilt::DFTU>::cal_occ(const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const double* dm_pointer, - std::vector& occ) +#endif + +// uramping_update() and u_converged() are now implemented in +// dftu_base.cpp as Plus_U_Base methods (inherited by Plus_U). + +#ifdef __LCAO + +void Plus_U::set_dmr(const elecstate::DensityMatrix, double>* dmr) { + this->dm_in_dftu_cd = dmr; + return; +} - // 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 = paraV->get_indexes_row(iat1); - auto col_indexes = paraV->get_indexes_col(iat2); - const int m_size = int(sqrt(occ.size()) / npol); - const int m_size2 = m_size * m_size; -#ifdef __DEBUG - assert(m_size * m_size == occ.size()); -#endif - // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 - std::vector step_trace(npol * npol, 0); - for (int is = 0; is < npol; is++) +void Plus_U::set_dmr(const elecstate::DensityMatrix* dmr) +{ + this->dm_in_dftu_d = dmr; + return; +} + +const hamilt::HContainer* Plus_U::get_dmr(int ispin) const +{ + if (this->dm_in_dftu_d != nullptr) { - for (int is2 = 0; is2 < npol; is2++) - { - step_trace[is * npol + is2] = paraV->get_ncol_atom(iat2) * is + is2; - } + return this->dm_in_dftu_d->get_DMR_pointer(ispin + 1); } - // calculate the local matrix - for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + else if (this->dm_in_dftu_cd != nullptr) { - 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 is1 = 0; is1 < npol; ++is1) - { - for (int is2 = 0; is2 < npol; ++is2) - { - for (int m1 = 0; m1 < m_size; ++m1) - { - for (int m2 = 0; m2 < m_size; ++m2) - { - occ[m1 * m_size + m2 + (is1 * npol + is2) * m_size2] - += nlm1[m1] * nlm2[m2] * dm_pointer[step_trace[is1 * npol + is2]]; - } - } - } - } - dm_pointer += npol; - } - dm_pointer += (npol - 1) * col_indexes.size(); + return this->dm_in_dftu_cd->get_DMR_pointer(ispin + 1); } -} - -template -void hamilt::DFTU>::transfer_vu(std::vector& vu_tmp, std::vector& vu) -{ -#ifdef __DEBUG - assert(vu.size() == vu_tmp.size()); -#endif - for (int i = 0; i < vu_tmp.size(); i++) + else { - vu[i] = vu_tmp[i]; + return nullptr; } } +//! dftu occupation matrix for gamma only using dm(double) template <> -void hamilt::DFTU, std::complex>>::transfer_vu( - std::vector& vu_tmp, - std::vector>& vu) +void dftu_cal_occup_m(const int iter, + const UnitCell& ucell, + const std::vector>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt* p_ham, + Plus_U &dftu) { -#ifdef __DEBUG - assert(vu.size() == vu_tmp.size()); -#endif - - // Pauli-to-spinor conversion for DFT+U potential: - // V = V_0*I + V_x*sigma_x + V_y*sigma_y + V_z*sigma_z - // sigma_y = [[0,-i],[i,0]], so: - // V_{up,up} = 0.5*(V_0 + V_z) - // V_{down,down} = 0.5*(V_0 - V_z) - // V_{up,down} = 0.5*(V_x - i*V_y) <- note: minus sign from sigma_y - // V_{down,up} = 0.5*(V_x + i*V_y) <- note: plus sign from sigma_y - // This is consistent with the convention in gint_common.cpp merge_hr_part_to_hR(). - const int m_size = int(sqrt(vu.size()) / 2); - const int m_size2 = m_size * m_size; - vu.resize(vu_tmp.size()); - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - int index[4]; - index[0] = m1 * m_size + m2; - index[1] = m1 * m_size + m2 + m_size2; - index[2] = m2 * m_size + m1 + m_size2 * 2; - index[3] = m2 * m_size + m1 + m_size2 * 3; - vu[index[0]] = 0.5 * (vu_tmp[index[0]] + vu_tmp[index[3]]); - vu[index[3]] = 0.5 * (vu_tmp[index[0]] - vu_tmp[index[3]]); - vu[index[1]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); - vu[index[2]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); - } - } + dftu.cal_occup_m_gamma(iter, ucell ,dm, mixing_beta, p_ham); } -template -void hamilt::DFTU>::cal_v_of_u(const std::vector& occ, - const int m_size, - const double u_value, - double* vu, - double& eu) +//! dftu occupation matrix for multiple k-points using dm(complex) +template <> +void dftu_cal_occup_m(const int iter, + const UnitCell& ucell, + const std::vector>>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt>* p_ham, + Plus_U &dftu) { - // calculate the local matrix - int spin_fold = occ.size() / m_size / m_size; - if (spin_fold < 4) { - for (int is = 0; is < spin_fold; ++is) - { - int start = is * m_size * m_size; - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - vu[start + m1 * m_size + m2] = u_value * (0.5 * (m1 == m2) - occ[start + m2 * m_size + m1]); - eu += u_value * 0.5 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; - } - } - } - } else - { - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - vu[m1 * m_size + m2] = u_value * (1.0 * (m1 == m2) - occ[m2 * m_size + m1]); - eu += u_value * 0.25 * occ[m2 * m_size + m1] * occ[m1 * m_size + m2]; - } - } - for (int is = 1; is < spin_fold; ++is) - { - int start = is * m_size * m_size; - for (int m1 = 0; m1 < m_size; m1++) - { - for (int m2 = 0; m2 < m_size; m2++) - { - vu[start + m1 * m_size + m2] = u_value * (0 - occ[start + m2 * m_size + m1]); - eu += u_value * 0.25 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; - } - } - } - } + dftu.cal_occup_m_k(iter,ucell, dm, kv, mixing_beta, p_ham); } -template class hamilt::DFTU>; -template class hamilt::DFTU, double>>; -template class hamilt::DFTU, std::complex>>; +#endif diff --git a/source/source_lcao/module_dftu/dftu_lcao.h b/source/source_lcao/module_dftu/dftu_lcao.h index 2ce7defb417..d617814a872 100644 --- a/source/source_lcao/module_dftu/dftu_lcao.h +++ b/source/source_lcao/module_dftu/dftu_lcao.h @@ -1,142 +1,280 @@ -#ifndef DFTPLUSU_H -#define DFTPLUSU_H -#include "dftu.hpp" -#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" +#ifndef DFTU_H +#define DFTU_H + +#include "source_cell/klist.h" #include "source_cell/unitcell.h" -#include "source_estate/module_dm/density_matrix.h" -#include "source_lcao/module_operator_lcao/operator_lcao.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_basis/module_ao/parallel_orbitals.h" +#include "source_estate/module_charge/charge_mixing.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" +#include "source_lcao/force_stress_arrays.h" // mohan add 2024-06-15 +#endif -#include +#include +#include -namespace hamilt -{ -/// DFTU class template specialization for OperatorLCAO base class -/// It is used to calculate the non-local pseudopotential matrix in real space and fold it to k-space -/// HR = D_{p1, p2} -/// HK = D_{p1, p2} = \sum_{R} e^{ikR} HR -/// Template parameters: -/// - TK: data type of k-space Hamiltonian -/// - TR: data type of real space Hamiltonian -template -class DFTU> : public OperatorLCAO +class Plus_U : public Plus_U_Base { + public: - DFTU>(HS_Matrix_K* hsk_in, - const std::vector>& kvec_d_in, - hamilt::HContainer* hR_in, - const UnitCell& ucell_in, - const Grid_Driver* gridD_in, - const TwoCenterIntegrator* intor, - const std::vector& orb_cutoff, - Plus_U* p_dftu); - ~DFTU>(); + Plus_U(); + ~Plus_U(); - /** - * @brief contributeHR() is used to calculate the HR matrix - * D_{p1, p2} - */ - virtual void contributeHR() override; + public: + // allocate relevant data strcutures + void init(UnitCell& cell, + const Parallel_Orbitals* pv, + const int npol, + 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, + const int nlocal, + const bool gamma_only_local, + const std::string& ks_solver, + const bool cal_force, + const bool cal_stress, + const std::string& device, + const int kpar, + const std::vector& hubbard_u, + const double uramping, + const int occ_mat_ctrl, + const int mixing_dftu +#ifdef __LCAO + , const LCAO_Orbitals* orb = nullptr +#endif + ); - /// calculate force and stress for DFT+U - void cal_force_stress(const bool cal_force, - const bool cal_stress, - ModuleBase::matrix& force, - ModuleBase::matrix& stress); + // calculate the energy correction + void cal_energy_correction(const UnitCell& ucell, const int istep); private: - const UnitCell* ucell = nullptr; - - Plus_U* dftu = nullptr; - hamilt::HContainer* HR = nullptr; + const Parallel_Orbitals* paraV = nullptr; - const TwoCenterIntegrator* intor_ = nullptr; + double yukawa_lambda = 0.0; + int npol = 1; + int nlocal = 0; + bool gamma_only_local = false; + std::string ks_solver; + bool cal_force = false; + bool cal_stress = false; +#ifdef __LCAO + const LCAO_Orbitals* ptr_orb_ = nullptr; std::vector orb_cutoff_; +#endif - /// @brief the number of spin components, 1 for no-spin, 2 for collinear spin case and 4 for non-collinear spin case - int nspin = 0; +#ifdef __LCAO + //============================================================= + // In dftu_hamilt.cpp + // For calculating contribution to Hamiltonian matrices + //============================================================= + public: + void cal_eff_pot_mat_complex(const int ik, + std::complex* eff_pot, + const std::vector& isk, + const std::complex* sk, + const int npol); - /** - * @brief search the nearest neighbor atoms and save them into this->adjs_all - * the size of HR will not change in DFTU, - * because I don't want to expand HR larger than Nonlocal operator caused by DFTU - */ - void initialize_HR(const Grid_Driver* gridD_in); + void cal_eff_pot_mat_real(const int ik, + double* eff_pot, + const std::vector& isk, + const double* sk, + const int npol); - /** - * @brief calculate the overlap values and save them in this->nlm_tot - * it will be reused in the calculation of calculate_HR() - */ - void cal_nlm_all(const Parallel_Orbitals* paraV); + void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR, const int npol); - /** - * @brief calculate the occ_mm' = \sum_R DMR* matrix for each atom to add U - */ - void cal_occ(const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const double* data_pointer, - std::vector& occupations); - - /// transfer VU format from pauli matrix to normal for non-collinear spin case - void transfer_vu(std::vector& vu_tmp, std::vector& vu); - /// VU_{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_v_of_u(const std::vector& occ, const int m_size, const double u_value, double* vu, double& eu); + void cal_eff_pot_mat_R_complex_double(const int ispin, + std::complex* SR, + std::complex* HR, + const int npol); +#endif + +#ifdef __LCAO + // calculate the local occupation number matrix + void cal_occup_m_k(const int iter, + const UnitCell& ucell, + const std::vector>>& dm_k, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt>* p_ham); + + void cal_occup_m_gamma(const int iter, + const UnitCell& ucell, + const std::vector>& dm_gamma, + const double& mixing_beta, + hamilt::Hamilt* p_ham); +#endif + +#ifdef __LCAO +private: + //============================================================= + // In dftu_tools.cpp + // For calculating onsite potential, which is used + // for both Hamiltonian and force/stress + //============================================================= + + void cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU, const int npol); + void cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU, const int npol); + + double get_onebody_eff_pot(const int T, + const int iat, + const int L, + const int N, + const int spin, + const int m0, + const int m1, + const bool newlocale); + + //============================================================= + // In dftu_folding.cpp + // Subroutines for folding S and dS matrix + //============================================================= + + void fold_dSR_gamma(const UnitCell& ucell, + const Parallel_Orbitals& pv, + const Grid_Driver* gd, + double* dsloc_x, + double* dsloc_y, + double* dsloc_z, + double* dh_r, + const int dim1, + const int dim2, + double* dSR_gamma); + + // dim = 0 : S, for Hamiltonian + // dim = 1-3 : dS, for force + // dim = 4-6 : dS * dR, for stress + + void folding_matrix_k(const UnitCell& ucell, + const Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const int dim1, + const int dim2, + std::complex* mat_k, + const ModuleBase::Vector3& kvec_d); /** - * @brief calculate the HR local matrix of atom pair - */ - void cal_HR_IJR(const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& vu_in, - TR* data_pointer); + * @brief new function of folding_S_matrix + * only for Hamiltonian now, for force and stress will be developed later + * use HContainer as input and output in mat_k + */ + void folding_matrix_k_new(const int ik, + hamilt::Hamilt>* p_ham); + + //============================================================= + // In dftu_force.cpp + // For calculating force and stress fomr DFT+U + //============================================================= + public: + void force_stress(const UnitCell& ucell, + const Grid_Driver& gd, + std::vector>* dmk_d, + std::vector>>* dmk_c, + const Parallel_Orbitals& pv, + ForceStressArrays& fsr, + ModuleBase::matrix& force_dftu, + ModuleBase::matrix& stress_dftu, + const K_Vectors& kv, + const int npol); + + private: + void cal_force_k(const UnitCell& ucell, + const Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const std::complex* rho_VU, + ModuleBase::matrix& force_dftu, + const ModuleBase::Vector3& kvec_d); + + void cal_stress_k(const UnitCell& ucell, + const Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const std::complex* rho_VU, + ModuleBase::matrix& stress_dftu, + const ModuleBase::Vector3& kvec_d); + + void cal_force_gamma(const UnitCell& ucell, + const double* rho_VU, + const Parallel_Orbitals& pv, + double* dsloc_x, + double* dsloc_y, + double* dsloc_z, + ModuleBase::matrix& force_dftu); + + void cal_stress_gamma(const UnitCell& ucell, + const Parallel_Orbitals& pv, + const Grid_Driver* gd, + double* dsloc_x, + double* dsloc_y, + double* dsloc_z, + double* dh_r, + const double* rho_VU, + ModuleBase::matrix& stress_dftu); +#endif + + //============================================================= + // In dftu_yukawa.cpp + // Relevant for calculating U using Yukawa potential + //============================================================= + + public: + void cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx); + private: + void cal_slater_Fk(const UnitCell& ucell,const int L, const int T); // L:angular momnet, T:atom type + void cal_yukawa_lambda(double** rho, const int& nrxx); + + double spherical_Bessel(const int k, const double r, const double lambda); + double spherical_Hankel(const int k, const double r, const double lambda); + +#ifdef __LCAO + public: /** - * @brief calculate the atomic Force of atom pair - */ - void cal_force_IJR(const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& vu_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - double* force1, - double* force2); + * @brief get the density matrix of target spin + * nspin = 1 and 4 : ispin should be 0 + * nspin = 2 : ispin should be 0/1 + */ + const hamilt::HContainer* get_dmr(int ispin) const; /** - * @brief calculate the Stress of atom pair - */ - void cal_stress_IJR(const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const std::vector& vu_in, - const hamilt::BaseMatrix** dmR_pointer, - const int nspin, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress); - - std::vector adjs_all; - /// @brief if the nlm_tot is calculated - bool precal_nlm_done = false; - /// @brief the overlap values for all [atoms][nerghbors][orb_index(iw) in NAOs][m of target_l in Projectors] - std::vector>>> nlm_tot; + * @brief set the density matrix for DFT+U calculation + * if the density matrix is not set or set to nullptr, the DFT+U calculation will not be performed + */ + void set_dmr(const elecstate::DensityMatrix* dm_in_dftu_d); + void set_dmr(const elecstate::DensityMatrix, double>* dm_in_dftu_cd); + + private: + const UnitCell* ucell = nullptr; + const elecstate::DensityMatrix* dm_in_dftu_d = nullptr; + const elecstate::DensityMatrix, double>* dm_in_dftu_cd = nullptr; +#endif }; -} // namespace hamilt + +#ifdef __LCAO +template +void dftu_cal_occup_m(const int iter, + const UnitCell& ucell, + const std::vector>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt* p_ham, + Plus_U &dftu); +#endif + + #endif diff --git a/source/source_lcao/module_dftu/dftu_lcao_op.cpp b/source/source_lcao/module_dftu/dftu_lcao_op.cpp new file mode 100644 index 00000000000..9497cacabc0 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_lcao_op.cpp @@ -0,0 +1,691 @@ +#include "dftu_lcao_op.h" + +#include "source_base/timer.h" +#include "source_base/tool_title.h" +#include "source_cell/module_neighbor/sltk_grid_driver.h" +#include "source_lcao/module_operator_lcao/operator_lcao.h" +#include "source_hamilt/module_hcontainer/hcontainer_funcs.h" +#include "source_io/module_parameter/parameter.h" +#ifdef _OPENMP +#include +#endif +#include "source_base/parallel_reduce.h" + +template +hamilt::DFTU>::DFTU(HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + hamilt::HContainer* hR_in, + const UnitCell& ucell_in, + const Grid_Driver* GridD_in, + const TwoCenterIntegrator* intor, + const std::vector& orb_cutoff, + Plus_U* p_dftu) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor), orb_cutoff_(orb_cutoff) +{ + this->cal_type = calculation_type::lcao_dftu; + this->ucell = &ucell_in; + this->dftu = p_dftu; +#ifdef __DEBUG + assert(this->ucell != nullptr); +#endif + // initialize HR to allocate sparse Nonlocal matrix memory + this->initialize_HR(GridD_in); + // set nspin + this->nspin = PARAM.inp.nspin; +} + +// destructor +template +hamilt::DFTU>::~DFTU() +{ +} + +// initialize_HR() +template +void hamilt::DFTU>::initialize_HR(const Grid_Driver* GridD) +{ + ModuleBase::TITLE("DFTU", "initialize_HR"); + ModuleBase::timer::start("DFTU", "initialize_HR"); + + this->adjs_all.clear(); + this->adjs_all.reserve(this->ucell->nat); + for (int iat0 = 0; iat0 < ucell->nat; iat0++) + { + 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); + + AdjacentAtomInfo adjs; + GridD->Find_atom(*ucell, tau0, T0, I0, &adjs); + std::vector is_adj(adjs.adj_num + 1, false); + 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); + const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; + const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + // choose the real adjacent atoms + // Note: the distance of atoms should less than the cutoff radius, + // When equal, the theoretical value of matrix element is zero, + // but the calculated value is not zero due to the numerical error, which would lead to result changes. + if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 + < orb_cutoff_[T1] + PARAM.inp.onsite_radius) + { + is_adj[ad1] = true; + } + } + filter_adjs(is_adj, adjs); + this->adjs_all.push_back(adjs); + } + + ModuleBase::timer::end("DFTU", "initialize_HR"); +} + +template +void hamilt::DFTU>::cal_nlm_all(const Parallel_Orbitals* paraV) +{ + ModuleBase::TITLE("DFTU", "cal_nlm_all"); + if (this->precal_nlm_done) + { + return; + } + + ModuleBase::timer::start("DFTU", "cal_nlm_all"); + nlm_tot.resize(this->ucell->nat); + const int npol = this->ucell->get_npol(); + int atom_index = 0; + for (int iat0 = 0; iat0 < ucell->nat; iat0++) + { + 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++]; + + // calculate and save the table of two-center integrals + nlm_tot[iat0].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 = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->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; + // only first zeta orbitals in target L of atom iat0 are needed + std::vector nlm_target(tlp1); + const int L1 = atom1->iw2l[iw1]; + const int N1 = atom1->iw2n[iw1]; + const int m1 = atom1->iw2m[iw1]; + 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) + + // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) + const 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, false /*cal_deri*/, nlm); + // select the elements of nlm with target_L + 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 < 2 * L0 + 1; m++) + { + nlm_target[m] = nlm[0][iw + m]; + } + break; + } + } + nlm_tot[iat0][ad].insert({all_indexes[iw1l], nlm_target}); + } + } + } + this->precal_nlm_done = true; + ModuleBase::timer::end("DFTU", "cal_nlm_all"); +} + +// contributeHR() +/** + * @brief Contribute DFT+U Hamiltonian to real-space HR matrix + * + * @details This function handles different scenarios based on: + * 1. Whether occ_mat (occupation matrix) is read from file (is_occ_mat_initialized) + * 2. Spin configuration (nspin=1, 2, or 4) + * 3. SCF iteration stage (first vs subsequent iterations) + * + * Case 1: Occ_mat NOT initialized (!is_occ_mat_initialized) + * - First electronic iteration: calculates occupation matrix from density matrix (DMR) + * * Uses get_dmr(current_spin) to get real-space density matrix + * * Accumulates contributions from all atom pairs via cal_occ() + * * Performs MPI reduction to sum occ across processes + * * Stores result via set_occ_mat_flat() for use in VU calculation + * * For nspin=1: occ is scaled by 0.5 (since only one spin channel computed) + * - Subsequent iterations: occ_mat is computed fresh each iteration from updated DMR + * + * Case 2: Occ_mat IS initialized (is_occ_mat_initialized, i.e., read from dm_onsite.txt file) + * - First electronic iteration: uses pre-read occ_mat directly without DMR calculation + * * Skips DMR-based occ calculation entirely + * * Reads locale from stored data via get_occ_mat() + * * Different indexing for nspin=4 vs nspin=1/2 (see below) + * - After first iteration: mark_occ_mat_dirty() is called to force recomputation + * + * Spin configurations: + * nspin=1 (non-spin-polarized): + * - Single spin channel, occ computed once + * - Energy correction doubled at end (set_double_energy) + * - current_spin always 0 + * + * nspin=2 (collinear spin-polarized): + * - Two separate spin channels (spin-up: 0, spin-down: 1) + * - current_spin toggles between 0 and 1 across iterations + * - mark_occ_mat_dirty() called when current_spin == 1 (last spin) + * - HR accumulated separately for each spin + * + * nspin=4 (non-collinear/SOC): + * - Single 4x4 Pauli matrix representation per atom + * - occ has 4*(2l+1)^2 elements (spin_fold=4) + * - get_occ_mat uses spin=0, ipol indices for Pauli blocks + * - mark_occ_mat_dirty() always called (current_spin check always true) + * - No current_spin toggling (all spins handled simultaneously) + * + * @warning THREAD SAFETY: cal_HR_IJR() updates shared HR matrix entries. + * Different iat0 may contribute to same HR(iat1, iat2, R), requiring + * critical section protection for multithreaded correctness. + * TODO: Consider refactoring to atom_row_list pattern (see nonlocal.cpp) + * for better parallel performance instead of critical section. + */ +template +void hamilt::DFTU>::contributeHR() +{ + ModuleBase::TITLE("DFTU", "contributeHR"); + // Early exit conditions: + // - get_dmr(0) == nullptr: DMR not available (typical in first iteration without file input) + // - !is_occ_mat_initialized(): occ_mat not read from file AND not yet computed from DMR + // When both true, skip DFT+U contribution entirely (first iteration, no file input) + const bool dmr_null = (this->dftu->get_dmr(0) == nullptr); + const bool occ_mat_not_init = !this->dftu->is_occ_mat_initialized(); + + if (dmr_null && occ_mat_not_init) + { + return; + } + else + { + // Reset DFT+U energy at start of each spin cycle + // For nspin=2: reset when current_spin==0 (start of spin-up calculation) + // For nspin=1/4: reset once (current_spin always 0) + if (this->current_spin == 0) + { + this->dftu->set_energy(0.0); + } + } + ModuleBase::timer::start("DFTU", "contributeHR"); + + const Parallel_Orbitals* paraV = this->hR->get_atom_pair(0).get_paraV(); + const int npol = this->ucell->get_npol(); + // 1. Calculate two-center integrals for all atom pairs + // This is reused in both occ and HR calculations + this->cal_nlm_all(paraV); + + // 2. Loop over all Hubbard-projector center atoms (iat0) + int atom_index = 0; + for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + { + auto tau0 = ucell->get_tau(iat0); + int T0, I0; + 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++]; + + ModuleBase::timer::start("DFTU", "cal_occ"); + // spin_fold: number of spin components in occ array + // nspin=4: 4 (Pauli matrix blocks), nspin=1/2: 1 (single spin channel) + const int spin_fold = (this->nspin == 4) ? 4 : 1; + std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); + + // ============================================================ + // BRANCH 1: Occ_mat NOT initialized (compute from DMR) + // ============================================================ + // This branch is taken when: + // - is_occ_mat_initialized() == false (no file read or omc != 0) + // - DMR is available (get_dmr() != nullptr) + // Typical scenario: normal SCF iterations after first step + if (!this->dftu->is_occ_mat_initialized()) + { + // TODO: UNSAFE - get_dmr(current_spin) assumes DMR has correct spin indexing. + // For nspin=2, current_spin must be correctly toggled (0 then 1). + // If current_spin is wrong, wrong spin channel's DMR is used. + const hamilt::HContainer* dmR_current = this->dftu->get_dmr(this->current_spin); + 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); + ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + const std::unordered_map>& nlm1 = nlm_tot[iat0][ad1]; + 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); + const std::unordered_map>& nlm2 = nlm_tot[iat0][ad2]; + ModuleBase::Vector3& R_index2 = adjs.box[ad2]; + 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 + = dmR_current->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); + if (tmp != nullptr) + { + this->cal_occ(iat1, iat2, paraV, nlm1, nlm2, tmp->get_pointer(), occ); + } + } + } +#ifdef __MPI + // CRITICAL: MPI reduction required for distributed DMR calculations. + // Each process computes partial occ from its local DMR blocks. + // Without this, occ would be incomplete and DFT+U potential wrong. + // TODO: Verify that occ size is consistent across processes before reduction. + // TODO: Consider using MPI_IN_PLACE to avoid extra buffer allocation. + Parallel_Reduce::reduce_all(occ.data(), occ.size()); +#endif + // For nspin=1: occ computed from single spin channel, but should represent + // total occupation (both spins). Scale by 0.5 to account for this. + if (this->nspin == 1) + { + for (auto& v : occ) { v *= 0.5; } + } + this->dftu->set_occ_mat_flat(iat0, target_L, this->current_spin, occ); + } + // ============================================================ + // BRANCH 2: Occ_mat IS initialized (use pre-read data) + // ============================================================ + // This branch is taken when: + // - is_occ_mat_initialized() == true (occ_mat read from dm_onsite.txt file) + // - OR omc != 0 (occupation matrix control with dm_onsite_ini.txt) + // Typical scenario: first SCF iteration with file input, or restart calculation + else + { + // nspin=4: Non-collinear case with Pauli matrix representation + // Occ_mat stored as single 4x4 block per atom, with spin indices embedded + // in the matrix indices (ipol0, ipol1 for Pauli block indices) + if (this->nspin == 4) + { + // 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); + } + // nspin=1 or nspin=2: Collinear spin case + // Occ_mat stored separately for each spin channel + else + { + for (int i = 0; i < static_cast(occ.size()); i++) + { + // 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, + i / (2 * target_L + 1), i % (2 * target_L + 1)); + } + } + } + ModuleBase::timer::end("DFTU", "cal_occ"); + + // 3. Calculate Hubbard potential VU from occupation matrix + // VU = U * (1/2 * delta(m,m') - occ(m,m')) for each spin channel + // Energy: EU = U * 1/2 * occ(m,m') * occ(m',m) + ModuleBase::timer::start("DFTU", "cal_vu"); + const double u_value = this->dftu->u_current[T0]; + std::vector VU_tmp(occ.size()); + + // mohan update 2025-11: get_energy/set_energy are now instance methods + // via this->dftu pointer, no longer global static state. + double u_energy = this->dftu->get_energy(); + this->cal_v_of_u(occ, tlp1, u_value, VU_tmp.data(), u_energy); + this->dftu->set_energy(u_energy); + + // 4. Convert VU to appropriate data type (real or complex) + // For nspin=4 with complex Hamiltonian, VU needs Pauli matrix transformation + std::vector VU(occ.size()); + this->transfer_vu(VU_tmp, VU); + + // 5. Second iteration: Calculate Hamiltonian matrix contribution + // HR += * VU(m,m') * + // for all atom pairs within cutoff + // Note: different iat0 may contribute to the same HR(iat1, iat2, R), so we need to protect the update + // to avoid race conditions in multithreading. Reference: nonlocal.cpp for the atom_row_list pattern. + // TODO: CRITICAL SECTION PERFORMANCE - This critical section serializes HR updates. + // For systems with many Hubbard atoms, this becomes a bottleneck. + // Consider refactoring to atom_row_list pattern (see nonlocal.cpp lines 127-220): + // 1. Use #pragma omp for to distribute iat0 across threads + // 2. Each thread records its assigned iat0 in thread-local atom_row_list + // 3. When updating HR(iat1, iat2, R), skip if iat1 not in thread's atom_row_list + // 4. This eliminates race conditions without critical section + 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); + ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + const std::unordered_map>& nlm1 = nlm_tot[iat0][ad1]; + 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); + const std::unordered_map>& nlm2 = nlm_tot[iat0][ad2]; + ModuleBase::Vector3& R_index2 = adjs.box[ad2]; + ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], + R_index2[1] - R_index1[1], + R_index2[2] - R_index1[2]); + hamilt::BaseMatrix* tmp = this->hR->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); + if (tmp != nullptr) + { +#ifdef _OPENMP +#pragma omp critical(dftu_hr_update) +#endif + { + this->cal_HR_IJR(iat1, iat2, paraV, nlm1, nlm2, VU, tmp->get_pointer()); + } + } + } + } + ModuleBase::timer::end("DFTU", "cal_vu"); + } + + // 6. Post-processing: Energy correction and occ_mat state management + // For nspin=1: DFT+U energy computed for single spin channel, but should count both spins + // set_double_energy() doubles the energy to account for degenerate spin-up/down + if (this->nspin == 1) + { + this->dftu->set_double_energy(); + } + + // 7. Mark occ_mat as dirty to force recomputation in next iteration + // This is called when: + // - nspin=4: Always (all spins handled simultaneously, current_spin==0==nspin-1) + // - nspin=2: When current_spin==1 (after spin-down calculation, last spin channel) + // - nspin=1: When current_spin==0==nspin-1 (always called) + // + // Purpose: Ensure occ_mat is recomputed from updated DMR in next SCF iteration, + // rather than using stale pre-read data from file. + // TODO: This logic is confusing. Consider explicit variable like `is_last_spin_channel`. + if (this->current_spin == this->nspin - 1 || this->nspin == 4) + { + this->dftu->mark_occ_mat_dirty(); + } + + // 8. Spin channel toggling for nspin=2 + // nspin=2 requires separate HR updates for spin-up (current_spin=0) and spin-down (current_spin=1) + // The HR matrix is updated twice per SCF iteration, once for each spin channel + // current_spin toggles: 0 -> 1 -> 0 -> 1 ... + // For nspin=1: current_spin always 0 (no toggling needed) + // For nspin=4: current_spin always 0 (all spins handled simultaneously via Pauli matrices) + // TODO: UNSAFE - This assumes contributeHR() is called in strict alternating order. + // If called out of order (e.g., due to parallel k-point distribution), current_spin may be wrong. + // TODO: Consider deriving current_spin from ik or explicit parameter instead of toggling. + if (this->nspin == 2) + { + this->current_spin = 1 - this->current_spin; + } + + ModuleBase::timer::end("DFTU", "contributeHR"); +} + +// cal_HR_IJR() +template +void hamilt::DFTU>::cal_HR_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& VU, + TR* data_pointer) +{ + + // 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 = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->get_indexes_col(iat2); + const int m_size = int(sqrt(VU.size()) / npol); + // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 + std::vector step_trace(npol * npol, 0); + for (int is = 0; is < npol; is++) + { + for (int is2 = 0; is2 < npol; is2++) + { + step_trace[is * npol + is2] = paraV->get_ncol_atom(iat2) * is + is2; + } + } + // calculate the local matrix + const TR* tmp_d = nullptr; + 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 is = 0; is < npol * npol; ++is) + { + int start = is * m_size * m_size; + TR nlm_tmp = TR(0); + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + nlm_tmp += nlm1[m1] * nlm2[m2] * VU[m1 * m_size + m2 + start]; + } + } + data_pointer[step_trace[is]] += nlm_tmp; + } + data_pointer += npol; + } + data_pointer += (npol - 1) * col_indexes.size(); + } +} + +template +void hamilt::DFTU>::cal_occ(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const double* dm_pointer, + std::vector& occ) +{ + + // 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 = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->get_indexes_col(iat2); + const int m_size = int(sqrt(occ.size()) / npol); + const int m_size2 = m_size * m_size; +#ifdef __DEBUG + assert(m_size * m_size == occ.size()); +#endif + // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 + std::vector step_trace(npol * npol, 0); + for (int is = 0; is < npol; is++) + { + for (int is2 = 0; is2 < npol; is2++) + { + step_trace[is * npol + is2] = paraV->get_ncol_atom(iat2) * is + is2; + } + } + // calculate the local matrix + 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 is1 = 0; is1 < npol; ++is1) + { + for (int is2 = 0; is2 < npol; ++is2) + { + for (int m1 = 0; m1 < m_size; ++m1) + { + for (int m2 = 0; m2 < m_size; ++m2) + { + occ[m1 * m_size + m2 + (is1 * npol + is2) * m_size2] + += nlm1[m1] * nlm2[m2] * dm_pointer[step_trace[is1 * npol + is2]]; + } + } + } + } + dm_pointer += npol; + } + dm_pointer += (npol - 1) * col_indexes.size(); + } +} + +template +void hamilt::DFTU>::transfer_vu(std::vector& vu_tmp, std::vector& vu) +{ +#ifdef __DEBUG + assert(vu.size() == vu_tmp.size()); +#endif + for (int i = 0; i < vu_tmp.size(); i++) + { + vu[i] = vu_tmp[i]; + } +} + +template <> +void hamilt::DFTU, std::complex>>::transfer_vu( + std::vector& vu_tmp, + std::vector>& vu) +{ +#ifdef __DEBUG + assert(vu.size() == vu_tmp.size()); +#endif + + // Pauli-to-spinor conversion for DFT+U potential: + // V = V_0*I + V_x*sigma_x + V_y*sigma_y + V_z*sigma_z + // sigma_y = [[0,-i],[i,0]], so: + // V_{up,up} = 0.5*(V_0 + V_z) + // V_{down,down} = 0.5*(V_0 - V_z) + // V_{up,down} = 0.5*(V_x - i*V_y) <- note: minus sign from sigma_y + // V_{down,up} = 0.5*(V_x + i*V_y) <- note: plus sign from sigma_y + // This is consistent with the convention in gint_common.cpp merge_hr_part_to_hR(). + const int m_size = int(sqrt(vu.size()) / 2); + const int m_size2 = m_size * m_size; + vu.resize(vu_tmp.size()); + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + int index[4]; + index[0] = m1 * m_size + m2; + index[1] = m1 * m_size + m2 + m_size2; + index[2] = m2 * m_size + m1 + m_size2 * 2; + index[3] = m2 * m_size + m1 + m_size2 * 3; + vu[index[0]] = 0.5 * (vu_tmp[index[0]] + vu_tmp[index[3]]); + vu[index[3]] = 0.5 * (vu_tmp[index[0]] - vu_tmp[index[3]]); + vu[index[1]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); + vu[index[2]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); + } + } +} + +template +void hamilt::DFTU>::cal_v_of_u(const std::vector& occ, + const int m_size, + const double u_value, + double* vu, + double& eu) +{ + // calculate the local matrix + int spin_fold = occ.size() / m_size / m_size; + if (spin_fold < 4) { + for (int is = 0; is < spin_fold; ++is) + { + int start = is * m_size * m_size; + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + vu[start + m1 * m_size + m2] = u_value * (0.5 * (m1 == m2) - occ[start + m2 * m_size + m1]); + eu += u_value * 0.5 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; + } + } + } + } else + { + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + vu[m1 * m_size + m2] = u_value * (1.0 * (m1 == m2) - occ[m2 * m_size + m1]); + eu += u_value * 0.25 * occ[m2 * m_size + m1] * occ[m1 * m_size + m2]; + } + } + for (int is = 1; is < spin_fold; ++is) + { + int start = is * m_size * m_size; + for (int m1 = 0; m1 < m_size; m1++) + { + for (int m2 = 0; m2 < m_size; m2++) + { + vu[start + m1 * m_size + m2] = u_value * (0 - occ[start + m2 * m_size + m1]); + eu += u_value * 0.25 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; + } + } + } + } +} + +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_lcao_op.h new file mode 100644 index 00000000000..23b361bbbe1 --- /dev/null +++ b/source/source_lcao/module_dftu/dftu_lcao_op.h @@ -0,0 +1,155 @@ +#ifndef DFTPLUSU_H +#define DFTPLUSU_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_estate/module_dm/density_matrix.h" +#include "source_lcao/module_operator_lcao/operator_lcao.h" +#include "source_lcao/module_dftu/dftu_lcao.h" +#include "source_hamilt/module_hcontainer/hcontainer.h" + +#include + +namespace hamilt +{ + +#ifndef __DFTUTEMPLATE +#define __DFTUTEMPLATE + +/// The DFTU class template inherits from class T +/// it is used to calculate the non-local pseudopotential of wavefunction basis +/// Template parameters: +/// - T: base class, it would be OperatorLCAO or OperatorPW +template +class DFTU : public T +{ +}; + +#endif + +/// DFTU class template specialization for OperatorLCAO base class +/// It is used to calculate the non-local pseudopotential matrix in real space and fold it to k-space +/// HR = D_{p1, p2} +/// HK = D_{p1, p2} = \sum_{R} e^{ikR} HR +/// Template parameters: +/// - TK: data type of k-space Hamiltonian +/// - TR: data type of real space Hamiltonian +template +class DFTU> : public OperatorLCAO +{ + public: + DFTU>(HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + hamilt::HContainer* hR_in, + const UnitCell& ucell_in, + const Grid_Driver* gridD_in, + const TwoCenterIntegrator* intor, + const std::vector& orb_cutoff, + Plus_U* p_dftu); + ~DFTU>(); + + /** + * @brief contributeHR() is used to calculate the HR matrix + * D_{p1, p2} + */ + virtual void contributeHR() override; + + /// calculate force and stress for DFT+U + void cal_force_stress(const bool cal_force, + const bool cal_stress, + ModuleBase::matrix& force, + ModuleBase::matrix& stress); + + private: + const UnitCell* ucell = nullptr; + + Plus_U* dftu = nullptr; + + hamilt::HContainer* HR = nullptr; + + const TwoCenterIntegrator* intor_ = nullptr; + + std::vector orb_cutoff_; + + /// @brief the number of spin components, 1 for no-spin, 2 for collinear spin case and 4 for non-collinear spin case + int nspin = 0; + + /** + * @brief search the nearest neighbor atoms and save them into this->adjs_all + * the size of HR will not change in DFTU, + * because I don't want to expand HR larger than Nonlocal operator caused by DFTU + */ + void initialize_HR(const Grid_Driver* gridD_in); + + /** + * @brief calculate the overlap values and save them in this->nlm_tot + * it will be reused in the calculation of calculate_HR() + */ + void cal_nlm_all(const Parallel_Orbitals* paraV); + + /** + * @brief calculate the occ_mm' = \sum_R DMR* matrix for each atom to add U + */ + void cal_occ(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const double* data_pointer, + std::vector& occupations); + + /// transfer VU format from pauli matrix to normal for non-collinear spin case + void transfer_vu(std::vector& vu_tmp, std::vector& vu); + /// VU_{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_v_of_u(const std::vector& occ, const int m_size, const double u_value, double* vu, double& eu); + + /** + * @brief calculate the HR local matrix of atom pair + */ + void cal_HR_IJR(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + TR* data_pointer); + + /** + * @brief calculate the atomic Force of atom pair + */ + void cal_force_IJR(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, + double* force2); + /** + * @brief calculate the Stress of atom pair + */ + void cal_stress_IJR(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + + std::vector adjs_all; + /// @brief if the nlm_tot is calculated + bool precal_nlm_done = false; + /// @brief the overlap values for all [atoms][nerghbors][orb_index(iw) in NAOs][m of target_l in Projectors] + std::vector>>> nlm_tot; +}; + +} // namespace hamilt +#endif diff --git a/source/source_lcao/module_dftu/op_dftu_lcao.cpp b/source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp similarity index 96% rename from source/source_lcao/module_dftu/op_dftu_lcao.cpp rename to source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp index 65198b46934..643a42f677e 100644 --- a/source/source_lcao/module_dftu/op_dftu_lcao.cpp +++ b/source/source_lcao/module_dftu/dftu_lcao_op_legacy.cpp @@ -1,7 +1,7 @@ -#include "op_dftu_lcao.h" +#include "dftu_lcao_op_legacy.h" #include "source_base/timer.h" #include "source_base/tool_title.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/module_dftu/dftu_lcao.h" namespace hamilt { diff --git a/source/source_lcao/module_dftu/op_dftu_lcao.h b/source/source_lcao/module_dftu/dftu_lcao_op_legacy.h similarity index 94% rename from source/source_lcao/module_dftu/op_dftu_lcao.h rename to source/source_lcao/module_dftu/dftu_lcao_op_legacy.h index 99f4fae8b6d..8bf8e8260d3 100644 --- a/source/source_lcao/module_dftu/op_dftu_lcao.h +++ b/source/source_lcao/module_dftu/dftu_lcao_op_legacy.h @@ -3,7 +3,7 @@ #include "source_base/timer.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 namespace hamilt { diff --git a/source/source_lcao/module_dftu/dftu_occup.cpp b/source/source_lcao/module_dftu/dftu_occup.cpp index 782b6b1cb14..ff9c5529eea 100644 --- a/source/source_lcao/module_dftu/dftu_occup.cpp +++ b/source/source_lcao/module_dftu/dftu_occup.cpp @@ -1,4 +1,4 @@ -#include "dftu.h" +#include "dftu_lcao.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" #ifdef __LCAO @@ -41,7 +41,7 @@ void Plus_U::cal_occup_m_k(const int iter, std::complex* s_k_pointer = nullptr; - if(Plus_U::nspin != 4) + if(this->nspin != 4) { s_k_pointer = dynamic_cast, double>*>(p_ham)->getSk(); } @@ -180,7 +180,7 @@ void Plus_U::cal_occup_m_k(const int iter, // set the local occupation mumber matrix of spin up and down zeros #ifdef __MPI - if (Plus_U::nspin == 1 || Plus_U::nspin == 4) + if (this->nspin == 1 || this->nspin == 4) { ModuleBase::matrix temp(occ_mat[iat][l][n][0]); MPI_Allreduce(&temp(0, 0), @@ -190,7 +190,7 @@ void Plus_U::cal_occup_m_k(const int iter, MPI_SUM, MPI_COMM_WORLD); } - else if (Plus_U::nspin == 2) + else if (this->nspin == 2) { ModuleBase::matrix temp0(occ_mat[iat][l][n][0]); MPI_Allreduce(&temp0(0, 0), @@ -210,7 +210,7 @@ void Plus_U::cal_occup_m_k(const int iter, } #endif - switch (Plus_U::nspin) + switch (this->nspin) { case 1: occ_mat[iat][l][n][0] += transpose(occ_mat[iat][l][n][0]); @@ -219,7 +219,7 @@ void Plus_U::cal_occup_m_k(const int iter, break; case 2: - for (int is = 0; is < Plus_U::nspin; is++) + for (int is = 0; is < this->nspin; is++) occ_mat[iat][l][n][is] += transpose(occ_mat[iat][l][n][is]); break; @@ -264,7 +264,7 @@ void Plus_U::cal_occup_m_gamma(const int iter, const double alpha = 1.0, beta = 0.0; std::vector srho(this->paraV->nloc); - for (int is = 0; is < Plus_U::nspin; is++) + for (int is = 0; is < this->nspin; is++) { double* s_gamma_pointer = dynamic_cast*>(p_ham)->getSk(); @@ -370,7 +370,7 @@ void Plus_U::cal_occup_m_gamma(const int iter, #endif // for the case spin independent calculation - switch (Plus_U::nspin) + switch (this->nspin) { case 1: occ_mat[iat][l][n][0] += transpose(occ_mat[iat][l][n][0]); diff --git a/source/source_lcao/module_dftu/dftu_tools.cpp b/source/source_lcao/module_dftu/dftu_tools.cpp index ab03942e64d..6e86daeccf7 100644 --- a/source/source_lcao/module_dftu/dftu_tools.cpp +++ b/source/source_lcao/module_dftu/dftu_tools.cpp @@ -1,4 +1,4 @@ -#include "dftu.h" +#include "dftu_lcao.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" @@ -10,7 +10,7 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool new_occ_mat, std: for (int it = 0; it < this->ucell->ntype; ++it) { - if (Plus_U::orbital_corr[it] == -1) + if (this->orbital_corr[it] == -1) { continue; } @@ -19,7 +19,7 @@ void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool new_occ_mat, std: const int iat = this->ucell->itia2iat(it, ia); for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != Plus_U::orbital_corr[it]) + if (L != this->orbital_corr[it]) { continue; } @@ -74,7 +74,7 @@ void Plus_U::cal_VU_pot_mat_real(const int spin, const bool new_occ_mat, double* for (int it = 0; it < this->ucell->ntype; ++it) { - if (Plus_U::orbital_corr[it] == -1) + if (this->orbital_corr[it] == -1) { continue; } @@ -83,7 +83,7 @@ void Plus_U::cal_VU_pot_mat_real(const int spin, const bool new_occ_mat, double* const int iat = this->ucell->itia2iat(it, ia); for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != Plus_U::orbital_corr[it]) + if (L != this->orbital_corr[it]) { continue; } @@ -158,7 +158,7 @@ double Plus_U::get_onebody_eff_pot(const int T, case 3: // simplified formalism and FLL double counting if (new_occ_mat) { - if (Yukawa) + if (use_yukawa) { if (m0 == m1) { @@ -171,15 +171,15 @@ double Plus_U::get_onebody_eff_pot(const int T, else { if (m0 == m1) { - VU = (this->U[T]) * (0.5 - this->occ_mat[iat][L][N][spin](m0, m1)); + VU = (this->u_current[T]) * (0.5 - this->occ_mat[iat][L][N][spin](m0, m1)); } else { - VU = -(this->U[T]) * this->occ_mat[iat][L][N][spin](m0, m1); + VU = -(this->u_current[T]) * this->occ_mat[iat][L][N][spin](m0, m1); } } } else { - if (Yukawa) + if (use_yukawa) { if (m0 == m1) { VU = (this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N]) @@ -192,9 +192,9 @@ double Plus_U::get_onebody_eff_pot(const int T, else { if (m0 == m1) { - VU = (this->U[T]) * (0.5 - this->occ_mat_save[iat][L][N][spin](m0, m1)); + VU = (this->u_current[T]) * (0.5 - this->occ_mat_save[iat][L][N][spin](m0, m1)); } else { - VU = -(this->U[T]) * this->occ_mat_save[iat][L][N][spin](m0, m1); + VU = -(this->u_current[T]) * this->occ_mat_save[iat][L][N][spin](m0, m1); } } } diff --git a/source/source_lcao/module_dftu/dftu_yukawa.cpp b/source/source_lcao/module_dftu/dftu_yukawa.cpp index a5909267b18..e743909243f 100644 --- a/source/source_lcao/module_dftu/dftu_yukawa.cpp +++ b/source/source_lcao/module_dftu/dftu_yukawa.cpp @@ -2,7 +2,7 @@ #include "source_io/module_parameter/parameter.h" #include "source_base/constants.h" #include "source_base/global_function.h" -#include "dftu.h" +#include "dftu_lcao.h" #include #include @@ -26,9 +26,9 @@ void Plus_U::cal_yukawa_lambda(double** rho, const int& nrxx) double sum_rho = 0.0; double sum_rho_lambda = 0.0; - for (int is = 0; is < Plus_U::nspin; is++) + for (int is = 0; is < this->nspin; is++) { - if(Plus_U::nspin == 4 && is > 0) + if(this->nspin == 4 && is > 0) { continue;// for non-collinear spin case, first spin contains the charge density } @@ -64,7 +64,7 @@ void Plus_U::cal_slater_Fk(const UnitCell& ucell, { ModuleBase::TITLE("Plus_U", "cal_slater_Fk"); - if (Yukawa) + if (use_yukawa) { for (int chi = 0; chi < ucell.atoms[T].l_nchi[L]; chi++) { @@ -111,7 +111,7 @@ void Plus_U::cal_slater_Fk(const UnitCell& ucell, void Plus_U::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) { ModuleBase::TITLE("Plus_U", "cal_slater_UJ"); - if (!Yukawa) + if (!use_yukawa) { return; } @@ -140,9 +140,9 @@ void Plus_U::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) { const int N = ucell.atoms[T].l_nchi[L]; - if (L >= Plus_U::get_orbital_corr(T) && Plus_U::get_orbital_corr(T) != -1) + if (L >= this->get_orbital_corr(T) && this->get_orbital_corr(T) != -1) { - if (L != Plus_U::get_orbital_corr(T)) + if (L != this->get_orbital_corr(T)) { continue; } @@ -171,7 +171,7 @@ void Plus_U::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) this->U_Yukawa[T][L][0] *= 2.0; this->J_Yukawa[T][L][0] *= 2.0; // update current U with calculated U-J from Slater integrals - this->U[T] = this->U_Yukawa[T][L][0] - this->J_Yukawa[T][L][0]; + this->u_current[T] = this->U_Yukawa[T][L][0] - this->J_Yukawa[T][L][0]; } // end if } // end L } // end T diff --git a/source/source_lcao/module_dftu/test/CMakeLists.txt b/source/source_lcao/module_dftu/test/CMakeLists.txt index e18eebb6982..506ee79523c 100644 --- a/source/source_lcao/module_dftu/test/CMakeLists.txt +++ b/source/source_lcao/module_dftu/test/CMakeLists.txt @@ -22,7 +22,7 @@ if(ENABLE_LCAO AND ENABLE_MPI) AddTest( TARGET dftu_lcao_test LIBS parameter psi base device container - SOURCES dftu_lcao_test.cpp ../dftu_lcao.cpp ../dftu_fs.cpp + SOURCES dftu_lcao_test.cpp ../dftu_lcao_op.cpp ../dftu_fs.cpp ../../../source_pw/module_pwdft/dftu_base.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.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 090a3998243..d7826090397 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.h" -#include "source_lcao/module_dftu/dftu.h" +#include "../dftu_lcao_op.h" +#include "source_lcao/module_dftu/dftu_lcao.h" Plus_U::Plus_U(){}; Plus_U::~Plus_U(){}; @@ -97,8 +97,8 @@ class DFTUTest : public ::testing::Test dftu.occ_mat[iat][l][0][1].create(2 * l + 1, 2 * l + 1); } } - Plus_U::U = {U_test}; - Plus_U::orbital_corr = {orbital_c_test}; + dftu.u_current = {U_test}; + dftu.orbital_corr = {orbital_c_test}; PARAM.input.onsite_radius = 1.0; } diff --git a/source/source_lcao/setup_dftu_lcao.cpp b/source/source_lcao/setup_dftu_lcao.cpp index a56449015a5..307fb213799 100644 --- a/source/source_lcao/setup_dftu_lcao.cpp +++ b/source/source_lcao/setup_dftu_lcao.cpp @@ -1,5 +1,5 @@ #include "setup_dftu_lcao.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/module_dftu/dftu_lcao.h" #include "source_pw/module_pwdft/dftu_output.h" // mohan add 2025-11-08 #include "source_estate/module_dm/density_matrix.h" #include "source_lcao/hamilt_lcao.h" @@ -61,7 +61,7 @@ void finish_dftu_lcao(const int iter, /// new DFT+U method calculates energy in Hamiltonian if (dft_plus_u == 2) { - if (dftu_ptr->omc != 2) + if (dftu_ptr->occ_mat_ctrl != 2) { dftu_cal_occup_m(iter, ucell, dm_vec, kv, mixing_beta, static_cast*>(hamilt_lcao_ptr), *dftu_ptr); diff --git a/source/source_lcao/spar_u.h b/source/source_lcao/spar_u.h index bdb71326ca4..b4456dda4e8 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.h" // mohan add 20251107 +#include "source_lcao/module_dftu/dftu_lcao.h" // mohan add 20251107 namespace sparse_format { diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index a0c99e7b439..abc2d918e82 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -17,24 +17,7 @@ // 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. // mohan refactored 2025-11-08 - -// static member definitions (mohan add 2025-11-06) -double Plus_U_Base::energy_u = 0.0; - -std::vector Plus_U_Base::U = {}; - -std::vector Plus_U_Base::U0 = {}; - -std::vector Plus_U_Base::orbital_corr = {}; - -double Plus_U_Base::uramping = 0.0; - -int Plus_U_Base::omc = 0; - -int Plus_U_Base::mixing_dftu = 0; -int Plus_U_Base::nspin = 0; - -bool Plus_U_Base::Yukawa = false; +// All members are now non-static; default values are in the header. Plus_U_Base::Plus_U_Base() @@ -52,17 +35,15 @@ 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, - const int nlocal, - const bool gamma_only_local, - const std::string& ks_solver, - const bool cal_force, - const bool cal_stress, const std::string& device, - const int kpar) + const int kpar, + const std::vector& hubbard_u, + const double uramping, + const int occ_mat_ctrl, + const int mixing_dftu) { ModuleBase::TITLE("Plus_U_Base", "init_base"); @@ -71,25 +52,24 @@ void Plus_U_Base::init_base(UnitCell& cell, exit(0); #endif - Plus_U_Base::nspin = nspin; - Plus_U_Base::orbital_corr = orbital_corr; - Plus_U_Base::Yukawa = yukawa_potential; - this->yukawa_lambda = yukawa_lambda; - - this->global_readin_dir = global_readin_dir; - this->global_out_dir = global_out_dir; - this->init_chg = init_chg; - this->npol = npol; - - this->nlocal = nlocal; - this->gamma_only_local = gamma_only_local; - this->ks_solver = ks_solver; - this->cal_force = cal_force; - this->cal_stress = cal_stress; + 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; + this->u_target = hubbard_u; + this->u_current = hubbard_u; + if (uramping > 0.01) + { + std::fill(this->u_current.begin(), + this->u_current.end(), + 0.0); + } this->device = device; this->kpar = kpar; - Plus_U_Base::energy_u = 0.0; + this->energy_u = 0.0; this->occ_mat.resize(cell.nat); this->occ_mat_save.resize(cell.nat); @@ -198,7 +178,7 @@ 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 (Yukawa) + if (use_yukawa) { this->Fk.resize(cell.ntype); @@ -229,11 +209,11 @@ void Plus_U_Base::init_base(UnitCell& cell, } } - if (omc != 0) + if (occ_mat_ctrl != 0) { std::stringstream sst; - sst << this->global_readin_dir << "dm_onsite_ini.txt"; - this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol); + sst << global_readin_dir << "dm_onsite_ini.txt"; + this->read_occup_m(cell, sst.str(), init_chg, nspin, npol); #ifdef __MPI this->local_occup_bcast(cell, nspin, npol); #endif @@ -243,11 +223,11 @@ void Plus_U_Base::init_base(UnitCell& cell, } else { - if (this->init_chg == "file") + if (init_chg == "file") { std::stringstream sst; - sst << this->global_readin_dir << "dm_onsite.txt"; - this->read_occup_m(cell, sst.str(), this->init_chg, nspin, npol); + sst << global_readin_dir << "dm_onsite.txt"; + this->read_occup_m(cell, sst.str(), init_chg, nspin, npol); #ifdef __MPI this->local_occup_bcast(cell, nspin, npol); #endif @@ -267,7 +247,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 (Yukawa) { + if (use_yukawa) { return; } // if uramping < 0.1, use the original U @@ -275,15 +255,15 @@ void Plus_U_Base::uramping_update() return; } // loop to change U - for (int i = 0; i < static_cast(this->U0.size()); i++) + for (int i = 0; i < static_cast(this->u_target.size()); i++) { - if (this->U[i] + this->uramping < this->U0[i]) + if (this->u_current[i] + this->uramping < this->u_target[i]) { - this->U[i] += this->uramping; + this->u_current[i] += this->uramping; } else { - this->U[i] = this->U0[i]; + this->u_current[i] = this->u_target[i]; } } } @@ -292,12 +272,12 @@ void Plus_U_Base::uramping_update() bool Plus_U_Base::u_converged() { // Yukawa calculates U directly every iteration, always considered converged - if (Yukawa) { + if (use_yukawa) { return true; } - for (int i = 0; i < static_cast(this->U0.size()); i++) + for (int i = 0; i < static_cast(this->u_target.size()); i++) { - if (this->U[i] != this->U0[i]) + if (this->u_current[i] != this->u_target[i]) { return false; } @@ -322,7 +302,7 @@ void Plus_U_Base::copy_occ_mat(const UnitCell& ucell) { const int iat = ucell.itia2iat(T, I); - if (Plus_U_Base::nspin == 4) + 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) @@ -334,7 +314,7 @@ void Plus_U_Base::copy_occ_mat(const UnitCell& ucell) } } } - else if (Plus_U_Base::nspin == 1 || Plus_U_Base::nspin == 2) + 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]; @@ -377,11 +357,11 @@ void Plus_U_Base::zero_occ_mat(const UnitCell& ucell) for (int n = 0; n < N; n++) { - if (Plus_U_Base::nspin == 4) + if (this->nspin == 4) { occ_mat[iat][l][n][0].zero_out(); } - else if (Plus_U_Base::nspin == 1 || Plus_U_Base::nspin == 2) + else if (this->nspin == 1 || this->nspin == 2) { occ_mat[iat][l][n][0].zero_out(); occ_mat[iat][l][n][1].zero_out(); @@ -412,7 +392,7 @@ void Plus_U_Base::mix_occ_mat(const UnitCell& ucell, { const int iat = ucell.itia2iat(T, I); - if (Plus_U_Base::nspin == 4) + 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++) @@ -427,7 +407,7 @@ void Plus_U_Base::mix_occ_mat(const UnitCell& ucell, } } } - else if (Plus_U_Base::nspin == 1 || Plus_U_Base::nspin == 2) + 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; @@ -463,18 +443,18 @@ void Plus_U_Base::set_occ_mat(const UnitCell& ucell) for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); - if (Plus_U_Base::nspin == 4) + 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[eff_pot_pw_index[iat] + mm]; } - else if (Plus_U_Base::nspin == 1 || Plus_U_Base::nspin == 2) + 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[eff_pot_pw_index[iat] + mm]; - if (Plus_U_Base::nspin == 2) + if (this->nspin == 2) { occ_mat[iat][l][0][1].c[mm] = this->uom_array[half_size + eff_pot_pw_index[iat] + mm]; } @@ -538,7 +518,7 @@ void Plus_U_Base::read_occup_m(const UnitCell& ucell, if (!ifdftu) { - if (omc > 0) + 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"); } diff --git a/source/source_pw/module_pwdft/dftu_base.h b/source/source_pw/module_pwdft/dftu_base.h index 278ec0ee2fb..7885f6f8fa5 100644 --- a/source/source_pw/module_pwdft/dftu_base.h +++ b/source/source_pw/module_pwdft/dftu_base.h @@ -22,62 +22,47 @@ 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, - const int nlocal, - const bool gamma_only_local, - const std::string& ks_solver, - const bool cal_force, - const bool cal_stress, const std::string& device, - const int kpar); + const int kpar, + const std::vector& hubbard_u, + const double uramping, + const int occ_mat_ctrl, + const int mixing_dftu); void uramping_update(); bool u_converged(); - // mohan change these parameters to static, 2025-11-07 - static std::vector U; - static std::vector U0; - static std::vector orbital_corr; - static double uramping; - static int omc; - static int mixing_dftu; - static int nspin; - - // --- Accessors for static data --- - - static double get_hubbard_u(int it) { return U[it]; } - static double get_hubbard_u0(int it) { return U0[it]; } - static int get_num_u_types() { return static_cast(U.size()); } - static int get_orbital_corr(int it) { return orbital_corr[it]; } - static bool has_correlated_orbital(int it) { return orbital_corr[it] != -1; } - static const int* get_orbital_corr_data() { return orbital_corr.data(); } - - // mohan add 2025-11-08 for dftu_io::output free function + std::vector u_current; + std::vector u_target; + std::vector orbital_corr; + double uramping = 0.0; + int occ_mat_ctrl = 0; + int mixing_dftu = 0; + int nspin = 0; + + // --- Accessors --- + + 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(); } + 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]; } - static double get_energy() { return energy_u; } - static void set_energy(const double &e) { energy_u = e; } - static void set_double_energy() { energy_u *= 2.0; } + double get_energy() const { return energy_u; } + void set_energy(const double &e) { energy_u = e; } + void set_double_energy() { energy_u *= 2.0; } protected: - static double energy_u; + double energy_u = 0.0; - protected: int cal_type = 3; - std::string global_readin_dir; - std::string global_out_dir; - double yukawa_lambda = 0.0; - std::string init_chg; - int npol = 1; - int nlocal = 0; - bool gamma_only_local = false; - std::string ks_solver; - bool cal_force = false; - bool cal_stress = false; std::string device; int kpar = 1; @@ -141,8 +126,8 @@ class Plus_U_Base void mark_occ_mat_initialized() { occ_mat_initialized = true; } void mark_occ_mat_dirty() { occ_mat_initialized = false; } - static bool is_mixing_enabled() { return mixing_dftu != 0; } - static void enable_mixing() { mixing_dftu = 1; } + bool is_mixing_enabled() const { return mixing_dftu != 0; } + void enable_mixing() { mixing_dftu = 1; } protected: void copy_occ_mat(const UnitCell& ucell); @@ -212,7 +197,7 @@ class Plus_U_Base //============================================================= public: - static bool Yukawa; + bool use_yukawa = false; }; diff --git a/source/source_pw/module_pwdft/dftu_output.cpp b/source/source_pw/module_pwdft/dftu_output.cpp index 13502a9c473..0a01b3ec5ad 100644 --- a/source/source_pw/module_pwdft/dftu_output.cpp +++ b/source/source_pw/module_pwdft/dftu_output.cpp @@ -106,17 +106,17 @@ void output(const Plus_U_Base& dftu, { const int N = ucell.atoms[T].l_nchi[L]; - if (L >= Plus_U_Base::get_orbital_corr(T) && Plus_U_Base::has_correlated_orbital(T)) + if (L >= dftu.get_orbital_corr(T) && dftu.has_correlated_orbital(T)) { - if (L != Plus_U_Base::get_orbital_corr(T)) + if (L != dftu.get_orbital_corr(T)) { continue; } - if (!Plus_U_Base::Yukawa) + if (!dftu.use_yukawa) { GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << 0 - << " U=" << dftu.get_hubbard_u(T) * ModuleBase::Ry_to_eV << " eV" << std::endl; + << " U=" << dftu.get_u_current(T) * ModuleBase::Ry_to_eV << " eV" << std::endl; } else { @@ -177,12 +177,12 @@ void write_occup_m(const Plus_U_Base& dftu, for (int T = 0; T < ucell.ntype; T++) { - if (!Plus_U_Base::has_correlated_orbital(T)) + if (!dftu.has_correlated_orbital(T)) { continue; } const int NL = ucell.atoms[T].nwl + 1; - const int LC = Plus_U_Base::get_orbital_corr(T); + const int LC = dftu.get_orbital_corr(T); for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -190,7 +190,7 @@ void write_occup_m(const Plus_U_Base& dftu, for (int l = 0; l < NL; l++) { - if (l != Plus_U_Base::get_orbital_corr(T)) + if (l != dftu.get_orbital_corr(T)) { continue; } diff --git a/source/source_pw/module_pwdft/dftu_pw.cpp b/source/source_pw/module_pwdft/dftu_pw.cpp index 08a8aafd24c..2a03d47169c 100644 --- a/source/source_pw/module_pwdft/dftu_pw.cpp +++ b/source/source_pw/module_pwdft/dftu_pw.cpp @@ -37,7 +37,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, const int npol = psi_p->get_npol(); for(int ik = 0; ik < psi_p->get_nk(); ik++) { - int is = (Plus_U_Base::nspin == 2) ? isk[ik] : 0; + int is = (this->nspin == 2) ? isk[ik] : 0; psi_p->fix_k(ik); onsite_p->tabulate_atomic(ik); @@ -59,7 +59,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, const int m_begin = target_l * target_l; const int tlp1 = 2 * target_l + 1; const int tlp1_2 = tlp1 * tlp1; - if(Plus_U_Base::nspin == 4) + if(this->nspin == 4) { for(int ib = 0;ibget_npol(); for(int ik = 0; ik < psi_p->get_nk(); ik++) { - int is = (Plus_U_Base::nspin == 2) ? isk[ik] : 0; + int is = (this->nspin == 2) ? isk[ik] : 0; psi_p->fix_k(ik); onsite_p->tabulate_atomic(ik); @@ -138,7 +138,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, const int m_begin = target_l * target_l; const int tlp1 = 2 * target_l + 1; const int tlp1_2 = tlp1 * tlp1; - if(Plus_U_Base::nspin == 4) + if(this->nspin == 4) { for(int ib = 0;ibnspin != 4) { Parallel_Reduce::reduce_double_allpool(this->kpar, GlobalV::NPROC_IN_POOL, this->occ_mat[iat][target_l][0][0].c, size); - if(Plus_U_Base::nspin == 2) + if(this->nspin == 2) { Parallel_Reduce::reduce_double_allpool(this->kpar, GlobalV::NPROC_IN_POOL, @@ -228,7 +228,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, { this->uom_array[eff_pot_pw_index[iat]+mm] = this->occ_mat[iat][target_l][0][0].c[mm]; } - if(Plus_U_Base::nspin == 2) + if(this->nspin == 2) { const int half_size = this->uom_array.size() / 2; for(int mm=0;mmset_occ_mat(cell); } - Plus_U_Base::energy_u = 0.0; - const double weight_eu = (Plus_U_Base::nspin == 1) ? 1.0 : (Plus_U_Base::nspin == 2) ? 0.5 : 0.25; - const double diag_coeff = (Plus_U_Base::nspin == 4) ? 1.0 : 0.5; + 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 VU and energy (occ_mat already reduced above) for(int iat = 0; iat < cell.nat; iat++) { @@ -261,11 +261,11 @@ void Plus_U_Base::cal_occ_pw(const int iter, const int size = (2 * target_l + 1) * (2 * target_l + 1); //update effective potential - const double u_value = this->U[it]; + const double u_value = this->u_current[it]; std::complex* vu_iat = &(this->eff_pot_pw[this->eff_pot_pw_index[iat]]); const int m_size = 2 * target_l + 1; - if(Plus_U_Base::nspin == 4) + if(this->nspin == 4) { for (int m1 = 0; m1 < m_size; m1++) { @@ -273,7 +273,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, { vu_iat[m1 * m_size + m2] = u_value * (diag_coeff * (m1 == m2) - this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1]); - Plus_U_Base::energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1] + this->energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1] * this->occ_mat[iat][target_l][0][0].c[m1 * m_size + m2]; } } @@ -286,7 +286,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, { vu_iat[start + m1 * m_size + m2] = u_value * (0 - this->occ_mat[iat][target_l][0][0].c[start + m2 * m_size + m1]); - Plus_U_Base::energy_u += u_value * weight_eu + this->energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][0].c[start + m2 * m_size + m1] * this->occ_mat[iat][target_l][0][0].c[start + m1 * m_size + m2]; } @@ -323,12 +323,12 @@ void Plus_U_Base::cal_occ_pw(const int iter, { vu_iat[m1 * m_size + m2] = u_value * (diag_coeff * (m1 == m2) - this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1]); - Plus_U_Base::energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1] + this->energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][0].c[m2 * m_size + m1] * this->occ_mat[iat][target_l][0][0].c[m1 * m_size + m2]; } } // spin-down channel for nspin=2 - if(Plus_U_Base::nspin == 2) + if(this->nspin == 2) { std::complex* vu_iat1 = &(this->eff_pot_pw[this->eff_pot_pw.size()/2 + this->eff_pot_pw_index[iat]]); for (int m1 = 0; m1 < m_size; m1++) @@ -337,7 +337,7 @@ void Plus_U_Base::cal_occ_pw(const int iter, { vu_iat1[m1 * m_size + m2] = u_value * (diag_coeff * (m1 == m2) - this->occ_mat[iat][target_l][0][1].c[m2 * m_size + m1]); - Plus_U_Base::energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][1].c[m2 * m_size + m1] + this->energy_u += u_value * weight_eu * this->occ_mat[iat][target_l][0][1].c[m2 * m_size + m1] * this->occ_mat[iat][target_l][0][1].c[m1 * m_size + m2]; } } diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index ad1fc3c0425..f02afc017b6 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -34,7 +34,7 @@ void Forces::cal_force(UnitCell& ucell, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, surchem& solvent, - const Plus_U *p_dftu, //mohan add 2025-11-06 + const Plus_U_Base* p_dftu, const pseudopot_cell_vl* locpp, const pseudopot_cell_vnl* p_nlpp, K_Vectors* pkv, diff --git a/source/source_pw/module_pwdft/force_pw.h b/source/source_pw/module_pwdft/force_pw.h index 3338fbee191..3a2f3c5fa9f 100644 --- a/source/source_pw/module_pwdft/force_pw.h +++ b/source/source_pw/module_pwdft/force_pw.h @@ -12,7 +12,7 @@ #include "source_base/kernels/math_kernel_op.h" #include "source_psi/psi.h" #include "stru_fac.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 +#include "source_pw/module_pwdft/dftu_base.h" class pseudopot_cell_vnl; @@ -53,7 +53,7 @@ class Forces ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, surchem& solvent, - const Plus_U *p_dftu, //mohan add 2025-11-06 + const Plus_U_Base* p_dftu, const pseudopot_cell_vl* locpp, const pseudopot_cell_vnl* nlpp = nullptr, K_Vectors* pkv = nullptr, @@ -107,7 +107,7 @@ class Forces const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, const UnitCell& ucell_in, - const Plus_U &dftu, // mohan add 2025-11-06 + const Plus_U_Base& dftu, const psi::Psi , Device>* psi_in = nullptr); void cal_force_scc(ModuleBase::matrix& forcescc, diff --git a/source/source_pw/module_pwdft/force_pw_onsite.cpp b/source/source_pw/module_pwdft/force_pw_onsite.cpp index 9b501e790b9..56c5943a50e 100644 --- a/source/source_pw/module_pwdft/force_pw_onsite.cpp +++ b/source/source_pw/module_pwdft/force_pw_onsite.cpp @@ -5,7 +5,7 @@ #include "source_pw/module_pwdft/onsite_proj.h" #include "source_pw/module_pwdft/kernels/force_op.h" #include "source_io/module_parameter/parameter.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "source_lcao/module_deltaspin/spin_constrain.h" template @@ -13,7 +13,7 @@ void Forces::cal_force_onsite(ModuleBase::matrix& force_onsite, const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, const UnitCell& ucell_in, - const Plus_U &dftu, + const Plus_U_Base& dftu, const psi::Psi , Device>* psi_in) { ModuleBase::TITLE("Forces", "cal_force_onsite"); diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index e0b309f29ea..9809e43ee46 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -20,7 +20,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl* nlpp, - Plus_U* p_dftu, // mohan add 2025-11-06 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info) : ucell(ucell) diff --git a/source/source_pw/module_pwdft/hamilt_pw.h b/source/source_pw/module_pwdft/hamilt_pw.h index ab5d6fa9611..732fc4833d9 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.h +++ b/source/source_pw/module_pwdft/hamilt_pw.h @@ -6,7 +6,7 @@ #include "source_cell/klist.h" #include "source_estate/module_pot/potential_new.h" #include "source_hamilt/hamilt.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 +#include "source_pw/module_pwdft/dftu_base.h" #include "source_pw/module_pwdft/exx_helper.h" #include "source_pw/module_pwdft/vnl_pw.h" @@ -31,7 +31,7 @@ class HamiltPW : public Hamilt ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, - Plus_U* p_dftu, // mohan add 2025-11-06 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info); diff --git a/source/source_pw/module_pwdft/onsite_proj.cpp b/source/source_pw/module_pwdft/onsite_proj.cpp index 67458d7e7bc..331d3331bd4 100644 --- a/source/source_pw/module_pwdft/onsite_proj.cpp +++ b/source/source_pw/module_pwdft/onsite_proj.cpp @@ -6,7 +6,7 @@ #include #include "source_pw/module_pwdft/onsite_proj.h" #include "source_pw/module_pwdft/onsite_proj_print.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_cell/cell_tools.h" #include "source_io/module_parameter/parameter.h" @@ -590,7 +590,7 @@ void projectors::OnsiteProjector::cal_occupations( template void projectors::OnsiteProjector::cal_force_onsite_dftu(int ik, int npm, T* force, - const Plus_U& dftu, int nks, + const Plus_U_Base& dftu, int nks, const double* wg_ik) const { const int isk_val = this->isk_ ? this->isk_[ik] : 0; @@ -602,7 +602,7 @@ void projectors::OnsiteProjector::cal_force_onsite_dftu(int ik, int n template double projectors::OnsiteProjector::cal_stress_onsite_dftu(int ik, int npm, - const Plus_U& dftu, int nks, + const Plus_U_Base& dftu, int nks, const double* wg_ik) const { const int isk_val = this->isk_ ? this->isk_[ik] : 0; diff --git a/source/source_pw/module_pwdft/onsite_proj.h b/source/source_pw/module_pwdft/onsite_proj.h index 8924ae79db8..a4181e4b4bb 100644 --- a/source/source_pw/module_pwdft/onsite_proj.h +++ b/source/source_pw/module_pwdft/onsite_proj.h @@ -7,7 +7,7 @@ #include "source_pw/module_pwdft/radial_proj.h" #include "source_psi/psi.h" #include "source_pw/module_pwdft/onsite_proj_tools.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/dftu_base.h" #include #include @@ -90,12 +90,12 @@ namespace projectors /// high-level: compute DFT+U force contribution for one k-point void cal_force_onsite_dftu(int ik, int npm, T* force, - const Plus_U& dftu, int nks, + const Plus_U_Base& dftu, int nks, const double* wg_ik) const; /// high-level: compute DFT+U stress contribution for one k-point double cal_stress_onsite_dftu(int ik, int npm, - const Plus_U& dftu, int nks, + const Plus_U_Base& dftu, int nks, const double* wg_ik) const; /// high-level: compute DeltaSpin force contribution for one k-point diff --git a/source/source_pw/module_pwdft/op_pw_proj.cpp b/source/source_pw/module_pwdft/op_pw_proj.cpp index 5294a2b7de6..7244b315f78 100644 --- a/source/source_pw/module_pwdft/op_pw_proj.cpp +++ b/source/source_pw/module_pwdft/op_pw_proj.cpp @@ -4,7 +4,7 @@ #include "source_base/parallel_reduce.h" #include "source_base/tool_quit.h" #include "source_lcao/module_deltaspin/spin_constrain.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "source_pw/module_pwdft/onsite_proj.h" #include "source_pw/module_pwdft/kernels/onsite_op.h" @@ -14,7 +14,7 @@ namespace hamilt { template OnsiteProj>::OnsiteProj(const int* isk_in, const UnitCell* ucell_in, - Plus_U *p_dftu, // mohan add 2025-11-06 + Plus_U_Base* p_dftu, const bool cal_delta_spin, const bool cal_dftu) { diff --git a/source/source_pw/module_pwdft/op_pw_proj.h b/source/source_pw/module_pwdft/op_pw_proj.h index bd8044724da..298cd8749c9 100644 --- a/source/source_pw/module_pwdft/op_pw_proj.h +++ b/source/source_pw/module_pwdft/op_pw_proj.h @@ -5,7 +5,7 @@ #include "source_cell/unitcell.h" #include "source_base/kernels/math_kernel_op.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 20251106 +#include "source_pw/module_pwdft/dftu_base.h" namespace hamilt { @@ -27,7 +27,7 @@ class OnsiteProj> : public OperatorPW public: OnsiteProj(const int* isk_in, const UnitCell* ucell_in, - Plus_U *p_dftu, // mohan add 2025-11-06 + Plus_U_Base* p_dftu, const bool cal_delta_spin, const bool cal_dftu); @@ -65,7 +65,7 @@ class OnsiteProj> : public OperatorPW const UnitCell* ucell = nullptr; - Plus_U *dftu = nullptr; // mohan add 2025-11-06 + Plus_U_Base* dftu = nullptr; mutable int* ip_iat = nullptr; mutable T* lambda_coeff = nullptr; diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index db5ad901b68..fbab39d9db3 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -25,7 +25,7 @@ void iter_init_dftu_pw(const int iter, return; } - if (dftu.omc != 2) + if (dftu.occ_mat_ctrl != 2) { dftu.cal_occ_pw(iter, psi, wg, ucell, p_chgmix, isk); } diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 2a44b0b8f0f..738664d700f 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -32,7 +32,7 @@ void pw::setup_pot(const int istep, //! 0) DFT-1/2 calculations, sep potential need to generate // before effective potential calculation //---------------------------------------------------------- - if (PARAM.inp.dfthalf_type > 0) + if (inp.dfthalf_type > 0) { vsep_cell->generate_vsep_r(pw_rhod[0], sf.strucFac, ucell.sep_cell); } @@ -70,16 +70,16 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- //! 4) Onsite projectors //---------------------------------------------------------- - if (PARAM.inp.onsite_radius > 0) + if (inp.onsite_radius > 0) { auto* onsite_p = projectors::OnsiteProjector::get_instance(); - onsite_p->init(PARAM.inp.orbital_dir, + onsite_p->init(inp.orbital_dir, &ucell, *(kspw_psi), kv, *(pw_wfc), sf, - PARAM.inp.onsite_radius, + inp.onsite_radius, PARAM.globalv.nqx, PARAM.globalv.dq, pelec->wg, @@ -89,20 +89,20 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- //! 5) Spin-constrained algorithms //---------------------------------------------------------- - if (PARAM.inp.sc_mag_switch) + if (inp.sc_mag_switch) { spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); - sc.init_sc(PARAM.inp.sc_thr, - PARAM.inp.nsc, - PARAM.inp.nsc_min, - PARAM.inp.alpha_trial, - PARAM.inp.sccut, - PARAM.inp.sc_drop_thr, + sc.init_sc(inp.sc_thr, + inp.nsc, + inp.nsc_min, + inp.alpha_trial, + inp.sccut, + inp.sc_drop_thr, ucell, - PARAM.inp.sc_direction_only, + inp.sc_direction_only, nullptr, // parallel orbitals - PARAM.inp.nspin, + inp.nspin, kv, p_hamilt, kspw_psi, @@ -119,22 +119,20 @@ void pw::setup_pot(const int istep, // called in before_all_runners (esolver), which should // be improved later. Mohan note 2025-11-06 //---------------------------------------------------------- - if (PARAM.inp.dft_plus_u) + if (inp.dft_plus_u) { - const int nlocal_dftu = 0; dftu.init_base(ucell, PARAM.globalv.npol, - PARAM.inp.nspin, PARAM.inp.orbital_corr, PARAM.inp.yukawa_potential, PARAM.inp.yukawa_lambda, + inp.nspin, inp.orbital_corr, inp.yukawa_potential, PARAM.globalv.global_readin_dir, PARAM.globalv.global_out_dir, - PARAM.inp.init_chg, - nlocal_dftu, - PARAM.globalv.gamma_only_local, - PARAM.inp.ks_solver, - PARAM.inp.cal_force, - PARAM.inp.cal_stress, - PARAM.inp.device, - PARAM.inp.kpar); + inp.init_chg, + inp.device, + inp.kpar, + PARAM.globalv.hubbard_u, + PARAM.globalv.uramping, + inp.occ_mat_ctrl, + inp.mixing_dftu); } return; diff --git a/source/source_pw/module_pwdft/stress_func.h b/source/source_pw/module_pwdft/stress_func.h index 23218eca54a..9ac596cb63c 100644 --- a/source/source_pw/module_pwdft/stress_func.h +++ b/source/source_pw/module_pwdft/stress_func.h @@ -14,7 +14,7 @@ #include "source_pw/module_pwdft/stru_fac.h" #include "source_base/kernels/math_kernel_op.h" #include "source_psi/psi.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 +#include "source_pw/module_pwdft/dftu_base.h" //------------------------------------------------------------------- // mohan reconstruction note: 2021-02-07 @@ -179,7 +179,7 @@ class Stress_Func const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, const UnitCell& ucell_in, - const Plus_U &dftu, // mohan add 2025-11-06 + const Plus_U_Base& dftu, const void* psi_in, ModuleSymmetry::Symmetry* p_symm); // nonlocal part in PW basis diff --git a/source/source_pw/module_pwdft/stress_onsite.cpp b/source/source_pw/module_pwdft/stress_onsite.cpp index 8223be49ef4..f0ac36f8e3c 100644 --- a/source/source_pw/module_pwdft/stress_onsite.cpp +++ b/source/source_pw/module_pwdft/stress_onsite.cpp @@ -3,7 +3,7 @@ #include "source_base/timer.h" #include "source_pw/module_pwdft/onsite_proj.h" #include "source_io/module_parameter/parameter.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "stress_func.h" @@ -27,7 +27,7 @@ void Stress_Func::stress_onsite( const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, const UnitCell& ucell_in, - const Plus_U &dftu, // mohan add 2025-11-06 + const Plus_U_Base& dftu, const void* psi_in, ModuleSymmetry::Symmetry* p_symm ) diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index fa81550434a..bcf6aff8db2 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -12,7 +12,7 @@ template void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, UnitCell& ucell, const vdw::VdwResult* vdw_result, - Plus_U &dftu, // mhan add 2025-11-07 + Plus_U_Base& dftu, const pseudopot_cell_vl& locpp, const pseudopot_cell_vnl& nlpp, ModulePW::PW_Basis* rho_basis, diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index 030a3414233..f203a642e71 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -4,7 +4,6 @@ #include "source_estate/elecstate.h" #include "source_pw/module_pwdft/vl_pw.h" #include "stress_func.h" -#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 #include "source_hamilt/module_xc/coulomb_config.h" namespace vdw @@ -25,7 +24,7 @@ class Stress_PW : public Stress_Func void cal_stress(ModuleBase::matrix& smearing_sigmatot, UnitCell& ucell, const vdw::VdwResult* vdw_result, - Plus_U &dftu, // mhan add 2025-11-07 + Plus_U_Base& dftu, const pseudopot_cell_vl& locpp, const pseudopot_cell_vnl& nlpp, ModulePW::PW_Basis* rho_basis, diff --git a/source/source_pw/module_stodft/test/test_hamilt_sto.cpp b/source/source_pw/module_stodft/test/test_hamilt_sto.cpp index 126ef923fb3..b749cc67b47 100644 --- a/source/source_pw/module_stodft/test/test_hamilt_sto.cpp +++ b/source/source_pw/module_stodft/test/test_hamilt_sto.cpp @@ -1,4 +1,5 @@ #include "../hamilt_sdft_pw.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "source_hamilt/operator.h" #include "gtest/gtest.h" @@ -14,7 +15,7 @@ hamilt::HamiltPW::HamiltPW( ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* ppcell, - Plus_U* p_dftu, // mohan add 20251108 + Plus_U_Base* p_dftu, const UnitCell* ucell, const General_Exx_Info* exx_info){} diff --git a/source/source_pw/module_stodft/test/test_sto_tool.cpp b/source/source_pw/module_stodft/test/test_sto_tool.cpp index 81081f869a5..e4f1c42868e 100644 --- a/source/source_pw/module_stodft/test/test_sto_tool.cpp +++ b/source/source_pw/module_stodft/test/test_sto_tool.cpp @@ -1,4 +1,5 @@ #include "../sto_tool.h" +#include "source_pw/module_pwdft/dftu_base.h" #include "mpi.h" #include @@ -12,7 +13,7 @@ hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl*, - Plus_U* p_dftu, // mohan add 20251108 + Plus_U_Base* p_dftu, const UnitCell*, const General_Exx_Info*){} diff --git a/tests/02_NAO_Gamma/CASES_CPU.txt b/tests/02_NAO_Gamma/CASES_CPU.txt index b7313ff1fe6..c002f71e58d 100644 --- a/tests/02_NAO_Gamma/CASES_CPU.txt +++ b/tests/02_NAO_Gamma/CASES_CPU.txt @@ -34,5 +34,6 @@ scf_out_wf #scf_out_wf_spin2 scf_solvation scf_u_spin2 +scf_u_spin2_old scf_upf100 test.sum diff --git a/tests/02_NAO_Gamma/scf_u_spin2_old/INPUT b/tests/02_NAO_Gamma/scf_u_spin2_old/INPUT new file mode 100644 index 00000000000..646bf5c2919 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_u_spin2_old/INPUT @@ -0,0 +1,37 @@ +INPUT_PARAMETERS +suffix autotest +nbands 28 + +calculation scf +ecutwfc 6 +scf_thr 1.0e-7 +scf_nmax 30 +out_chg 0 + +smearing_method gaussian +smearing_sigma 0.07 + +cal_force 0 +cal_stress 0 + +#force_thr_ev 0.01 +#relax_method cg +#relax_bfgs_init 0.5 + +mixing_type broyden +mixing_beta 0.2 +mixing_gg0 1.0 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 1 +symmetry 0 +nspin 2 + +#Parameter DFT+U +dft_plus_u 2 # 1: new 2: old +orbital_corr 2 2 +hubbard_u 5.0 5.0 +onsite_radius 5.0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB diff --git a/tests/02_NAO_Gamma/scf_u_spin2_old/README b/tests/02_NAO_Gamma/scf_u_spin2_old/README new file mode 100644 index 00000000000..665790f4420 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_u_spin2_old/README @@ -0,0 +1 @@ +DFTU + NSPIN2 but use old DFTU method, Fe2O2, gamma-only case diff --git a/tests/02_NAO_Gamma/scf_u_spin2_old/STRU b/tests/02_NAO_Gamma/scf_u_spin2_old/STRU new file mode 100644 index 00000000000..b0f79de153c --- /dev/null +++ b/tests/02_NAO_Gamma/scf_u_spin2_old/STRU @@ -0,0 +1,29 @@ +ATOMIC_SPECIES +Fe 1.000 Fe_ONCV_PBE-1.0.upf +O 1.000 O_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_2s1p1d.orb +O_gga_7au_60Ry_2s2p1d.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 +0.50 0.50 0.50 mag -1.0 + +O +0.0 +2 +0.25 0.25 0.25 1 1 1 +0.75 0.75 0.75 1 1 1 diff --git a/tests/02_NAO_Gamma/scf_u_spin2_old/result.ref b/tests/02_NAO_Gamma/scf_u_spin2_old/result.ref new file mode 100644 index 00000000000..afd5b985d03 --- /dev/null +++ b/tests/02_NAO_Gamma/scf_u_spin2_old/result.ref @@ -0,0 +1,3 @@ +etotref -6304.2287463338798261 +etotperatomref -1576.0571865835 +totaltimeref 1.17