Skip to content

Update ESolver, delete pelec->charge, instead using chr defined in ESolver_fp directly #5963

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
merged 6 commits into from
Mar 7, 2025
Merged

Update ESolver, delete pelec->charge, instead using chr defined in ESolver_fp directly #5963

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
285e9fb
reorganize esolver_ks_lcao.h
mohanchen Mar 6, 2025
56b7ca6
delete some useless files in esolver_ks.h
mohanchen Mar 6, 2025
b10c6b8
change pelec->charge to this->chr
mohanchen Mar 6, 2025
9f0f689
will delete pelec->charge in near future, we just directly use chr de…
mohanchen Mar 6, 2025
27eae7c
fix problem in ML_KEDF
mohanchen Mar 6, 2025
82ba0b6
update this->chr
mohanchen Mar 6, 2025
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
50 changes: 26 additions & 24 deletions source/module_esolver/esolver_fp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,15 +23,16 @@ namespace ModuleESolver

ESolver_FP::ESolver_FP()
{
// pw_rho = new ModuleBase::PW_Basis();
// LCAO basis doesn't support GPU acceleration on FFT currently
std::string fft_device = PARAM.inp.device;

// LCAO basis doesn't support GPU acceleration on FFT currently
if(PARAM.inp.basis_type == "lcao")
{
fft_device = "cpu";
}

pw_rho = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
if ( PARAM.globalv.double_grid)
if (PARAM.globalv.double_grid)
{
pw_rhod = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
}
Expand All @@ -44,6 +45,7 @@ ESolver_FP::ESolver_FP()
pw_big = static_cast<ModulePW::PW_Basis_Big*>(pw_rhod);
pw_big->setbxyz(PARAM.inp.bx, PARAM.inp.by, PARAM.inp.bz);
sf.set(pw_rhod, PARAM.inp.nbspline);

}

ESolver_FP::~ESolver_FP()
Expand Down Expand Up @@ -140,7 +142,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
// 2) write fermi energy
ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);

// 3) update delta rho for charge extrapolation
// 3) update delta_rho for charge extrapolation
CE.update_delta_rho(ucell, &(this->chr), &(this->sf));

