Skip to content

refactor(ofdft): reduce solver global state dependencies #8157

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

Closed
zzlinpku wants to merge 1 commit into from
Closed

refactor(ofdft): reduce solver global state dependencies #8157

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
49 changes: 30 additions & 19 deletions source/source_esolver/esolver_of.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,14 @@ void ESolver_OF::before_all_runners(BaseCell& basecell, const Input_para& inp)

ESolver_FP::before_all_runners(ucell, inp);

// Snapshot the global settings needed by OFDFT at the solver boundary.
this->of_runtime_config_.global_readin_dir = PARAM.globalv.global_readin_dir;

@ABACUS-PR-BOT ABACUS-PR-BOT Oct 17, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Could you add a focused regression test for the captured settings and an injected std::ostringstream warning sink? MODULE_PW_OfDFT only compiles KEDF code, and the calculation cases do not distinguish the supplied sink from GlobalV::ofs_warning or a retained snapshot from later global reads. The existing integration coverage is useful; this additional coverage can follow after merge.

this->of_runtime_config_.global_out_dir = PARAM.globalv.global_out_dir;
this->of_runtime_config_.domag = PARAM.globalv.domag;
this->of_runtime_config_.domag_z = PARAM.globalv.domag_z;
this->of_runtime_config_.npol = PARAM.globalv.npol;
this->of_runtime_config_.gamma_only_pw = PARAM.globalv.gamma_only_pw;

// save necessary parameters
this->of_kinetic_ = inp.of_kinetic;
this->of_method_ = inp.of_method;
Expand All @@ -73,21 +81,21 @@ void ESolver_OF::before_all_runners(BaseCell& basecell, const Input_para& inp)
init_rho_cfg.init_chg = inp.init_chg;
init_rho_cfg.suffix = inp.suffix;
init_rho_cfg.esolver_type = inp.esolver_type;
init_rho_cfg.global_readin_dir = PARAM.globalv.global_readin_dir;
init_rho_cfg.global_readin_dir = this->of_runtime_config_.global_readin_dir;
init_rho_cfg.nelec = inp.nelec;
init_rho_cfg.nbands = inp.nbands;
init_rho_cfg.test_charge = inp.test_charge;
init_rho_cfg.domag = PARAM.globalv.domag;
init_rho_cfg.domag_z = PARAM.globalv.domag_z;
init_rho_cfg.npol = PARAM.globalv.npol;
init_rho_cfg.domag = this->of_runtime_config_.domag;
init_rho_cfg.domag_z = this->of_runtime_config_.domag_z;
init_rho_cfg.npol = this->of_runtime_config_.npol;
init_rho_cfg.meta_gga = XC_Functional::get_ked_flag();
this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, nullptr, init_rho_cfg);
module_charge::check_rho(this->chr.rho, this->chr.nspin, this->chr.rhopw->nrxx, ucell.omega,
this->chr.rhopw->nxyz, inp.nelec); // check the rho

// initialize local pseudopotential
this->locpp.init_vloc(ucell,pw_rho);
this->locpp.print_vloc(ucell, pw_rho, inp.out_element_info, PARAM.globalv.global_out_dir);
this->locpp.print_vloc(ucell, pw_rho, inp.out_element_info, this->of_runtime_config_.global_out_dir);
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");


Expand All @@ -96,8 +104,8 @@ void ESolver_OF::before_all_runners(BaseCell& basecell, const Input_para& inp)

// calculate the total local pseudopotential in real space
const int istep=0;
elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep,
PARAM.globalv.global_out_dir, *this->inp_, this->pelec);
elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep,
this->of_runtime_config_.global_out_dir, *this->inp_, this->pelec);

