Skip to content

Precommit: Format source_esolver #8189

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

Open
Growl1234 wants to merge 1 commit into
base: develop
Choose a base branch
from
+2,749 −1,849
Open

Precommit: Format source_esolver #8189

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
4 changes: 3 additions & 1 deletion source/source_esolver/esolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,9 @@ class ESolver

//! deal with exx and other calculation than scf/md/relax/cell-relax:
//! such as nscf, get_wf and get_pchg
virtual void others(BaseCell&, const int) {}
virtual void others(BaseCell&, const int)
{
}

//! calculate total energy of a given system
virtual double cal_energy() = 0;
Expand Down
79 changes: 46 additions & 33 deletions source/source_esolver/esolver_dfpt_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@

#include "esolver_dfpt_pw.h"

#include "source_estate/module_charge/charge.h"
#include "source_estate/module_pot/pot_xc_fdm.h"
#include "source_base/macros.h"
#include "source_base/math_polyint.h"
#include "source_base/tool_quit.h"
#include "source_estate/module_charge/charge.h"
#include "source_estate/module_pot/pot_xc_fdm.h"
#include "source_hamilt/module_xc/libxc_abacus.h"
#include "source_hamilt/module_xc/xc_functional.h"
#include "source_io/module_parameter/parameter.h"
Expand All @@ -30,7 +30,8 @@
#include <sstream>
#include <vector>

