Skip to content

Commit 203c66d

Browse files
YuLiu98YuLiu98
authored
Refactor: remove set_matrix_grid (#5558)
1 parent a6d0ba1 commit 203c66d

File tree

5 files changed

+277
-154
lines changed

5 files changed

+277
-154
lines changed

‎source/Makefile.Objects‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -251,7 +251,6 @@ OBJS_ESOLVER=esolver.o\
251251
OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
252252
esolver_ks_lcao_tddft.o\
253253
dpks_cal_e_delta_band.o\
254-
set_matrix_grid.o\
255254
lcao_before_scf.o\
256255
esolver_gets.o\
257256
lcao_others.o\

‎source/module_esolver/CMakeLists.txt‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,6 @@ if(ENABLE_LCAO)
1818
esolver_ks_lcao.cpp
1919
esolver_ks_lcao_tddft.cpp
2020
dpks_cal_e_delta_band.cpp
21-
set_matrix_grid.cpp
2221
lcao_before_scf.cpp
2322
esolver_gets.cpp
2423
lcao_others.cpp

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 97 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,94 @@ namespace ModuleESolver
3939
{
4040

4141
template <typename TK, typename TR>
42-
void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
42+
void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
4343
{
44-
ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver");
45-
ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver");
44+
ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");
45+
46+
//! 1) call before_scf() of ESolver_FP
47+
ESolver_FP::before_scf(istep);
48+
49+
if (GlobalC::ucell.ionic_position_updated)
50+
{
51+
this->CE.update_all_dis(GlobalC::ucell);
52+
this->CE.extrapolate_charge(
53+
#ifdef __MPI
54+
&(GlobalC::Pgrid),
55+
#endif
56+
GlobalC::ucell,
57+
this->pelec->charge,
58+
&(this->sf),
59+
GlobalV::ofs_running,
60+
GlobalV::ofs_warning);
61+
}
62+
63+
//----------------------------------------------------------
64+
// about vdw, jiyy add vdwd3 and linpz add vdwd2
65+
//----------------------------------------------------------
66+
auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
67+
if (vdw_solver != nullptr)
68+
{
69+
this->pelec->f_en.evdw = vdw_solver->get_energy();
70+
}
4671

4772
// 1. prepare HS matrices, prepare grid integral
48-
this->set_matrix_grid(this->RA);
73+
// (1) Find adjacent atoms for each atom.
74+
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
75+
PARAM.inp.out_level,
76+
orb_.get_rcutmax_Phi(),
77+
GlobalC::ucell.infoNL.get_rcutmax_Beta(),
78+
PARAM.globalv.gamma_only_local);
79+
80+
atom_arrange::search(PARAM.inp.search_pbc,
81+
GlobalV::ofs_running,
82+
GlobalC::GridD,
83+
GlobalC::ucell,
84+
search_radius,
85+
PARAM.inp.test_atom_input);
86+
87+
// (3) Periodic condition search for each grid.
88+
double dr_uniform = 0.001;
89+
std::vector<double> rcuts;
90+
std::vector<std::vector<double>> psi_u;
91+
std::vector<std::vector<double>> dpsi_u;
92+
std::vector<std::vector<double>> d2psi_u;
93+
94+
Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u);
95+
96+
this->GridT.set_pbc_grid(this->pw_rho->nx,
97+
this->pw_rho->ny,
98+
this->pw_rho->nz,
99+
this->pw_big->bx,
100+
this->pw_big->by,
101+
this->pw_big->bz,
102+
this->pw_big->nbx,
103+
this->pw_big->nby,
104+
this->pw_big->nbz,
105+
this->pw_big->nbxx,
106+
this->pw_big->nbzp_start,
107+
this->pw_big->nbzp,
108+
this->pw_rho->ny,
109+
this->pw_rho->nplane,
110+
this->pw_rho->startz_current,
111+
GlobalC::ucell,
112+
GlobalC::GridD,
113+
dr_uniform,
114+
rcuts,
115+
psi_u,
116+
dpsi_u,
117+
d2psi_u,
118+
PARAM.inp.nstream);
119+
psi_u.clear();
120+
psi_u.shrink_to_fit();
121+
dpsi_u.clear();
122+
dpsi_u.shrink_to_fit();
123+
d2psi_u.clear();
124+
d2psi_u.shrink_to_fit();
125+
126+
// (2)For each atom, calculate the adjacent atoms in different cells
127+
// and allocate the space for H(R) and S(R).
128+
// If k point is used here, allocate HlocR after atom_arrange.
129+
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
49130

50131
// 2. density matrix extrapolation
51132

@@ -63,12 +144,10 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
63144
{
64145
nsk = PARAM.inp.nspin;
65146
ncol = this->pv.ncol_bands;
66-
if (PARAM.inp.ks_solver == "genelpa"
67-
|| PARAM.inp.ks_solver == "elpa"
68-
|| PARAM.inp.ks_solver == "lapack"
69-
|| PARAM.inp.ks_solver == "pexsi"
70-
|| PARAM.inp.ks_solver == "cusolver"
71-
|| PARAM.inp.ks_solver == "cusolvermp") {
147+
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
148+
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
149+
|| PARAM.inp.ks_solver == "cusolvermp")
150+
{
72151
ncol = this->pv.ncol;
73152
}
74153
}
@@ -85,12 +164,11 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
85164
}
86165

