Skip to content

Commit 5b3f106

Browse files
YuLiu98YuLiu98
committed
Refactor: before_scf of esolver_fp
1 parent 1fa5e3a commit 5b3f106

File tree

4 files changed

+86
-182
lines changed

4 files changed

+86
-182
lines changed

‎source/module_esolver/esolver_fp.cpp‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,11 @@
11
#include "esolver_fp.h"
22

33
#include "module_base/global_variable.h"
4+
#include "module_elecstate/cal_ux.h"
45
#include "module_elecstate/module_charge/symmetry_rho.h"
56
#include "module_elecstate/read_pseudo.h"
7+
#include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
8+
#include "module_hamilt_general/module_vdw/vdw.h"
69
#include "module_hamilt_pw/hamilt_pwdft/global.h"
710
#include "module_io/cif_io.h"
811
#include "module_io/cube_io.h"
@@ -283,6 +286,88 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
283286
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
284287
}
285288

289+
// charge extrapolation
290+
if (ucell.ionic_position_updated)
291+
{
292+
this->CE.update_all_dis(ucell);
293+
this->CE.extrapolate_charge(&(this->Pgrid),
294+
ucell,
295+
this->pelec->charge,
296+
&(this->sf),
297+
GlobalV::ofs_running,
298+
GlobalV::ofs_warning);
299+
}
300+
301+
//----------------------------------------------------------
302+
// about vdw, jiyy add vdwd3 and linpz add vdwd2
303+
//----------------------------------------------------------
304+
auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
305+
if (vdw_solver != nullptr)
306+
{
307+
this->pelec->f_en.evdw = vdw_solver->get_energy();
308+
}
309+
310+
// calculate ewald energy
311+
if (!PARAM.inp.test_skip_ewald)
312+
{
313+
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
314+
}
315+
316+
//=========================================================
317+
// cal_ux should be called before init_scf because
318+
// the direction of ux is used in noncoline_rho
319+
//=========================================================
320+
elecstate::cal_ux(ucell);
321+
322+
this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
323+
324+
//! Symmetry_rho should behind init_scf, because charge should be
325+
//! initialized first. liuyu comment: Symmetry_rho should be located between
326+
//! init_rho and v_of_rho?
327+
Symmetry_rho srho;
328+
for (int is = 0; is < PARAM.inp.nspin; is++)
329+
{
330+
srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
331+
}
332+
333+
//! output the initial charge density
334+
if (PARAM.inp.out_chg[0] == 2)
335+
{
336+
for (int is = 0; is < PARAM.inp.nspin; is++)
337+
{
338+
std::stringstream ss;
339+
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
340+
ModuleIO::write_vdata_palgrid(this->Pgrid,
341+
this->pelec->charge->rho[is],
342+
is,
343+
PARAM.inp.nspin,
344+
istep,
345+
ss.str(),
346+
this->pelec->eferm.ef,
347+
&(ucell));
348+
}
349+
}
350+
351+
//! output total local potential of the initial charge density
352+
if (PARAM.inp.out_pot == 3)
353+
{
354+
for (int is = 0; is < PARAM.inp.nspin; is++)
355+
{
356+
std::stringstream ss;
357+
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube";
358+
ModuleIO::write_vdata_palgrid(this->Pgrid,
359+
this->pelec->pot->get_effective_v(is),
360+
is,
361+
PARAM.inp.nspin,
362+
istep,
363+
ss.str(),
364+
0.0, // efermi
365+
&(ucell),
366+
11, // precsion
367+
0); // out_fermi
368+
}
369+
}
370+
286371
return;
287372
}
288373

‎source/module_esolver/esolver_ks_pw.cpp‎