namespace {
namespace
{

/**
* @brief Re/Im-split finite-difference adapter over elecstate::PotXC_FDM.
Expand All @@ -44,12 +45,12 @@ class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order
{
public:
XC_First_Order_FDM(const int nspin,
const bool domag,
const bool domag_z,
const int gga_grad,
const bool out_elf,
const int test_charge,
ModulePW::PW_Basis* rho_basis,
const bool domag,
const bool domag_z,
const int gga_grad,
const bool out_elf,
const int test_charge,
ModulePW::PW_Basis* rho_basis,
const Charge* chg0,
const UnitCell* ucell)
: ucell_(ucell)
Expand Down Expand Up @@ -130,12 +131,10 @@ class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order
imax = ir;
}
}
std::cout << "XCDBG call=" << call
<< " <drho|dvxc>/|drho|^2=" << (dp_dv / dp_dp)
std::cout << "XCDBG call=" << call << " <drho|dvxc>/|drho|^2=" << (dp_dv / dp_dp)
<< " |dvxc|/|drho|=" << (std::sqrt(dv_dv) / std::sqrt(dp_dp))
<< " ratio@max=" << (dvxc_r[imax].real() / drho_r[imax].real())
<< " drho@max=" << drho_r[imax].real()
<< std::endl;
<< " drho@max=" << drho_r[imax].real() << std::endl;
}
++call;
}
Expand All @@ -158,8 +157,7 @@ class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order
fdm_->cal_v_eff(chg1_, ucell_, veff_1_);
for (int ir = 0; ir < nrxx; ++ir)
{
dvxc_r[ir] += std::complex<double>(0.0, 1.0)
* (v_plus[ir] - veff_1_(0, ir)) / (2.0 * eta);
dvxc_r[ir] += std::complex<double>(0.0, 1.0) * (v_plus[ir] - veff_1_(0, ir)) / (2.0 * eta);
}
}

Expand Down Expand Up @@ -287,8 +285,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
}
if (this->pelec->charge->nrxx != this->pw_rho->nrxx)
{
ModuleBase::WARNING_QUIT("ESolver_DFPT_PW::init_dfpt",
"charge is not on the rho grid (DFPT supports NCPP)");
ModuleBase::WARNING_QUIT("ESolver_DFPT_PW::init_dfpt", "charge is not on the rho grid (DFPT supports NCPP)");
}

// converged effective potential on the shared rho grid (row 0: nspin = 1)
Expand All @@ -312,9 +309,8 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
if (xc_func_type != 1 && xc_func_type != 2)
{
std::stringstream msg;
msg << "DFPT supports LDA/GGA functionals only; the ground-state functional '"
<< ucell.atoms[0].ncpp.xc_func << "' (func_type " << xc_func_type
<< ") needs response channels that are not implemented";
msg << "DFPT supports LDA/GGA functionals only; the ground-state functional '" << ucell.atoms[0].ncpp.xc_func
<< "' (func_type " << xc_func_type << ") needs response channels that are not implemented";
ModuleBase::WARNING_QUIT("ESolver_DFPT_PW::init_dfpt", msg.str());
}

Expand Down Expand Up @@ -343,7 +339,9 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
// factors (e.g. BLYP_LR subtracts its short-range components)
const std::vector<int> xc_func_id = XC_Functional::get_func_id();
const std::vector<std::vector<double>> xc_ext_params
= XC_Functional_Libxc::ext_params_for(xc_func_id, XC_Functional::get_hybrid_alpha(), XC_Functional::get_hse_omega());
= XC_Functional_Libxc::ext_params_for(xc_func_id,
XC_Functional::get_hybrid_alpha(),
XC_Functional::get_hse_omega());
std::vector<double> xc_weights(xc_func_id.size(), 1.0);
for (size_t i = 0; i < xc_func_id.size(); ++i)
{
Expand All @@ -360,8 +358,15 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
"built without libxc: the DFPT XC kernel falls back to the finite-difference "
"amplitude path (LDA-only support; the GGA response at q != 0 lacks the "
"i(q.f_grad)A carrier term)");
xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad,
this->inp_->out_elf[0] > 0, this->inp_->test_charge, this->pw_rho, this->pelec->charge, &ucell);
xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin,
PARAM.globalv.domag,
PARAM.globalv.domag_z,
this->inp_->gga_grad,
this->inp_->out_elf[0] > 0,
this->inp_->test_charge,
this->pw_rho,
this->pelec->charge,
&ucell);
#endif

if (getenv("DFPT_VKB") != nullptr)
Expand All @@ -383,8 +388,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
pert.build_vkb(0, 1, gk, vkb1);
const int nh = ucell.atoms[0].ncpp.nh;
const std::complex<double>* gsvkb = this->ppcell.get_vkb_data<double>();
std::cout << "VKBCHK npw=" << npw << " nh=" << nh
<< " nkb=" << this->ppcell.nkb << std::endl;
std::cout << "VKBCHK npw=" << npw << " nh=" << nh << " nkb=" << this->ppcell.nkb << std::endl;
for (int mu = 0; mu < nh; ++mu)
{
std::complex<double> dot(0.0, 0.0);
Expand All @@ -397,8 +401,8 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
nrm_mine += std::norm(vkb[mu][ig]);
nrm_gs += std::norm(gs);
}
std::cout << "VKBCHK mu=" << mu << " <mine|gs>=" << dot
<< " |mine|^2=" << nrm_mine << " |gs|^2=" << nrm_gs << std::endl;
std::cout << "VKBCHK mu=" << mu << " <mine|gs>=" << dot << " |mine|^2=" << nrm_mine << " |gs|^2=" << nrm_gs
<< std::endl;
}
for (int mu = 0; mu < nh; ++mu)
{
Expand All @@ -412,20 +416,29 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
nrm_mine += std::norm(vkb1[mu][ig]);
nrm_gs += std::norm(gs);
}
std::cout << "VKBCHK a1 mu=" << mu << " <mine|gs>=" << dot
<< " |mine|^2=" << nrm_mine << " |gs|^2=" << nrm_gs << std::endl;
std::cout << "VKBCHK a1 mu=" << mu << " <mine|gs>=" << dot << " |mine|^2=" << nrm_mine
<< " |gs|^2=" << nrm_gs << std::endl;
}
for (int ig = 0; ig < 6; ++ig)
{
std::cout << "VKBEL a1 mu=0 ig=" << ig << " mine=" << vkb1[0][ig]
<< " gs=" << gsvkb[4 * this->pw_wfc->npwk_max + ig]
<< " gcar=" << gk[ig].x << "," << gk[ig].y << "," << gk[ig].z << std::endl;
<< " gs=" << gsvkb[4 * this->pw_wfc->npwk_max + ig] << " gcar=" << gk[ig].x << "," << gk[ig].y
<< "," << gk[ig].z << std::endl;
}
}

this->stp.sync_cpu();
dfpt_->init(ucell, *this->stp.psi_cpu, this->pw_rho, this->pw_wfc, &this->sf, veff_r,
this->pelec->wg, this->pelec->ekb, xc_adapter_, nelec_, ecutwfc_,
dfpt_->init(ucell,
*this->stp.psi_cpu,
this->pw_rho,
this->pw_wfc,
&this->sf,
veff_r,
this->pelec->wg,
this->pelec->ekb,
xc_adapter_,
nelec_,
ecutwfc_,
dft_plus_u_ ? this->dftu_.get() : nullptr);
}

Expand Down
3 changes: 2 additions & 1 deletion source/source_esolver/esolver_dfpt_pw.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@
#include "esolver_ks_pw.h"
#include "source_pw/module_dfpt/dfpt_pw.h"

namespace ModuleDFPT {
namespace ModuleDFPT
{
class XC_First_Order;
}

Expand Down
9 changes: 7 additions & 2 deletions source/source_esolver/esolver_dm2rho.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
#include "source_estate/elecstate_lcao.h"
#include "source_io/module_ml/io_npz.h"
#include "source_io/module_output/cube_io.h"
#include "source_lcao/lcao_domain.h"
#include "source_lcao/hamilt_lcao.h"
#include "source_lcao/lcao_domain.h"
#include "source_lcao/module_operator_lcao/operator_lcao.h"
#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24

Expand Down Expand Up @@ -66,7 +66,12 @@ void ESolver_DM2rho<TK, TR>::runner(BaseCell& basecell, const int istep)

// it's dangerous to design psiToRho function like this, mohan note 20251024
// this->pelec->psiToRho(*this->psi);
LCAO_domain::dm2rho(this->dmat.dm->get_dmr_vec(), this->inp_->nspin, &this->chr, this->inp_->nelec, ucell.omega, false);
LCAO_domain::dm2rho(this->dmat.dm->get_dmr_vec(),
this->inp_->nspin,
&this->chr,
this->inp_->nelec,
ucell.omega,
false);

// Write every density channel (rho[0..3] for nspin=4) so the file set
// matches ctrl_output_fp: the init_chg=file reader reads chgs1..chgs<nspin>
Expand Down
49 changes: 29 additions & 20 deletions source/source_esolver/esolver_double_xc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,9 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(BaseCell& basecell, const Inpu
{
nsk = this->inp_->nspin;
ncol = this->pv.ncol_bands;
if (this->inp_->ks_solver == "genelpa" || this->inp_->ks_solver == "elpa" || this->inp_->ks_solver == "lapack"
|| this->inp_->ks_solver == "pexsi" || this->inp_->ks_solver == "cusolver"
|| this->inp_->ks_solver == "cusolvermp")
if (this->inp_->ks_solver == "genelpa" || this->inp_->ks_solver == "elpa"
|| this->inp_->ks_solver == "lapack" || this->inp_->ks_solver == "pexsi"
|| this->inp_->ks_solver == "cusolver" || this->inp_->ks_solver == "cusolvermp")
{
ncol = this->pv.ncol;
}
Expand Down Expand Up @@ -104,14 +104,16 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(BaseCell& basecell, const Inpu
init_rho_cfg.domag_z = PARAM.globalv.domag_z;
init_rho_cfg.npol = PARAM.globalv.npol;
init_rho_cfg.meta_gga = XC_Functional::get_ked_flag();
this->chr_base.set_rhopw(this->pw_rhod); // mohan add 20251130
this->chr_base.set_rhopw(this->pw_rhod); // mohan add 20251130
const bool kin_den = XC_Functional::get_ked_flag() || (this->inp_->out_elf[0] > 0); // mohan add 20251202
this->chr_base.allocate(this->inp_->nspin, kin_den, XC_Functional::get_ked_flag(),
this->inp_->test_charge);
this->chr_base.allocate(this->inp_->nspin, kin_den, XC_Functional::get_ked_flag(), this->inp_->test_charge);
this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, nullptr, init_rho_cfg);
module_charge::check_rho(this->chr_base.rho, this->chr_base.nspin,
this->chr_base.rhopw->nrxx, ucell.omega,
this->chr_base.rhopw->nxyz, this->inp_->nelec);
module_charge::check_rho(this->chr_base.rho,
this->chr_base.nspin,
this->chr_base.rhopw->nrxx,
ucell.omega,
this->chr_base.rhopw->nxyz,
this->inp_->nelec);

// 11) initialize the potential
if (this->pelec_base->pot == nullptr)
Expand Down Expand Up @@ -158,8 +160,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
}
if (this->p_hamilt_base == nullptr)
{
const bool load_exx_flag = !GlobalC::restart.info_load.restart_exx
&& GlobalC::restart.info_load.load_H;
const bool load_exx_flag = !GlobalC::restart.info_load.restart_exx && GlobalC::restart.info_load.load_H;
this->p_hamilt_base = new hamilt::HamiltLCAO<TK, TR>(ucell,
this->gd,
&this->pv,
Expand Down Expand Up @@ -305,10 +306,10 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
this->inp_->imp_sol,
this->inp_->sc_mag_switch,
this->inp_->dft_plus_u,
this->inp_->assume_isolated); // 2 means Kohn-Sham functional
// std::cout<<"in double_xc------"<<std::endl;
// this->pelec_base->f_en.print_all();
// std::cout<<"in double_xc------"<<std::endl;
this->inp_->assume_isolated); // 2 means Kohn-Sham functional
// std::cout<<"in double_xc------"<<std::endl;
// this->pelec_base->f_en.print_all();
// std::cout<<"in double_xc------"<<std::endl;
// GlobalV::ofs_running << std::setprecision(15) << " etot of base functional (Ry) " << pelec_base->f_en.etot <<
// std::endl;

Expand Down Expand Up @@ -424,10 +425,16 @@ void ESolver_DoubleXC<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&

this->deepks.dpks_out_type = "base"; // for deepks method

const FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h,
this->inp_->sc_mag_switch, this->inp_->device,
PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad,
PARAM.globalv.gamma_only_pw, this->inp_->imp_sol};
const FSCalcConfig fs_cfg{this->inp_->nspin,
this->inp_->nbands,
this->inp_->t_in_h,
this->inp_->sc_mag_switch,
this->inp_->device,
PARAM.globalv.domag,
PARAM.globalv.domag_z,
this->inp_->gga_grad,
PARAM.globalv.gamma_only_pw,
this->inp_->imp_sol};

fsl.getForceStress(ucell,
this->get_vdw_result(),
Expand All @@ -454,7 +461,9 @@ void ESolver_DoubleXC<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
this->exx_nao,
&ucell.symm,
this->exx_info_,
fs_cfg, nullptr, nullptr);
fs_cfg,
nullptr,
nullptr);

// restore to original xc
XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
Expand Down
Loading
Loading