diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 082271ee2e..1e80391dfc 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -146,7 +146,8 @@ void ESolver_DoubleXC::before_scf(UnitCell& ucell, const int istep) this->kv, this->two_center_bundle_, this->orb_, - this->dmat_base.dm, + this->dmat_base.dm, + &this->dftu, this->deepks, istep, this->exx_nao); @@ -396,6 +397,7 @@ void ESolver_DoubleXC::cal_force(UnitCell& ucell, ModuleBase::matrix& fo this->kv, this->pw_rho, this->solvent, + this->dftu, this->deepks, this->exx_nao, &ucell.symm); diff --git a/source/source_esolver/esolver_gets.h b/source/source_esolver/esolver_gets.h index a0d163ef4d..564fd55035 100644 --- a/source/source_esolver/esolver_gets.h +++ b/source/source_esolver/esolver_gets.h @@ -37,7 +37,7 @@ class ESolver_GetS : public ESolver_KS> TwoCenterBundle two_center_bundle_; - // temporary introduced during removing GlobalC::ORB + // temporary introduced LCAO_Orbitals orb_; }; } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index 2818152b22..18b08d0c64 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -15,6 +15,7 @@ #include "source_estate/elecstate_print.h" // print_etot #include "source_io/print_info.h" // print_parameters #include "source_psi/setup_psi.h" // mohan add 20251009 +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 namespace ModuleESolver { @@ -281,9 +282,19 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i this->pelec->nelec_spin.data()); // 2.2) charge mixing + // SCF will continue if U is not converged for uramping calculation + bool converged_u = true; + // to avoid unnecessary dependence on dft+u, refactor is needed +#ifdef __LCAO + if (PARAM.inp.dft_plus_u) + { + converged_u = this->dftu.u_converged(); + } +#endif + module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr, this->p_chgmix, this->pw_rhod->nrxx, this->drho, this->oscillate_esolver, conv_esolver, hsolver_error, - this->scf_thr, this->scf_ene_thr, PARAM.inp); + this->scf_thr, this->scf_ene_thr, converged_u, PARAM.inp); // 2.3) Update potentials (should be done every SF iter) elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index bd47446a68..787b58ba74 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -7,6 +7,7 @@ #include "source_estate/module_charge/charge_mixing.h" // use charge mixing #include "source_psi/psi.h" // use electronic wave functions #include "source_hamilt/hamilt.h" // use Hamiltonian +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace ModuleESolver { @@ -61,6 +62,9 @@ class ESolver_KS : public ESolver_FP //! Electronic wavefunctions psi::Psi* psi = nullptr; + //! DFT+U method, mohan add 2025-11-07 + Plus_U dftu; + std::string basisname; //! esolver_ks_lcao.cpp double esolver_KS_ne = 0.0; //! number of electrons double diag_ethr; //! the threshold for diagonalization diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index b3f83c42b8..d92bf2ab3b 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -4,7 +4,6 @@ #include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp #include "source_estate/module_charge/symmetry_rho.h" #include "source_lcao/LCAO_domain.h" // need DeePKS_init -#include "source_lcao/module_dftu/dftu.h" #include "source_lcao/FORCE_STRESS.h" #include "source_estate/elecstate_lcao.h" #include "source_lcao/hamilt_lcao.h" @@ -99,8 +98,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // 9) initialize DFT+U if (inp.dft_plus_u) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_); + this->dftu.init(ucell, &this->pv, this->kv.get_nks(), &orb_); } // 10) init local pseudopotentials @@ -189,7 +187,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { this->p_hamilt = new hamilt::HamiltLCAO( ucell, this->gd, &this->pv, this->pelec->pot, this->kv, - two_center_bundle_, orb_, this->dmat.dm, this->deepks, istep, exx_nao); + two_center_bundle_, orb_, this->dmat.dm, &this->dftu, this->deepks, istep, exx_nao); } // 9) for each ionic step, the overlap must be rebuilt @@ -274,7 +272,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for this->gd, this->pv, this->pelec, this->dmat, this->psi, two_center_bundle_, orb_, force, this->scs, this->locpp, this->sf, this->kv, - this->pw_rho, this->solvent, this->deepks, + this->pw_rho, this->solvent, this->dftu, this->deepks, this->exx_nao, &ucell.symm); // delete RA after cal_force @@ -338,7 +336,8 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // call iter_init() of ESolver_KS ESolver_KS::iter_init(ucell, istep, iter); - module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); + module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dftu, + this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); // mohan update 2012-06-05 this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell); @@ -377,10 +376,10 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const { if (istep != 0 || iter != 1) { - GlobalC::dftu.set_dmr(this->dmat.dm); + this->dftu.set_dmr(this->dmat.dm); } // Calculate U and J if Yukawa potential is used - GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx); + this->dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx); } #ifdef __MLALGO @@ -489,18 +488,18 @@ 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 if (PARAM.inp.dft_plus_u) { - // only old DFT+U method should calculated energy correction in esolver, - // new DFT+U method will calculate energy in calculating Hamiltonian + // old DFT+U method calculates energy correction in esolver, + // new DFT+U method calculates energy in Hamiltonian if (PARAM.inp.dft_plus_u == 2) { - if (GlobalC::dftu.omc != 2) + if (this->dftu.omc != 2) { - ModuleDFTU::dftu_cal_occup_m(iter, ucell, dm_vec, this->kv, - this->p_chgmix->get_mixing_beta(), hamilt_lcao); + dftu_cal_occup_m(iter, ucell, dm_vec, this->kv, + this->p_chgmix->get_mixing_beta(), hamilt_lcao, this->dftu); } - GlobalC::dftu.cal_energy_correction(ucell, istep); + this->dftu.cal_energy_correction(ucell, istep); } - GlobalC::dftu.output(ucell); + this->dftu.output(ucell); } // 2) for deepks, calculate delta_e, output labels during electronic steps @@ -532,7 +531,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // use the converged occupation matrix for next MD/Relax SCF calculation if (PARAM.inp.dft_plus_u && conv_esolver) { - GlobalC::dftu.initialed_locale = true; + this->dftu.initialed_locale = true; } // control the output related to the finished iteration @@ -566,7 +565,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const //! 2) output of lcao every few ionic steps ModuleIO::ctrl_scf_lcao(ucell, PARAM.inp, this->kv, this->pelec, this->dmat.dm, this->pv, - this->gd, this->psi, hamilt_lcao, this->two_center_bundle_, + this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_, this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf, this->rdmft_solver, this->deepks, this->exx_nao, this->conv_esolver, this->scf_nmax_flag, istep); diff --git a/source/source_esolver/esolver_ks_lcao.h b/source/source_esolver/esolver_ks_lcao.h index 2bcf0c672c..4191306788 100644 --- a/source/source_esolver/esolver_ks_lcao.h +++ b/source/source_esolver/esolver_ks_lcao.h @@ -69,32 +69,34 @@ class ESolver_KS_LCAO : public ESolver_KS //! GintInfo: used to store some basic infomation about module_gint std::unique_ptr gint_info_; + //! NAO: store related information + LCAO_Orbitals orb_; + //! NAO orbitals: two-center integrations TwoCenterBundle two_center_bundle_; //! Add density matrix class, mohan add 2025-10-30 LCAO_domain::Setup_DM dmat; + + // For deepks method, mohan add 2025-10-08 + Setup_DeePKS deepks; + + // For exact-exchange energy, mohan add 2025-10-08 + Exx_NAO exx_nao; + //! For RDMFT calculations, added by jghan, 2024-03-16 rdmft::RDMFT rdmft_solver; - //! NAO: store related information - LCAO_Orbitals orb_; + //! For linear-response TDDFT + friend class LR::ESolver_LR; + friend class LR::ESolver_LR, double>; // Temporarily store the stress to unify the interface with PW, // because it's hard to seperate force and stress calculation in LCAO. ModuleBase::matrix scs; bool have_force = false; - // deepks method, mohan add 2025-10-08 - Setup_DeePKS deepks; - - // exact-exchange energy, mohan add 2025-10-08 - Exx_NAO exx_nao; - - friend class LR::ESolver_LR; - friend class LR::ESolver_LR, double>; - public: const Record_adj & get_RA() const { return RA; } diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index f1886eb004..6c19a6c5a1 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -55,7 +55,13 @@ ESolver_KS_PW::~ESolver_KS_PW() template void ESolver_KS_PW::allocate_hamilt(const UnitCell& ucell) { - this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell); + this->p_hamilt = new hamilt::HamiltPW( + this->pelec->pot, + this->pw_wfc, + &this->kv, + &this->ppcell, + &this->dftu, + &ucell); } template @@ -133,8 +139,10 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->allocate_hamilt(ucell); //! Setup potentials (local, non-local, sc, +U, DFT-1/2) + // note: init DFT+U is done here for pw basis for every scf iteration, however, + // init DFT+U is done in "before_all_runners" in LCAO basis. This should be refactored, mohan note 2025-11-06 pw::setup_pot(istep, ucell, this->kv, this->sf, this->pelec, this->Pgrid, - this->chr, this->locpp, this->ppcell, this->vsep_cell, + this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell, this->stp.psi_t, this->p_hamilt, this->pw_wfc, this->pw_rhod, PARAM.inp); // setup psi (electronic wave functions) @@ -162,7 +170,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const ESolver_KS::iter_init(ucell, istep, iter); // 2) perform charge mixing for KSDFT using pw basis - module_charge::chgmixing_ks_pw(iter, this->p_chgmix, PARAM.inp); + module_charge::chgmixing_ks_pw(iter, this->p_chgmix, this->dftu, PARAM.inp); // 3) mohan move harris functional here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) @@ -172,14 +180,13 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const // should before lambda loop in DeltaSpin if (PARAM.inp.dft_plus_u && (iter != 1 || istep != 0)) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); // only old DFT+U method should calculate energy correction in esolver, // new DFT+U method will calculate energy when evaluating the Hamiltonian - if (dftu->omc != 2) + if (this->dftu.omc != 2) { - dftu->cal_occ_pw(iter, this->stp.psi_t, this->pelec->wg, ucell, PARAM.inp.mixing_beta); + this->dftu.cal_occ_pw(iter, this->stp.psi_t, this->pelec->wg, ucell, PARAM.inp.mixing_beta); } - dftu->output(ucell); + this->dftu.output(ucell); } } @@ -389,7 +396,7 @@ void ESolver_KS_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& fo // Calculate forces ff.cal_force(ucell, force, *this->pelec, this->pw_rhod, &ucell.symm, - &this->sf, this->solvent, &this->locpp, &this->ppcell, + &this->sf, this->solvent, &this->dftu, &this->locpp, &this->ppcell, &this->kv, this->pw_wfc, this->stp.psi_d); } @@ -401,7 +408,7 @@ void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& s // mohan add 2025-10-12 this->stp.update_psi_d(); - ss.cal_stress(stress, ucell, this->locpp, this->ppcell, this->pw_rhod, + ss.cal_stress(stress, ucell, this->dftu, this->locpp, this->ppcell, this->pw_rhod, &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->stp.psi_d); // external stress diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index e3fad24322..f609b3e2c1 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -575,7 +575,10 @@ double ESolver_OF::cal_energy() void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { Forces ff(ucell.nat); - ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, &this->locpp); + + // here nullptr is for DFT+U, which may cause bugs, mohan note 2025-11-07 + // solvent can be used? mohan ask 2025-11-07 + ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, nullptr, &this->locpp); } /** diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 9c413836d4..0c404dbbcf 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -182,6 +182,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) two_center_bundle_, orb_, this->dmat.dm, + &this->dftu, this->deepks, istep, this->exx_nao); diff --git a/source/source_estate/elecstate_energy_terms.cpp b/source/source_estate/elecstate_energy_terms.cpp index 73abe5f17b..651fd494a2 100644 --- a/source/source_estate/elecstate_energy_terms.cpp +++ b/source/source_estate/elecstate_energy_terms.cpp @@ -4,7 +4,7 @@ #include "source_estate/module_pot/gatefield.h" #include "source_lcao/module_deepks/LCAO_deepks.h" #include "source_lcao/module_deltaspin/spin_constrain.h" -#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 namespace elecstate { @@ -36,7 +36,7 @@ double ElecState::get_solvent_model_Acav() double ElecState::get_dftu_energy() { - return GlobalC::dftu.get_energy(); + return Plus_U::get_energy(); } double ElecState::get_local_pp_energy() @@ -52,4 +52,4 @@ double ElecState::get_local_pp_energy() return local_pseudopot_energy; } -} // namespace elecstate \ No newline at end of file +} // namespace elecstate diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index 4a5918f75e..5bdc9e733b 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -192,28 +192,52 @@ void print_etot(const Magnetism& magnet, if( (iter % PARAM.inp.out_freq_elec == 0) || converged || iter == PARAM.inp.scf_nmax ) { int n_order = std::max(0, Occupy::gaussian_type); + + //! Kohn-Sham functional energy titles.push_back("E_KohnSham"); energies_Ry.push_back(elec.f_en.etot); + + //! Kohn-Sham energy with sigma->0 titles.push_back("E_KS(sigma->0)"); energies_Ry.push_back(elec.f_en.etot - elec.f_en.demet / (2 + n_order)); + + //! Harris functional energy titles.push_back("E_Harris"); energies_Ry.push_back(elec.f_en.etot_harris); + + //! band energy titles.push_back("E_band"); energies_Ry.push_back(elec.f_en.eband); + + //! one-electron energy titles.push_back("E_one_elec"); energies_Ry.push_back(elec.f_en.eband + elec.f_en.deband); + + //! Hartree energy titles.push_back("E_Hartree"); energies_Ry.push_back(elec.f_en.hartree_energy); + + //! exchange-correlation energy titles.push_back("E_xc"); energies_Ry.push_back(elec.f_en.etxc - elec.f_en.etxcc); + + //! Ewald energy titles.push_back("E_Ewald"); energies_Ry.push_back(elec.f_en.ewald_energy); + + //! entropy energy titles.push_back("E_entropy(-TS)"); energies_Ry.push_back(elec.f_en.demet); + + //! correction energy for scf titles.push_back("E_descf"); energies_Ry.push_back(elec.f_en.descf); - titles.push_back("E_LocalPP"); + + //! local potential energy + titles.push_back("E_localpp"); energies_Ry.push_back(elec.f_en.e_local_pp); + + //! vdw energy std::string vdw_method = PARAM.inp.vdw_method; if (vdw_method == "d2") // Peize Lin add 2014-04, update 2021-03-09 { @@ -225,8 +249,19 @@ void print_etot(const Magnetism& magnet, titles.push_back("E_vdwD3"); energies_Ry.push_back(elec.f_en.evdw); } + + // mohan add 20251108 + if (PARAM.inp.dft_plus_u) + { + titles.push_back("E_plusU"); + energies_Ry.push_back(elec.f_en.edftu); + } + + //! hybrid functional energy titles.push_back("E_exx"); energies_Ry.push_back(elec.f_en.exx); + + //! solvation energy if (PARAM.inp.imp_sol) { titles.push_back("E_sol_el"); @@ -234,17 +269,22 @@ void print_etot(const Magnetism& magnet, titles.push_back("E_sol_cav"); energies_Ry.push_back(elec.f_en.esol_cav); } + + //! electric field energy if (PARAM.inp.efield_flag) { titles.push_back("E_efield"); energies_Ry.push_back(elecstate::Efield::etotefield); } + + //! gate energy if (PARAM.inp.gate_flag) { titles.push_back("E_gatefield"); energies_Ry.push_back(elecstate::Gatefield::etotgatefield); } + //! deepks energy #ifdef __MLALGO if (PARAM.inp.deepks_scf) { diff --git a/source/source_estate/module_charge/chgmixing.cpp b/source/source_estate/module_charge/chgmixing.cpp index 59467fc010..45e5c5b350 100644 --- a/source/source_estate/module_charge/chgmixing.cpp +++ b/source/source_estate/module_charge/chgmixing.cpp @@ -15,6 +15,7 @@ void module_charge::chgmixing_ks(const int iter, // scf iteration number const double &hsolver_error, const double &scf_thr, const double &scf_ene_thr, + const bool converged_u, // mohan add 2025-11-06 const Input_para& inp) // input parameters { @@ -36,17 +37,8 @@ void module_charge::chgmixing_ks(const int iter, // scf iteration number // drho will be 0 at p_chgmix->mixing_restart step, which is // not ground state bool not_restart_step = !(iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0); - // SCF will continue if U is not converged for uramping calculation - bool is_U_converged = true; - // to avoid unnecessary dependence on dft+u, refactor is needed -#ifdef __LCAO - if (inp.dft_plus_u) - { - is_U_converged = GlobalC::dftu.u_converged(); - } -#endif - conv_esolver = (drho < scf_thr && not_restart_step && is_U_converged); + conv_esolver = (drho < scf_thr && not_restart_step && converged_u); // add energy threshold for SCF convergence if (scf_ene_thr > 0.0) @@ -67,8 +59,6 @@ void module_charge::chgmixing_ks(const int iter, // scf iteration number } } - - // If drho < hsolver_error in the first iter or drho < scf_thr, we // do not change rho. if (drho < hsolver_error || conv_esolver || inp.calculation == "nscf") @@ -127,9 +117,9 @@ 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 const Input_para& inp) // input parameters { ModuleBase::TITLE("module_charge", "chgmixing_ks_pw"); @@ -148,12 +138,11 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number if (inp.dft_plus_u) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - if (dftu->uramping > 0.01 && !dftu->u_converged()) + if (dftu.uramping > 0.01 && !dftu.u_converged()) { p_chgmix->mixing_restart_step = inp.scf_nmax + 1; } - if (dftu->uramping > 0.01) + if (dftu.uramping > 0.01) { bool do_uramping = true; if (inp.sc_mag_switch) @@ -167,11 +156,11 @@ void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number } if (do_uramping) { - dftu->uramping_update(); // update U by uramping if uramping > 0.01 + dftu.uramping_update(); // update U by uramping if uramping > 0.01 std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < dftu->U0.size(); i++) + for (int i = 0; i < dftu.U0.size(); i++) { - std::cout << dftu->U[i] * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.U[i] * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } @@ -184,6 +173,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 const int nnr, // dimension of density matrix const Input_para& inp) // input parameters { @@ -195,12 +185,12 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number p_chgmix->mixing_restart_step = inp.scf_nmax + 1; p_chgmix->mixing_restart_count = 0; // this output will be removed once the feeature is stable - if (GlobalC::dftu.uramping > 0.01) + if (dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + for (int i = 0; i < dftu.U0.size(); i++) { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.U[i] * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } @@ -213,17 +203,17 @@ void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number p_chgmix->mixing_restart_count++; if (inp.dft_plus_u) { - GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 - if (GlobalC::dftu.uramping > 0.01) + dftu.uramping_update(); // update U by uramping if uramping > 0.01 + if (dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + for (int i = 0; i < dftu.U0.size(); i++) { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + std::cout << dftu.U[i] * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } - if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) + if (dftu.uramping > 0.01 && !dftu.u_converged()) { p_chgmix->mixing_restart_step = inp.scf_nmax + 1; } diff --git a/source/source_estate/module_charge/chgmixing.h b/source/source_estate/module_charge/chgmixing.h index 3be8660013..4a04a0880d 100644 --- a/source/source_estate/module_charge/chgmixing.h +++ b/source/source_estate/module_charge/chgmixing.h @@ -6,6 +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 namespace module_charge { @@ -22,14 +23,17 @@ void chgmixing_ks(const int iter, // scf iteration number const double &hsolver_error, const double &scf_thr, const double &scf_ene_thr, + const bool converged_u, // mohan add 2025-11-06 const Input_para& inp); // input parameters void chgmixing_ks_pw(const int iter, Charge_Mixing* p_chgmix, + Plus_U &dftu, // mohan add 2025-11-06 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 const int nnr, // dimension of density matrix const Input_para& inp); // input parameters diff --git a/source/source_hsolver/test/diago_bpcg_test.cpp b/source/source_hsolver/test/diago_bpcg_test.cpp index 93e1147ccf..f36a09dafe 100644 --- a/source/source_hsolver/test/diago_bpcg_test.cpp +++ b/source/source_hsolver/test/diago_bpcg_test.cpp @@ -99,7 +99,7 @@ class DiagoBPCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/source_hsolver/test/diago_cg_float_test.cpp b/source/source_hsolver/test/diago_cg_float_test.cpp index b446bcb912..9af83f174c 100644 --- a/source/source_hsolver/test/diago_cg_float_test.cpp +++ b/source/source_hsolver/test/diago_cg_float_test.cpp @@ -108,7 +108,7 @@ class DiagoCGPrepare float *en = new float[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); psi::Psi> psi; psi.resize(ik,nband,npw); //psi.fix_k(0); diff --git a/source/source_hsolver/test/diago_cg_test.cpp b/source/source_hsolver/test/diago_cg_test.cpp index 65c283d94b..d360e47e53 100644 --- a/source/source_hsolver/test/diago_cg_test.cpp +++ b/source/source_hsolver/test/diago_cg_test.cpp @@ -104,7 +104,7 @@ class DiagoCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); psi::Psi> psi; psi.resize(ik,nband,npw); //psi.fix_k(0); diff --git a/source/source_hsolver/test/diago_david_float_test.cpp b/source/source_hsolver/test/diago_david_float_test.cpp index f907f939e4..76c5226ce3 100644 --- a/source/source_hsolver/test/diago_david_float_test.cpp +++ b/source/source_hsolver/test/diago_david_float_test.cpp @@ -82,7 +82,7 @@ class DiagoDavPrepare //do Diago_David::diag() float* en = new float[npw]; hamilt::Hamilt> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/source_hsolver/test/diago_david_real_test.cpp b/source/source_hsolver/test/diago_david_real_test.cpp index b8670bc74d..77a64dcd01 100644 --- a/source/source_hsolver/test/diago_david_real_test.cpp +++ b/source/source_hsolver/test/diago_david_real_test.cpp @@ -81,7 +81,7 @@ class DiagoDavPrepare //do Diago_David::diag() double* en = new double[npw]; hamilt::Hamilt* phm; - phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr,nullptr); + phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/source_hsolver/test/diago_david_test.cpp b/source/source_hsolver/test/diago_david_test.cpp index 643eeed4bf..6239be7a2e 100644 --- a/source/source_hsolver/test/diago_david_test.cpp +++ b/source/source_hsolver/test/diago_david_test.cpp @@ -77,14 +77,16 @@ class DiagoDavPrepare { //calculate eigenvalues by LAPACK; double* e_lapack = new double[npw]; - double* ev; - if(mypnum == 0) { lapackEigen(npw, DIAGOTEST::hmatrix, e_lapack,DETAILINFO); -} + double* ev = nullptr; + if(mypnum == 0) + { + lapackEigen(npw, DIAGOTEST::hmatrix, e_lapack,DETAILINFO); + } //do Diago_David::diag() double* en = new double[npw]; hamilt::Hamilt> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/source_hsolver/test/diago_mock.h b/source/source_hsolver/test/diago_mock.h index ea2c2affaa..16e1359cae 100644 --- a/source/source_hsolver/test/diago_mock.h +++ b/source/source_hsolver/test/diago_mock.h @@ -572,7 +572,13 @@ template<> void hamilt::HamiltPW::updateHk(const int ik) return; } -template<> hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) +template<> hamilt::HamiltPW::HamiltPW( + elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* pkv, + pseudopot_cell_vnl* ppcell, + Plus_U* p_dftu, // mohan add 20251108 + const UnitCell* ucell) { this->ops = new OperatorMock; } @@ -587,7 +593,13 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) +template<> hamilt::HamiltPW>::HamiltPW( + elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* pkv, + pseudopot_cell_vnl* ppcell, + Plus_U* p_dftu, // mohan add 20251108 + const UnitCell* ucell) { this->ops = new OperatorMock>; } @@ -602,7 +614,13 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) +template<> hamilt::HamiltPW>::HamiltPW( + elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* pkv, + pseudopot_cell_vnl* ppcell, + Plus_U* p_dftu, // mohan add 20251108 + const UnitCell* ucell) { this->ops = new OperatorMock>; } diff --git a/source/source_io/ctrl_scf_lcao.cpp b/source/source_io/ctrl_scf_lcao.cpp index f66fd2720b..590451b32d 100644 --- a/source/source_io/ctrl_scf_lcao.cpp +++ b/source/source_io/ctrl_scf_lcao.cpp @@ -42,6 +42,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, Grid_Driver& gd, psi::Psi* psi, hamilt::HamiltLCAO* p_hamilt, + Plus_U &dftu, // mohan add 2025-11-07 TwoCenterBundle &two_center_bundle, LCAO_Orbitals &orb, const ModulePW::PW_Basis_K* pw_wfc, // for berryphase @@ -222,7 +223,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, ucell, gd, kv, - p_ham_tk); + p_ham_tk, + &dftu); //------------------------------------------------------------------ //! 8) Output kinetic matrix @@ -474,6 +476,7 @@ template void ModuleIO::ctrl_scf_lcao( Grid_Driver& gd, psi::Psi* psi, hamilt::HamiltLCAO* p_hamilt, + Plus_U &dftu, // mohan add 2025-11-07 TwoCenterBundle &two_center_bundle, LCAO_Orbitals &orb, const ModulePW::PW_Basis_K* pw_wfc, // for berryphase @@ -498,6 +501,7 @@ template void ModuleIO::ctrl_scf_lcao, double>( Grid_Driver& gd, psi::Psi>* psi, hamilt::HamiltLCAO, double>* p_hamilt, + Plus_U &dftu, // mohan add 2025-11-07 TwoCenterBundle &two_center_bundle, LCAO_Orbitals &orb, const ModulePW::PW_Basis_K* pw_wfc, // for berryphase @@ -521,6 +525,7 @@ template void ModuleIO::ctrl_scf_lcao, std::complex Grid_Driver& gd, psi::Psi>* psi, hamilt::HamiltLCAO, std::complex>* p_hamilt, + Plus_U &dftu, // mohan add 2025-11-07 TwoCenterBundle &two_center_bundle, LCAO_Orbitals &orb, const ModulePW::PW_Basis_K* pw_wfc, // for berryphase diff --git a/source/source_io/ctrl_scf_lcao.h b/source/source_io/ctrl_scf_lcao.h index 9dd381cd2a..d6780b52ed 100644 --- a/source/source_io/ctrl_scf_lcao.h +++ b/source/source_io/ctrl_scf_lcao.h @@ -16,6 +16,7 @@ #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008 #include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04 +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace ModuleIO { @@ -30,6 +31,7 @@ namespace ModuleIO Grid_Driver& gd, psi::Psi* psi, hamilt::HamiltLCAO* p_hamilt, + Plus_U &dftu, // mohan add 2025-11-07 TwoCenterBundle &two_center_bundle, LCAO_Orbitals &orb, const ModulePW::PW_Basis_K* pw_wfc, // for berryphase diff --git a/source/source_io/input_conv.cpp b/source/source_io/input_conv.cpp index 279c09c405..f4edf1da6e 100644 --- a/source/source_io/input_conv.cpp +++ b/source/source_io/input_conv.cpp @@ -43,6 +43,8 @@ #include "source_hsolver/hsolver_pw.h" #include "source_md/md_func.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 + #ifdef __LCAO std::vector Input_Conv::convert_units(std::string params, double c) { std::vector params_ori; @@ -213,16 +215,16 @@ void Input_Conv::Convert() if (PARAM.inp.dft_plus_u) { - GlobalC::dftu.Yukawa = PARAM.inp.yukawa_potential; - GlobalC::dftu.omc = PARAM.inp.omc; - GlobalC::dftu.orbital_corr = PARAM.inp.orbital_corr; - GlobalC::dftu.uramping = PARAM.globalv.uramping; - GlobalC::dftu.mixing_dftu = PARAM.inp.mixing_dftu; - GlobalC::dftu.U = PARAM.globalv.hubbard_u; - GlobalC::dftu.U0 = PARAM.globalv.hubbard_u; + Plus_U::Yukawa = PARAM.inp.yukawa_potential; + 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(GlobalC::dftu.U.data(), PARAM.inp.ntype); + ModuleBase::GlobalFunc::ZEROS(Plus_U::U.data(), PARAM.inp.ntype); } } diff --git a/source/source_io/output_mat_sparse.cpp b/source/source_io/output_mat_sparse.cpp index 70992b1766..a98e6c8297 100644 --- a/source/source_io/output_mat_sparse.cpp +++ b/source/source_io/output_mat_sparse.cpp @@ -19,14 +19,15 @@ void output_mat_sparse(const bool& out_mat_hsR, UnitCell& ucell, const Grid_Driver& grid, const K_Vectors& kv, - hamilt::Hamilt* p_ham) + hamilt::Hamilt* p_ham, + Plus_U* p_dftu) { LCAO_HS_Arrays HS_Arrays; // store sparse arrays //! generate a file containing the Hamiltonian and S(overlap) matrices if (out_mat_hsR) { - output_HSR(ucell, istep, v_eff, pv, HS_Arrays, grid, kv, p_ham); + output_HSR(ucell, istep, v_eff, pv, HS_Arrays, grid, kv, *p_dftu, p_ham); } //! generate a file containing the kinetic energy matrix @@ -92,7 +93,9 @@ template void output_mat_sparse(const bool& out_mat_hsR, UnitCell& ucell, const Grid_Driver& grid, const K_Vectors& kv, - hamilt::Hamilt* p_ham); + hamilt::Hamilt* p_ham, + Plus_U* p_dftu); + template void output_mat_sparse>(const bool& out_mat_hsR, const bool& out_mat_dh, const bool& out_mat_ds, @@ -106,6 +109,7 @@ template void output_mat_sparse>(const bool& out_mat_hsR, UnitCell& ucell, const Grid_Driver& grid, const K_Vectors& kv, - hamilt::Hamilt>* p_ham); + hamilt::Hamilt>* p_ham, + Plus_U* p_dftu); } // namespace ModuleIO diff --git a/source/source_io/output_mat_sparse.h b/source/source_io/output_mat_sparse.h index bce47f7fb0..ec9af7af72 100644 --- a/source/source_io/output_mat_sparse.h +++ b/source/source_io/output_mat_sparse.h @@ -6,6 +6,8 @@ #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 + namespace ModuleIO { /// @brief the output interface to write the sparse matrix of H, S, T, and r @@ -23,7 +25,8 @@ void output_mat_sparse(const bool& out_mat_hsR, UnitCell& ucell, const Grid_Driver& grid, // mohan add 2024-04-06 const K_Vectors& kv, - hamilt::Hamilt* p_ham); + hamilt::Hamilt* p_ham, + Plus_U* p_dftu); // mohan add 20251107 } // namespace ModuleIO #endif // OUTPUT_MAT_SPARSE_H diff --git a/source/source_io/write_HS_R.cpp b/source/source_io/write_HS_R.cpp index 149e61f948..51ad6c5487 100644 --- a/source/source_io/write_HS_R.cpp +++ b/source/source_io/write_HS_R.cpp @@ -19,8 +19,9 @@ void ModuleIO::output_HSR(const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, // mohan add 2024-04-06 - const K_Vectors& kv, - hamilt::Hamilt* p_ham, + const K_Vectors& kv, + Plus_U &dftu, // mohan add 20251107 + hamilt::Hamilt* p_ham, #ifdef __EXX const std::vector>>>* Hexxd, const std::vector>>>>* Hexxc, @@ -48,18 +49,19 @@ void ModuleIO::output_HSR(const UnitCell& ucell, { const int spin_now = 0; // jingan add 2021-6-4, modify 2021-12-2 - sparse_format::cal_HSR(ucell, - pv, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham + sparse_format::cal_HSR(ucell, + dftu, // mohan add 20251107 + pv, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham #ifdef __EXX - , - Hexxd, - Hexxc + , + Hexxd, + Hexxc #endif ); } @@ -68,20 +70,21 @@ void ModuleIO::output_HSR(const UnitCell& ucell, int spin_now = 1; // save HR of spin down first (the current spin always be down) - sparse_format::cal_HSR(ucell, - pv, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham + sparse_format::cal_HSR(ucell, + dftu, + pv, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham #ifdef __EXX - , - Hexxd, - Hexxc + , + Hexxd, + Hexxc #endif - ); + ); // cal HR of the spin up if (PARAM.inp.vl_in_h) @@ -92,18 +95,19 @@ void ModuleIO::output_HSR(const UnitCell& ucell, spin_now = 0; } - sparse_format::cal_HSR(ucell, - pv, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham + sparse_format::cal_HSR(ucell, + dftu, + pv, + HS_Arrays, + grid, + spin_now, + sparse_thr, + kv.nmp, + p_ham #ifdef __EXX - , - Hexxd, - Hexxc + , + Hexxd, + Hexxc #endif ); } @@ -305,6 +309,7 @@ template void ModuleIO::output_HSR( LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, const K_Vectors& kv, + Plus_U &dftu, // mohan add 20251107 hamilt::Hamilt* p_ham, #ifdef __EXX const std::vector>>>* Hexxd, @@ -324,6 +329,7 @@ template void ModuleIO::output_HSR>( LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, const K_Vectors& kv, + Plus_U &dftu, // mohan add 20251107 hamilt::Hamilt>* p_ham, #ifdef __EXX const std::vector>>>* Hexxd, @@ -346,4 +352,4 @@ template void ModuleIO::output_SR>(Parallel_Orbitals& pv, hamilt::Hamilt>* p_ham, const std::string& SR_filename, const bool& binary, - const double& sparse_thr); \ No newline at end of file + const double& sparse_thr); diff --git a/source/source_io/write_HS_R.h b/source/source_io/write_HS_R.h index 9d290d5343..d92c40e59c 100644 --- a/source/source_io/write_HS_R.h +++ b/source/source_io/write_HS_R.h @@ -7,6 +7,7 @@ #include "source_hamilt/hamilt.h" #include "source_lcao/LCAO_HS_arrays.hpp" #include "source_pw/module_pwdft/global.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace ModuleIO { @@ -19,6 +20,7 @@ void output_HSR(const UnitCell& ucell, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, // mohan add 2024-04-06 const K_Vectors& kv, + Plus_U &dftu, // mohan add 20251107 hamilt::Hamilt* p_ham, #ifdef __EXX const std::vector>>>* Hexxd = nullptr, diff --git a/source/source_io/write_vxc.hpp b/source/source_io/write_vxc.hpp index ad503265c0..f579fc8b60 100644 --- a/source/source_io/write_vxc.hpp +++ b/source/source_io/write_vxc.hpp @@ -203,7 +203,7 @@ void write_Vxc(const int nspin, &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); std::vector> e_orb_exx; // orbital energy (EXX) #endif - hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, kv.isk); + hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, nullptr, kv.isk); // 4. calculate and write the MO-matrix Exc Parallel_2D p2d; diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index 54fe92ed9a..878f54264c 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -77,7 +77,8 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, const Structure_Factor& sf, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, - surchem& solvent, + surchem& solvent, + Plus_U &dftu, // mohan add 2025-11-07 Setup_DeePKS& deepks, Exx_NAO &exx_nao, ModuleSymmetry::Symmetry* symm) @@ -259,41 +260,41 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, } //! atomic forces from DFT+U (Quxin version) - ModuleBase::matrix force_dftu; - ModuleBase::matrix stress_dftu; + ModuleBase::matrix force_u; + ModuleBase::matrix stress_u; if (PARAM.inp.dft_plus_u) // Quxin add for DFT+U on 20201029 { if (isforce) { - force_dftu.create(nat, 3); + force_u.create(nat, 3); } if (isstress) { - stress_dftu.create(3, 3); + stress_u.create(3, 3); } if (PARAM.inp.dft_plus_u == 2) { - // GlobalC::dftu.force_stress(ucell, gd, pelec, pv, fsr, force_dftu, stress_dftu, kv); + // dftu.force_stress(ucell, gd, pelec, pv, fsr, force_u, stress_u, kv); // mohan modify 2025-11-03 std::vector>* dmk_d = nullptr; std::vector>>* dmk_c = nullptr; // add a new template function assign_dmk_ptr(dmat.dm, dmk_d, dmk_c, PARAM.globalv.gamma_only_local); - GlobalC::dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr, force_dftu, stress_dftu, kv); + dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr, force_u, stress_u, kv); } else { - hamilt::DFTU> tmp_dftu(nullptr, // HK and SK are not used for force&stress + hamilt::DFTU> tmpu(nullptr, // HK and SK are not used for force&stress kv.kvec_d, nullptr, // HR are not used for force&stress ucell, &gd, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), - &GlobalC::dftu); + &dftu); - tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu); + tmpu.cal_force_stress(isforce, isstress, force_u, stress_u); } } @@ -393,7 +394,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // Force contribution from DFT+U, Quxin add on 20201029 if (PARAM.inp.dft_plus_u) { - fcs(iat, i) += force_dftu(iat, i); + fcs(iat, i) += force_u(iat, i); } if (PARAM.inp.sc_mag_switch) { @@ -536,7 +537,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, } if (PARAM.inp.dft_plus_u) { - ModuleIO::print_force(GlobalV::ofs_running, ucell, "DFT+U FORCE", force_dftu, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "DFT+U FORCE", force_u, false); } if (PARAM.inp.sc_mag_switch) { @@ -604,7 +605,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // DFT plus U stress from qux if (PARAM.inp.dft_plus_u) { - scs(i, j) += stress_dftu(i, j); + scs(i, j) += stress_u(i, j); } if (PARAM.inp.sc_mag_switch) { @@ -673,7 +674,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, } if (PARAM.inp.dft_plus_u) { - ModuleIO::print_stress("DFTU STRESS", stress_dftu, screen, ry, GlobalV::ofs_running); + ModuleIO::print_stress("DFTU STRESS", stress_u, screen, ry, GlobalV::ofs_running); } if (PARAM.inp.sc_mag_switch) { diff --git a/source/source_lcao/FORCE_STRESS.h b/source/source_lcao/FORCE_STRESS.h index b62dd1460b..1c2ebb790e 100644 --- a/source/source_lcao/FORCE_STRESS.h +++ b/source/source_lcao/FORCE_STRESS.h @@ -17,6 +17,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 template @@ -49,7 +50,8 @@ class Force_Stress_LCAO const Structure_Factor& sf, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, - surchem& solvent, + surchem& solvent, + Plus_U &dftu, // mohan add 2025-11-07 Setup_DeePKS &deepks, Exx_NAO &exx_nao, ModuleSymmetry::Symmetry* symm); diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index 3946c62b2a..201dcadee9 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -78,8 +78,9 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, const K_Vectors& kv_in, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, - elecstate::DensityMatrix* DM_in, - Setup_DeePKS &deepks, + elecstate::DensityMatrix* DM_in, + Plus_U* p_dftu, // mohan add 2025-11-05 + Setup_DeePKS &deepks, const int istep, Exx_NAO &exx_nao) { @@ -212,26 +213,27 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, // end node should be OperatorDFTU if (PARAM.inp.dft_plus_u) { - Operator* dftu = nullptr; + Operator* plus_u = nullptr; if (PARAM.inp.dft_plus_u == 2) { - dftu = new OperatorDFTU>(this->hsk, + plus_u = new OperatorDFTU>(this->hsk, this->kv->kvec_d, - this->hR, // no explicit call yet - this->kv->isk); + this->hR, // no explicit call yet + p_dftu, // mohan add 2025-11-07 + this->kv->isk); } else { - dftu = new DFTU>(this->hsk, + plus_u = new DFTU>(this->hsk, this->kv->kvec_d, this->hR, ucell, &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), - &GlobalC::dftu); + p_dftu); } - this->getOperator()->add(dftu); + this->getOperator()->add(plus_u); } } // multi-k-points case to initialize HamiltLCAO, ops will be used @@ -367,26 +369,27 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, } if (PARAM.inp.dft_plus_u) { - Operator* dftu = nullptr; + Operator* plus_u = nullptr; if (PARAM.inp.dft_plus_u == 2) { - dftu = new OperatorDFTU>(this->hsk, + plus_u = new OperatorDFTU>(this->hsk, this->kv->kvec_d, - this->hR, // no explicit call yet + this->hR, // no explicit call yet + p_dftu, // mohan add 2025-11-07 this->kv->isk); } else { - dftu = new DFTU>(this->hsk, + plus_u = new DFTU>(this->hsk, this->kv->kvec_d, this->hR, ucell, &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), - &GlobalC::dftu); + p_dftu); } - this->getOperator()->add(dftu); + this->getOperator()->add(plus_u); } if (PARAM.inp.sc_mag_switch) { diff --git a/source/source_lcao/hamilt_lcao.h b/source/source_lcao/hamilt_lcao.h index dbebffe161..5f0d854544 100644 --- a/source/source_lcao/hamilt_lcao.h +++ b/source/source_lcao/hamilt_lcao.h @@ -19,6 +19,7 @@ #endif #include "source_lcao/setup_exx.h" // for exx, mohan add 20251022 +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-05 namespace hamilt { @@ -48,6 +49,7 @@ class HamiltLCAO : public Hamilt const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, elecstate::DensityMatrix* DM_in, + Plus_U* p_dftu, // mohan add 2025-11-05 Setup_DeePKS &deepks, const int istep, Exx_NAO &exx_nao); diff --git a/source/source_lcao/module_dftu/dftu.cpp b/source/source_lcao/module_dftu/dftu.cpp index 1bf05e98d9..8c990379fd 100644 --- a/source/source_lcao/module_dftu/dftu.cpp +++ b/source/source_lcao/module_dftu/dftu.cpp @@ -20,22 +20,30 @@ #include #include -namespace GlobalC -{ -ModuleDFTU::DFTU dftu; -} + // mohan add 2025-11-06 +double Plus_U::energy_u = 0.0; -namespace ModuleDFTU -{ -DFTU::DFTU() -{ -} +std::vector Plus_U::U = {}; // U (Hubbard parameter U) -DFTU::~DFTU() -{ -} +std::vector Plus_U::U0 = {}; // U0 (target Hubbard parameter U0) + +std::vector Plus_U::orbital_corr = {}; // + +double Plus_U::uramping = 0.0; // increase U by uramping, default is -1.0 -void DFTU::init(UnitCell& cell, // unitcell class +int Plus_U::omc=0; // occupation matrix control + +int Plus_U::mixing_dftu=0; //whether to mix locale + +bool Plus_U::Yukawa=false; // whether to use Yukawa potential + +Plus_U::Plus_U() +{} + +Plus_U::~Plus_U() +{} + +void Plus_U::init(UnitCell& cell, // unitcell class const Parallel_Orbitals* pv, const int nks #ifdef __LCAO @@ -43,7 +51,7 @@ void DFTU::init(UnitCell& cell, // unitcell class #endif ) { - ModuleBase::TITLE("DFTU", "init"); + ModuleBase::TITLE("Plus_U", "init"); #ifndef __MPI std::cout << "DFT+U module is only accessible in mpi versioin" << std::endl; @@ -67,7 +75,8 @@ void DFTU::init(UnitCell& cell, // unitcell class const int nlocal = PARAM.globalv.nlocal; // number of total local orbitals const int nspin = PARAM.inp.nspin; // number of spins - this->EU = 0.0; + // mohan update 2025-11-06 + Plus_U::energy_u = 0.0; this->locale.resize(cell.nat); this->locale_save.resize(cell.nat); @@ -221,24 +230,27 @@ void DFTU::init(UnitCell& cell, // unitcell class } } - ModuleBase::Memory::record("DFTU::locale", sizeof(double) * num_locale); + ModuleBase::Memory::record("Plus_U::locale", sizeof(double) * num_locale); return; } #ifdef __LCAO -void DFTU::cal_energy_correction(const UnitCell& ucell, +void Plus_U::cal_energy_correction(const UnitCell& ucell, const int istep) { - ModuleBase::TITLE("DFTU", "cal_energy_correction"); - ModuleBase::timer::tick("DFTU", "cal_energy_correction"); + ModuleBase::TITLE("Plus_U", "cal_energy_correction"); + ModuleBase::timer::tick("Plus_U", "cal_energy_correction"); if (!initialed_locale) { - ModuleBase::timer::tick("DFTU", "cal_energy_correction"); + ModuleBase::timer::tick("Plus_U", "cal_energy_correction"); return; } - this->EU = 0.0; - double EU_dc = 0.0; + + // mohan update 20251106 + Plus_U::energy_u = 0.0; + + double energy_dc = 0.0; for (int T = 0; T < ucell.ntype; T++) { @@ -291,12 +303,12 @@ void DFTU::cal_energy_correction(const UnitCell& ucell, } if (Yukawa) { - this->EU += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) + Plus_U::energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) * (nm_trace - nm2_trace); } else { - this->EU += 0.5 * this->U[T] * (nm_trace - nm2_trace); + Plus_U::energy_u += 0.5 * this->U[T] * (nm_trace - nm2_trace); } } } @@ -326,12 +338,12 @@ void DFTU::cal_energy_correction(const UnitCell& ucell, } if (Yukawa) { - this->EU - += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) * (nm_trace - nm2_trace); + Plus_U::energy_u += 0.5 * (this->U_Yukawa[T][l][n] - this->J_Yukawa[T][l][n]) + * (nm_trace - nm2_trace); } else { - this->EU += 0.5 * this->U[T] * (nm_trace - nm2_trace); + Plus_U::energy_u += 0.5 * this->U[T] * (nm_trace - nm2_trace); } } @@ -353,14 +365,14 @@ void DFTU::cal_energy_correction(const UnitCell& ucell, { double VU = 0.0; VU = get_onebody_eff_pot(T, iat, l, n, is, m1_all, m2_all, false); - EU_dc += VU * this->locale[iat][l][n][is](m1_all, m2_all); + energy_dc += VU * this->locale[iat][l][n][is](m1_all, m2_all); } } else if (PARAM.inp.nspin == 4) // SOC { double VU = 0.0; VU = get_onebody_eff_pot(T, iat, l, n, 0, m1_all, m2_all, false); - EU_dc += VU * this->locale[iat][l][n][0](m1_all, m2_all); + energy_dc += VU * this->locale[iat][l][n][0](m1_all, m2_all); } } } @@ -371,16 +383,16 @@ void DFTU::cal_energy_correction(const UnitCell& ucell, } // end I } // end T - // substract the double counting EU_dc included in band energy eband - this->EU -= EU_dc; + // substract the double counting energy_dc included in band energy eband + Plus_U::energy_u -= energy_dc; - ModuleBase::timer::tick("DFTU", "cal_energy_correction"); + ModuleBase::timer::tick("Plus_U", "cal_energy_correction"); return; } #endif -void DFTU::uramping_update() +void Plus_U::uramping_update() { // if uramping < 0.1, use the original U if (this->uramping < 0.01) { @@ -400,7 +412,7 @@ void DFTU::uramping_update() } } -bool DFTU::u_converged() +bool Plus_U::u_converged() { for (int i = 0; i < this->U0.size(); i++) { @@ -414,19 +426,19 @@ bool DFTU::u_converged() #ifdef __LCAO -void DFTU::set_dmr(const elecstate::DensityMatrix, double>* dmr) +void Plus_U::set_dmr(const elecstate::DensityMatrix, double>* dmr) { this->dm_in_dftu_cd = dmr; return; } -void DFTU::set_dmr(const elecstate::DensityMatrix* dmr) +void Plus_U::set_dmr(const elecstate::DensityMatrix* dmr) { this->dm_in_dftu_d = dmr; return; } -const hamilt::HContainer* DFTU::get_dmr(int ispin) const +const hamilt::HContainer* Plus_U::get_dmr(int ispin) const { if (this->dm_in_dftu_d != nullptr) { @@ -449,9 +461,10 @@ void dftu_cal_occup_m(const int iter, const std::vector>& dm, const K_Vectors& kv, const double& mixing_beta, - hamilt::Hamilt* p_ham) + hamilt::Hamilt* p_ham, + Plus_U &dftu) { - GlobalC::dftu.cal_occup_m_gamma(iter, ucell ,dm, mixing_beta, p_ham); + dftu.cal_occup_m_gamma(iter, ucell ,dm, mixing_beta, p_ham); } //! dftu occupation matrix for multiple k-points using dm(complex) @@ -461,11 +474,10 @@ void dftu_cal_occup_m(const int iter, const std::vector>>& dm, const K_Vectors& kv, const double& mixing_beta, - hamilt::Hamilt>* p_ham) + hamilt::Hamilt>* p_ham, + Plus_U &dftu) { - GlobalC::dftu.cal_occup_m_k(iter,ucell, dm, kv, mixing_beta, p_ham); + dftu.cal_occup_m_k(iter,ucell, dm, kv, mixing_beta, p_ham); } #endif - -} // namespace ModuleDFTU diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h index 457debc8bd..56d213386d 100644 --- a/source/source_lcao/module_dftu/dftu.h +++ b/source/source_lcao/module_dftu/dftu.h @@ -7,7 +7,6 @@ #ifdef __LCAO #include "source_estate/module_charge/charge_mixing.h" #include "source_hamilt/hamilt.h" -// #include "source_estate/elecstate.h" // mohan update 2025-11-02 #include "source_lcao/module_hcontainer/hcontainer.h" #include "source_estate/module_dm/density_matrix.h" #include "source_lcao/force_stress_arrays.h" // mohan add 2024-06-15 @@ -16,24 +15,14 @@ #include #include -//========================================================== -// CLASS : -// NAME : DTFU (DFT+U) -//========================================================== -namespace ModuleDFTU -{ -class DFTU +class Plus_U { public: - DFTU(); // constructor - ~DFTU(); // deconstructor + Plus_U(); + ~Plus_U(); - //============================================================= - // In dftu.cpp - // Initialization & Calculating energy - //============================================================= public: // allocate relevant data strcutures void init(UnitCell& cell, // unitcell class @@ -44,23 +33,41 @@ class DFTU #endif ); - static DFTU* get_instance(); - // calculate the energy correction void cal_energy_correction(const UnitCell& ucell, const int istep); - double get_energy(){return EU;} + + // mohan change the function to static, 20251106 + static double get_energy() + { + return Plus_U::energy_u; + } + + static void set_energy(const double &e) + { + Plus_U::energy_u = e; + } + + static void set_double_energy() // mohan add 20251107 + { + Plus_U::energy_u *= 2.0; + } + void uramping_update(); // update U by uramping bool u_converged(); // check if U is converged - std::vector U = {}; // U (Hubbard parameter U) - std::vector U0; // U0 (target Hubbard parameter U0) - std::vector orbital_corr = {}; // - double uramping; // increase U by uramping, default is -1.0 - int omc; // occupation matrix control - int mixing_dftu; //whether to mix locale + // mohan change these parameters to static, 2025-11-07 + static std::vector U; // U (Hubbard parameter U) + static std::vector U0; // U0 (target Hubbard parameter U0) + static std::vector orbital_corr; // + static double uramping; // increase U by uramping, default is -1.0 + static int omc; // occupation matrix control + static int mixing_dftu; //whether to mix locale - double EU; //+U energy private: + + // mohan change the variable to static, 20251106 + static double energy_u; //+U energy, mohan update 2025-11-06, change this to private + const Parallel_Orbitals* paraV = nullptr; int cal_type = 3; // 1:dftu_tpye=1, dc=1; 2:dftu_type=1, dc=2; 3:dftu_tpye=2, dc=1; 4:dftu_tpye=2, dc=2; @@ -82,10 +89,21 @@ class DFTU // 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); - void cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk); + void cal_eff_pot_mat_complex(const int ik, + std::complex* eff_pot, + const std::vector& isk, + const std::complex* sk); + + void cal_eff_pot_mat_real(const int ik, + double* eff_pot, + const std::vector& isk, + const double* sk); + void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR); - void cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR); + + void cal_eff_pot_mat_R_complex_double(const int ispin, + std::complex* SR, + std::complex* HR); #endif //============================================================= @@ -95,13 +113,26 @@ class DFTU //============================================================= public: /// interface for PW base - /// calculate the local occupation number matrix for PW based wave functions - void cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matrix& wg_in, const UnitCell& cell, const double& mixing_beta); + /// calculate the local occupation number matrix for PW based wave functions + void cal_occ_pw(const int iter, + const void* psi_in, + const ModuleBase::matrix& wg_in, + const UnitCell& cell, + const double& mixing_beta); + /// calculate the local DFT+U effective potential matrix for PW base. void cal_VU_pot_pw(const int spin); + /// get effective potential matrix for PW base - const std::complex* get_eff_pot_pw(const int iat) const { return &(eff_pot_pw[this->eff_pot_pw_index[iat]]); } - int get_size_eff_pot_pw() const { return eff_pot_pw.size(); } + const std::complex* get_eff_pot_pw(const int iat) const + { + return &(eff_pot_pw[this->eff_pot_pw_index[iat]]); + } + + int get_size_eff_pot_pw() const + { + return eff_pot_pw.size(); + } #ifdef __LCAO // calculate the local occupation number matrix @@ -111,6 +142,7 @@ class DFTU 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, @@ -122,6 +154,7 @@ class DFTU bool initialed_locale = false; private: + void copy_locale(const UnitCell& ucell); void zero_locale(const UnitCell& ucell); void mix_locale(const UnitCell& ucell,const double& mixing_beta); @@ -129,12 +162,13 @@ class DFTU std::vector> eff_pot_pw; std::vector eff_pot_pw_index; -public: - // local occupancy matrix of the correlated subspace + public: + // local occupancy matrix of the correlated subspace // locale: the out put local occupation number matrix of correlated electrons in the current electronic step // locale_save: the input local occupation number matrix of correlated electrons in the current electronic step std::vector>>> locale; // locale[iat][l][n][spin](m1,m2) std::vector>>> locale_save; // locale_save[iat][l][n][spin](m1,m2) + #ifdef __LCAO private: //============================================================= @@ -267,7 +301,7 @@ class DFTU //============================================================= public: - bool Yukawa; // 1:use Yukawa potential; 0: do not use Yukawa potential + static bool Yukawa; // 1:use Yukawa potential; 0: do not use Yukawa potential void cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx); private: @@ -304,6 +338,7 @@ class DFTU #endif }; + #ifdef __LCAO template void dftu_cal_occup_m(const int iter, @@ -311,13 +346,9 @@ void dftu_cal_occup_m(const int iter, const std::vector>& dm, const K_Vectors& kv, const double& mixing_beta, - hamilt::Hamilt* p_ham); + hamilt::Hamilt* p_ham, + Plus_U &dftu); #endif -} // namespace ModuleDFTU -namespace GlobalC -{ - extern ModuleDFTU::DFTU dftu; -} #endif diff --git a/source/source_lcao/module_dftu/dftu_folding.cpp b/source/source_lcao/module_dftu/dftu_folding.cpp index 26e70d0d55..bd2c8ad44a 100644 --- a/source/source_lcao/module_dftu/dftu_folding.cpp +++ b/source/source_lcao/module_dftu/dftu_folding.cpp @@ -8,10 +8,7 @@ #include "source_lcao/module_hcontainer/hcontainer.h" #include "source_lcao/module_hcontainer/hcontainer_funcs.h" -namespace ModuleDFTU -{ - -void DFTU::fold_dSR_gamma(const UnitCell& ucell, +void Plus_U::fold_dSR_gamma(const UnitCell& ucell, const Parallel_Orbitals& pv, const Grid_Driver* gd, double* dsloc_x, @@ -22,7 +19,7 @@ void DFTU::fold_dSR_gamma(const UnitCell& ucell, const int dim2, double* dSR_gamma) { - ModuleBase::TITLE("DFTU", "fold_dSR_gamma"); + ModuleBase::TITLE("Plus_U", "fold_dSR_gamma"); ModuleBase::GlobalFunc::ZEROS(dSR_gamma, pv.nloc); @@ -125,7 +122,7 @@ void DFTU::fold_dSR_gamma(const UnitCell& ucell, return; } -void DFTU::folding_matrix_k(const UnitCell& ucell, +void Plus_U::folding_matrix_k(const UnitCell& ucell, const Grid_Driver& gd, ForceStressArrays& fsr, const Parallel_Orbitals& pv, @@ -135,11 +132,11 @@ void DFTU::folding_matrix_k(const UnitCell& ucell, std::complex* mat_k, const ModuleBase::Vector3& kvec_d) { - ModuleBase::TITLE("DFTU", "folding_matrix_k"); - ModuleBase::timer::tick("DFTU", "folding_matrix_k"); + ModuleBase::TITLE("Plus_U", "folding_matrix_k"); + ModuleBase::timer::tick("Plus_U", "folding_matrix_k"); ModuleBase::GlobalFunc::ZEROS(mat_k, pv.nloc); - double* mat_ptr; + double* mat_ptr = nullptr; if (dim1 == 1 || dim1 == 4) { mat_ptr = fsr.DSloc_Rx; @@ -277,16 +274,16 @@ void DFTU::folding_matrix_k(const UnitCell& ucell, } // ad } // I1 } // T1 - ModuleBase::timer::tick("DFTU", "folding_matrix_k"); + ModuleBase::timer::tick("Plus_U", "folding_matrix_k"); return; } -void DFTU::folding_matrix_k_new(const int ik, +void Plus_U::folding_matrix_k_new(const int ik, hamilt::Hamilt>* p_ham) { - ModuleBase::TITLE("DFTU", "folding_matrix_k_new"); - ModuleBase::timer::tick("DFTU", "folding_matrix_k_new"); + ModuleBase::TITLE("Plus_U", "folding_matrix_k_new"); + ModuleBase::timer::tick("Plus_U", "folding_matrix_k_new"); int hk_type = 0; if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) @@ -313,7 +310,8 @@ void DFTU::folding_matrix_k_new(const int ik, ->updateSk(ik, hk_type); } } + + ModuleBase::timer::tick("Plus_U", "folding_matrix_k_new"); } -} // namespace ModuleDFTU #endif // __LCAO diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_force.cpp index 5d76611470..b54026663a 100644 --- a/source/source_lcao/module_dftu/dftu_force.cpp +++ b/source/source_lcao/module_dftu/dftu_force.cpp @@ -67,10 +67,8 @@ extern "C" const int* DESCC); } -namespace ModuleDFTU -{ -void DFTU::force_stress(const UnitCell& ucell, +void Plus_U::force_stress(const UnitCell& ucell, const Grid_Driver& gd, std::vector>* dmk_d, // mohan modify 2025-11-02 std::vector>>* dmk_c, // dmat.get_dm()->get_DMK_vector(); @@ -80,8 +78,8 @@ void DFTU::force_stress(const UnitCell& ucell, ModuleBase::matrix& stress_dftu, const K_Vectors& kv) { - ModuleBase::TITLE("DFTU", "force_stress"); - ModuleBase::timer::tick("DFTU", "force_stress"); + ModuleBase::TITLE("Plus_U", "force_stress"); + ModuleBase::timer::tick("Plus_U", "force_stress"); const int nlocal = PARAM.globalv.nlocal; @@ -207,12 +205,12 @@ void DFTU::force_stress(const UnitCell& ucell, } } } - ModuleBase::timer::tick("DFTU", "force_stress"); + ModuleBase::timer::tick("Plus_U", "force_stress"); return; } -void DFTU::cal_force_k(const UnitCell& ucell, +void Plus_U::cal_force_k(const UnitCell& ucell, const Grid_Driver& gd, ForceStressArrays& fsr, const Parallel_Orbitals& pv, @@ -221,8 +219,8 @@ void DFTU::cal_force_k(const UnitCell& ucell, ModuleBase::matrix& force_dftu, const ModuleBase::Vector3& kvec_d) { - ModuleBase::TITLE("DFTU", "cal_force_k"); - ModuleBase::timer::tick("DFTU", "cal_force_k"); + ModuleBase::TITLE("Plus_U", "cal_force_k"); + ModuleBase::timer::tick("Plus_U", "cal_force_k"); const char transN = 'N'; const char transC = 'C'; @@ -337,12 +335,12 @@ void DFTU::cal_force_k(const UnitCell& ucell, } // ia } // it } // end dim - ModuleBase::timer::tick("DFTU", "cal_force_k"); + ModuleBase::timer::tick("Plus_U", "cal_force_k"); return; } -void DFTU::cal_stress_k(const UnitCell& ucell, +void Plus_U::cal_stress_k(const UnitCell& ucell, const Grid_Driver& gd, ForceStressArrays& fsr, const Parallel_Orbitals& pv, @@ -351,8 +349,8 @@ void DFTU::cal_stress_k(const UnitCell& ucell, ModuleBase::matrix& stress_dftu, const ModuleBase::Vector3& kvec_d) { - ModuleBase::TITLE("DFTU", "cal_stress_k"); - ModuleBase::timer::tick("DFTU", "cal_stress_k"); + ModuleBase::TITLE("Plus_U", "cal_stress_k"); + ModuleBase::timer::tick("Plus_U", "cal_stress_k"); const int nlocal = PARAM.globalv.nlocal; @@ -408,12 +406,12 @@ void DFTU::cal_stress_k(const UnitCell& ucell, } // end dim2 } // end dim1 - ModuleBase::timer::tick("DFTU", "cal_stress_k"); + ModuleBase::timer::tick("Plus_U", "cal_stress_k"); return; } -void DFTU::cal_force_gamma(const UnitCell& ucell, +void Plus_U::cal_force_gamma(const UnitCell& ucell, const double* rho_VU, const Parallel_Orbitals& pv, double* dsloc_x, @@ -421,8 +419,8 @@ void DFTU::cal_force_gamma(const UnitCell& ucell, double* dsloc_z, ModuleBase::matrix& force_dftu) { - ModuleBase::TITLE("DFTU", "cal_force_gamma"); - ModuleBase::timer::tick("DFTU", "cal_force_gamma"); + ModuleBase::TITLE("Plus_U", "cal_force_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_force_gamma"); const char transN = 'N', transT = 'T'; const int one_int = 1; const double one = 1.0, zero = 0.0, minus_one = -1.0; @@ -548,12 +546,12 @@ void DFTU::cal_force_gamma(const UnitCell& ucell, } // it } // end dim - ModuleBase::timer::tick("DFTU", "cal_force_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_force_gamma"); return; } -void DFTU::cal_stress_gamma(const UnitCell& ucell, +void Plus_U::cal_stress_gamma(const UnitCell& ucell, const Parallel_Orbitals& pv, const Grid_Driver* gd, double* dsloc_x, @@ -563,8 +561,8 @@ void DFTU::cal_stress_gamma(const UnitCell& ucell, const double* rho_VU, ModuleBase::matrix& stress_dftu) { - ModuleBase::TITLE("DFTU", "cal_stress_gamma"); - ModuleBase::timer::tick("DFTU", "cal_stress_gamma"); + ModuleBase::TITLE("Plus_U", "cal_stress_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_stress_gamma"); const char transN = 'N'; const int one_int = 1; @@ -621,8 +619,7 @@ void DFTU::cal_stress_gamma(const UnitCell& ucell, } // end dim2 } // end dim1 - ModuleBase::timer::tick("DFTU", "cal_stress_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_stress_gamma"); return; } -} // namespace ModuleDFTU #endif diff --git a/source/source_lcao/module_dftu/dftu_hamilt.cpp b/source/source_lcao/module_dftu/dftu_hamilt.cpp index 9ba7fe1b14..7731aafdb7 100644 --- a/source/source_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/source_lcao/module_dftu/dftu_hamilt.cpp @@ -4,20 +4,21 @@ #include "source_base/timer.h" #include "source_pw/module_pwdft/global.h" -namespace ModuleDFTU -{ #ifdef __LCAO -void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::vector& isk, const std::complex* sk) +void Plus_U::cal_eff_pot_mat_complex(const int ik, + std::complex* eff_pot, + const std::vector& isk, + const std::complex* sk) { - ModuleBase::TITLE("DFTU", "cal_eff_pot_mat"); - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); + ModuleBase::TITLE("Plus_U", "cal_eff_pot_c"); if (!this->initialed_locale) { - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); return; } + ModuleBase::timer::tick("Plus_U", "cal_eff_pot_c"); + int spin = isk[ik]; ModuleBase::GlobalFunc::ZEROS(eff_pot, this->paraV->nloc); @@ -44,8 +45,10 @@ void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, eff_pot, &one_int, &one_int, this->paraV->desc); #endif - for (int irc = 0; irc < this->paraV->nloc; irc++) - VU[irc] = eff_pot[irc]; + for (int irc = 0; irc < this->paraV->nloc; irc++) + { + VU[irc] = eff_pot[irc]; + } #ifdef __MPI pztranc_(&PARAM.globalv.nlocal, &PARAM.globalv.nlocal, @@ -55,19 +58,18 @@ void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, eff_pot, &one_int, &one_int, this->paraV->desc); #endif - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); + ModuleBase::timer::tick("Plus_U", "cal_eff_pot_c"); return; } -void DFTU::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk) +void Plus_U::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk) { - ModuleBase::TITLE("DFTU", "cal_eff_pot_mat"); - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); + ModuleBase::TITLE("Plus_U", "cal_eff_pot_r"); if (!this->initialed_locale) { - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); return; } + ModuleBase::timer::tick("Plus_U", "cal_eff_pot_r"); int spin = isk[ik]; @@ -104,11 +106,11 @@ void DFTU::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector eff_pot, &one_int, &one_int, const_cast(this->paraV->desc)); #endif - ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); + ModuleBase::timer::tick("Plus_U", "cal_eff_pot_r"); return; } -void DFTU::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) +void Plus_U::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) { const char transN = 'N', transT = 'T'; const int one_int = 1; @@ -138,7 +140,7 @@ void DFTU::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) return; } -void DFTU::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR) +void Plus_U::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR) { const char transN = 'N', transT = 'T'; const int one_int = 1; @@ -169,4 +171,3 @@ void DFTU::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex= orbital_corr[T] && orbital_corr[T] != -1) { - if (L != orbital_corr[T]) { - continue; -} + if (L != orbital_corr[T]) + { + continue; + } if (!Yukawa) { @@ -35,10 +34,11 @@ void DFTU::output(const UnitCell &ucell) { for (int n = 0; n < N; n++) { - if (n != 0) { - continue; -} - double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV; + if (n != 0) + { + continue; + } + double Ueff = (this->U_Yukawa[T][L][n] - this->J_Yukawa[T][L][n]) * ModuleBase::Ry_to_eV; GlobalV::ofs_running << "atom_type=" << T << " L=" << L << " chi=" << n << " U=" << this->U_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV " << "J=" << this->J_Yukawa[T][L][n] * ModuleBase::Ry_to_eV << "eV" @@ -61,7 +61,7 @@ void DFTU::output(const UnitCell &ucell) } } if(!ofdftu){ - std::cout << "DFTU::write_occup_m. Can't create file onsite.dm!" << std::endl; + std::cout << "Plus_U::write_occup_m. Can't create file onsite.dm!" << std::endl; exit(0); } this->write_occup_m(ucell,ofdftu); @@ -73,21 +73,24 @@ void DFTU::output(const UnitCell &ucell) // define the function calculate the eigenvalues of a matrix std::vector CalculateEigenvalues(std::vector>& A, int n); -void DFTU::write_occup_m(const UnitCell& ucell, +void Plus_U::write_occup_m(const UnitCell& ucell, std::ofstream &ofs, bool diag) { - ModuleBase::TITLE("DFTU", "write_occup_m"); + ModuleBase::TITLE("Plus_U", "write_occup_m"); - if(GlobalV::MY_RANK != 0) { return; -} + if(GlobalV::MY_RANK != 0) + { + return; + } for (int T = 0; T < ucell.ntype; T++) { - if (orbital_corr[T] == -1) { - continue; -} - const int NL = ucell.atoms[T].nwl + 1; + if (orbital_corr[T] == -1) + { + continue; + } + const int NL = ucell.atoms[T].nwl + 1; const int LC = orbital_corr[T]; for (int I = 0; I < ucell.atoms[T].na; I++) @@ -98,9 +101,10 @@ void DFTU::write_occup_m(const UnitCell& ucell, for (int l = 0; l < NL; l++) { - if (l != orbital_corr[T]) { - continue; -} + if (l != orbital_corr[T]) + { + continue; + } const int N = ucell.atoms[T].l_nchi[l]; ofs << "L" @@ -109,9 +113,10 @@ void DFTU::write_occup_m(const UnitCell& ucell, for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } ofs << "zeta" << " " << n << std::endl; @@ -158,7 +163,8 @@ void DFTU::write_occup_m(const UnitCell& ucell, } if(diag) { - ofs << std::setw(12) << std::setprecision(8) << std::fixed<< "atomic mag: "< 0) { std::cout - << "DFTU::read_occup_m. Can not find the file initial_onsite.dm . Please check your initial_onsite.dm" + << "Plus_U::read_occup_m. Can not find the file initial_onsite.dm . Please check your initial_onsite.dm" << std::endl; } else { if (PARAM.inp.init_chg == "file") { - std::cout << "DFTU::read_occup_m. Can not find the file onsite.dm . Please do scf calculation first" + std::cout << "Plus_U::read_occup_m. Can not find the file onsite.dm . Please do scf calculation first" << std::endl; } } @@ -257,16 +267,21 @@ void DFTU::read_occup_m(const UnitCell& ucell, char word[10]; - int T, iat, spin, L, zeta; + int T=0; + int iat=0; + int spin=0; + int L=0; + int zeta=0; ifdftu.rdstate(); while (ifdftu.good()) { ifdftu >> word; - if (ifdftu.eof()) { - break; -} + if (ifdftu.eof()) + { + break; + } if (strcmp("atoms", word) == 0) { @@ -279,9 +294,10 @@ void DFTU::read_occup_m(const UnitCell& ucell, for (int l = 0; l < NL; l++) { - if (l != orbital_corr[T]) { - continue; -} + if (l != orbital_corr[T]) + { + continue; + } ifdftu >> word; @@ -294,9 +310,10 @@ void DFTU::read_occup_m(const UnitCell& ucell, for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } ifdftu >> word; if (strcmp("zeta", word) == 0) @@ -327,7 +344,7 @@ void DFTU::read_occup_m(const UnitCell& ucell, } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM DFTU FILE" + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; exit(0); } @@ -358,21 +375,21 @@ void DFTU::read_occup_m(const UnitCell& ucell, } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM DFTU FILE" << std::endl; + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; exit(0); } } } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM DFTU FILE" << std::endl; + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; exit(0); } } } else { - std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM DFTU FILE" << std::endl; + std::cout << "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE" << std::endl; exit(0); } @@ -387,15 +404,16 @@ void DFTU::read_occup_m(const UnitCell& ucell, return; } -void DFTU::local_occup_bcast(const UnitCell& ucell) +void Plus_U::local_occup_bcast(const UnitCell& ucell) { - ModuleBase::TITLE("DFTU", "local_occup_bcast"); + ModuleBase::TITLE("Plus_U", "local_occup_bcast"); for (int T = 0; T < ucell.ntype; T++) { - if (orbital_corr[T] == -1) { - continue; -} + if (orbital_corr[T] == -1) + { + continue; + } for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -404,16 +422,18 @@ void DFTU::local_occup_bcast(const UnitCell& ucell) for (int l = 0; l <= ucell.atoms[T].nwl; l++) { - if (l != orbital_corr[T]) { - continue; -} + if (l != orbital_corr[T]) + { + continue; + } for (int n = 0; n < ucell.atoms[T].l_nchi[l]; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { @@ -462,9 +482,11 @@ void DFTU::local_occup_bcast(const UnitCell& ucell) return; } -inline void JacobiRotate(std::vector>& A, int p, int q, int n) { - if (std::abs(A[p][q]) > 1e-10) { - double r = (A[q][q] - A[p][p]) / (2.0 * A[p][q]); +inline void JacobiRotate(std::vector>& A, int p, int q, int n) +{ + if (std::abs(A[p][q]) > 1e-10) + { + double r = (A[q][q] - A[p][p]) / (2.0 * A[p][q]); double t; if (r >= 0) { t = 1.0 / (r + sqrt(1.0 + r * r)); @@ -489,7 +511,8 @@ inline void JacobiRotate(std::vector>& A, int p, int q, int } } -inline std::vector CalculateEigenvalues(std::vector>& A, int n) { +inline std::vector CalculateEigenvalues(std::vector>& A, int n) +{ std::vector eigenvalues(n); while (true) { int p = 0, q = 1; @@ -513,4 +536,3 @@ inline std::vector CalculateEigenvalues(std::vector> } return eigenvalues; } -} // namespace ModuleDFTU \ No newline at end of file diff --git a/source/source_lcao/module_dftu/dftu_occup.cpp b/source/source_lcao/module_dftu/dftu_occup.cpp index 1bc771d42f..aaabbb81ac 100644 --- a/source/source_lcao/module_dftu/dftu_occup.cpp +++ b/source/source_lcao/module_dftu/dftu_occup.cpp @@ -7,12 +7,10 @@ #endif #include "source_base/module_external/scalapack_connector.h" -namespace ModuleDFTU +void Plus_U::copy_locale(const UnitCell& ucell) { -void DFTU::copy_locale(const UnitCell& ucell) -{ - ModuleBase::TITLE("DFTU", "copy_locale"); - ModuleBase::timer::tick("DFTU", "copy_locale"); + ModuleBase::TITLE("Plus_U", "copy_locale"); + ModuleBase::timer::tick("Plus_U", "copy_locale"); for (int T = 0; T < ucell.ntype; T++) { @@ -44,13 +42,13 @@ void DFTU::copy_locale(const UnitCell& ucell) } } } - ModuleBase::timer::tick("DFTU", "copy_locale"); + ModuleBase::timer::tick("Plus_U", "copy_locale"); } -void DFTU::zero_locale(const UnitCell& ucell) +void Plus_U::zero_locale(const UnitCell& ucell) { - ModuleBase::TITLE("DFTU", "zero_locale"); - ModuleBase::timer::tick("DFTU", "zero_locale"); + ModuleBase::TITLE("Plus_U", "zero_locale"); + ModuleBase::timer::tick("Plus_U", "zero_locale"); for (int T = 0; T < ucell.ntype; T++) { @@ -82,14 +80,14 @@ void DFTU::zero_locale(const UnitCell& ucell) } } } - ModuleBase::timer::tick("DFTU", "zero_locale"); + ModuleBase::timer::tick("Plus_U", "zero_locale"); } -void DFTU::mix_locale(const UnitCell& ucell, +void Plus_U::mix_locale(const UnitCell& ucell, const double& mixing_beta) { - ModuleBase::TITLE("DFTU", "mix_locale"); - ModuleBase::timer::tick("DFTU", "mix_locale"); + ModuleBase::TITLE("Plus_U", "mix_locale"); + ModuleBase::timer::tick("Plus_U", "mix_locale"); double beta = mixing_beta; @@ -123,20 +121,20 @@ void DFTU::mix_locale(const UnitCell& ucell, } } } - ModuleBase::timer::tick("DFTU", "mix_locale"); + ModuleBase::timer::tick("Plus_U", "mix_locale"); } #ifdef __LCAO -void DFTU::cal_occup_m_k(const int iter, +void Plus_U::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) { - ModuleBase::TITLE("DFTU", "cal_occup_m_k"); - ModuleBase::timer::tick("DFTU", "cal_occup_m_k"); + ModuleBase::TITLE("Plus_U", "cal_occup_m_k"); + ModuleBase::timer::tick("Plus_U", "cal_occup_m_k"); this->copy_locale(ucell); this->zero_locale(ucell); @@ -380,18 +378,18 @@ void DFTU::cal_occup_m_k(const int iter, } this->initialed_locale = true; - ModuleBase::timer::tick("DFTU", "cal_occup_m_k"); + ModuleBase::timer::tick("Plus_U", "cal_occup_m_k"); return; } -void DFTU::cal_occup_m_gamma(const int iter, +void Plus_U::cal_occup_m_gamma(const int iter, const UnitCell &ucell, const std::vector> &dm_gamma, const double& mixing_beta, hamilt::Hamilt* p_ham) { - ModuleBase::TITLE("DFTU", "cal_occup_m_gamma"); - ModuleBase::timer::tick("DFTU", "cal_occup_m_gamma"); + ModuleBase::TITLE("Plus_U", "cal_occup_m_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_occup_m_gamma"); this->copy_locale(ucell); this->zero_locale(ucell); @@ -428,26 +426,6 @@ void DFTU::cal_occup_m_gamma(const int iter, one_int, one_int, &this->paraV->desc[0]); - /*pdgemm_(&transN, - &transT, - &PARAM.globalv.nlocal, - &PARAM.globalv.nlocal, - &PARAM.globalv.nlocal, - &alpha, - s_gamma_pointer, - &one_int, - &one_int, - this->paraV->desc, - dm_gamma[is].data(), - //dm_gamma[is].c, - &one_int, - &one_int, - this->paraV->desc, - &beta, - &srho[0], - &one_int, - &one_int, - this->paraV->desc);*/ #endif for (int it = 0; it < ucell.ntype; it++) @@ -558,8 +536,7 @@ void DFTU::cal_occup_m_gamma(const int iter, } this->initialed_locale = true; - ModuleBase::timer::tick("DFTU", "cal_occup_m_gamma"); + ModuleBase::timer::tick("Plus_U", "cal_occup_m_gamma"); return; } #endif -} // namespace ModuleDFTU diff --git a/source/source_lcao/module_dftu/dftu_pw.cpp b/source/source_lcao/module_dftu/dftu_pw.cpp index 03968cabee..e3107a313a 100644 --- a/source/source_lcao/module_dftu/dftu_pw.cpp +++ b/source/source_lcao/module_dftu/dftu_pw.cpp @@ -5,16 +5,14 @@ #include "source_base/timer.h" -namespace ModuleDFTU -{ -DFTU* DFTU::get_instance() -{ - return &GlobalC::dftu; -} /// calculate occupation matrix for DFT+U -void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matrix& wg_in, const UnitCell& cell, const double& mixing_beta) +void Plus_U::cal_occ_pw(const int iter, + const void* psi_in, + const ModuleBase::matrix& wg_in, + const UnitCell& cell, + const double& mixing_beta) { - ModuleBase::timer::tick("DFTU", "cal_occ_pw"); + ModuleBase::timer::tick("Plus_U", "cal_occ_pw"); this->copy_locale(cell); this->zero_locale(cell); @@ -137,7 +135,7 @@ void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matr } #endif - this->EU = 0.0; + Plus_U::energy_u = 0.0; // reduce mag from all k-pools for(int iat = 0; iat < cell.nat; iat++) { @@ -148,7 +146,12 @@ void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matr continue; } const int size = (2 * target_l + 1) * (2 * target_l + 1); - Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, PARAM.globalv.nproc_in_pool, this->locale[iat][target_l][0][0].c, size * PARAM.inp.nspin); + + Parallel_Reduce::reduce_double_allpool(PARAM.inp.kpar, + PARAM.globalv.nproc_in_pool, + this->locale[iat][target_l][0][0].c, + size * PARAM.inp.nspin); + //update effective potential const double u_value = this->U[it]; std::complex* vu_iat = &(this->eff_pot_pw[this->eff_pot_pw_index[iat]]); @@ -157,8 +160,10 @@ void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matr { for (int m2 = 0; m2 < m_size; m2++) { - vu_iat[m1 * m_size + m2] = u_value * (1.0 * (m1 == m2) - this->locale[iat][target_l][0][0].c[m2 * m_size + m1]); - this->EU += u_value * 0.25 * this->locale[iat][target_l][0][0].c[m2 * m_size + m1] * this->locale[iat][target_l][0][0].c[m1 * m_size + m2]; + vu_iat[m1 * m_size + m2] = u_value * + (1.0 * (m1 == m2) - this->locale[iat][target_l][0][0].c[m2 * m_size + m1]); + Plus_U::energy_u += u_value * 0.25 * this->locale[iat][target_l][0][0].c[m2 * m_size + m1] + * this->locale[iat][target_l][0][0].c[m1 * m_size + m2]; } } for (int is = 1; is < 4; ++is) @@ -168,8 +173,11 @@ void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matr { for (int m2 = 0; m2 < m_size; m2++) { - vu_iat[start + m1 * m_size + m2] = u_value * (0 - this->locale[iat][target_l][0][0].c[start + m2 * m_size + m1]); - this->EU += u_value * 0.25 * this->locale[iat][target_l][0][0].c[start + m2 * m_size + m1] * this->locale[iat][target_l][0][0].c[start + m1 * m_size + m2]; + vu_iat[start + m1 * m_size + m2] = u_value * + (0 - this->locale[iat][target_l][0][0].c[start + m2 * m_size + m1]); + Plus_U::energy_u += u_value * 0.25 + * this->locale[iat][target_l][0][0].c[start + m2 * m_size + m1] + * this->locale[iat][target_l][0][0].c[start + m1 * m_size + m2]; } } } @@ -201,12 +209,11 @@ void DFTU::cal_occ_pw(const int iter, const void* psi_in, const ModuleBase::matr this->mix_locale(cell, mixing_beta); } // update effective potential - ModuleBase::timer::tick("DFTU", "cal_occ_pw"); + ModuleBase::timer::tick("Plus_U", "cal_occ_pw"); } /// calculate the local DFT+U effective potential matrix for PW base. -void DFTU::cal_VU_pot_pw(const int spin) +void Plus_U::cal_VU_pot_pw(const int spin) { } -} // namespace ModuleDFTU \ No newline at end of file diff --git a/source/source_lcao/module_dftu/dftu_tools.cpp b/source/source_lcao/module_dftu/dftu_tools.cpp index dddb347e71..3ceece7ae0 100644 --- a/source/source_lcao/module_dftu/dftu_tools.cpp +++ b/source/source_lcao/module_dftu/dftu_tools.cpp @@ -3,13 +3,11 @@ #include "source_io/module_parameter/parameter.h" #include "source_pw/module_pwdft/global.h" -namespace ModuleDFTU -{ #ifdef __LCAO -void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU) +void Plus_U::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU) { - ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_complex"); + ModuleBase::TITLE("Plus_U", "cal_VU_pot_mat_complex"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); for (int it = 0; it < this->ucell->ntype; ++it) @@ -71,9 +69,9 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com return; } -void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) +void Plus_U::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) { - ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_real"); + ModuleBase::TITLE("Plus_U", "cal_VU_pot_mat_real"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); for (int it = 0; it < this->ucell->ntype; ++it) @@ -136,7 +134,7 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) return; } -double DFTU::get_onebody_eff_pot(const int T, +double Plus_U::get_onebody_eff_pot(const int T, const int iat, const int L, const int N, @@ -145,7 +143,7 @@ double DFTU::get_onebody_eff_pot(const int T, const int m1, const bool newlocale) { - ModuleBase::TITLE("DFTU", "get_onebody_eff_pot"); + ModuleBase::TITLE("Plus_U", "get_onebody_eff_pot"); double VU = 0.0; @@ -213,4 +211,3 @@ double DFTU::get_onebody_eff_pot(const int T, return VU; } #endif -} // namespace ModuleDFTU \ No newline at end of file diff --git a/source/source_lcao/module_dftu/dftu_yukawa.cpp b/source/source_lcao/module_dftu/dftu_yukawa.cpp index 73a9910052..4f124e9509 100644 --- a/source/source_lcao/module_dftu/dftu_yukawa.cpp +++ b/source/source_lcao/module_dftu/dftu_yukawa.cpp @@ -1,9 +1,5 @@ -//========================================================== -// Author:Xin Qu #ifdef __LCAO #include "source_io/module_parameter/parameter.h" -// DATE : 2019-12-10 -//========================================================== #include "source_base/constants.h" #include "source_base/global_function.h" #include "source_pw/module_pwdft/global.h" @@ -18,12 +14,10 @@ #include #include -namespace ModuleDFTU -{ -void DFTU::cal_yukawa_lambda(double** rho, const int& nrxx) +void Plus_U::cal_yukawa_lambda(double** rho, const int& nrxx) { - ModuleBase::TITLE("DFTU", "cal_yukawa_lambda"); + ModuleBase::TITLE("Plus_U", "cal_yukawa_lambda"); if (PARAM.inp.yukawa_lambda > 0) { @@ -35,8 +29,10 @@ void DFTU::cal_yukawa_lambda(double** rho, const int& nrxx) double sum_rho_lambda = 0.0; for (int is = 0; is < PARAM.inp.nspin; is++) { - if(PARAM.inp.nspin == 4 && is > 0) { continue;// for non-collinear spin case, first spin contains the charge density -} + if(PARAM.inp.nspin == 4 && is > 0) + { + continue;// for non-collinear spin case, first spin contains the charge density + } for (int ir = 0; ir < nrxx; ir++) { double rho_ir = rho[is][ir]; @@ -63,11 +59,11 @@ void DFTU::cal_yukawa_lambda(double** rho, const int& nrxx) return; } -void DFTU::cal_slater_Fk(const UnitCell& ucell, +void Plus_U::cal_slater_Fk(const UnitCell& ucell, const int L, const int T) { - ModuleBase::TITLE("DFTU", "cal_slater_Fk"); + ModuleBase::TITLE("Plus_U", "cal_slater_Fk"); if (Yukawa) { @@ -113,12 +109,13 @@ void DFTU::cal_slater_Fk(const UnitCell& ucell, return; } -void DFTU::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) +void Plus_U::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) { - ModuleBase::TITLE("DFTU", "cal_slater_UJ"); - if (!Yukawa) { - return; -} + ModuleBase::TITLE("Plus_U", "cal_slater_UJ"); + if (!Yukawa) + { + return; + } this->cal_yukawa_lambda(rho, nrxx); @@ -146,16 +143,18 @@ void DFTU::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) if (L >= PARAM.inp.orbital_corr[T] && PARAM.inp.orbital_corr[T] != -1) { - if (L != PARAM.inp.orbital_corr[T]) { - continue; -} - this->cal_slater_Fk(ucell,L, T); + if (L != PARAM.inp.orbital_corr[T]) + { + continue; + } + this->cal_slater_Fk(ucell,L, T); for (int n = 0; n < N; n++) { - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } switch (L) { @@ -170,9 +169,10 @@ void DFTU::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) break; case 3: // f electrons - if (Yukawa) { - this->U_Yukawa[T][L][n] = this->Fk[T][L][n][0]; -} + if (Yukawa) + { + this->U_Yukawa[T][L][n] = this->Fk[T][L][n][0]; + } this->J_Yukawa[T][L][n] = (286.0 * this->Fk[T][L][n][1] + 195.0 * this->Fk[T][L][n][2] + 250.0 * this->Fk[T][L][n][3]) / 6435.0; @@ -192,9 +192,9 @@ void DFTU::cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx) return; } -double DFTU::spherical_Bessel(const int k, const double r, const double lambda) +double Plus_U::spherical_Bessel(const int k, const double r, const double lambda) { - ModuleBase::TITLE("DFTU", "spherical_Bessel"); + ModuleBase::TITLE("Plus_U", "spherical_Bessel"); double val=0.0; double x = r * lambda; @@ -247,9 +247,9 @@ double DFTU::spherical_Bessel(const int k, const double r, const double lambda) return val; } -double DFTU::spherical_Hankel(const int k, const double r, const double lambda) +double Plus_U::spherical_Hankel(const int k, const double r, const double lambda) { - ModuleBase::TITLE("DFTU", "spherical_Bessel"); + ModuleBase::TITLE("Plus_U", "spherical_Bessel"); double val=0.0; double x = r * lambda; @@ -305,6 +305,4 @@ double DFTU::spherical_Hankel(const int k, const double r, const double lambda) return val; } -} // namespace ModuleDFTU - -#endif \ No newline at end of file +#endif diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp index c35aa51bef..327076fab5 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp @@ -19,12 +19,12 @@ hamilt::DFTU>::DFTU(HS_Matrix_K* hsk_in, const Grid_Driver* GridD_in, const TwoCenterIntegrator* intor, const std::vector& orb_cutoff, - ModuleDFTU::DFTU* dftu_in) + 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 = dftu_in; + this->dftu = p_dftu; #ifdef __DEBUG assert(this->ucell != nullptr); #endif @@ -186,7 +186,7 @@ void hamilt::DFTU>::contributeHR() // will update this->dftu->locale and this->dftu->EU if (this->current_spin == 0) { - this->dftu->EU = 0.0; + this->dftu->set_energy(0.0); } } ModuleBase::timer::tick("DFTU", "contributeHR"); @@ -272,7 +272,12 @@ void hamilt::DFTU>::contributeHR() ModuleBase::timer::tick("DFTU", "cal_vu"); const double u_value = this->dftu->U[T0]; std::vector VU_tmp(occ.size()); - this->cal_v_of_u(occ, tlp1, u_value, VU_tmp.data(), this->dftu->EU); + + // mohan add 2025-11-08 + 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); + // transfer occ from pauli matrix format to normal format std::vector VU(occ.size()); this->transfer_vu(VU_tmp, VU); @@ -311,7 +316,7 @@ void hamilt::DFTU>::contributeHR() // energy correction for NSPIN=1 if (this->nspin == 1) { - this->dftu->EU *= 2.0; + this->dftu->set_double_energy(); } // for readin onsite_dm, set initialed_locale to false to avoid using readin locale in next iteration if (this->current_spin == this->nspin - 1 || this->nspin == 4) diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.h b/source/source_lcao/module_operator_lcao/dftu_lcao.h index 4ef1890f76..d24c03bd8d 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.h +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.h @@ -33,7 +33,7 @@ class DFTU> : public OperatorLCAO const Grid_Driver* gridD_in, const TwoCenterIntegrator* intor, const std::vector& orb_cutoff, - ModuleDFTU::DFTU* dftu_in); + Plus_U* p_dftu); ~DFTU>(); /** @@ -51,7 +51,7 @@ class DFTU> : public OperatorLCAO private: const UnitCell* ucell = nullptr; - ModuleDFTU::DFTU* dftu = nullptr; + Plus_U* dftu = nullptr; hamilt::HContainer* HR = nullptr; diff --git a/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp index 9fc4418bb2..c3ed6bb51e 100644 --- a/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/op_dftu_lcao.cpp @@ -27,7 +27,9 @@ void OperatorDFTU>::contributeHk(int ik) ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_real(ik, &eff_pot[0], isk, this->hsk->get_sk()); + + this->dftu->cal_eff_pot_mat_real(ik, &eff_pot[0], isk, this->hsk->get_sk()); + double* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) @@ -43,9 +45,12 @@ void OperatorDFTU, double>>::contributeHk(int { ModuleBase::TITLE("OperatorDFTU", "contributeHk"); ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); + // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + + this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + std::complex* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) @@ -63,7 +68,8 @@ void OperatorDFTU, std::complex>>::con ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + + this->dftu->cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); std::complex* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) @@ -74,4 +80,4 @@ void OperatorDFTU, std::complex>>::con ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); } -} \ No newline at end of file +} diff --git a/source/source_lcao/module_operator_lcao/op_dftu_lcao.h b/source/source_lcao/module_operator_lcao/op_dftu_lcao.h index e15f4ae9dd..a8b275f045 100644 --- a/source/source_lcao/module_operator_lcao/op_dftu_lcao.h +++ b/source/source_lcao/module_operator_lcao/op_dftu_lcao.h @@ -1,7 +1,9 @@ #ifndef OPDFTULCAO_H #define OPDFTULCAO_H + #include "source_base/timer.h" #include "operator_lcao.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace hamilt { @@ -22,11 +24,13 @@ class OperatorDFTU> : public OperatorLCAO public: OperatorDFTU>(HS_Matrix_K* hsk_in, const std::vector>& kvec_d_in, - hamilt::HContainer* hR_in, - const std::vector& isk_in) + hamilt::HContainer* hR_in, + Plus_U* dftu_in, // mohan add 2025-11-05 + const std::vector& isk_in) : isk(isk_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_dftu; + this->dftu = dftu_in; // mohan add 2025-11-07 } virtual void contributeHR() override; @@ -35,9 +39,11 @@ class OperatorDFTU> : public OperatorLCAO private: + Plus_U *dftu; // mohan add 20251107 + bool HR_fixed_done = false; const std::vector& isk; }; } // namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/source_lcao/module_operator_lcao/test/test_dftu.cpp b/source/source_lcao/module_operator_lcao/test/test_dftu.cpp index 12be312278..31adb426ad 100644 --- a/source/source_lcao/module_operator_lcao/test/test_dftu.cpp +++ b/source/source_lcao/module_operator_lcao/test/test_dftu.cpp @@ -7,23 +7,27 @@ #undef private #include "../dftu_lcao.h" #include "source_lcao/module_dftu/dftu.h" -ModuleDFTU::DFTU::DFTU(){}; -ModuleDFTU::DFTU::~DFTU(){}; -namespace GlobalC -{ -ModuleDFTU::DFTU dftu; -} + +Plus_U::Plus_U(){}; +Plus_U::~Plus_U(){}; + +Plus_U dftu; +double Plus_U::energy_u = 0.0; +std::vector Plus_U::U = {}; // U (Hubbard parameter U) +std::vector Plus_U::orbital_corr = {}; + const hamilt::HContainer* tmp_DMR; -const hamilt::HContainer* ModuleDFTU::DFTU::get_dmr(int ispin) const + +const hamilt::HContainer* Plus_U::get_dmr(int ispin) const { return tmp_DMR; } //--------------------------------------- -// Unit test of DFTU class -// DFTU is a derivative class of Operator, it is used to calculate the kinetic matrix +// Unit test of Plus_U class +// Plus_U is a derivative class of Operator, it is used to calculate the kinetic matrix // It use HContainer to store the real space HR matrix -// In this test, we test the correctness and time consuming of 3 functions in DFTU class +// In this test, we test the correctness and time consuming of 3 functions in Plus_U class // - initialize_HR() called in constructor // - contributeHR() // - contributeHk() @@ -35,6 +39,7 @@ const hamilt::HContainer* ModuleDFTU::DFTU::get_dmr(int ispin) const // modify test_size to test different size of unitcell int test_size = 10; int test_nw = 10; // please larger than 5 + class DFTUTest : public ::testing::Test { protected: @@ -82,20 +87,20 @@ class DFTUTest : public ::testing::Test tmp_DMR = DMR; // setting of DFTU - GlobalC::dftu.locale.resize(test_size); + dftu.locale.resize(test_size); for (int iat = 0; iat < test_size; iat++) { - GlobalC::dftu.locale[iat].resize(3); + dftu.locale[iat].resize(3); for (int l = 0; l < 3; l++) { - GlobalC::dftu.locale[iat][l].resize(1); - GlobalC::dftu.locale[iat][l][0].resize(2); - GlobalC::dftu.locale[iat][l][0][0].create(2 * l + 1, 2 * l + 1); - GlobalC::dftu.locale[iat][l][0][1].create(2 * l + 1, 2 * l + 1); + dftu.locale[iat][l].resize(1); + dftu.locale[iat][l][0].resize(2); + dftu.locale[iat][l][0][0].create(2 * l + 1, 2 * l + 1); + dftu.locale[iat][l][0][1].create(2 * l + 1, 2 * l + 1); } } - GlobalC::dftu.U = {U_test}; - GlobalC::dftu.orbital_corr = {orbital_c_test}; + Plus_U::U = {U_test}; + Plus_U::orbital_corr = {orbital_c_test}; PARAM.input.onsite_radius = 1.0; } @@ -155,7 +160,7 @@ TEST_F(DFTUTest, constructHRd2d) } std::chrono::high_resolution_clock::time_point start_time = std::chrono::high_resolution_clock::now(); hamilt::DFTU> - op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &GlobalC::dftu); + op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &dftu); std::chrono::high_resolution_clock::time_point end_time = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed_time = std::chrono::duration_cast>(end_time - start_time); @@ -169,7 +174,7 @@ TEST_F(DFTUTest, constructHRd2d) { for (int icc = 0; icc < 25; icc++) { - EXPECT_NEAR(GlobalC::dftu.locale[iat][2][0][0].c[icc], 0.5, 1e-10); + EXPECT_NEAR(dftu.locale[iat][2][0][0].c[icc], 0.5, 1e-10); } } // check the value of HR @@ -220,14 +225,14 @@ TEST_F(DFTUTest, constructHRd2cd) HR->get_wrapper()[i] = 0.0; } hamilt::DFTU, double>> - op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &GlobalC::dftu); + op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &dftu); op.contributeHR(); // check the occupations of dftu for spin-up for (int iat = 0; iat < test_size; iat++) { for (int icc = 0; icc < 25; icc++) { - EXPECT_NEAR(GlobalC::dftu.locale[iat][2][0][0].c[icc], 0.5, 1e-10); + EXPECT_NEAR(dftu.locale[iat][2][0][0].c[icc], 0.5, 1e-10); } } // check the value of HR @@ -260,7 +265,7 @@ TEST_F(DFTUTest, constructHRd2cd) { for (int icc = 0; icc < 25; icc++) { - EXPECT_NEAR(GlobalC::dftu.locale[iat][2][0][1].c[icc], 0.5, 1e-10); + EXPECT_NEAR(dftu.locale[iat][2][0][1].c[icc], 0.5, 1e-10); } } } diff --git a/source/source_lcao/spar_hsr.cpp b/source/source_lcao/spar_hsr.cpp index 4fc970749a..b0c1ce2139 100644 --- a/source/source_lcao/spar_hsr.cpp +++ b/source/source_lcao/spar_hsr.cpp @@ -69,19 +69,20 @@ void sparse_format::sync_all_R_coor(std::set>& all_R_co template void sparse_format::cal_HSR(const UnitCell& ucell, - const Parallel_Orbitals& pv, - LCAO_HS_Arrays& HS_Arrays, - const Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt* p_ham + Plus_U &dftu, // mohan add 20251107 + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + const Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int (&nmp)[3], + hamilt::Hamilt* p_ham #ifdef __EXX - , - const std::vector>>>* Hexxd, - const std::vector>>>>* Hexxc + , + const std::vector>>>* Hexxd, + const std::vector>>>>* Hexxc #endif -) + ) { ModuleBase::TITLE("sparse_format", "cal_HSR"); @@ -143,17 +144,15 @@ void sparse_format::cal_HSR(const UnitCell& ucell, { if (nspin == 1 || nspin == 2) { - cal_HR_dftu(pv, HS_Arrays.all_R_coor, HS_Arrays.SR_sparse, HS_Arrays.HR_sparse, current_spin, sparse_thr); + cal_HR_dftu(dftu, pv, HS_Arrays.all_R_coor, + HS_Arrays.SR_sparse, HS_Arrays.HR_sparse, current_spin, sparse_thr); } else if (nspin == 4) { - cal_HR_dftu_soc(pv, - HS_Arrays.all_R_coor, - HS_Arrays.SR_soc_sparse, - HS_Arrays.HR_soc_sparse, - current_spin, - sparse_thr); - } + cal_HR_dftu_soc(dftu, pv, HS_Arrays.all_R_coor, + HS_Arrays.SR_soc_sparse, HS_Arrays.HR_soc_sparse, + current_spin, sparse_thr); + } else { ModuleBase::WARNING_QUIT("cal_HSR", "check the value of nspin."); @@ -374,47 +373,53 @@ void sparse_format::destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays) } template void sparse_format::cal_HSR( - const UnitCell& ucell, - const Parallel_Orbitals& pv, - LCAO_HS_Arrays& HS_Arrays, - const Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt* p_ham + const UnitCell& ucell, + Plus_U &dftu, // mohan add 20251107 + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + const Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int (&nmp)[3], + hamilt::Hamilt* p_ham #ifdef __EXX - , - const std::vector>>>* Hexxd, - const std::vector>>>>* Hexxc + , + const std::vector>>>* Hexxd, + const std::vector>>>>* Hexxc #endif -); + ); + template void sparse_format::cal_HSR>( - const UnitCell& ucell, - const Parallel_Orbitals& pv, - LCAO_HS_Arrays& HS_Arrays, - const Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt>* p_ham + const UnitCell& ucell, + Plus_U &dftu, // mohan add 20251107 + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + const Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int (&nmp)[3], + hamilt::Hamilt>* p_ham #ifdef __EXX - , - const std::vector>>>* Hexxd, - const std::vector>>>>* Hexxc + , + const std::vector>>>* Hexxd, + const std::vector>>>>* Hexxc #endif -); + ); + template void sparse_format::cal_HContainer( const Parallel_Orbitals& pv, const double& sparse_thr, const hamilt::HContainer& hR, std::map, std::map>>& target); + template void sparse_format::cal_HContainer>( const Parallel_Orbitals& pv, const double& sparse_thr, const hamilt::HContainer& hR, std::map, std::map>>>& target); + template void sparse_format::cal_HContainer>( const Parallel_Orbitals& pv, const double& sparse_thr, const hamilt::HContainer>& hR, - std::map, std::map>>>& target); \ No newline at end of file + std::map, std::map>>>& target); diff --git a/source/source_lcao/spar_hsr.h b/source/source_lcao/spar_hsr.h index f622317e7a..b9c1334914 100644 --- a/source/source_lcao/spar_hsr.h +++ b/source/source_lcao/spar_hsr.h @@ -3,6 +3,7 @@ #include "source_lcao/LCAO_HS_arrays.hpp" #include "source_lcao/hamilt_lcao.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace sparse_format { @@ -39,19 +40,20 @@ std::set> get_R_range(const hamilt::HContainer& hR) using TAC = std::pair>; template void cal_HSR(const UnitCell& ucell, - const Parallel_Orbitals& pv, - LCAO_HS_Arrays& HS_Arrays, - const Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt* p_ham + Plus_U &dftu, // mohan add 2025-11-07 + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + const Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int (&nmp)[3], + hamilt::Hamilt* p_ham #ifdef __EXX - , - const std::vector>>>* Hexxd = nullptr, - const std::vector>>>>* Hexxc = nullptr + , + const std::vector>>>* Hexxd = nullptr, + const std::vector>>>>* Hexxc = nullptr #endif -); + ); template void cal_HContainer(const Parallel_Orbitals& pv, @@ -65,4 +67,4 @@ void destroy_HS_R_sparse(LCAO_HS_Arrays& HS_Arrays); } // namespace sparse_format -#endif \ No newline at end of file +#endif diff --git a/source/source_lcao/spar_u.cpp b/source/source_lcao/spar_u.cpp index af5f4982ac..6a7d48709a 100644 --- a/source/source_lcao/spar_u.cpp +++ b/source/source_lcao/spar_u.cpp @@ -3,9 +3,9 @@ #include "source_io/module_parameter/parameter.h" #include "source_pw/module_pwdft/global.h" #include "source_base/timer.h" -#include "source_lcao/module_dftu/dftu.h" void sparse_format::cal_HR_dftu( + Plus_U &dftu, // mohan add 2025-11-07 const Parallel_Orbitals &pv, std::set> &all_R_coor, std::map, std::map>> &SR_sparse, @@ -74,7 +74,7 @@ void sparse_format::cal_HR_dftu( } } - GlobalC::dftu.cal_eff_pot_mat_R_double(current_spin, SR_tmp, HR_tmp); + dftu.cal_eff_pot_mat_R_double(current_spin, SR_tmp, HR_tmp); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { @@ -128,6 +128,7 @@ void sparse_format::cal_HR_dftu( void sparse_format::cal_HR_dftu_soc( + Plus_U &dftu, // mohan add 2025-11-07 const Parallel_Orbitals &pv, std::set> &all_R_coor, std::map, std::map>>> &SR_soc_sparse, @@ -194,7 +195,7 @@ void sparse_format::cal_HR_dftu_soc( } } - GlobalC::dftu.cal_eff_pot_mat_R_complex_double(current_spin, SR_soc_tmp, HR_soc_tmp); + dftu.cal_eff_pot_mat_R_complex_double(current_spin, SR_soc_tmp, HR_soc_tmp); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { diff --git a/source/source_lcao/spar_u.h b/source/source_lcao/spar_u.h index ced9781d79..972033b47d 100644 --- a/source/source_lcao/spar_u.h +++ b/source/source_lcao/spar_u.h @@ -5,12 +5,14 @@ #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_pw/module_pwdft/global.h" #include "source_lcao/hamilt_lcao.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 namespace sparse_format { void cal_HR_dftu( + Plus_U &dftu, // mohan add 2025-11-07 const Parallel_Orbitals &pv, std::set> &all_R_coor, std::map, std::map>> &SR_sparse, @@ -19,6 +21,7 @@ namespace sparse_format const double &sparse_thr); void cal_HR_dftu_soc( + Plus_U &dftu, // mohan add 2025-11-07 const Parallel_Orbitals &pv, std::set> &all_R_coor, std::map, std::map>>> &SR_soc_sparse, diff --git a/source/source_pw/module_pwdft/forces.cpp b/source/source_pw/module_pwdft/forces.cpp index 113217657f..819f0cdf23 100644 --- a/source/source_pw/module_pwdft/forces.cpp +++ b/source/source_pw/module_pwdft/forces.cpp @@ -29,6 +29,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 pseudopot_cell_vl* locpp, const pseudopot_cell_vnl* p_nlpp, K_Vectors* pkv, @@ -70,9 +71,10 @@ void Forces::cal_force(UnitCell& ucell, } // DFT+U and DeltaSpin + // here maybe a bug when OFDFT calls +U, mohan add 20251107 if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) { - this->cal_force_onsite(forceonsite, wg, wfc_basis, ucell, psi_in); + this->cal_force_onsite(forceonsite, wg, wfc_basis, ucell, *p_dftu, psi_in); } } diff --git a/source/source_pw/module_pwdft/forces.h b/source/source_pw/module_pwdft/forces.h index e67cac0f24..61466fcdd6 100644 --- a/source/source_pw/module_pwdft/forces.h +++ b/source/source_pw/module_pwdft/forces.h @@ -14,6 +14,7 @@ #include "source_base/kernels/math_kernel_op.h" #include "source_psi/psi.h" #include "structure_factor.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 template class Forces @@ -40,8 +41,9 @@ class Forces const ModulePW::PW_Basis* const rho_basis, ModuleSymmetry::Symmetry* p_symm, Structure_Factor* p_sf, - surchem& solvent, - const pseudopot_cell_vl* locpp, + surchem& solvent, + const Plus_U *p_dftu, //mohan add 2025-11-06 + const pseudopot_cell_vl* locpp, const pseudopot_cell_vnl* nlpp = nullptr, K_Vectors* pkv = nullptr, ModulePW::PW_Basis_K* psi_basis = nullptr, @@ -94,7 +96,9 @@ class Forces const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, const UnitCell& ucell_in, - const psi::Psi , Device>* psi_in = nullptr); + const Plus_U &dftu, // mohan add 2025-11-06 + const psi::Psi , Device>* psi_in = nullptr); + void cal_force_scc(ModuleBase::matrix& forcescc, const ModulePW::PW_Basis* const rho_basis, const ModuleBase::matrix& v_current, diff --git a/source/source_pw/module_pwdft/forces_onsite.cpp b/source/source_pw/module_pwdft/forces_onsite.cpp index e0650ccc36..9d93e9ef54 100644 --- a/source/source_pw/module_pwdft/forces_onsite.cpp +++ b/source/source_pw/module_pwdft/forces_onsite.cpp @@ -11,8 +11,9 @@ template 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 psi::Psi , Device>* psi_in) + const UnitCell& ucell_in, + const Plus_U &dftu, // mohan add 2025-11-06 + const psi::Psi , Device>* psi_in) { ModuleBase::TITLE("Forces", "cal_force_onsite"); if(psi_in == nullptr || wfc_basis == nullptr) @@ -53,12 +54,13 @@ void Forces::cal_force_onsite(ModuleBase::matrix& force_onsite, // force for DFT+U if(PARAM.inp.dft_plus_u) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - onsite_p->get_fs_tools()->cal_force_dftu(ik, npm, force, dftu->orbital_corr.data(), dftu->get_eff_pot_pw(0), dftu->get_size_eff_pot_pw(), wg.c); + onsite_p->get_fs_tools()->cal_force_dftu(ik, npm, force, + dftu.orbital_corr.data(), dftu.get_eff_pot_pw(0), dftu.get_size_eff_pot_pw(), wg.c); } if(PARAM.inp.sc_mag_switch) { - spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); + spinconstrain::SpinConstrain>& sc = + spinconstrain::SpinConstrain>::getScInstance(); const std::vector>& lambda = sc.get_sc_lambda(); onsite_p->get_fs_tools()->cal_force_dspin(ik, npm, force, lambda.data(), wg.c); } @@ -76,4 +78,4 @@ void Forces::cal_force_onsite(ModuleBase::matrix& force_onsite, template class Forces; #if ((defined __CUDA) || (defined __ROCM)) template class Forces; -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/hamilt_lcaopw.h b/source/source_pw/module_pwdft/hamilt_lcaopw.h index 63d6c7c5f3..a5b2cedb6e 100644 --- a/source/source_pw/module_pwdft/hamilt_lcaopw.h +++ b/source/source_pw/module_pwdft/hamilt_lcaopw.h @@ -18,7 +18,7 @@ namespace hamilt K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, const UnitCell* ucell) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell){}; + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp, nullptr, ucell){}; #ifdef __EXX HamiltLIP(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, @@ -26,11 +26,12 @@ namespace hamilt pseudopot_cell_vnl* nlpp, const UnitCell* ucell, Exx_Lip& exx_lip_in) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell), exx_lip(exx_lip_in){}; + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp, nullptr, ucell), + exx_lip(exx_lip_in){}; Exx_Lip& exx_lip; #endif }; } // namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/hamilt_pw.cpp b/source/source_pw/module_pwdft/hamilt_pw.cpp index a48f230bdc..9370c50b68 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.cpp +++ b/source/source_pw/module_pwdft/hamilt_pw.cpp @@ -22,6 +22,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 const UnitCell* ucell): ucell(ucell) { this->classname = "HamiltPW"; @@ -120,7 +121,8 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, if(PARAM.inp.sc_mag_switch || PARAM.inp.dft_plus_u) { Operator* onsite_proj - = new OnsiteProj>(isk, ucell, PARAM.inp.sc_mag_switch, (PARAM.inp.dft_plus_u>0)); + = new OnsiteProj>(isk, ucell, p_dftu, + PARAM.inp.sc_mag_switch, (PARAM.inp.dft_plus_u>0)); this->ops->add(onsite_proj); } if (GlobalC::exx_info.info_global.cal_exx) diff --git a/source/source_pw/module_pwdft/hamilt_pw.h b/source/source_pw/module_pwdft/hamilt_pw.h index 34c51c78fd..ca67b6f424 100644 --- a/source/source_pw/module_pwdft/hamilt_pw.h +++ b/source/source_pw/module_pwdft/hamilt_pw.h @@ -9,6 +9,7 @@ #include "source_pw/module_pwdft/VNL_in_pw.h" #include "source_base/kernels/math_kernel_op.h" #include "source_pw/module_pwdft/module_exx_helper/exx_helper.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 namespace hamilt { @@ -21,10 +22,19 @@ class HamiltPW : public Hamilt // return T if T is real type(float, double), // otherwise return the real type of T(complex, std::complex) using Real = typename GetTypeReal::type; + public: - HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp,const UnitCell* ucell); + + HamiltPW(elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* p_kv, + pseudopot_cell_vnl* nlpp, + Plus_U *p_dftu, // mohan add 2025-11-06 + const UnitCell* ucell); + template explicit HamiltPW(const HamiltPW* hamilt); + ~HamiltPW(); // for target K point, update consequence of hPsi() and matrix() @@ -58,4 +68,4 @@ class HamiltPW : public Hamilt } // namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.cpp b/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.cpp index 56b6d3756a..249f78e4a7 100644 --- a/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.cpp +++ b/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.cpp @@ -13,9 +13,10 @@ namespace hamilt { template OnsiteProj>::OnsiteProj(const int* isk_in, - const UnitCell* ucell_in, - const bool cal_delta_spin, - const bool cal_dftu) + const UnitCell* ucell_in, + Plus_U *p_dftu, // mohan add 2025-11-06 + const bool cal_delta_spin, + const bool cal_dftu) { this->classname = "OnsiteProj"; this->cal_type = calculation_type::pw_onsite; @@ -23,6 +24,7 @@ OnsiteProj>::OnsiteProj(const int* isk_in, this->ucell = ucell_in; this->has_delta_spin = cal_delta_spin; this->has_dftu = cal_dftu; + this->dftu = p_dftu; // mohan add 2025-11-08 } template @@ -209,15 +211,18 @@ void OnsiteProj>::cal_ps_delta_spin(const int npol, const } template -void OnsiteProj>::cal_ps_dftu(const int npol, const int m) const +void OnsiteProj>::cal_ps_dftu( + const int npol, + const int m) const { - if(!this->has_dftu) return; + if(!this->has_dftu) + { + return; + } auto* onsite_p = projectors::OnsiteProjector::get_instance(); const std::complex* becp = onsite_p->get_becp(); - auto* dftu = ModuleDFTU::DFTU::get_instance(); - // T *ps = new T[tnp * m]; // ModuleBase::GlobalFunc::ZEROS(ps, m * tnp); if (this->nkb_m < m * tnp) { @@ -247,7 +252,7 @@ void OnsiteProj>::cal_ps_dftu(const int npol, const int m) for(int iat=0;iatucell->nat;iat++) { const int it = this->ucell->iat2it[iat]; - const int target_l = dftu->orbital_corr[it]; + const int target_l = this->dftu->orbital_corr[it]; orb_l_iat0[iat] = target_l; const int nproj = onsite_p->get_nh(iat); if(target_l == -1) @@ -354,30 +359,56 @@ void OnsiteProj>::cal_ps_dftu(const int npol, const int m) } template<> -void OnsiteProj, base_device::DEVICE_CPU>>::add_onsite_proj(std::complex *hpsi_in, const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_CPU>>::add_onsite_proj( + std::complex *hpsi_in, + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_CPU>>::update_becp(const std::complex *psi_in, const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_CPU>>::update_becp( + const std::complex *psi_in, + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_CPU>>::cal_ps_delta_spin(const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_CPU>>::cal_ps_delta_spin( + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_CPU>>::cal_ps_dftu(const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_CPU>>::cal_ps_dftu( + const int npol, + const int m) const {} #if ((defined __CUDA) || (defined __ROCM)) template<> -void OnsiteProj, base_device::DEVICE_GPU>>::add_onsite_proj(std::complex *hpsi_in, const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_GPU>>::add_onsite_proj( + std::complex *hpsi_in, + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_GPU>>::update_becp(const std::complex *psi_in, const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_GPU>>::update_becp( + const std::complex *psi_in, + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_GPU>>::cal_ps_delta_spin(const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_GPU>>::cal_ps_delta_spin( + const int npol, + const int m) const {} + template<> -void OnsiteProj, base_device::DEVICE_GPU>>::cal_ps_dftu(const int npol, const int m) const +void OnsiteProj, base_device::DEVICE_GPU>>::cal_ps_dftu( + const int npol, + const int m) const {} #endif @@ -405,7 +436,6 @@ hamilt::OnsiteProj>::OnsiteProj(const OnsiteProjclassname = "OnsiteProj"; this->cal_type = calculation_type::pw_nonlocal; - // FIXME: } template class OnsiteProj, base_device::DEVICE_CPU>>; @@ -415,4 +445,4 @@ template class OnsiteProj, base_device::DEVICE_C template class OnsiteProj, base_device::DEVICE_GPU>>; template class OnsiteProj, base_device::DEVICE_GPU>>; #endif -} // namespace hamilt \ No newline at end of file +} // namespace hamilt diff --git a/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.h b/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.h index b8606aa087..3eca4d99d4 100644 --- a/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.h +++ b/source/source_pw/module_pwdft/operator_pw/onsite_proj_pw.h @@ -5,6 +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 namespace hamilt { @@ -22,9 +23,11 @@ class OnsiteProj> : public OperatorPW { private: using Real = typename GetTypeReal::type; + public: OnsiteProj(const int* isk_in, const UnitCell* ucell_in, + Plus_U *p_dftu, // mohan add 2025-11-06 const bool cal_delta_spin, const bool cal_dftu); @@ -40,22 +43,27 @@ class OnsiteProj> : public OperatorPW const int npol, const T* tmpsi_in, T* tmhpsi, - const int ngk = 0, - const bool is_first_node = false)const override; + const int ngk = 0, + const bool is_first_node = false)const override; const int *get_isk() const {return this->isk;} const UnitCell *get_ucell() const {return this->ucell;} private: void cal_ps_delta_spin(const int npol, const int m) const; + void cal_ps_dftu(const int npol, const int m) const; + void update_becp(const T* psi_in, const int npol, const int m) const; + void add_onsite_proj(T *hpsi_in, const int npol, const int m) const; const int* isk = nullptr; const UnitCell* ucell = nullptr; + Plus_U *dftu = nullptr; // mohan add 2025-11-06 + mutable int* ip_iat = nullptr; mutable T* lambda_coeff = nullptr; mutable int* orb_l_iat = nullptr; @@ -95,4 +103,4 @@ class OnsiteProj> : public OperatorPW } // namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index fb20b547e8..e4bc159c80 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -16,6 +16,7 @@ void pw::setup_pot(const int istep, const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi* kspw_psi, // electronic wave functions hamilt::Hamilt* p_hamilt, // hamiltonian @@ -116,11 +117,13 @@ void pw::setup_pot(const int istep, //---------------------------------------------------------- //! 6) DFT+U algorithm + // This should not called in before_scf (esolver), it should be + // called in before_all_runners (esolver), which should + // be improved later. Mohan note 2025-11-06 //---------------------------------------------------------- if (PARAM.inp.dft_plus_u) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - dftu->init(ucell, nullptr, kv.get_nks()); + dftu.init(ucell, nullptr, kv.get_nks()); } return; @@ -136,6 +139,7 @@ template void pw::setup_pot, base_device::DEVICE_CPU>( const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions hamilt::Hamilt, base_device::DEVICE_CPU>* p_hamilt, // hamiltonian @@ -154,6 +158,7 @@ template void pw::setup_pot, base_device::DEVICE_CPU>( const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions hamilt::Hamilt, base_device::DEVICE_CPU>* p_hamilt, // hamiltonian @@ -173,6 +178,7 @@ template void pw::setup_pot, base_device::DEVICE_GPU>( const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions hamilt::Hamilt, base_device::DEVICE_GPU>* p_hamilt, // hamiltonian @@ -190,6 +196,7 @@ template void pw::setup_pot, base_device::DEVICE_GPU>( const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions hamilt::Hamilt, base_device::DEVICE_GPU>* p_hamilt, // hamiltonian diff --git a/source/source_pw/module_pwdft/setup_pot.h b/source/source_pw/module_pwdft/setup_pot.h index b89e33c21e..3849ee5087 100644 --- a/source/source_pw/module_pwdft/setup_pot.h +++ b/source/source_pw/module_pwdft/setup_pot.h @@ -8,6 +8,7 @@ #include "source_estate/elecstate.h" #include "source_pw/module_pwdft/VL_in_pw.h" #include "source_hamilt/hamilt.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-06 namespace pw { @@ -22,6 +23,7 @@ void setup_pot(const int istep, const Charge &chr, // charge density pseudopot_cell_vl &locpp, // local pseudopotentials pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + Plus_U &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi* kspw_psi, // electronic wave functions hamilt::Hamilt* p_hamilt, // hamiltonian diff --git a/source/source_pw/module_pwdft/stress_func.h b/source/source_pw/module_pwdft/stress_func.h index 3d4ed38c89..3f0a79b6d3 100644 --- a/source/source_pw/module_pwdft/stress_func.h +++ b/source/source_pw/module_pwdft/stress_func.h @@ -16,6 +16,7 @@ #include "source_pw/module_pwdft/structure_factor.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 //------------------------------------------------------------------- // mohan reconstruction note: 2021-02-07 @@ -179,8 +180,9 @@ class Stress_Func void stress_onsite(ModuleBase::matrix& sigma, const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, - const UnitCell& ucell_in, - const psi::Psi , Device>* psi_in, + const UnitCell& ucell_in, + const Plus_U &dftu, // mohan add 2025-11-06 + const psi::Psi , Device>* psi_in, ModuleSymmetry::Symmetry* p_symm); // nonlocal part in PW basis void get_dvnl1(ModuleBase::ComplexMatrix& vkb, diff --git a/source/source_pw/module_pwdft/stress_func_onsite.cpp b/source/source_pw/module_pwdft/stress_func_onsite.cpp index 1be9831556..bcd37d900f 100644 --- a/source/source_pw/module_pwdft/stress_func_onsite.cpp +++ b/source/source_pw/module_pwdft/stress_func_onsite.cpp @@ -10,7 +10,8 @@ template void Stress_Func::stress_onsite(ModuleBase::matrix& sigma, const ModuleBase::matrix& wg, const ModulePW::PW_Basis_K* wfc_basis, - const UnitCell& ucell_in, + const UnitCell& ucell_in, + const Plus_U &dftu, // mohan add 2025-11-06 const psi::Psi , Device>* psi_in, ModuleSymmetry::Symmetry* p_symm) { @@ -56,13 +57,12 @@ void Stress_Func::stress_onsite(ModuleBase::matrix& sigma, onsite_p->get_fs_tools()->cal_dbecp_s(ik, npm, ipol, jpol); if(PARAM.inp.dft_plus_u) { - auto* dftu = ModuleDFTU::DFTU::get_instance(); onsite_p->get_fs_tools()->cal_stress_dftu(ik, npm, stress_device_tmp, - dftu->orbital_corr.data(), - dftu->get_eff_pot_pw(0), - dftu->get_size_eff_pot_pw(), + dftu.orbital_corr.data(), + dftu.get_eff_pot_pw(0), + dftu.get_size_eff_pot_pw(), wg.c); } if(PARAM.inp.sc_mag_switch) diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index fb5bfd026b..ee069ad898 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -8,6 +8,7 @@ template void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, UnitCell& ucell, + Plus_U &dftu, // mhan add 2025-11-07 const pseudopot_cell_vl& locpp, const pseudopot_cell_vnl& nlpp, ModulePW::PW_Basis* rho_basis, @@ -119,7 +120,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, // DFT+U and DeltaSpin stress if (PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) { - this->stress_onsite(sigmaonsite, this->pelec->wg, wfc_basis, ucell, d_psi_in, p_symm); + this->stress_onsite(sigmaonsite, this->pelec->wg, wfc_basis, ucell, dftu, d_psi_in, p_symm); } // EXX PW stress diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index a9ac9decb5..0fb31fecaa 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -4,6 +4,7 @@ #include "source_estate/elecstate.h" #include "source_pw/module_pwdft/VL_in_pw.h" #include "stress_func.h" +#include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 template class Stress_PW : public Stress_Func @@ -13,15 +14,16 @@ class Stress_PW : public Stress_Func // calculate the stress in PW basis void cal_stress(ModuleBase::matrix& smearing_sigmatot, - UnitCell& ucell, - const pseudopot_cell_vl& locpp, - const pseudopot_cell_vnl& nlpp, - ModulePW::PW_Basis* rho_basis, - ModuleSymmetry::Symmetry* p_symm, - Structure_Factor* p_sf, - K_Vectors* p_kv, - ModulePW::PW_Basis_K* wfc_basis, - const psi::Psi , Device>* d_psi_in = nullptr); + UnitCell& ucell, + Plus_U &dftu, // mhan add 2025-11-07 + const pseudopot_cell_vl& locpp, + const pseudopot_cell_vnl& nlpp, + ModulePW::PW_Basis* rho_basis, + ModuleSymmetry::Symmetry* p_symm, + Structure_Factor* p_sf, + K_Vectors* p_kv, + ModulePW::PW_Basis_K* wfc_basis, + const psi::Psi , Device>* d_psi_in = nullptr); protected: // call the vdw stress diff --git a/source/source_pw/module_stodft/hamilt_sdft_pw.cpp b/source/source_pw/module_stodft/hamilt_sdft_pw.cpp index fac4e7f867..b1e752c50e 100644 --- a/source/source_pw/module_stodft/hamilt_sdft_pw.cpp +++ b/source/source_pw/module_stodft/hamilt_sdft_pw.cpp @@ -10,11 +10,11 @@ HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, - const UnitCell* ucell, - const int& npol, + const UnitCell* ucell, + const int& npol, Real* emin_in, Real* emax_in) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell), ngk(p_kv->ngk) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp, nullptr, ucell), ngk(p_kv->ngk) { this->classname = "HamiltSdftPW"; this->npwk_max = wfc_basis->npwk_max; @@ -69,4 +69,4 @@ template class HamiltSdftPW, base_device::DEVICE_GPU>; template class HamiltSdftPW, base_device::DEVICE_GPU>; #endif -} // namespace hamilt \ No newline at end of file +} // namespace hamilt diff --git a/source/source_pw/module_stodft/hamilt_sdft_pw.h b/source/source_pw/module_stodft/hamilt_sdft_pw.h index 39bc2bc13e..282ebec424 100644 --- a/source/source_pw/module_stodft/hamilt_sdft_pw.h +++ b/source/source_pw/module_stodft/hamilt_sdft_pw.h @@ -26,7 +26,7 @@ class HamiltSdftPW : public HamiltPW K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, const UnitCell* ucell, - const int& npol, + const int& npol, Real* emin_in, Real* emax_in); /** @@ -64,4 +64,4 @@ class HamiltSdftPW : public HamiltPW } // namespace hamilt -#endif \ No newline at end of file +#endif 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 24ff21ffae..a1937caf05 100644 --- a/source/source_pw/module_stodft/test/test_hamilt_sto.cpp +++ b/source/source_pw/module_stodft/test/test_hamilt_sto.cpp @@ -9,7 +9,14 @@ void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase: void elecstate::Potential::cal_fixed_v(double*){} template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl*,const UnitCell*){} +hamilt::HamiltPW::HamiltPW( + elecstate::Potential* pot_in, + ModulePW::PW_Basis_K* wfc_basis, + K_Vectors* p_kv, + pseudopot_cell_vnl* ppcell, + Plus_U* p_dftu, // mohan add 20251108 + const UnitCell* ucell){} + template hamilt::HamiltPW::~HamiltPW(){ delete this->ops; 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 a0654e1e7f..79d4cdcf22 100644 --- a/source/source_pw/module_stodft/test/test_sto_tool.cpp +++ b/source/source_pw/module_stodft/test/test_sto_tool.cpp @@ -8,7 +8,13 @@ ***********************************************/ template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv,pseudopot_cell_vnl*,const UnitCell*){} +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 + const UnitCell*){} + template hamilt::HamiltPW::~HamiltPW(){}; template @@ -25,9 +31,10 @@ hamilt::HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, const int& npol, Real* emin_in, Real* emax_in) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp, ucell), ngk(p_kv->ngk) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp, nullptr, ucell), ngk(p_kv->ngk) { } + template void hamilt::HamiltSdftPW::hPsi_norm(const T* psi_in, T* hpsi, const int& nbands){} @@ -35,6 +42,7 @@ template class hamilt::HamiltPW, base_device::DEVICE_CPU>; template class hamilt::HamiltSdftPW, base_device::DEVICE_CPU>; template class hamilt::HamiltPW, base_device::DEVICE_CPU>; template class hamilt::HamiltSdftPW, base_device::DEVICE_CPU>; + #if ((defined __CUDA) || (defined __ROCM)) template class hamilt::HamiltPW, base_device::DEVICE_GPU>; template class hamilt::HamiltSdftPW, base_device::DEVICE_GPU>;