// liuyu move here 2023-10-09
// D in uspp need vloc, thus behind init_scf()
Expand Down Expand Up @@ -136,6 +144,7 @@ void ESolver_OF::runner(BaseCell& basecell, const int istep)
{
basecell.require_kind(BaseCell::Kind::unitcell, __FUNCTION__);
UnitCell& ucell = static_cast<UnitCell&>(basecell);
std::ostream& warning_stream = GlobalV::ofs_warning;

ModuleBase::timer::start("ESolver_OF", "runner");
// get Ewald energy, initial rho and phi if necessary
Expand All @@ -157,13 +166,13 @@ void ESolver_OF::runner(BaseCell& basecell, const int istep)


// check if the job is done
if (this->check_exit(conv_esolver))
if (this->check_exit(conv_esolver, warning_stream))
{
break;
}

// find the optimization direction and step lenghth theta according to the potential
this->optimize(ucell);
this->optimize(ucell, warning_stream);

// update the rho and phi based on the direction and theta
this->update_rho();
Expand Down Expand Up @@ -227,7 +236,8 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
this->resize_spin_arrays();
}

elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, *this->inp_, this->pelec);
elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep,
this->of_runtime_config_.global_out_dir, *this->inp_, this->pelec);

const int nspin = this->inp_->nspin;
if (this->inp_->init_chg == "file")
Expand Down Expand Up @@ -323,7 +333,7 @@ void ESolver_OF::update_potential(UnitCell& ucell)
*
* @param ucell
*/
void ESolver_OF::optimize(UnitCell& ucell)
void ESolver_OF::optimize(UnitCell& ucell, std::ostream& warning_stream)
{
// (1) get |d0> with optimization algorithm
this->get_direction(ucell);
Expand All @@ -347,11 +357,11 @@ void ESolver_OF::optimize(UnitCell& ucell)
// (3) make sure that dEdtheta<0 at theta = 0
std::vector<double> dEdtheta(this->inp_->nspin, 0.0); // dE/dtheta of tempPhi

this->check_direction(dEdtheta.data(), ptemp_phi.data(), ucell);
// this->test_direction(dEdtheta, ptemp_phi, ucell);
this->check_direction(dEdtheta.data(), ptemp_phi.data(), ucell, warning_stream);
// this->test_direction(dEdtheta, ptemp_phi, ucell, warning_stream);

// (4) call line search to find the best theta (step length)
this->get_step_length(dEdtheta.data(), ptemp_phi.data(), ucell);
this->get_step_length(dEdtheta.data(), ptemp_phi.data(), ucell, warning_stream);
}

/**
Expand Down Expand Up @@ -379,7 +389,7 @@ void ESolver_OF::update_rho()
*
* @return exit or not
*/
bool ESolver_OF::check_exit(bool& conv_esolver)
bool ESolver_OF::check_exit(bool& conv_esolver, std::ostream& warning_stream)
{
conv_esolver = false;
bool potConv = false;
Expand Down Expand Up @@ -415,7 +425,7 @@ bool ESolver_OF::check_exit(bool& conv_esolver)
// ============ temporary solution of potential convergence ===========
else if (this->of_conv_ == "potential" && potHold)
{
GlobalV::ofs_warning << "ESolver_OF WARNING: "
warning_stream << "ESolver_OF WARNING: "
<< "The convergence of potential has not been reached, but the norm of potential nearly "
"remains unchanged, set of_full_pw = 1 may work."
<< std::endl;
Expand Down Expand Up @@ -540,7 +550,8 @@ void ESolver_OF::cal_force(BaseCell& basecell, ModuleBase::matrix& force)

// 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(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad,
ff.cal_force(this->inp_->nspin, this->of_runtime_config_.domag, this->of_runtime_config_.domag_z,
this->inp_->gga_grad,
this->inp_->dft_plus_u || this->inp_->sc_mag_switch,
ucell, force, this->get_vdw_result(), *pelec, this->pw_rho, &ucell.symm, &sf,
this->solvent, nullptr, &this->locpp);
Expand All @@ -562,8 +573,8 @@ void ESolver_OF::cal_stress(BaseCell& basecell, ModuleBase::matrix& stress)
this->pphi_, this->pw_rho, kinetic_stress_); // kinetic stress

OF_Stress_PW ss(this->pelec, this->pw_rho);
ss.cal_stress(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z,
this->inp_->gga_grad, PARAM.globalv.gamma_only_pw,
ss.cal_stress(this->inp_->nspin, this->of_runtime_config_.domag, this->of_runtime_config_.domag_z,
this->inp_->gga_grad, this->of_runtime_config_.gamma_only_pw,
stress, kinetic_stress_, ucell, this->get_vdw_result(), &ucell.symm, this->locpp, &sf, &kv);
}
} // namespace ModuleESolver
28 changes: 23 additions & 5 deletions source/source_esolver/esolver_of.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define ESOLVER_OF_H

