Skip to content

Commit 53d83b2

Browse files
YuLiu98YuLiu98
authored
Refactor: before_scf of esolver_fp (#5943)
* Refactor: before_scf of esolver_fp * reset tool_title.cpp
1 parent 83f966a commit 53d83b2

File tree

4 files changed

+106
-193
lines changed

4 files changed

+106
-193
lines changed

‎source/module_esolver/esolver_fp.cpp‎

Lines changed: 74 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,77 @@ 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+
//! output the initial charge density
323+
if (PARAM.inp.out_chg[0] == 2)
324+
{
325+
for (int is = 0; is < PARAM.inp.nspin; is++)
326+
{
327+
std::stringstream ss;
328+
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
329+
ModuleIO::write_vdata_palgrid(this->Pgrid,
330+
this->pelec->charge->rho[is],
331+
is,
332+
PARAM.inp.nspin,
333+
istep,
334+
ss.str(),
335+
this->pelec->eferm.ef,
336+
&(ucell));
337+
}
338+
}
339+
340+
//! output total local potential of the initial charge density
341+
if (PARAM.inp.out_pot == 3)
342+
{
343+
for (int is = 0; is < PARAM.inp.nspin; is++)
344+
{
345+
std::stringstream ss;
346+
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube";
347+
ModuleIO::write_vdata_palgrid(this->Pgrid,
348+
this->pelec->pot->get_effective_v(is),
349+
is,
350+
PARAM.inp.nspin,
351+
istep,
352+
ss.str(),
353+
0.0, // efermi
354+
&(ucell),
355+
11, // precsion
356+
0); // out_fermi
357+
}
358+
}
359+
286360
return;
287361
}
288362

‎source/module_esolver/esolver_ks_pw.cpp‎

Lines changed: 32 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -6,25 +6,23 @@
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"
10-
#include "module_hamilt_pw/hamilt_pwdft/global.h"
11-
#include "module_io/print_info.h"
129
#include "module_hamilt_pw/hamilt_pwdft/forces.h"
10+
#include "module_hamilt_pw/hamilt_pwdft/global.h"
1311
#include "module_hamilt_pw/hamilt_pwdft/stress_pw.h"
12+
#include "module_io/print_info.h"
1413
//---------------------------------------------------
1514
#include "module_base/formatter.h"
1615
#include "module_base/global_variable.h"
16+
#include "module_base/kernels/math_kernel_op.h"
1717
#include "module_base/memory.h"
1818
#include "module_base/module_device/device.h"
1919
#include "module_elecstate/elecstate_pw.h"
2020
#include "module_elecstate/elecstate_pw_sdft.h"
21-
#include "module_hamilt_general/module_vdw/vdw.h"
2221
#include "module_hamilt_pw/hamilt_pwdft/elecond.h"
2322
#include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h"
2423
#include "module_hsolver/diago_iter_assist.h"
2524
#include "module_hsolver/hsolver_pw.h"
2625
#include "module_hsolver/kernels/dngvd_op.h"
27-
#include "module_base/kernels/math_kernel_op.h"
2826
#include "module_io/berryphase.h"
2927
#include "module_io/cube_io.h"
3028
#include "module_io/get_pchg_pw.h"
@@ -109,7 +107,7 @@ ESolver_KS_PW<T, Device>::~ESolver_KS_PW()
109107
std::cout << " ** Closing DSP Hardware..." << std::endl;
110108
mtfunc::dspDestoryHandle(GlobalV::MY_RANK);
111109
#endif
112-
if(PARAM.inp.device == "gpu" || PARAM.inp.precision == "single")
110+
if (PARAM.inp.device == "gpu" || PARAM.inp.precision == "single")
113111
{
114112
delete this->kspw_psi;
115113
}
@@ -191,7 +189,6 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
191189
&(this->pelec->f_en.vtxc));
192190
}
193191

194-
195192
//! initalize local pseudopotential
196193
this->locpp.init_vloc(ucell, this->pw_rhod);
197194
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
@@ -253,16 +250,6 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
253250

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

267254
// init Hamilt, this should be allocated before each scf loop
268255
// Operators in HamiltPW should be reallocated once cell changed
@@ -272,77 +259,12 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
272259
// allocate HamiltPW
273260
this->allocate_hamilt(ucell);
274261

275-
//----------------------------------------------------------
276-
// about vdw, jiyy add vdwd3 and linpz add vdwd2
277-
//----------------------------------------------------------
278-
auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
279-
if (vdw_solver != nullptr)
280-
{
281-
this->pelec->f_en.evdw = vdw_solver->get_energy();
282-
}
283-
284-
//----------------------------------------------------------
285-
// calculate ewald energy
286-
//----------------------------------------------------------
287-
if (!PARAM.inp.test_skip_ewald)
288-
{
289-
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
290-
}
291-
292-
//----------------------------------------------------------
293-
//! cal_ux should be called before init_scf because
294-
//! the direction of ux is used in noncoline_rho
295-
//----------------------------------------------------------
296-
elecstate::cal_ux(ucell);
297-
298262
//----------------------------------------------------------
299263
//! calculate the total local pseudopotential in real space
300264
//----------------------------------------------------------
301265
this->pelec
302266
->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm, (void*)this->pw_wfc);
303267

