Skip to content

Refactor: move psi init to before_all_runners #5992

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 1 commit into from
Mar 11, 2025
Merged

Refactor: move psi init to before_all_runners #5992

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
61 changes: 49 additions & 12 deletions source/module_esolver/esolver_ks_lcao.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "module_base/formatter.h"
#include "module_base/global_variable.h"
#include "module_base/tool_title.h"
#include "module_elecstate/elecstate_tools.h"
#include "module_elecstate/module_dm/cal_dm_psi.h"
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
#include "module_hamilt_lcao/module_dftu/dftu.h"
Expand All @@ -17,6 +18,7 @@
#include "module_io/output_mat_sparse.h"
#include "module_io/output_mulliken.h"
#include "module_io/output_sk.h"
#include "module_io/read_wfc_nao.h"
#include "module_io/to_qo.h"
#include "module_io/to_wannier90_lcao.h"
#include "module_io/to_wannier90_lcao_in_pw.h"
Expand All @@ -27,7 +29,6 @@
#include "module_io/write_proj_band_lcao.h"
#include "module_io/write_wfc_nao.h"
#include "module_parameter/parameter.h"
#include "module_elecstate/elecstate_tools.h"

//be careful of hpp, there may be multiple definitions of functions, 20250302, mohan
#include "module_io/write_eband_terms.hpp"
Expand Down Expand Up @@ -133,12 +134,49 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
two_center_bundle_,
orb_);

// 4) initialize the density matrix
// 4) initialize electronic wave function psi
if (this->psi == nullptr)
{
int nsk = 0;
int ncol = 0;
if (PARAM.globalv.gamma_only_local)
{
nsk = PARAM.inp.nspin;
ncol = this->pv.ncol_bands;
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
|| PARAM.inp.ks_solver == "cusolvermp")
{
ncol = this->pv.ncol;
}
}
else
{
nsk = this->kv.get_nks();
#ifdef __MPI
ncol = this->pv.ncol_bands;
#else
ncol = PARAM.inp.nbands;
#endif
}
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
}

// 5) read psi from file
if (PARAM.inp.init_wfc == "file")
{
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
{
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
}
}

// 6) initialize the density matrix
// DensityMatrix is allocated here, DMK is also initialized here
// DMR is not initialized here, it will be constructed in each before_scf
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);

// 5) initialize exact exchange calculations
// 7) initialize exact exchange calculations
#ifdef __EXX
if (PARAM.inp.calculation == "scf"
|| PARAM.inp.calculation == "relax"
Expand All @@ -163,22 +201,22 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
}
#endif

// 6) initialize DFT+U
// 8) initialize DFT+U
if (PARAM.inp.dft_plus_u)
{
auto* dftu = ModuleDFTU::DFTU::get_instance();
dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
}

// 7) initialize local pseudopotentials
// 9) initialize local pseudopotentials
this->locpp.init_vloc(ucell, this->pw_rho);
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");

// 8) inititlize the charge density
// 10) inititlize the charge density
this->chr.allocate(PARAM.inp.nspin);
this->pelec->omega = ucell.omega;

// 9) initialize the potential
// 11) initialize the potential
if (this->pelec->pot == nullptr)
{
this->pelec->pot = new elecstate::Potential(this->pw_rhod,
Expand All @@ -191,8 +229,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
&(this->pelec->f_en.vtxc));
}


// 10) initialize deepks
// 12) initialize deepks
#ifdef __DEEPKS
LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld, GlobalV::ofs_running);
if (PARAM.inp.deepks_scf)
Expand All @@ -212,7 +249,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
}
#endif

// 11) set occupations
// 13) set occupations
// tddft does not need to set occupations in the first scf
if (PARAM.inp.ocp && inp.esolver_type != "tddft")
{
Expand All @@ -224,7 +261,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
this->pelec->skip_weights);
}