if (istep % PARAM.inp.out_interval == 0)
Expand All @@ -153,13 +155,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
double* data = nullptr;
if (PARAM.inp.dm_to_rho)
{
data = this->pelec->charge->rho[is];
this->pw_rhod->real2recip(this->pelec->charge->rho[is], this->pelec->charge->rhog[is]);
data = this->chr.rho[is];
this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]);
}
else
{
data = this->pelec->charge->rho_save[is];
this->pw_rhod->real2recip(this->pelec->charge->rho_save[is], this->pelec->charge->rhog_save[is]);
data = this->chr.rho_save[is];
this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
}
std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
ModuleIO::write_vdata_palgrid(Pgrid,
Expand All @@ -176,7 +178,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
{
fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube";
ModuleIO::write_vdata_palgrid(Pgrid,
this->pelec->charge->kin_r_save[is],
this->chr.kin_r_save[is],
is,
PARAM.inp.nspin,
istep,
Expand Down Expand Up @@ -217,7 +219,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
fn,
istep,
this->pw_rhod,
this->pelec->charge,
&this->chr,
&(ucell),
this->pelec->pot->get_fixed_v(),
this->solvent);
Expand All @@ -226,11 +228,11 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
// 6) write ELF
if (PARAM.inp.out_elf[0] > 0)
{
this->pelec->charge->cal_elf = true;
this->chr.cal_elf = true;
Symmetry_rho srho;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
}

std::string out_dir =PARAM.globalv.global_out_dir;
Expand All @@ -242,8 +244,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
out_dir,
istep,
PARAM.inp.nspin,
this->pelec->charge->rho,
this->pelec->charge->kin_r,
this->chr.rho,
this->chr.kin_r,
this->pw_rhod,
this->Pgrid,
&(ucell),
Expand All @@ -260,9 +262,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
{
// only G-vector and K-vector are changed due to the change of lattice
// vector FFT grids do not change!!
pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
pw_rho->collect_local_pw();
pw_rho->collect_uniqgg();
this->pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
this->pw_rho->collect_local_pw();
this->pw_rho->collect_uniqgg();

if (PARAM.globalv.double_grid)
{
Expand Down Expand Up @@ -292,7 +294,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
this->CE.update_all_dis(ucell);
this->CE.extrapolate_charge(&(this->Pgrid),
ucell,
this->pelec->charge,
&this->chr,
&(this->sf),
GlobalV::ofs_running,
GlobalV::ofs_warning);
Expand Down Expand Up @@ -327,7 +329,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
std::stringstream ss;
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
ModuleIO::write_vdata_palgrid(this->Pgrid,
this->pelec->charge->rho[is],
this->chr.rho[is],
is,
PARAM.inp.nspin,
istep,
Expand Down Expand Up @@ -368,8 +370,8 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
{
std::complex<double>** rhog_tot
= (PARAM.inp.dm_to_rho) ? this->pelec->charge->rhog : this->pelec->charge->rhog_save;
double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rho : this->pelec->charge->rho_save;
= (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save;
double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]);
Expand All @@ -386,12 +388,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&

if (XC_Functional::get_ked_flag())
{
std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->pelec->charge->ngmc, {0.0, 0.0});
std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->chr.ngmc, {0.0, 0.0});
std::vector<std::complex<double>*> kin_g;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
kin_g.push_back(kin_g_space.data() + is * this->pelec->charge->ngmc);
this->pw_rhod->real2recip(this->pelec->charge->kin_r_save[is], kin_g[is]);
kin_g.push_back(kin_g_space.data() + is * this->chr.ngmc);
this->pw_rhod->real2recip(this->chr.kin_r_save[is], kin_g[is]);
}
ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-TAU-DENSITY.restart",
PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local,
Expand Down
52 changes: 31 additions & 21 deletions source/module_esolver/esolver_fp.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,31 @@
#define ESOLVER_FP_H

#include "esolver.h"

//! plane wave basis
#include "module_basis/module_pw/pw_basis.h"

//! symmetry analysis
#include "module_cell/module_symmetry/symmetry.h"

//! electronic states
#include "module_elecstate/elecstate.h"

//! charge extrapolation
#include "module_elecstate/module_charge/charge_extra.h"

//! solvation model
#include "module_hamilt_general/module_surchem/surchem.h"

//! local pseudopotential
#include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"

//! structure factor related to plane wave basis
#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"

#include <fstream>


//! The First-Principles (FP) Energy Solver Class
/**
* This class represents components that needed in
Expand All @@ -22,7 +37,7 @@

namespace ModuleESolver
{
class ESolver_FP : public ESolver
class ESolver_FP: public ESolver
{
public:
//! Constructor
Expand All @@ -49,39 +64,34 @@ class ESolver_FP : public ESolver
//! ------------------------------------------------------------------------------
elecstate::ElecState* pelec = nullptr; ///< Electronic states

//! ------------------------------------------------------------------------------
//! K points in Brillouin zone
K_Vectors kv;

//! Electorn charge density
Charge chr;

//! Structure factors that used with plane-wave basis set
Structure_Factor sf;

//! K points in Brillouin zone
K_Vectors kv;

//! Plane-wave basis set for charge density
//! pw_rho: Plane-wave basis set for charge density
//! pw_rhod: same as pw_rho for NCPP. Here 'd' stands for 'dense',
//! dense grid for for uspp, used for ultrasoft augmented charge density.
//! charge density and potential are defined on dense grids,
//! but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
ModulePW::PW_Basis* pw_rho;
ModulePW::PW_Basis* pw_rhod; //! dense grid for USPP
ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class

//! parallel for rho grid
Parallel_Grid Pgrid;

//! pointer to local pseudopotential
pseudopot_cell_vl locpp;
//! Structure factors that used with plane-wave basis set
Structure_Factor sf;

/**
* @brief same as pw_rho for ncpp. Here 'd' stands for 'dense'
* dense grid for for uspp, used for ultrasoft augmented charge density.
* charge density and potential are defined on dense grids,
* but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
*/
ModulePW::PW_Basis* pw_rhod;
ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class
//! local pseudopotentials
pseudopot_cell_vl locpp;

//! Charge extrapolation
//! charge extrapolation method
Charge_Extra CE;

// solvent model
//! solvent model
surchem solvent;
};
} // namespace ModuleESolver
Expand Down
21 changes: 11 additions & 10 deletions source/module_esolver/esolver_ks.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "module_io/write_istate_info.h"
#include "module_parameter/parameter.h"
#include "module_elecstate/elecstate_print.h"
#include "module_hsolver/hsolver.h"

#include <ctime>
#include <iostream>
Expand Down Expand Up @@ -367,7 +368,7 @@ void ESolver_KS<T, Device>::hamilt2density(UnitCell& ucell, const int istep, con
// double drho = this->estate.caldr2();
// EState should be used after it is constructed.

drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec);
drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
hsolver_error = 0.0;
if (iter == 1 && PARAM.inp.calculation != "nscf")
{
Expand All @@ -389,7 +390,7 @@ void ESolver_KS<T, Device>::hamilt2density(UnitCell& ucell, const int istep, con

this->hamilt2density_single(ucell, istep, iter, diag_ethr);

drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec);
drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);

hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type,
PARAM.inp.esolver_type,
Expand Down Expand Up @@ -521,7 +522,7 @@ void ESolver_KS<T, Device>::iter_init(UnitCell& ucell, const int istep, const in
}

// 1) save input rho
this->pelec->charge->save_rho_before_sum_band();
this->chr.save_rho_before_sum_band();
}

template <typename T, typename Device>
Expand Down Expand Up @@ -551,9 +552,9 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i

// compute magnetization, only for LSDA(spin==2)
ucell.magnet.compute_magnetization(ucell.omega,
this->pelec->charge->nrxx,
this->pelec->charge->nxyz,
this->pelec->charge->rho,
this->chr.nrxx,
this->chr.nxyz,
this->chr.rho,
this->pelec->nelec_spin.data());

if (PARAM.globalv.ks_run)
Expand Down Expand Up @@ -628,11 +629,11 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
}
else
{
p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing
p_chgmix->mix_rho(&this->chr); // update chr->rho by mixing
}
if (PARAM.inp.scf_thr_type == 2)
{
pelec->charge->renormalize_rho(); // renormalize rho in R-space would
this->chr.renormalize_rho(); // renormalize rho in R-space would
// induce a error in K-space
}
//----------charge mixing done-----------
Expand All @@ -644,7 +645,7 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i

// be careful! conv_esolver is bool, not double !! Maybe a bug 20250302 by mohan
MPI_Bcast(&conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD);
MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD);
MPI_Bcast(this->chr.rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD);
#endif

// update potential
Expand Down Expand Up @@ -676,7 +677,7 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
double dkin = 0.0; // for meta-GGA
if (XC_Functional::get_ked_flag())
{
dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec);
dkin = p_chgmix->get_dkin(&this->chr, PARAM.inp.nelec);
}


Expand Down
Loading
Loading