Skip to content

Commit e3bff8f

Browse files
mohanchenmohanchen
committed
reduce the lines of esolver_ks_lcao
1 parent b5b66d6 commit e3bff8f

File tree

1 file changed

+48
-140
lines changed

1 file changed

+48
-140
lines changed

‎source/source_esolver/esolver_ks_lcao.cpp‎

Lines changed: 48 additions & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -112,32 +112,19 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
112112
ESolver_KS<TK>::before_all_runners(ucell, inp);
113113

114114
// 2) init ElecState
115-
// autoset nbands in ElecState before basis_init (for Psi 2d division)
115+
// autoset nbands in ElecState before init_basis (for Psi 2d division)
116116
if (this->pelec == nullptr)
117117
{
118118
// TK stands for double and std::complex<double>?
119-
this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), // use which parameter?
120-
&(this->kv),
121-
this->kv.get_nks(),
122-
&(this->GG),
123-
&(this->GK),
124-
this->pw_rho,
125-
this->pw_big);
126-
}
127-
128-
// 3) init LCAO basis
129-
// reading the localized orbitals/projectors
130-
// construct the interpolation tables.
131-
LCAO_domain::init_basis_lcao(this->pv,
132-
inp.onsite_radius,
133-
inp.lcao_ecut,
134-
inp.lcao_dk,
135-
inp.lcao_dr,
136-
inp.lcao_rmax,
137-
ucell,
138-
two_center_bundle_,
139-
orb_);
119+
this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), &(this->kv),
120+
this->kv.get_nks(), &(this->GG), &(this->GK), this->pw_rho, this->pw_big);
121+
}
122+
123+
// 3) read the LCAO orbitals/projectors and construct the interpolation tables.
124+
LCAO_domain::init_basis_lcao(this->pv, inp.onsite_radius, inp.lcao_ecut,
125+
inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax, ucell, two_center_bundle_, orb_);
140126

127+
// 4) setup EXX calculations
141128
if (PARAM.inp.calculation == "gen_opt_abfs")
142129
{
143130
#ifdef __EXX
@@ -149,7 +136,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
149136
return;
150137
}
151138

152-
// 4) initialize electronic wave function psi
139+
// 5) initialize electronic wave function psi
153140
if (this->psi == nullptr)
154141
{
155142
int nsk = 0;
@@ -177,7 +164,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
177164
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, this->kv.ngk, true);
178165
}
179166

180-
// 5) read psi from file
167+
// 6) read psi from file
181168
if (inp.init_wfc == "file" && inp.esolver_type != "tddft")
182169
{
183170
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir,
@@ -192,12 +179,11 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
192179
}
193180
}
194181

195-
// 6) initialize the density matrix
196-
// DensityMatrix is allocated here, DMK is also initialized here
197-
// DMR is not initialized here, it will be constructed in each before_scf
182+
// 7) initialize the density matrix
183+
// DMK are allocated here, but DMR is constructed in before_scf()
198184
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin);
199185

200-
// 7) initialize exact exchange calculations
186+
// 8) initialize exact exchange calculations
201187
#ifdef __EXX
202188
if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax"
203189
|| inp.calculation == "md")
@@ -224,67 +210,52 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
224210
}
225211
#endif
226212

227-
// 8) initialize DFT+U
213+
// 9) initialize DFT+U
228214
if (inp.dft_plus_u)
229215
{
230216
auto* dftu = ModuleDFTU::DFTU::get_instance();
231217
dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
232218
}
233219

234-
// 9) initialize local pseudopotentials
220+
// 10) initialize local pseudopotentials
235221
this->locpp.init_vloc(ucell, this->pw_rho);
236222
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
237223

238-
// 10) inititlize the charge density
224+
// 11) inititlize the charge density
239225
this->chr.allocate(inp.nspin);
240226
this->pelec->omega = ucell.omega;
241227