87166
// init wfc from file
88-
if(istep == 0 && PARAM.inp.init_wfc == "file")
167+
if (istep == 0 && PARAM.inp.init_wfc == "file")
89168
{
90-
if (! ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
169+
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
91170
{
92-
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver",
93-
"read wfc nao failed");
171+
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver", "read wfc nao failed");
94172
}
95173
}
96174

@@ -116,10 +194,11 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
116194
orb_,
117195
DM
118196
#ifdef __EXX
119-
, istep
120-
, GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step
121-
, GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr
122-
, GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
197+
,
198+
istep,
199+
GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step,
200+
GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr,
201+
GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
123202
#endif
124203
);
125204
}
@@ -169,41 +248,6 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
169248
{
170249
GlobalC::ucell.cal_ux();
171250
}
172-
ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver");
173-
}
174-
175-
template <typename TK, typename TR>
176-
void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
177-
{
178-
ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");
179-
180-
//! 1) call before_scf() of ESolver_FP
181-
ESolver_FP::before_scf(istep);
182-
183-
if (GlobalC::ucell.ionic_position_updated)
184-
{
185-
this->CE.update_all_dis(GlobalC::ucell);
186-
this->CE.extrapolate_charge(
187-
#ifdef __MPI
188-
&(GlobalC::Pgrid),
189-
#endif
190-
GlobalC::ucell,
191-
this->pelec->charge,
192-
&(this->sf),
193-
GlobalV::ofs_running,
194-
GlobalV::ofs_warning);
195-
}
196-
197-
//----------------------------------------------------------
198-
// about vdw, jiyy add vdwd3 and linpz add vdwd2
199-
//----------------------------------------------------------
200-
auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running));
201-
if (vdw_solver != nullptr)
202-
{
203-
this->pelec->f_en.evdw = vdw_solver->get_energy();
204-
}
205-
206-
this->beforesolver(istep);
207251

208252
// Peize Lin add 2016-12-03
209253
#ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf

‎source/module_esolver/lcao_others.cpp‎

Lines changed: 180 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,186 @@ void ESolver_KS_LCAO<TK, TR>::others(const int istep)
7070
return;
7171
}
7272