Lines changed: 1 addition & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,6 @@
66
#include "module_elecstate/cal_ux.h"
77
#include "module_elecstate/module_charge/symmetry_rho.h"
88
#include "module_elecstate/occupy.h"
9-
#include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
109
#include "module_hamilt_pw/hamilt_pwdft/global.h"
1110
#include "module_io/print_info.h"
1211
//-----force-------------------
@@ -16,17 +15,16 @@
1615
//---------------------------------------------------
1716
#include "module_base/formatter.h"
1817
#include "module_base/global_variable.h"
18+
#include "module_base/kernels/math_kernel_op.h"
1919
#include "module_base/memory.h"
2020
#include "module_base/module_device/device.h"
2121
#include "module_elecstate/elecstate_pw.h"
2222
#include "module_elecstate/elecstate_pw_sdft.h"
23-
#include "module_hamilt_general/module_vdw/vdw.h"
2423
#include "module_hamilt_pw/hamilt_pwdft/elecond.h"
2524
#include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h"
2625
#include "module_hsolver/diago_iter_assist.h"
2726
#include "module_hsolver/hsolver_pw.h"
2827
#include "module_hsolver/kernels/dngvd_op.h"
29-
#include "module_base/kernels/math_kernel_op.h"
3028
#include "module_io/berryphase.h"
3129
#include "module_io/cube_io.h"
3230
#include "module_io/get_pchg_pw.h"
@@ -254,16 +252,6 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
254252

255253
this->p_psi_init->prepare_init(PARAM.inp.pw_seed);
256254
}
257-
if (ucell.ionic_position_updated)
258-
{
259-
this->CE.update_all_dis(ucell);
260-
this->CE.extrapolate_charge(&this->Pgrid,
261-
ucell,
262-
this->pelec->charge,
263-
&this->sf,
264-
GlobalV::ofs_running,
265-
GlobalV::ofs_warning);
266-
}
267255

268256
// init Hamilt, this should be allocated before each scf loop
269257
// Operators in HamiltPW should be reallocated once cell changed
@@ -273,76 +261,6 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
273261
// allocate HamiltPW
274262
this->allocate_hamilt(ucell);
275263

276-
//----------------------------------------------------------
277-
// about vdw, jiyy add vdwd3 and linpz add vdwd2
278-
//----------------------------------------------------------
279-
auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
280-
if (vdw_solver != nullptr)
281-
{
282-
this->pelec->f_en.evdw = vdw_solver->get_energy();
283-
}
284-
285-
// calculate ewald energy
286-
if (!PARAM.inp.test_skip_ewald)
287-
{
288-
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
289-
}
290-
291-
//! cal_ux should be called before init_scf because
292-
//! the direction of ux is used in noncoline_rho
293-
elecstate::cal_ux(ucell);
294-
295-
//! calculate the total local pseudopotential in real space
296-
this->pelec
297-
->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm, (void*)this->pw_wfc);
298-
299-
//! output the initial charge density
300-
if (PARAM.inp.out_chg[0] == 2)
301-
{
302-
for (int is = 0; is < PARAM.inp.nspin; is++)
303-
{
304-
std::stringstream ss;
305-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
306-
ModuleIO::write_vdata_palgrid(this->Pgrid,
307-
this->pelec->charge->rho[is],
308-
is,
309-
PARAM.inp.nspin,
310-
istep,
311-
ss.str(),
312-
this->pelec->eferm.ef,
313-
&(ucell));
314-
}
315-
}
316-
317-
//! output total local potential of the initial charge density
318-
if (PARAM.inp.out_pot == 3)
319-
{
320-
for (int is = 0; is < PARAM.inp.nspin; is++)
321-
{
322-
std::stringstream ss;
323-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube";
324-
ModuleIO::write_vdata_palgrid(this->Pgrid,
325-
this->pelec->pot->get_effective_v(is),
326-
is,
327-
PARAM.inp.nspin,
328-
istep,
329-
ss.str(),
330-
0.0, // efermi
331-
&(ucell),
332-
11, // precsion
333-
0); // out_fermi
334-
}
335-
}
336-
337-
//! Symmetry_rho should behind init_scf, because charge should be
338-
//! initialized first. liuyu comment: Symmetry_rho should be located between
339-
//! init_rho and v_of_rho?
340-
Symmetry_rho srho;
341-
for (int is = 0; is < PARAM.inp.nspin; is++)
342-
{
343-
srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
344-
}
345-
346264
// liuyu move here 2023-10-09
347265
// D in uspp need vloc, thus behind init_scf()
348266
// calculate the effective coefficient matrix for non-local pseudopotential

‎source/module_esolver/esolver_of.cpp‎

Lines changed: 0 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -254,22 +254,6 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
254254
this->precip_dir_[is] = new std::complex<double>[pw_rho->npw];
255255
}
256256
}
257-
if (ucell.ionic_position_updated)
258-
{
259-
CE.update_all_dis(ucell);
260-
CE.extrapolate_charge(&Pgrid, ucell, pelec->charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
261-
}
262-
263-
this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
264-
265-
// calculate ewald energy
266-
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac);
267-
268-
Symmetry_rho srho;
269-
for (int is = 0; is < PARAM.inp.nspin; is++)
270-
{
271-
srho.begin(is, *(pelec->charge), this->pw_rho, ucell.symm);
272-
}
273257