304-
//----------------------------------------------------------
305-
//! output the initial charge density
306-
//----------------------------------------------------------
307-
if (PARAM.inp.out_chg[0] == 2)
308-
{
309-
for (int is = 0; is < PARAM.inp.nspin; is++)
310-
{
311-
std::stringstream ss;
312-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
313-
ModuleIO::write_vdata_palgrid(this->Pgrid,
314-
this->pelec->charge->rho[is],
315-
is,
316-
PARAM.inp.nspin,
317-
istep,
318-
ss.str(),
319-
this->pelec->eferm.ef,
320-
&(ucell));
321-
}
322-
}
323-
324-
//----------------------------------------------------------
325-
//! output total local potential of the initial charge density
326-
//----------------------------------------------------------
327-
if (PARAM.inp.out_pot == 3)
328-
{
329-
for (int is = 0; is < PARAM.inp.nspin; is++)
330-
{
331-
std::stringstream ss;
332-
ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube";
333-
ModuleIO::write_vdata_palgrid(this->Pgrid,
334-
this->pelec->pot->get_effective_v(is),
335-
is,
336-
PARAM.inp.nspin,
337-
istep,
338-
ss.str(),
339-
0.0, // efermi
340-
&(ucell),
341-
11, // precsion
342-
0); // out_fermi
343-
}
344-
}
345-
346268
//----------------------------------------------------------
347269
//! Symmetry_rho should behind init_scf, because charge should be
348270
//! initialized first. liuyu comment: Symmetry_rho should be located between
@@ -593,7 +515,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
593515

594516
// 2) Update USPP-related quantities
595517
// D in USPP needs vloc, thus needs update when veff updated
596-
// calculate the effective coefficient matrix for non-local pp projectors
518+
// calculate the effective coefficient matrix for non-local pp projectors
597519
// liuyu 2023-10-24
598520
if (PARAM.globalv.use_uspp)
599521
{
@@ -604,9 +526,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
604526
// 3) Print out electronic wavefunctions in pw basis
605527
if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2)
606528
{
607-
if (iter % PARAM.inp.out_freq_elec == 0 ||
608-
iter == PARAM.inp.scf_nmax ||
609-
conv_esolver)
529+
if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
610530
{
611531
std::stringstream ssw;
612532
ssw << PARAM.globalv.global_out_dir << "WAVEFUNC";
@@ -941,35 +861,35 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
941861

942862
#ifdef __MLKEDF
943863
// generate training data for ML-KEDF
944-
if(PARAM.inp.of_ml_gene_data == 1)
864+
if (PARAM.inp.of_ml_gene_data == 1)
945865
{
946866
this->pelec->pot->update_from_charge(this->pelec->charge, &ucell);
947867

948-
ML_data ml_data;
949-
ml_data.set_para(this->pelec->charge->nrxx,
950-
PARAM.inp.nelec,
951-
PARAM.inp.of_tf_weight,
952-
PARAM.inp.of_vw_weight,
953-
PARAM.inp.of_ml_chi_p,
954-
PARAM.inp.of_ml_chi_q,
955-
PARAM.inp.of_ml_chi_xi,
956-
PARAM.inp.of_ml_chi_pnl,
957-
PARAM.inp.of_ml_chi_qnl,
958-
PARAM.inp.of_ml_nkernel,
959-
PARAM.inp.of_ml_kernel,
960-
PARAM.inp.of_ml_kernel_scaling,
961-
PARAM.inp.of_ml_yukawa_alpha,
962-
PARAM.inp.of_ml_kernel_file,
963-
ucell.omega,
964-
this->pw_rho);
965-
966-
ml_data.generateTrainData_KS(this->kspw_psi,
967-
this->pelec,
968-
this->pw_wfc,
969-
this->pw_rho,
970-
ucell,
971-
this->pelec->pot->get_effective_v(0));
972-
}
868+
ML_data ml_data;
869+
ml_data.set_para(this->pelec->charge->nrxx,
870+
PARAM.inp.nelec,
871+
PARAM.inp.of_tf_weight,
872+
PARAM.inp.of_vw_weight,
873+
PARAM.inp.of_ml_chi_p,
874+
PARAM.inp.of_ml_chi_q,
875+
PARAM.inp.of_ml_chi_xi,
876+
PARAM.inp.of_ml_chi_pnl,
877+
PARAM.inp.of_ml_chi_qnl,
878+
PARAM.inp.of_ml_nkernel,
879+
PARAM.inp.of_ml_kernel,
880+
PARAM.inp.of_ml_kernel_scaling,
881+
PARAM.inp.of_ml_yukawa_alpha,
882+
PARAM.inp.of_ml_kernel_file,
883+
ucell.omega,
884+
this->pw_rho);
885+
886+
ml_data.generateTrainData_KS(this->kspw_psi,
887+
this->pelec,
888+
this->pw_wfc,
889+
this->pw_rho,
890+
ucell,
891+
this->pelec->pot->get_effective_v(0));
892+
}
973893
#endif
974894
}
975895

‎source/module_esolver/esolver_of.cpp‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -258,17 +258,9 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
258258
this->precip_dir_[is] = new std::complex<double>[pw_rho->npw];
259259
}
260260
}
261-
if (ucell.ionic_position_updated)
262-
{
263-
CE.update_all_dis(ucell);
264-
CE.extrapolate_charge(&Pgrid, ucell, pelec->charge, &sf, GlobalV::ofs_running, GlobalV::ofs_warning);
265-
}
266261

267262
this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
268263

269-
// calculate ewald energy
270-
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac);
271-
272264
Symmetry_rho srho;
273265
for (int is = 0; is < PARAM.inp.nspin; is++)
274266
{

0 commit comments

Comments
 (0)