Skip to content

Commit f3c6cfc

Browse files
YuLiu98YuLiu98
authored
Refactor: move psi init to before_all_runners (#5992)
1 parent 33d0b4f commit f3c6cfc

File tree

2 files changed

+60
-69
lines changed

2 files changed

+60
-69
lines changed

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 49 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "module_base/formatter.h"
44
#include "module_base/global_variable.h"
55
#include "module_base/tool_title.h"
6+
#include "module_elecstate/elecstate_tools.h"
67
#include "module_elecstate/module_dm/cal_dm_psi.h"
78
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
89
#include "module_hamilt_lcao/module_dftu/dftu.h"
@@ -17,6 +18,7 @@
1718
#include "module_io/output_mat_sparse.h"
1819
#include "module_io/output_mulliken.h"
1920
#include "module_io/output_sk.h"
21+
#include "module_io/read_wfc_nao.h"
2022
#include "module_io/to_qo.h"
2123
#include "module_io/to_wannier90_lcao.h"
2224
#include "module_io/to_wannier90_lcao_in_pw.h"
@@ -27,7 +29,6 @@
2729
#include "module_io/write_proj_band_lcao.h"
2830
#include "module_io/write_wfc_nao.h"
2931
#include "module_parameter/parameter.h"
30-
#include "module_elecstate/elecstate_tools.h"
3132

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

136-
// 4) initialize the density matrix
137+
// 4) initialize electronic wave function psi
138+
if (this->psi == nullptr)
139+
{
140+
int nsk = 0;
141+
int ncol = 0;
142+
if (PARAM.globalv.gamma_only_local)
143+
{
144+
nsk = PARAM.inp.nspin;
145+
ncol = this->pv.ncol_bands;
146+
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
147+
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
148+
|| PARAM.inp.ks_solver == "cusolvermp")
149+
{
150+
ncol = this->pv.ncol;
151+
}
152+
}
153+
else
154+
{
155+
nsk = this->kv.get_nks();
156+
#ifdef __MPI
157+
ncol = this->pv.ncol_bands;
158+
#else
159+
ncol = PARAM.inp.nbands;
160+
#endif
161+
}
162+
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
163+
}
164+
165+
// 5) read psi from file
166+
if (PARAM.inp.init_wfc == "file")
167+
{
168+
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
169+
{
170+
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
171+
}
172+
}
173+
174+
// 6) initialize the density matrix
137175
// DensityMatrix is allocated here, DMK is also initialized here
138176
// DMR is not initialized here, it will be constructed in each before_scf
139177
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
140178

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

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

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

177-
// 8) inititlize the charge density
215+
// 10) inititlize the charge density
178216
this->chr.allocate(PARAM.inp.nspin);
179217
this->pelec->omega = ucell.omega;
180218

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

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

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

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

248-
// 13) initialize rdmft, added by jghan
285+
// 15) initialize rdmft, added by jghan
249286
if (PARAM.inp.rdmft == true)
250287
{
251288
rdmft_solver.init(this->GG,

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 11 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,6 @@
2727
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
2828
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
2929
#include "module_io/cube_io.h"
30-
#include "module_io/read_wfc_nao.h"
3130
#include "module_io/write_elecstat_pot.h"
3231
#include "module_io/write_wfc_nao.h"
3332
#ifdef __EXX
@@ -128,52 +127,10 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
128127
// If k point is used here, allocate HlocR after atom_arrange.
129128
this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
130129

131-
132-
// 8) initialize electronic wave function psi
133-
// this code belongs to before_all_runners? mohan add 2025-03-10
134-
if (this->psi == nullptr)
135-
{
136-
int nsk = 0;
137-
int ncol = 0;
138-
if (PARAM.globalv.gamma_only_local)
139-
{
140-
nsk = PARAM.inp.nspin;
141-
ncol = this->pv.ncol_bands;
142-
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
143-
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
144-
|| PARAM.inp.ks_solver == "cusolvermp")
145-
{
146-
ncol = this->pv.ncol;
147-
}
148-
}
149-
else
150-
{
151-
nsk = this->kv.get_nks();
152-
#ifdef __MPI
153-
ncol = this->pv.ncol_bands;
154-
#else
155-
ncol = PARAM.inp.nbands;
156-
#endif
157-
}
158-
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
159-
}
160-
161-
// 9) read psi from file
162-
// this code belongs to before_all_runners? mohan add 2025-03-10
163-
if (istep == 0 && PARAM.inp.init_wfc == "file")
164-
{
165-
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
166-
{
167-
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed");
168-
}
169-
}
170-
171-
172-
// 10) after ions move, prepare grid in Gint
130+
// 8) after ions move, prepare grid in Gint
173131
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
174132

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

214171

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

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

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

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

283-
284-
// 16) initalize DMR
240+
// 14) initalize DMR
285241
// DMR should be same size with Hamiltonian(R)
286242
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
287243
->get_DM()
288244
->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt)->getHR()));
289245

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

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

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

361-
362-
// 20) update of RDMFT, added by jghan
316+
// 18) update of RDMFT, added by jghan
363317
if (PARAM.inp.rdmft == true)
364318
{
365319
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp

0 commit comments

Comments
 (0)