#include <string>
#include <ostream>

@ABACUS-PR-BOT ABACUS-PR-BOT Oct 17, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Please use <iosfwd> here, since this header only declares std::ostream& parameters and does not need the complete stream definition. The output operations are implemented in the .cpp files. This can be handled in a follow-up.

#include <vector>

#include "esolver_fp.h"
Expand All @@ -12,6 +13,22 @@

namespace ModuleESolver
{
/**
* @brief Runtime settings owned by PARAM but needed by the OFDFT solver.
*
* ESolver_OF snapshots these values at its initialization boundary so later
* OFDFT stages do not read PARAM directly.
*/
struct OFDFTRuntimeConfig
{
std::string global_readin_dir;
std::string global_out_dir;
bool domag = false;
bool domag_z = false;
int npol = 1;
bool gamma_only_pw = false;
};

class ESolver_OF : public ESolver_FP
{
public:
Expand Down Expand Up @@ -50,6 +67,7 @@ class ESolver_OF : public ESolver_FP
int max_iter_ = 50; // scf_nmax

// ------------------ parameters from other module --------------
OFDFTRuntimeConfig of_runtime_config_;
double dV_ = 0; // volume of one grid point in real space
std::vector<double> nelec_; // number of electrons with each spin

Expand Down Expand Up @@ -83,9 +101,9 @@ class ESolver_OF : public ESolver_FP
// ==================== main process of OFDFT ======================
void before_opt(const int istep, UnitCell& ucell);
void update_potential(UnitCell& ucell);
void optimize(UnitCell& ucell);
void optimize(UnitCell& ucell, std::ostream& warning_stream);
void update_rho();
bool check_exit(bool& conv_esolver);
bool check_exit(bool& conv_esolver, std::ostream& warning_stream);
void after_opt(const int istep, UnitCell& ucell, const bool conv_esolver);

// ============================ tools ===============================
Expand All @@ -103,8 +121,8 @@ class ESolver_OF : public ESolver_FP

// --------------------- determine the optimization direction -------
void adjust_direction();
void check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
void check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell, std::ostream& warning_stream);
void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell, std::ostream& warning_stream);

// --------------------- interface to blas --------------------------
double inner_product(double* pa, double* pb, int length, double dV = 1) const
Expand All @@ -117,7 +135,7 @@ class ESolver_OF : public ESolver_FP
// ---------------------- interfaces to optimization methods --------
void init_opt();
void get_direction(UnitCell& ucell);
void get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
void get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& ucell, std::ostream& warning_stream);
};
} // namespace ModuleESolver