// 12) if kpar is not divisible by nks, print a warning
// 14) if kpar is not divisible by nks, print a warning
if (PARAM.globalv.kpar_lcao > 1)
{
if (this->kv.get_nks() % PARAM.globalv.kpar_lcao != 0)
Expand All @@ -245,7 +282,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
}
}

// 13) initialize rdmft, added by jghan
// 15) initialize rdmft, added by jghan
if (PARAM.inp.rdmft == true)
{
rdmft_solver.init(this->GG,
Expand Down
68 changes: 11 additions & 57 deletions source/module_esolver/lcao_before_scf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
#include "module_io/cube_io.h"
#include "module_io/read_wfc_nao.h"
#include "module_io/write_elecstat_pot.h"
#include "module_io/write_wfc_nao.h"
#ifdef __EXX
Expand Down Expand Up @@ -128,52 +127,10 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
// If k point is used here, allocate HlocR after atom_arrange.
this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());


// 8) initialize electronic wave function psi
// this code belongs to before_all_runners? mohan add 2025-03-10
if (this->psi == nullptr)
{
int nsk = 0;
int ncol = 0;
if (PARAM.globalv.gamma_only_local)
{
nsk = PARAM.inp.nspin;
ncol = this->pv.ncol_bands;
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
|| PARAM.inp.ks_solver == "cusolvermp")
{
ncol = this->pv.ncol;
}
}
else
{
nsk = this->kv.get_nks();
#ifdef __MPI
ncol = this->pv.ncol_bands;
#else
ncol = PARAM.inp.nbands;
#endif
}
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
}

// 9) read psi from file
// this code belongs to before_all_runners? mohan add 2025-03-10
if (istep == 0 && PARAM.inp.init_wfc == "file")
{
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
{
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
}
}


// 10) after ions move, prepare grid in Gint
// 8) after ions move, prepare grid in Gint
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);


// 11) initialize the Hamiltonian operators
// 9) initialize the Hamiltonian operators
// if atom moves, then delete old pointer and add a new one
if (this->p_hamilt != nullptr)
{
Expand Down Expand Up @@ -213,7 +170,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)


#ifdef __DEEPKS
// 12) for each ionic step, the overlap <phi|alpha> must be rebuilt
// 10) for each ionic step, the overlap <phi|alpha> must be rebuilt
// since it depends on ionic positions
if (PARAM.globalv.deepks_setorb)
{
Expand Down Expand Up @@ -242,7 +199,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
}
#endif

// 13) prepare sc calculation
// 11) prepare sc calculation
if (PARAM.inp.sc_mag_switch)
{
spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
Expand All @@ -261,7 +218,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
this->pelec);
}

// 14) set xc type before the first cal of xc in pelec->init_scf
// 12) set xc type before the first cal of xc in pelec->init_scf
// Peize Lin add 2016-12-03
#ifdef __EXX
if (PARAM.inp.calculation != "nscf")
Expand All @@ -277,18 +234,16 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
}
#endif

// 15) init_scf, should be before_scf? mohan add 2025-03-10
// 13) init_scf, should be before_scf? mohan add 2025-03-10
this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);


// 16) initalize DMR
// 14) initalize DMR
// DMR should be same size with Hamiltonian(R)
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
->get_DM()
->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt)->getHR()));


// 17) two cases are considered:
// 15) two cases are considered:
// 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
// 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
if (istep > 0)
Expand Down Expand Up @@ -347,19 +302,18 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
return;
}

// 18) the electron charge density should be symmetrized,
// 16) the electron charge density should be symmetrized,
// here is the initialization
Symmetry_rho srho;
for (int is = 0; is < PARAM.inp.nspin; is++)
{
srho.begin(is, this->chr, this->pw_rho, ucell.symm);
}

// 19) why we need to set this sentence? mohan add 2025-03-10
// 17) why we need to set this sentence? mohan add 2025-03-10
this->p_hamilt->non_first_scf = istep;


// 20) update of RDMFT, added by jghan
// 18) update of RDMFT, added by jghan
if (PARAM.inp.rdmft == true)
{
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
Expand Down
Loading