242-
// 11) initialize the potential
228+
// 12) initialize the potential
243229
if (this->pelec->pot == nullptr)
244230
{
245-
this->pelec->pot = new elecstate::Potential(this->pw_rhod,
246-
this->pw_rho,
247-
&ucell,
248-
&(this->locpp.vloc),
249-
&(this->sf),
250-
&(this->solvent),
251-
&(this->pelec->f_en.etxc),
252-
&(this->pelec->f_en.vtxc));
231+
this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho,
232+
&ucell, &(this->locpp.vloc), &(this->sf), &(this->solvent),
233+
&(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc));
253234
}
254235

255-
// 12) initialize deepks
236+
// 13) initialize deepks
256237
#ifdef __MLALGO
257238
LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld, GlobalV::ofs_running);
258239
if (inp.deepks_scf)
259240
{
260241
// load the DeePKS model from deep neural network
261242
DeePKS_domain::load_model(inp.deepks_model, ld.model_deepks);
262243
// read pdm from file for NSCF or SCF-restart, do it only once in whole calculation
263-
DeePKS_domain::read_pdm((inp.init_chg == "file"),
264-
inp.deepks_equiv,
265-
ld.init_pdm,
266-
ucell.nat,
267-
orb_.Alpha[0].getTotal_nchi() * ucell.nat,
268-
ld.lmaxd,
269-
ld.inl2l,
270-
*orb_.Alpha,
271-
ld.pdm);
244+
DeePKS_domain::read_pdm((inp.init_chg == "file"), inp.deepks_equiv,
245+
ld.init_pdm, ucell.nat, orb_.Alpha[0].getTotal_nchi() * ucell.nat,
246+
ld.lmaxd, ld.inl2l, *orb_.Alpha, ld.pdm);
272247
}
273248
#endif
274249

275-
// 13) set occupations
250+
// 14) set occupations
276251
// tddft does not need to set occupations in the first scf
277252
if (inp.ocp && inp.esolver_type != "tddft")
278253
{
279-
elecstate::fixed_weights(inp.ocp_kb,
280-
inp.nbands,
281-
inp.nelec,
282-
this->pelec->klist,
283-
this->pelec->wg,
284-
this->pelec->skip_weights);
254+
elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec,
255+
this->pelec->klist, this->pelec->wg, this->pelec->skip_weights);
285256
}
286257

287-
// 14) if kpar is not divisible by nks, print a warning
258+
// 15) if kpar is not divisible by nks, print a warning
288259
if (PARAM.globalv.kpar_lcao > 1)
289260
{
290261
if (this->kv.get_nks() % PARAM.globalv.kpar_lcao != 0)
@@ -305,20 +276,12 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
305276
}
306277
}
307278

308-
// 15) initialize rdmft, added by jghan
279+
// 16) initialize rdmft, added by jghan
309280
if (inp.rdmft == true)
310281
{
311-
rdmft_solver.init(this->GG,
312-
this->GK,
313-
this->pv,
314-
ucell,
315-
this->gd,
316-
this->kv,
317-
*(this->pelec),
318-
this->orb_,
319-
two_center_bundle_,
320-
inp.dft_functional,
321-
inp.rdmft_power_alpha);
282+
rdmft_solver.init(this->GG, this->GK, this->pv, ucell,
283+
this->gd, this->kv, *(this->pelec), this->orb_,
284+
two_center_bundle_, inp.dft_functional, inp.rdmft_power_alpha);
322285
}
323286

324287
ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners");
@@ -568,28 +531,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
568531
this->pelec->psiToRho(*this->psi);
569532
this->pelec->skip_weights = false;
570533

571-
// calculate the local potential(rho) again.
572-
// the grid integration will do in later grid integration.
573-
574-
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
575-
// a puzzle remains here.
576-
// if I don't renew potential,
577-
// The scf_thr is very small.
578-
// OneElectron, Hartree and
579-
// Exc energy are all correct
580-
// except the band energy.
581-
//
582-
// solved by mohan 2010-09-10
583-
// there are there rho here:
584-
// rho1: formed by read in orbitals.
585-
// rho2: atomic rho, used to construct H
586-
// rho3: generated by after diagonalize
587-
// here converged because rho3 and rho1
588-
// are very close.
589-
// so be careful here, make sure
590-
// rho1 and rho2 are the same rho.
591-
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592-
593534
elecstate::cal_ux(ucell);
594535