73-
this->beforesolver(istep);
73+
// 1. prepare HS matrices, prepare grid integral
74+
// (1) Find adjacent atoms for each atom.
75+
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
76+
PARAM.inp.out_level,
77+
orb_.get_rcutmax_Phi(),
78+
GlobalC::ucell.infoNL.get_rcutmax_Beta(),
79+
PARAM.globalv.gamma_only_local);
80+
81+
atom_arrange::search(PARAM.inp.search_pbc,
82+
GlobalV::ofs_running,
83+
GlobalC::GridD,
84+
GlobalC::ucell,
85+
search_radius,
86+
PARAM.inp.test_atom_input);
87+
88+
// (3) Periodic condition search for each grid.
89+
double dr_uniform = 0.001;
90+
std::vector<double> rcuts;
91+
std::vector<std::vector<double>> psi_u;
92+
std::vector<std::vector<double>> dpsi_u;
93+
std::vector<std::vector<double>> d2psi_u;
94+
95+
Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u);
96+
97+
this->GridT.set_pbc_grid(this->pw_rho->nx,
98+
this->pw_rho->ny,
99+
this->pw_rho->nz,
100+
this->pw_big->bx,
101+
this->pw_big->by,
102+
this->pw_big->bz,
103+
this->pw_big->nbx,
104+
this->pw_big->nby,
105+
this->pw_big->nbz,
106+
this->pw_big->nbxx,
107+
this->pw_big->nbzp_start,
108+
this->pw_big->nbzp,
109+
this->pw_rho->ny,
110+
this->pw_rho->nplane,
111+
this->pw_rho->startz_current,
112+
GlobalC::ucell,
113+
GlobalC::GridD,
114+
dr_uniform,
115+
rcuts,
116+
psi_u,
117+
dpsi_u,
118+
d2psi_u,
119+
PARAM.inp.nstream);
120+
psi_u.clear();
121+
psi_u.shrink_to_fit();
122+
dpsi_u.clear();
123+
dpsi_u.shrink_to_fit();
124+
d2psi_u.clear();
125+
d2psi_u.shrink_to_fit();
126+
127+
// (2)For each atom, calculate the adjacent atoms in different cells
128+
// and allocate the space for H(R) and S(R).
129+
// If k point is used here, allocate HlocR after atom_arrange.
130+
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
131+
132+
// 2. density matrix extrapolation
133+
134+
// set the augmented orbitals index.
135+
// after ParaO and GridT,
136+
// this information is used to calculate
137+
// the force.
138+
139+
// init psi
140+
if (this->psi == nullptr)
141+
{
142+
int nsk = 0;
143+
int ncol = 0;
144+
if (PARAM.globalv.gamma_only_local)
145+
{
146+
nsk = PARAM.inp.nspin;
147+
ncol = this->pv.ncol_bands;
148+
if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack"
149+
|| PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver"
150+
|| PARAM.inp.ks_solver == "cusolvermp")
151+
{
152+
ncol = this->pv.ncol;
153+
}
154+
}
155+
else
156+
{
157+
nsk = this->kv.get_nks();
158+
#ifdef __MPI
159+
ncol = this->pv.ncol_bands;
160+
#else
161+
ncol = PARAM.inp.nbands;
162+
#endif
163+
}
164+
this->psi = new psi::Psi<TK>(nsk, ncol, this->pv.nrow, nullptr);
165+
}
166+
167+
// init wfc from file
168+
if (istep == 0 && PARAM.inp.init_wfc == "file")
169+
{
170+
if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec))
171+
{
172+
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO<TK, TR>::beforesolver", "read wfc nao failed");
173+
}
174+
}
175+
176+
// prepare grid in Gint
177+
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);
178+
179+
// init Hamiltonian
180+
if (this->p_hamilt != nullptr)
181+
{
182+
delete this->p_hamilt;
183+
this->p_hamilt = nullptr;
184+
}
185+
if (this->p_hamilt == nullptr)
186+
{
187+
elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
188+
this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
189+
PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
190+
PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
191+
&this->pv,
192+
this->pelec->pot,
193+
this->kv,
194+
two_center_bundle_,
195+
orb_,
196+
DM
197+
#ifdef __EXX
198+
,
199+
istep,
200+
GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step,
201+
GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr,
202+
GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
203+
#endif
204+
);
205+
}
206+
207+
#ifdef __DEEPKS
208+
// for each ionic step, the overlap <psi|alpha> must be rebuilt
209+
// since it depends on ionic positions
210+
if (PARAM.globalv.deepks_setorb)
211+
{
212+
const Parallel_Orbitals* pv = &this->pv;
213+
// build and save <psi(0)|alpha(R)> at beginning
214+
GlobalC::ld.build_psialpha(PARAM.inp.cal_force,
215+
GlobalC::ucell,
216+
orb_,
217+
GlobalC::GridD,
218+
*(two_center_bundle_.overlap_orb_alpha));
219+
220+
if (PARAM.inp.deepks_out_unittest)
221+
{
222+
GlobalC::ld.check_psialpha(PARAM.inp.cal_force, GlobalC::ucell, orb_, GlobalC::GridD);
223+
}
224+
}
225+
#endif
226+
if (PARAM.inp.sc_mag_switch)
227+
{
228+
spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
229+
sc.init_sc(PARAM.inp.sc_thr,
230+
PARAM.inp.nsc,
231+
PARAM.inp.nsc_min,
232+
PARAM.inp.alpha_trial,
233+
PARAM.inp.sccut,
234+
PARAM.inp.sc_drop_thr,
235+
GlobalC::ucell,
236+
&(this->pv),
237+
PARAM.inp.nspin,
238+
this->kv,
239+
PARAM.inp.ks_solver,
240+
this->p_hamilt,
241+
this->psi,
242+
this->pelec);
243+
}
244+
//=========================================================
245+
// cal_ux should be called before init_scf because
246+
// the direction of ux is used in noncoline_rho
247+
//=========================================================
248+
if (PARAM.inp.nspin == 4)
249+
{
250+
GlobalC::ucell.cal_ux();
251+
}
252+
74253
// pelec should be initialized before these calculations
75254
this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm);
76255
// self consistent calculations for electronic ground state

0 commit comments

Comments
 (0)