From 6398a05121681a7976c6bbfdd0cb4710542916c5 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Thu, 22 Oct 2026 11:17:27 +0800 Subject: [PATCH] Fix PW DFT+U initialization lifecycle --- source/source_esolver/esolver_ks_pw.cpp | 26 ++++- source/source_pw/module_pwdft/dftu_base.cpp | 78 ++++--------- source/source_pw/module_pwdft/dftu_base.h | 13 +-- source/source_pw/module_pwdft/setup_pot.cpp | 29 ----- source/source_pw/module_pwdft/setup_pot.h | 2 - .../module_pwdft/test/dftu_base_test.cpp | 104 +++++++++--------- 6 files changed, 102 insertions(+), 150 deletions(-) diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 1a10a37c316..f1e6c5fa47c 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -20,6 +20,7 @@ #include "source_estate/setup_estate_pw.h" // mohan add 20251005 #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type used via general_exx_info_ #include "source_io/module_ctrl/ctrl_output_pw.h" // mohan add 20250927 +#include "source_pw/module_pwdft/dftu_base.h" #include "source_pw/module_pwdft/dftu_base_io.h" // append_ion_step_snapshot #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/deltaspin_pw.h" // mohan add 20250309 @@ -101,6 +102,26 @@ void ESolver_KS_PW::before_all_runners(BaseCell& basecell, const Inpu this->solvent, inp); + if (inp.dft_plus_u) + { + const System_para& sys = PARAM.globalv; + this->dftu_->init_base(ucell, + sys.npol, + inp.nspin, + inp.l_channel, + inp.yukawa_potential, + inp.yukawa_lambda, + sys.global_readin_dir, + sys.global_out_dir, + inp.init_chg, + inp.device, + sys.hubbard_u, + sys.uramping, + inp.init_occ_mat, + inp.mixing_dftu, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + this->stp.before_runner(ucell, this->kv, this->sf, *this->pw_wfc, this->ppcell.lmaxkb, *this->inp_); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); @@ -166,9 +187,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) //! Allocate HamiltPW 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 + //! Setup potentials (local, non-local, spin constraint, DFT-1/2) pw::setup_pot(istep, ucell, this->kv, @@ -178,7 +197,6 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->chr, this->locpp, this->ppcell, - *this->dftu_, this->vsep_cell, this->stp.template get_psi_t(), this->p_hamilt, diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 35e33c45816..b11c1487181 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -66,17 +66,6 @@ void Plus_U_Base::init_base(UnitCell& cell, this->energy_u = 0.0; - // The occupation-matrix file is read only once for the whole run. - // occmat_.init() below reallocates and clears the storage on every - // ionic step, so preserve the current in-memory matrix first and - // restore it after init() when the file has already been loaded. - OccupationMatrix occmat_kept; - const bool restore_occmat = (init_occ_mat != 0 && this->occmat_file_loaded_); - if (restore_occmat) - { - occmat_kept = this->occmat_; - } - this->occmat_.init(cell, l_channel, nspin, npol); this->uterm_mat_index.resize(cell.nat); @@ -141,60 +130,33 @@ void Plus_U_Base::init_base(UnitCell& cell, this->yukawa_.reset(new YukawaScreening()); this->yukawa_->init(cell, l_channel, yukawa_lambda); } - else - { - // Clear any stale object from a previous init_base() call with - // yukawa_potential == true, preserving the old explicit-flag semantics. - this->yukawa_.reset(); - } - if (init_occ_mat != 0) { - if (this->occmat_file_loaded_) + // Try the user-prepared dm_onsite_ini.txt first, then fall back to + // the compact occupation-matrix snapshot written by a previous + // run (occ_mat.txt). Legacy dm_onsite.txt input is not supported. + const std::vector candidates = {"dm_onsite_ini.txt", + "occ_mat.txt"}; + const std::string readin_fn = DFTU_BASE::find_first_existing_file(global_readin_dir, candidates); + if (readin_fn.empty()) { - // The file was read at the first ionic step; keep using the - // in-memory occupation matrix instead of reading the file - // again. For init_occ_mat=2 this matrix never changes; for - // init_occ_mat=1 it carries the result of the previous - // ionic step. - this->occmat_ = occmat_kept; - this->set_occmat_ready(); - this->occmat_.copy_to_save(cell, this->l_channel); - if (this->has_occ_mixer()) - { - // seed the freshly rebuilt mixer with the current matrix - this->occ_mixer().seed_save(this->occmat_); - } + ModuleBase::WARNING_QUIT("Plus_U_Base::init_base", + "init_occ_mat is set but no occupation-matrix file found in " + + global_readin_dir + + ". Tried: dm_onsite_ini.txt, occ_mat.txt"); } - else - { - // Try the user-prepared dm_onsite_ini.txt first, then fall back to - // the compact occupation-matrix snapshot written by a previous - // run (occ_mat.txt). Legacy dm_onsite.txt input is not supported. - const std::vector candidates = {"dm_onsite_ini.txt", - "occ_mat.txt"}; - const std::string readin_fn = DFTU_BASE::find_first_existing_file(global_readin_dir, candidates); - if (readin_fn.empty()) - { - ModuleBase::WARNING_QUIT("Plus_U_Base::init_base", - "init_occ_mat is set but no occupation-matrix file found in " - + global_readin_dir - + ". Tried: dm_onsite_ini.txt, occ_mat.txt"); - } - DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, - readin_fn, nspin, npol, soc_layout); + DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, + readin_fn, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); #endif - this->set_occmat_ready(); - this->occmat_.copy_to_save(cell, this->l_channel); - if (this->has_occ_mixer()) - { - // seed the mixing history with the file-loaded occupation matrix - this->occ_mixer().seed_save(this->occmat_); - } - this->occmat_file_loaded_ = true; + this->set_occmat_ready(); + this->occmat_.copy_to_save(cell, this->l_channel); + if (this->has_occ_mixer()) + { + // seed the mixing history with the file-loaded occupation matrix + this->occ_mixer().seed_save(this->occmat_); } } else diff --git a/source/source_pw/module_pwdft/dftu_base.h b/source/source_pw/module_pwdft/dftu_base.h index af590f53305..3203d60e114 100644 --- a/source/source_pw/module_pwdft/dftu_base.h +++ b/source/source_pw/module_pwdft/dftu_base.h @@ -35,7 +35,11 @@ class Plus_U_Base Plus_U_Base(); ~Plus_U_Base(); - /// allocate relevant data structures (base part, no LCAO types) + /// Initialize the run-lifetime DFT+U state (base part, no LCAO types). + /// + /// Call this exactly once on a freshly constructed object for each solver + /// run, before entering the electronic or ionic loops. A new calculation + /// configuration requires a new Plus_U_Base object. void init_base(UnitCell& cell, const int npol, const int nspin, @@ -152,13 +156,6 @@ class Plus_U_Base // --- Occupation matrices --- OccupationMatrix occmat_; - /// Whether the initial occupation matrix has already been read from - /// file during this run. The file is read exactly once; on later - /// ionic steps the in-memory matrix is preserved (it stays fixed for - /// init_occ_mat=2 and evolves from the previous ionic step for - /// init_occ_mat=1). - bool occmat_file_loaded_ = false; - // Occupation-matrix mixer; constructed only when mixing_dftu != 0. // Owns the flat uom/uom_save buffers and the mixing orchestration. std::unique_ptr occ_mixer_; diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index eec2cf13ad6..52761bb5270 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -4,7 +4,6 @@ #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_proj/onsite_proj.h" #include "source_pw/module_pwdft/vnl_pw.h" -#include "source_pw/module_pwdft/dftu_base.h" // mohan add 2025-11-06 #include "source_pw/module_pwdft/vsep_pw.h" template @@ -17,7 +16,6 @@ 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_Base &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi* kspw_psi, // electronic wave functions hamilt::HamiltBase* p_hamilt, // hamiltonian @@ -112,29 +110,6 @@ void pw::setup_pot(const int istep, pw_wfc); } - //---------------------------------------------------------- - //! 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 (inp.dft_plus_u) - { - dftu.init_base(ucell, - PARAM.globalv.npol, - inp.nspin, inp.l_channel, inp.yukawa_potential, - inp.yukawa_lambda, - PARAM.globalv.global_readin_dir, - PARAM.globalv.global_out_dir, - inp.init_chg, - inp.device, - PARAM.globalv.hubbard_u, - PARAM.globalv.uramping, - inp.init_occ_mat, - inp.mixing_dftu, - DFTU_BASE::SOC_LAYOUT_PAULI); - } - return; } @@ -148,7 +123,6 @@ 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_Base &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::HamiltBase* p_hamilt, // hamiltonian @@ -168,7 +142,6 @@ 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_Base &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::HamiltBase* p_hamilt, // hamiltonian @@ -189,7 +162,6 @@ 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_Base &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::HamiltBase* p_hamilt, // hamiltonian @@ -208,7 +180,6 @@ 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_Base &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::HamiltBase* p_hamilt, // hamiltonian diff --git a/source/source_pw/module_pwdft/setup_pot.h b/source/source_pw/module_pwdft/setup_pot.h index fa604da3561..6842159e261 100644 --- a/source/source_pw/module_pwdft/setup_pot.h +++ b/source/source_pw/module_pwdft/setup_pot.h @@ -7,7 +7,6 @@ #include "source_estate/elecstate.h" #include "source_pw/module_pwdft/vl_pw.h" #include "source_hamilt/hamilt.h" -#include "source_pw/module_pwdft/dftu_base.h" // mohan add 2025-11-06 class pseudopot_cell_vnl; @@ -24,7 +23,6 @@ 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_Base &dftu, // mohan add 2025-11-06 VSep* vsep_cell, // U-1/2 method psi::Psi* kspw_psi, // electronic wave functions hamilt::HamiltBase* p_hamilt, // hamiltonian diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 596da0b7cab..233af229513 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -1,11 +1,8 @@ /********************************************** * Unit tests for Plus_U_Base::init_base. * - * Focus: the Yukawa state must follow the - * yukawa_potential argument on every call, - * including a true -> false re-initialization - * (before_scf() -> setup_pot() may call - * init_base() repeatedly on the same object). + * Focus: initialize base DFT+U state before + * electronic and ionic iterations begin. ***********************************************/ #include "source_pw/module_pwdft/dftu_base.h" @@ -169,29 +166,54 @@ void write_compact_collinear_atom(std::ostream& os, } // namespace -/// After a true -> false re-initialization the Yukawa object must be -/// released so that use_yukawa() reflects the latest argument. -TEST_F(DFTUBaseTest, InitBaseYukawaTrueThenFalseClearsState) +/// Enabling Yukawa screening constructs its run-lifetime state. +TEST_F(DFTUBaseTest, InitBaseYukawaEnabledCreatesObject) { Plus_U_Base dftu; init_dftu(dftu, true); EXPECT_TRUE(dftu.use_yukawa()); +} + +/// Disabling Yukawa screening leaves its optional state empty. +TEST_F(DFTUBaseTest, InitBaseYukawaDisabledLeavesObjectEmpty) +{ + Plus_U_Base dftu; init_dftu(dftu, false); EXPECT_FALSE(dftu.use_yukawa()); } -/// A false -> true re-initialization must create the Yukawa object. -TEST_F(DFTUBaseTest, InitBaseYukawaFalseThenTrueCreatesObject) +/// U-ramping starts from zero and advances toward the configured target. +TEST_F(DFTUBaseTest, InitBaseSeedsURampingState) { Plus_U_Base dftu; + const std::vector l_channel = {2}; + const std::vector hubbard_u = {5.0}; - init_dftu(dftu, false); - EXPECT_FALSE(dftu.use_yukawa()); + dftu.init_base(ucell, + 1, + 2, + l_channel, + false, + 0.5, + "", + "", + "none", + "cpu", + hubbard_u, + 0.5, + 0, + 0, + DFTU_BASE::SOC_LAYOUT_PAULI); - init_dftu(dftu, true); - EXPECT_TRUE(dftu.use_yukawa()); + EXPECT_DOUBLE_EQ(dftu.get_u_current(0), 0.0); + dftu.uramping_update(); + EXPECT_DOUBLE_EQ(dftu.get_u_current(0), 0.5); + EXPECT_FALSE(dftu.u_converged()); + dftu.uramping_update(); + EXPECT_DOUBLE_EQ(dftu.get_u_current(0), 1.0); + EXPECT_FALSE(dftu.u_converged()); } // ===================================================================== @@ -639,11 +661,8 @@ TEST(FindFirstExistingFileTest, ReturnsFirstExistingCandidate) } } -/// For init_occ_mat=2 the occupation-matrix file is read exactly once. -/// Later init_base() calls (one per ionic step in a relax run) must keep -/// the in-memory matrix instead of reading the file again, so pointing -/// the second call at a non-existent readin dir must not abort the run. -TEST_F(DFTUBaseTest, InitBaseReadsOccMatFileOnlyOnce) +/// init_occ_mat=2 loads the fixed occupation matrix during run initialization. +TEST_F(DFTUBaseTest, InitBaseReadsFixedOccMat) { // Use the gtest-managed temporary directory so the test never creates // or deletes directories itself (see AGENTS.md rule 17). @@ -651,8 +670,7 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatFileOnlyOnce) const std::string fn = dir + "occ_mat.txt"; // One Fe atom, L=2, two spin channels of 5x5, filled with distinct - // constants so a re-read would be easy to distinguish from a - // preserved in-memory matrix. + // constants. { std::ofstream ofs(fn); ASSERT_TRUE(ofs.is_open()); @@ -681,35 +699,23 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatFileOnlyOnce) Plus_U_Base dftu; const std::vector l_channel = {2}; const std::vector hubbard_u = {0.0}; - auto call_init = [&](const std::string& readin_dir) - { - dftu.init_base(ucell, - 1, // npol - 2, // nspin - l_channel, - false, // yukawa_potential - 0.5, // yukawa_lambda - readin_dir, // global_readin_dir - "", // global_out_dir - "none", // init_chg - "cpu", // device - hubbard_u, - 0.0, // uramping - 2, // init_occ_mat - 0, // mixing_dftu - DFTU_BASE::SOC_LAYOUT_PAULI); - }; - - // First ionic step: the file is read. - call_init(dir); - ASSERT_TRUE(dftu.is_occmat_ready()); - EXPECT_NEAR(dftu.occmat().get(0, 2, 0, 0, 0), 0.1, 1e-12); - EXPECT_NEAR(dftu.occmat().get(0, 2, 1, 4, 4), 0.2, 1e-12); + dftu.init_base(ucell, + 1, // npol + 2, // nspin + l_channel, + false, // yukawa_potential + 0.5, // yukawa_lambda + dir, // global_readin_dir + "", // global_out_dir + "none", // init_chg + "cpu", // device + hubbard_u, + 0.0, // uramping + 2, // init_occ_mat + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); - // Simulate a later ionic step: the in-memory matrix is preserved and - // the file is not consulted again, so a non-existent readin dir is fine. - call_init(dir + "does_not_exist/"); - EXPECT_TRUE(dftu.is_occmat_ready()); + ASSERT_TRUE(dftu.is_occmat_ready()); EXPECT_NEAR(dftu.occmat().get(0, 2, 0, 0, 0), 0.1, 1e-12); EXPECT_NEAR(dftu.occmat().get(0, 2, 1, 4, 4), 0.2, 1e-12); }