595536
//! update the potentials by using new electron charge density
@@ -606,19 +547,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
606547
{
607548
if (GlobalC::exx_info.info_ri.real_number)
608549
{
609-
this->exd->exx_eachiterinit(istep,
610-
ucell,
611-
*dm,
612-
this->kv,
613-
iter);
550+
this->exd->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
614551
}
615552
else
616553
{
617-
this->exc->exx_eachiterinit(istep,
618-
ucell,
619-
*dm,
620-
this->kv,
621-
iter);
554+
this->exc->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
622555
}
623556
}
624557
#endif
@@ -747,8 +680,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
747680
{
748681
ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
749682

750-
// 1) calculate the local occupation number matrix and energy correction
751-
// in DFT+U
683+
// 1) calculate the local occupation number matrix and energy correction in DFT+U
752684
if (PARAM.inp.dft_plus_u)
753685
{
754686
// only old DFT+U method should calculated energy correction in esolver,
@@ -801,7 +733,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
801733
{
802734
if (PARAM.inp.mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && PARAM.inp.mixing_dmr)
803735
{
804-
elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
736+
elecstate::DensityMatrix<TK, double>* dm =
737+
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
805738
this->p_chgmix->mix_dmr(dm);
806739
}
807740
}
@@ -821,24 +754,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
821754
{
822755
if (GlobalC::exx_info.info_global.cal_exx)
823756
{
824-
GlobalC::exx_info.info_ri.real_number ? this->exd->exx_iter_finish(this->kv,
825-
ucell,
826-
*this->p_hamilt,
827-
*this->pelec,
828-
*this->p_chgmix,
829-
this->scf_ene_thr,
830-
iter,
831-
istep,
832-
conv_esolver)
833-
: this->exc->exx_iter_finish(this->kv,
834-
ucell,
835-
*this->p_hamilt,
836-
*this->pelec,
837-
*this->p_chgmix,
838-
this->scf_ene_thr,
839-
iter,
840-
istep,
841-
conv_esolver);
757+
GlobalC::exx_info.info_ri.real_number ?
758+
this->exd->exx_iter_finish(this->kv, ucell, *this->p_hamilt, *this->pelec,
759+
*this->p_chgmix, this->scf_ene_thr, iter, istep, conv_esolver) :
760+
this->exc->exx_iter_finish(this->kv, ucell, *this->p_hamilt, *this->pelec,
761+
*this->p_chgmix, this->scf_ene_thr, iter, istep, conv_esolver);
842762
}
843763
}
844764
#endif
@@ -859,23 +779,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
859779
std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&ld, [](LCAO_Deepks<TK>*) {});
860780
LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
861781

862-
deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
863-
this->kv.get_nks(),
864-
ucell.nat,
865-
PARAM.globalv.nlocal,
866-
this->pelec->ekb,
867-
this->kv.kvec_d,
868-
ucell,
869-
orb_,
870-
this->gd,
871-
&(this->pv),
872-
*(this->psi),
873-
dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
874-
p_ham_deepks,
875-
iter,
876-
conv_esolver,
877-
GlobalV::MY_RANK,
878-
GlobalV::ofs_running);
782+
deepks_interface.out_deepks_labels(this->pelec->f_en.etot, this->kv.get_nks(),
783+
ucell.nat, PARAM.globalv.nlocal, this->pelec->ekb, this->kv.kvec_d,
784+
ucell, orb_, this->gd, &(this->pv), *(this->psi),
785+
dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
786+
p_ham_deepks, iter, conv_esolver, GlobalV::MY_RANK, GlobalV::ofs_running);
879787
}
880788
}
881789
#endif

0 commit comments

Comments
 (0)