Expand Down
11 changes: 7 additions & 4 deletions source/source_esolver/esolver_of_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,10 @@ void ESolver_OF::get_direction(UnitCell& ucell)
* @param ptemp_phi
* @param ucell
*/
void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& ucell)
void ESolver_OF::get_step_length(double* dEdtheta,
double** ptemp_phi,
UnitCell& ucell,
std::ostream& warning_stream)
{
double temp_energy = 0.0; // energy of temp_phi and temp_rho
double kinetic_energy = 0.0; // kinetic energy
Expand Down Expand Up @@ -145,7 +148,7 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell&

if (numDC > this->max_dcsrch_)
{
GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING "
warning_stream << "ESolver_OF linesearch: WARNING "
<< "excedd the max iter number." << std::endl;
break;
}
Expand All @@ -156,13 +159,13 @@ void ESolver_OF::get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell&
}
else if (this->task_.compare(0, 2, "WA") == 0) // warning of line search
{
GlobalV::ofs_warning << "ESolver_OF linesearch: WARNING " << this->task_ << std::endl;
warning_stream << "ESolver_OF linesearch: WARNING " << this->task_ << std::endl;
std::cout << this->task_ << std::endl;
break;
}
else if (this->task_.compare(0, 2, "ER") == 0) // ERROR in line search
{
GlobalV::ofs_warning << "ESolver_OF linesearch: ERROR " << this->task_ << std::endl;
warning_stream << "ESolver_OF linesearch: ERROR " << this->task_ << std::endl;
std::cout << this->task_ << std::endl;
break;
}
Expand Down
5 changes: 3 additions & 2 deletions source/source_esolver/esolver_of_tddft.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ void ESolver_OF_TDDFT::runner(BaseCell& basecell, const int istep)
{
basecell.require_kind(BaseCell::Kind::unitcell, __FUNCTION__);
UnitCell& ucell = static_cast<UnitCell&>(basecell);
std::ostream& warning_stream = GlobalV::ofs_warning;

ModuleBase::timer::start("ESolver_OF_TDDFT", "runner");
// get Ewald energy, initial rho and phi if necessary
Expand Down Expand Up @@ -61,13 +62,13 @@ void ESolver_OF_TDDFT::runner(BaseCell& basecell, const int istep)


// check if the job is done
if (this->check_exit(conv_esolver))
if (this->check_exit(conv_esolver, warning_stream))
{
break;
}

// find the optimization direction and step lenghth theta according to the potential
this->optimize(ucell);
this->optimize(ucell, warning_stream);

// update the rho and phi based on the direction and theta
this->update_rho();
Expand Down
16 changes: 11 additions & 5 deletions source/source_esolver/esolver_of_tool.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,7 +314,10 @@ void ESolver_OF::adjust_direction(void)
* @param ptemp_phi
* @param ucell
*/
void ESolver_OF::check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell)
void ESolver_OF::check_direction(double* dEdtheta,
double** ptemp_phi,
UnitCell& ucell,
std::ostream& warning_stream)
{
assert(this->inp_->nspin > 0);
std::vector<double> temp_theta(this->inp_->nspin, 0.0);
Expand All @@ -332,7 +335,7 @@ void ESolver_OF::check_direction(double* dEdtheta, double** ptemp_phi, UnitCell&
}
else if (dEdtheta[is] > 0)
{
GlobalV::ofs_warning << "ESolver_OF: WARNING "
warning_stream << "ESolver_OF: WARNING "
<< "dEdphi > 0, replace direct with steepest descent method." << std::endl;
for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
{
Expand All @@ -347,7 +350,7 @@ void ESolver_OF::check_direction(double* dEdtheta, double** ptemp_phi, UnitCell&
}
else if (dEdtheta[is] > 0)
{
GlobalV::ofs_warning << "ESolver_OF: WARNING "
warning_stream << "ESolver_OF: WARNING "
<< "when use steepest dencent method, "
"dEdphi > 0, so we might get minimum."
<< std::endl;
Expand All @@ -364,7 +367,10 @@ void ESolver_OF::check_direction(double* dEdtheta, double** ptemp_phi, UnitCell&
* @param ptemp_phi
* @param ucell
*/
void ESolver_OF::test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell)
void ESolver_OF::test_direction(double* dEdtheta,
double** ptemp_phi,
UnitCell& ucell,
std::ostream& warning_stream)
{
double temp_energy = 0.;
if (this->iter_ == 0)
Expand Down Expand Up @@ -394,7 +400,7 @@ void ESolver_OF::test_direction(double* dEdtheta, double** ptemp_phi, UnitCell&
this->dV_);
Parallel_Reduce::reduce_all(pseudopot_energy);
temp_energy += kinetic_energy + pseudopot_energy;
GlobalV::ofs_warning << i << " " << dEdtheta[0] << " " << temp_energy << std::endl;
warning_stream << i << " " << dEdtheta[0] << " " << temp_energy << std::endl;
if (this->theta_[0] == 0)
{
std::cout << "dEdtheta " << dEdtheta[0] << std::endl;
Expand Down
Loading