Skip to content

Commit 408b661

Browse files
19helloFei Yang
19hello
and
Fei Yang
authored
Refactor init_esolver factory interface (deepmodeling#7702)
Co-authored-by: Fei Yang <2501213217@stu.pku.edu.cn>
1 parent cab0ecd commit 408b661

5 files changed

Lines changed: 52 additions & 63 deletions

File tree

‎source/source_esolver/esolver.cpp‎

Lines changed: 5 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ std::string determine_type()
123123
}
124124

125125
// Some API to operate E_Solver
126-
ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
126+
ESolver* init_esolver(const Input_para& inp)
127127
{
128128
// determine type of esolver based on INPUT information
129129
const std::string esolver_type = determine_type();
@@ -273,58 +273,25 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
273273
}
274274
else if (esolver_type == "lr_lcao")
275275
{
276-
// use constructor rather than Init function to initialize reference (instead of pointers) to ucell
277276
if (PARAM.globalv.gamma_only_local)
278277
{
279-
return new LR::ESolver_LR<double, double>(inp, ucell);
278+
return new LR::ESolver_LR<double, double>(inp);
280279
}
281280
else
282281
{
283-
return new LR::ESolver_LR<std::complex<double>, double>(inp, ucell);
282+
return new LR::ESolver_LR<std::complex<double>, double>(inp);
284283
}
285284
}
286285
else if (esolver_type == "ksdft_lr_lcao")
287286
{
288-
// initialize the 1st ESolver_KS
289-
ModuleESolver::ESolver* p_esolver = nullptr;
290287
if (PARAM.globalv.gamma_only_local)
291288
{
292-
p_esolver = new ESolver_KS_LCAO<double, double>();
293-
}
294-
else if (PARAM.inp.nspin < 4)
295-
{
296-
p_esolver = new ESolver_KS_LCAO<std::complex<double>, double>();
289+
return new LR::ESolver_LR<double, double>(inp);
297290
}
298291
else
299292
{
300-
p_esolver = new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
293+
return new LR::ESolver_LR<std::complex<double>, double>(inp);
301294
}
302-
p_esolver->before_all_runners(ucell, inp);
303-
p_esolver->runner(ucell, 0); // scf-only
304-
305-
// force and stress is not needed currently,
306-
// they will be supported after the analytical gradient
307-
// of LR-TDDFT is implemented.
308-
std::cout << " PREPARING FOR EXCITED STATES." << std::endl;
309-
// initialize the 2nd ESolver_LR at the temporary pointer
310-
ModuleESolver::ESolver* p_esolver_lr = nullptr;
311-
if (PARAM.globalv.gamma_only_local)
312-
{
313-
p_esolver_lr = new LR::ESolver_LR<double, double>(
314-
std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<double, double>*>(p_esolver)),
315-
inp,
316-
ucell);
317-
}
318-
else
319-
{
320-
p_esolver_lr = new LR::ESolver_LR<std::complex<double>, double>(
321-
std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<std::complex<double>, double>*>(p_esolver)),
322-
inp,
323-
ucell);
324-
}
325-
// clean the 1st ESolver_KS and swap the pointer
326-
delete p_esolver;
327-
return p_esolver_lr;
328295
}
329296
#endif
330297
else if (esolver_type == "ofdft")

‎source/source_esolver/esolver.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ std::string determine_type();
6868
*
6969
* @return [out] A pointer to an ESolver object that will be initialized.
7070
*/
71-
ESolver* init_esolver(const Input_para& inp, UnitCell& ucell);
71+
ESolver* init_esolver(const Input_para& inp);
7272

7373

7474

‎source/source_lcao/module_lr/esolver_lrtd_lcao.cpp‎

