Skip to content

Fix PW DFT+U initialization lifecycle #8203

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mohanchen merged 1 commit into from
Oct 23, 2026
+102 −150
Merged

Fix PW DFT+U initialization lifecycle #8203

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 22 additions & 4 deletions source/source_esolver/esolver_ks_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -101,6 +102,26 @@ void ESolver_KS_PW<T, Device>::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");
Expand Down Expand Up @@ -166,9 +187,7 @@ void ESolver_KS_PW<T, Device>::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,
Expand All @@ -178,7 +197,6 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
this->chr,
this->locpp,
this->ppcell,
*this->dftu_,
this->vsep_cell,
this->stp.template get_psi_t<T, Device>(),
this->p_hamilt,
Expand Down
78 changes: 20 additions & 58 deletions source/source_pw/module_pwdft/dftu_base.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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<std::string> 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<std::string> 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
Expand Down
13 changes: 5 additions & 8 deletions source/source_pw/module_pwdft/dftu_base.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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<OccMatMixer> occ_mixer_;
Expand Down
29 changes: 0 additions & 29 deletions source/source_pw/module_pwdft/setup_pot.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <typename T, typename Device>
Expand All @@ -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<T, Device>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand Down Expand Up @@ -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;
}

Expand All @@ -148,7 +123,6 @@ template void pw::setup_pot<std::complex<float>, 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<std::complex<float>, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand All @@ -168,7 +142,6 @@ template void pw::setup_pot<std::complex<double>, 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<std::complex<double>, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand All @@ -189,7 +162,6 @@ template void pw::setup_pot<std::complex<float>, 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<std::complex<float>, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand All @@ -208,7 +180,6 @@ template void pw::setup_pot<std::complex<double>, 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<std::complex<double>, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand Down
2 changes: 0 additions & 2 deletions source/source_pw/module_pwdft/setup_pot.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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<T, Device>* kspw_psi, // electronic wave functions
hamilt::HamiltBase* p_hamilt, // hamiltonian
Expand Down
Loading
Loading