274258
for (int is = 0; is < PARAM.inp.nspin; ++is)
275259
{

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 0 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,3 @@
1-
#include "module_elecstate/cal_ux.h"
21
#include "module_elecstate/module_charge/symmetry_rho.h"
32
#include "module_esolver/esolver_ks_lcao.h"
43
#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
@@ -22,8 +21,6 @@
2221
#include "module_base/formatter.h"
2322
#include "module_elecstate/elecstate_lcao.h"
2423
#include "module_elecstate/module_dm/cal_dm_psi.h"
25-
#include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
26-
#include "module_hamilt_general/module_vdw/vdw.h"
2724
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
2825
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h"
2926
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
@@ -48,26 +45,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
4845
//! 1) call before_scf() of ESolver_KS
4946
ESolver_KS<TK>::before_scf(ucell, istep);
5047

51-
if (ucell.ionic_position_updated)
52-
{
53-
this->CE.update_all_dis(ucell);
54-
this->CE.extrapolate_charge(&(this->Pgrid),
55-
ucell,
56-
this->pelec->charge,
57-
&(this->sf),
58-
GlobalV::ofs_running,
59-
GlobalV::ofs_warning);
60-
}
61-
62-
//----------------------------------------------------------
63-
// about vdw, jiyy add vdwd3 and linpz add vdwd2
64-
//----------------------------------------------------------
65-
auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
66-
if (vdw_solver != nullptr)
67-
{
68-
this->pelec->f_en.evdw = vdw_solver->get_energy();
69-
}
70-
7148
// 1. prepare HS matrices, prepare grid integral
7249
// (1) Find adjacent atoms for each atom.
7350
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
@@ -274,11 +251,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
274251
this->psi,
275252
this->pelec);
276253
}
277-
//=========================================================
278-
// cal_ux should be called before init_scf because
279-
// the direction of ux is used in noncoline_rho
280-
//=========================================================
281-
elecstate::cal_ux(ucell);
282254

283255
// Peize Lin add 2016-12-03
284256
#ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf
@@ -295,46 +267,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
295267
}
296268
#endif // __EXX
297269

298-
this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
299-
300-
//! output the initial charge density
301-
if (PARAM.inp.out_chg[0] == 2)
302-
{
303-
for (int is = 0; is < PARAM.inp.nspin; is++)
304-
{
305-
std::stringstream ss;
306-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
307-
ModuleIO::write_vdata_palgrid(this->Pgrid,
308-
this->pelec->charge->rho[is],
309-
is,
310-
PARAM.inp.nspin,
311-
istep,
312-
ss.str(),
313-
this->pelec->eferm.ef,
314-
&(ucell));
315-
}
316-
}
317-
318-
//! output total local potential of the initial charge density
319-
if (PARAM.inp.out_pot == 3)
320-
{
321-
for (int is = 0; is < PARAM.inp.nspin; is++)
322-
{
323-
std::stringstream ss;
324-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube";
325-
ModuleIO::write_vdata_palgrid(this->Pgrid,
326-
this->pelec->pot->get_effective_v(is),
327-
is,
328-
PARAM.inp.nspin,
329-
istep,
330-
ss.str(),
331-
0.0, // efermi
332-
&(ucell),
333-
11, // precsion
334-
0); // out_fermi
335-
}
336-
}
337-
338270
// initalize DMR
339271
// DMR should be same size with Hamiltonian(R)
340272
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
@@ -383,21 +315,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
383315
return;
384316
}
385317

386-
// the electron charge density should be symmetrized,
387-
// here is the initialization
388-
Symmetry_rho srho;
389-
for (int is = 0; is < PARAM.inp.nspin; is++)
390-
{
391-
srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm);
392-
}
393-
394-
// 1. calculate ewald energy.
395-
// mohan update 2021-02-25
396-
if (!PARAM.inp.test_skip_ewald)
397-
{
398-
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, this->sf.strucFac);
399-
}
400-
401318
this->p_hamilt->non_first_scf = istep;
402319

403320
// update in ion-step

0 commit comments

Comments
 (0)