Lines changed: 37 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ void LR::ESolver_LR<T, TR>::set_dimension()
107107
this->nbasis = PARAM.globalv.nlocal;
108108
// calculate the number of occupied and unoccupied states
109109
// which determines the basis size of the excited states
110-
this->nocc_max = LR_Util::cal_nocc(LR_Util::cal_nelec(ucell));
110+
this->nocc_max = LR_Util::cal_nocc(LR_Util::cal_nelec(*this->ucell_));
111111
this->nocc_in = std::max(1, std::min(input.nocc, this->nocc_max));
112112
this->nvirt_in = PARAM.inp.nbands - this->nocc_max; //nbands-nocc
113113
if (input.nvirt > this->nvirt_in) { GlobalV::ofs_warning << "ESolver_LR: input nvirt is too large to cover by nbands, set nvirt = nbands - nocc = " << this->nvirt_in << std::endl; }
@@ -128,7 +128,7 @@ void LR::ESolver_LR<T, TR>::set_dimension()
128128
// calculate total number of basis funcs, see https://en.cppreference.com/w/cpp/algorithm/inner_product
129129
this->nbasis = std::inner_product(input.aims_nbasis.begin(), /* iterator1.begin */
130130
input.aims_nbasis.end(), /* iterator1.end */
131-
ucell.atoms, /* iterator2.begin */
131+
this->ucell_->atoms, /* iterator2.begin */
132132
0, /* init value */
133133
std::plus<int>(), /* iter op1 */
134134
[](const int& a, const Atom& b) { return a * b.na; }); /* iter op2 */
@@ -173,12 +173,35 @@ void LR::ESolver_LR<T, TR>::reset_dim_spin2()
173173
}
174174

175175
template <typename T, typename TR>
176-
LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol,
177-
const Input_para& inp, UnitCell& ucell)
178-
: input(inp), ucell(ucell)
176+
LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp)
177+
: input(inp)
179178
#ifdef __EXX
180179
, exx_info(GlobalC::exx_info)
181180
#endif
181+
{
182+
}
183+
184+
template <typename T, typename TR>
185+
void LR::ESolver_LR<T, TR>::before_all_runners(UnitCell& ucell, const Input_para& inp)
186+
{
187+
this->ucell_ = &ucell;
188+
if (inp.esolver_type == "ks-lr")
189+
{
190+
ModuleESolver::ESolver_KS_LCAO<T, TR> ks_solver;
191+
ks_solver.before_all_runners(ucell, inp);
192+
ks_solver.runner(ucell, 0);
193+
this->initialize_from_ks_(std::move(ks_solver), ucell, inp);
194+
}
195+
else
196+
{
197+
this->initialize_from_unitcell_(ucell, inp);
198+
}
199+
}
200+
201+
template <typename T, typename TR>
202+
void LR::ESolver_LR<T, TR>::initialize_from_ks_(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol,
203+
UnitCell& ucell,
204+
const Input_para& inp)
182205
{
183206
ModuleBase::TITLE("ESolver_LR", "ESolver_LR(KS)");
184207

@@ -289,10 +312,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol
289312
}
290313

291314
template <typename T, typename TR>
292-
LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : input(inp), ucell(ucell)
293-
#ifdef __EXX
294-
, exx_info(GlobalC::exx_info)
295-
#endif
315+
void LR::ESolver_LR<T, TR>::initialize_from_unitcell_(UnitCell& ucell, const Input_para& inp)
296316
{
297317
ModuleBase::TITLE("ESolver_LR", "ESolver_LR(from scratch)");
298318
// xc kernel
@@ -392,7 +412,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
392412
atom_arrange::search(PARAM.globalv.search_pbc,
393413
GlobalV::ofs_running,
394414
this->gd,
395-
this->ucell,
415+
*this->ucell_,
396416
search_radius,
397417
PARAM.inp.test_atom_input);
398418
gint_info_.reset(
@@ -463,7 +483,7 @@ void LR::ESolver_LR<T, TR>::runner(UnitCell& ucell, const int istep)
463483
this->nbasis,
464484
this->nocc,
465485
this->nvirt,
466-
this->ucell,
486+
*this->ucell_,
467487
orb_cutoff_,
468488
this->gd,
469489
*this->psi_ks,
@@ -493,7 +513,7 @@ void LR::ESolver_LR<T, TR>::runner(UnitCell& ucell, const int istep)
493513
this->nbasis,
494514
this->nocc,
495515
this->nvirt,
496-
this->ucell,
516+
*this->ucell_,
497517
orb_cutoff_,
498518
this->gd,
499519
*this->psi_ks,
@@ -559,7 +579,7 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
559579
for (int is = 0;is < this->X.size();++is)
560580
{
561581
LR_Spectrum<T> spectrum(nspin, this->nbasis, this->nocc, this->nvirt, *this->pw_rho, *this->psi_ks,
562-
this->ucell, this->kv, this->gd, this->orb_cutoff_, this->two_center_bundle_,
582+
*this->ucell_, this->kv, this->gd, this->orb_cutoff_, this->two_center_bundle_,
563583
this->paraX_, this->paraC_, this->paraMat_,
564584
&this->pelec->ekb.c[is * nstates], this->X[is].template data<T>(), nstates, openshell,
565585
LR_Util::tolower(input.abs_gauge));
@@ -656,11 +676,11 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
656676
{
657677
using ST = PotHxcLR::SpinType;
658678
case 1:
659-
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, Pgrid, ST::S1, input.lr_init_xc_kernel);
679+
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, *this->ucell_, chg_gs, Pgrid, ST::S1, input.lr_init_xc_kernel);
660680
break;
661681
case 2:
662-
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.lr_init_xc_kernel);
663-
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, ucell, chg_gs, Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel);
682+
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, *this->ucell_, chg_gs, Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.lr_init_xc_kernel);
683+
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, *this->pw_rho, *this->ucell_, chg_gs, Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.lr_init_xc_kernel);
664684
break;
665685
default:
666686
throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");
@@ -717,7 +737,7 @@ void LR::ESolver_LR<T, TR>::read_ks_chg(Charge& chg_gs)
717737
GlobalV::ofs_running,
718738
ssc.str(),
719739
chg_gs.rho[is],
720-
ucell.nat)) {
740+
this->ucell_->nat)) {
721741
GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
722742
} else { // prenspin for nspin=4 is not supported currently
723743
ModuleBase::WARNING_QUIT(

‎source/source_lcao/module_lr/esolver_lrtd_lcao.h‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -26,17 +26,14 @@ namespace LR
2626
class ESolver_LR : public ModuleESolver::ESolver_FP
2727
{
2828
public:
29-
/// @brief a move constructor from ESolver_KS_LCAO
30-
ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol, const Input_para& inp, UnitCell& ucell);
31-
/// @brief a from-scratch constructor
32-
ESolver_LR(const Input_para& inp, UnitCell& ucell);
29+
explicit ESolver_LR(const Input_para& inp);
3330
~ESolver_LR() {
3431
delete this->psi_ks;
3532
}
3633

3734
///input: input, call, basis(LCAO), psi(ground state), elecstate
3835
// initialize sth. independent of the ground state
39-
virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override {};
36+
virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override;
4037
virtual void runner(UnitCell& ucell, int istep) override;
4138
virtual void after_all_runners(UnitCell& ucell) override;
4239

@@ -46,7 +43,7 @@ namespace LR
4643

4744
protected:
4845
const Input_para& input;
49-
const UnitCell& ucell;
46+
const UnitCell* ucell_ = nullptr;
5047
Grid_Driver gd;
5148
std::vector<double> orb_cutoff_;
5249

@@ -87,6 +84,11 @@ namespace LR
8784
bool openshell = false;
8885
std::string xc_kernel;
8986

87+
void initialize_from_unitcell_(UnitCell& ucell, const Input_para& inp);
88+
void initialize_from_ks_(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol,
89+
UnitCell& ucell,
90+
const Input_para& inp);
91+
9092
std::unique_ptr<ModuleGint::GintInfo> gint_info_ = nullptr;
9193
void set_gint();
9294

‎source/source_main/driver_run.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -64,7 +64,7 @@ void Driver::driver_run()
6464
//! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`)
6565
this->init_hardware();
6666

67-
ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell);
67+
ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp);
6868

6969
//! 3: initialize Esolver and fill json-structure
7070
p_esolver->before_all_runners(ucell, PARAM.inp);

0 commit comments

Comments
 (0)