From 285e9fb9d6dc7db97f2b8431327f38ff39097d54 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 09:35:13 +0800 Subject: [PATCH 1/6] reorganize esolver_ks_lcao.h --- source/module_esolver/esolver_ks_lcao.cpp | 7 +- source/module_esolver/esolver_ks_lcao.h | 42 ++++++---- source/module_esolver/lcao_others.cpp | 2 +- .../hamilt_lcaodft/FORCE_k.cpp | 77 ------------------- 4 files changed, 31 insertions(+), 97 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 665a4fbde5..b9f865b593 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -331,7 +331,6 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for &ucell.symm); // delete RA after cal_force - this->RA.delete_grid(); this->have_force = true; @@ -349,12 +348,16 @@ void ESolver_KS_LCAO::cal_stress(UnitCell& ucell, ModuleBase::matrix& st ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress"); + // if the users do not want to calculate forces but want stress, + // we call cal_force if (!this->have_force) { ModuleBase::matrix fcs; this->cal_force(ucell, fcs); } - stress = this->scs; // copy the stress + + // the 'scs' stress has already been calculated in 'cal_force' + stress = this->scs; this->have_force = false; ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress"); diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 7cc1fd2375..3a7d0d0e4c 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -1,35 +1,48 @@ #ifndef ESOLVER_KS_LCAO_H #define ESOLVER_KS_LCAO_H + #include "esolver_ks.h" #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" + // for grid integration #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_gint/temp_gint/gint.h" #include "module_hamilt_lcao/module_gint/temp_gint/gint_info.h" + +// for DeePKS #ifdef __DEEPKS #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #endif + +// for EXX #ifdef __EXX #include "module_ri/Exx_LRI_interface.h" #include "module_ri/Mix_DMk_2D.h" #endif + #include "module_basis/module_nao/two_center_bundle.h" #include "module_io/output_mat_sparse.h" -// added by jghan for rdmft calculation +// for RDMFT #include "module_rdmft/rdmft.h" - #include + +// for Linear Response namespace LR { template class ESolver_LR; } + +//----------------------------------- +// ESolver for LCAO +//----------------------------------- namespace ModuleESolver { + template class ESolver_KS_LCAO : public ESolver_KS { @@ -62,43 +75,38 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void others(UnitCell& ucell, const int istep) override; - // we will get rid of this class soon, don't use it, mohan 2024-03-28 + //! Store information about Adjacent Atoms Record_adj RA; + //! Store information about Adjacent Atoms Grid_Driver gd; - // 2d block-cyclic distribution info + //! NAO orbitals: 2d block-cyclic distribution info Parallel_Orbitals pv; - // used for k-dependent grid integration. + //! Grid integration: used for k-point-dependent algorithm Gint_k GK; - // used for gamma only algorithms. + //! Grid integration: used for gamma only algorithms. Gint_Gamma GG; + //! Grid integration: used to store some basic information Grid_Technique GridT; + //! NAO orbitals: two-center integrations TwoCenterBundle two_center_bundle_; - rdmft::RDMFT rdmft_solver; // added by jghan for rdmft calculation, 2024-03-16 + //! For RDMFT calculations, added by jghan, 2024-03-16 + rdmft::RDMFT rdmft_solver; - // temporary introduced during removing GlobalC::ORB + //! NAO: store related information LCAO_Orbitals orb_; // Temporarily store the stress to unify the interface with PW, // because it's hard to seperate force and stress calculation in LCAO. - // The copy costs memory and time ! - // Are there any better way to solve this problem ? ModuleBase::matrix scs; bool have_force = false; - //--------------common for all calculation, not only scf------------- - // set matrix and grid integral - void set_matrix_grid(Record_adj& ra); - - void beforesolver(const int istep); - //--------------------------------------------------------------------- - #ifdef __DEEPKS LCAO_Deepks ld; #endif diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index af96e30f6d..2e24e0e02b 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -173,7 +173,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) { if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec)) { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", "read wfc nao failed"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "read wfc nao failed"); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index 2fcd83dda6..3fc8390b7c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -187,83 +187,6 @@ void Force_LCAO>::finish_ftable(ForceStressArrays& fsr) return; } -// template <> -// void Force_LCAO>::test(Parallel_Orbitals& pv, double* mmm, const std::string& name) -//{ -// // mohan remove 'const' for pv, 2024-03-31 -// if (GlobalV::NPROC != 1) -// { -// return; -// } -// -// std::cout << "test!" << std::endl; -// -// int irr = 0; -// int ca = 0; -// -// GlobalV::ofs_running << " Calculate the test in Force_LCAO_k" << std::endl; -// Record_adj RA; -// -// // mohan update 2024-03-31 -// RA.for_2d(pv, GlobalV::GAMMA_ONLY_LOCAL, GlobalC::ORB.cutoffs()); -// -// double* test; -// test = new double[PARAM.globalv.nlocal * PARAM.globalv.nlocal]; -// ModuleBase::GlobalFunc::ZEROS(test, PARAM.globalv.nlocal * PARAM.globalv.nlocal); -// -// for (int T1 = 0; T1 < ucell.ntype; T1++) -// { -// Atom* atom1 = &ucell.atoms[T1]; -// for (int I1 = 0; I1 < atom1->na; I1++) -// { -// // const int iat = ucell.itia2iat(T1,I1); -// const int start1 = ucell.itiaiw2iwt(T1, I1, 0); -// for (int cb = 0; cb < RA.na_each[ca]; cb++) -// { -// const int T2 = RA.info[ca][cb][3]; -// const int I2 = RA.info[ca][cb][4]; -// Atom* atom2 = &ucell.atoms[T2]; -// const int start2 = ucell.itiaiw2iwt(T2, I2, 0); -// -// for (int jj = 0; jj < atom1->nw; jj++) -// { -// const int iw1_all = start1 + jj; -// for (int kk = 0; kk < atom2->nw; kk++) -// { -// const int iw2_all = start2 + kk; -// assert(irr < pv.nnr); -// test[iw1_all * PARAM.globalv.nlocal + iw2_all] += mmm[irr]; -// ++irr; -// } -// } -// } -// ++ca; -// } -// } -// -// std::cout << "\n " << name << std::endl; -// std::cout << std::setprecision(4); -// for (int i = 0; i < PARAM.globalv.nlocal; i++) -// { -// for (int j = 0; j < PARAM.globalv.nlocal; j++) -// { -// if (std::abs(test[i * PARAM.globalv.nlocal + j]) > 1.0e-5) -// { -// std::cout << std::setw(12) << test[i * PARAM.globalv.nlocal + j]; -// } -// else -// { -// std::cout << std::setw(12) << "0"; -// } -// } -// std::cout << std::endl; -// } -// delete[] test; -// -// RA.delete_grid(); // xiaohui add 2015-02-04 -// return; -// } - // be called in Force_LCAO::start_force_calculation template <> void Force_LCAO>::ftable(const bool isforce, From 56b7ca67d152e6845fae706c540babe04fa5ac0e Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 09:59:06 +0800 Subject: [PATCH 2/6] delete some useless files in esolver_ks.h --- source/module_esolver/esolver_ks.cpp | 1 + source/module_esolver/esolver_ks.h | 43 +++++++++++++------------ source/module_esolver/esolver_ks_lcao.h | 8 +++-- source/module_esolver/lcao_others.cpp | 3 ++ source/module_esolver/pw_others.cpp | 3 ++ 5 files changed, 35 insertions(+), 23 deletions(-) diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 55511a1977..9cb971a798 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -10,6 +10,7 @@ #include "module_io/write_istate_info.h" #include "module_parameter/parameter.h" #include "module_elecstate/elecstate_print.h" +#include "module_hsolver/hsolver.h" #include #include diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 652c193d5c..7dd356430a 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -1,14 +1,23 @@ #ifndef ESOLVER_KS_H #define ESOLVER_KS_H + #include "esolver_fp.h" + +// for plane wave basis set #include "module_basis/module_pw/pw_basis_k.h" + +// for k-points in Brillouin zone #include "module_cell/klist.h" + +// for charge mixing #include "module_elecstate/module_charge/charge_mixing.h" -#include "module_hamilt_general/hamilt.h" -#include "module_hsolver/hsolver.h" -#include "module_io/cal_test.h" + +// for electronic wave functions #include "module_psi/psi.h" +// for Hamiltonian +#include "module_hamilt_general/hamilt.h" + #ifdef __MPI #include #else @@ -66,21 +75,12 @@ class ESolver_KS : public ESolver_FP //! Charge mixing method, only used in KDSFT, not in OFDFT Charge_Mixing* p_chgmix = nullptr; - //! nonlocal pseudo potential + //! nonlocal pseudopotentials pseudopot_cell_vnl ppcell; //! Electronic wavefunctions psi::Psi* psi = nullptr; - //! plane wave or LCAO - std::string basisname; - - //! number of electrons - double esolver_KS_ne = 0.0; - - //! whether esolver is oscillated - bool oscillate_esolver = false; - //! the start time of scf iteration #ifdef __MPI double iter_time; @@ -88,13 +88,16 @@ class ESolver_KS : public ESolver_FP std::chrono::system_clock::time_point iter_time; #endif - double diag_ethr; //! the threshold for diagonalization - double scf_thr; //! scf density threshold - double scf_ene_thr; //! scf energy threshold - double drho; //! the difference between rho_in (before HSolver) and rho_out (After HSolver) - double hsolver_error; //! the error of HSolver - int maxniter; //! maximum iter steps for scf - int niter; //! iter steps actually used in scf + std::string basisname; //! esolver_ks_lcao.cpp + double esolver_KS_ne = 0.0; //! number of electrons + double diag_ethr; //! the threshold for diagonalization + double scf_thr; //! scf density threshold + double scf_ene_thr; //! scf energy threshold + double drho; //! the difference between rho_in (before HSolver) and rho_out (After HSolver) + double hsolver_error; //! the error of HSolver + int maxniter; //! maximum iter steps for scf + int niter; //! iter steps actually used in scf + bool oscillate_esolver = false; // whether esolver is oscillated }; } // namespace ModuleESolver #endif diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 3a7d0d0e4c..031400b2cc 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -2,8 +2,13 @@ #define ESOLVER_KS_LCAO_H #include "esolver_ks.h" + +// for adjacent atoms #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" +// for NAO basis +#include "module_basis/module_nao/two_center_bundle.h" + // for grid integration #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" @@ -21,9 +26,6 @@ #include "module_ri/Mix_DMk_2D.h" #endif -#include "module_basis/module_nao/two_center_bundle.h" -#include "module_io/output_mat_sparse.h" - // for RDMFT #include "module_rdmft/rdmft.h" diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 2e24e0e02b..8b211063cb 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -34,6 +34,9 @@ #include "module_io/restart_exx_csr.h" #endif +// mohan add 2025-03-06 +#include "module_io/cal_test.h" + namespace ModuleESolver { diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index 0f2be0a998..c36ddae6e7 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -46,6 +46,9 @@ #include #include "module_base/formatter.h" +// mohan add 2025-03-06 +#include "module_io/cal_test.h" + namespace ModuleESolver { template From b10c6b8da9280a8a6711f23c1c27bf4bdac73a3b Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 10:28:23 +0800 Subject: [PATCH 3/6] change pelec->charge to this->chr --- source/module_esolver/esolver_fp.cpp | 50 +++++++++++++------------- source/module_esolver/esolver_fp.h | 52 +++++++++++++++++----------- source/module_esolver/esolver_ks.h | 2 +- 3 files changed, 58 insertions(+), 46 deletions(-) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 297309bed5..ada41292ca 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -23,15 +23,16 @@ namespace ModuleESolver ESolver_FP::ESolver_FP() { - // pw_rho = new ModuleBase::PW_Basis(); - // LCAO basis doesn't support GPU acceleration on FFT currently std::string fft_device = PARAM.inp.device; + + // LCAO basis doesn't support GPU acceleration on FFT currently if(PARAM.inp.basis_type == "lcao") { fft_device = "cpu"; } + pw_rho = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision); - if ( PARAM.globalv.double_grid) + if (PARAM.globalv.double_grid) { pw_rhod = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision); } @@ -44,6 +45,7 @@ ESolver_FP::ESolver_FP() pw_big = static_cast(pw_rhod); pw_big->setbxyz(PARAM.inp.bx, PARAM.inp.by, PARAM.inp.bz); sf.set(pw_rhod, PARAM.inp.nbspline); + } ESolver_FP::~ESolver_FP() @@ -140,7 +142,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso // 2) write fermi energy ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef); - // 3) update delta rho for charge extrapolation + // 3) update delta_rho for charge extrapolation CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); if (istep % PARAM.inp.out_interval == 0) @@ -153,13 +155,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso double* data = nullptr; if (PARAM.inp.dm_to_rho) { - data = this->pelec->charge->rho[is]; - this->pw_rhod->real2recip(this->pelec->charge->rho[is], this->pelec->charge->rhog[is]); + data = this->chr->rho[is]; + this->pw_rhod->real2recip(this->chr->rho[is], this->chr->rhog[is]); } else { - data = this->pelec->charge->rho_save[is]; - this->pw_rhod->real2recip(this->pelec->charge->rho_save[is], this->pelec->charge->rhog_save[is]); + data = this->chr->rho_save[is]; + this->pw_rhod->real2recip(this->chr->rho_save[is], this->chr->rhog_save[is]); } std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(Pgrid, @@ -176,7 +178,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso { fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(Pgrid, - this->pelec->charge->kin_r_save[is], + this->chr->kin_r_save[is], is, PARAM.inp.nspin, istep, @@ -217,7 +219,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso fn, istep, this->pw_rhod, - this->pelec->charge, + this->chr, &(ucell), this->pelec->pot->get_fixed_v(), this->solvent); @@ -226,11 +228,11 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso // 6) write ELF if (PARAM.inp.out_elf[0] > 0) { - this->pelec->charge->cal_elf = true; + this->chr->cal_elf = true; Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); + srho.begin(is, *(this->chr), this->pw_rhod, ucell.symm); } std::string out_dir =PARAM.globalv.global_out_dir; @@ -242,8 +244,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso out_dir, istep, PARAM.inp.nspin, - this->pelec->charge->rho, - this->pelec->charge->kin_r, + this->chr->rho, + this->chr->kin_r, this->pw_rhod, this->Pgrid, &(ucell), @@ -260,9 +262,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) { // only G-vector and K-vector are changed due to the change of lattice // vector FFT grids do not change!! - pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz); - pw_rho->collect_local_pw(); - pw_rho->collect_uniqgg(); + this->pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz); + this->pw_rho->collect_local_pw(); + this->pw_rho->collect_uniqgg(); if (PARAM.globalv.double_grid) { @@ -292,7 +294,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->CE.update_all_dis(ucell); this->CE.extrapolate_charge(&(this->Pgrid), ucell, - this->pelec->charge, + this->chr, &(this->sf), GlobalV::ofs_running, GlobalV::ofs_warning); @@ -327,7 +329,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; ModuleIO::write_vdata_palgrid(this->Pgrid, - this->pelec->charge->rho[is], + this->chr->rho[is], is, PARAM.inp.nspin, istep, @@ -368,8 +370,8 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) { std::complex** rhog_tot - = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rhog : this->pelec->charge->rhog_save; - double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rho : this->pelec->charge->rho_save; + = (PARAM.inp.dm_to_rho) ? this->chr->rhog : this->chr->rhog_save; + double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr->rho : this->chr->rho_save; for (int is = 0; is < PARAM.inp.nspin; is++) { this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]); @@ -386,12 +388,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& if (XC_Functional::get_ked_flag()) { - std::vector> kin_g_space(PARAM.inp.nspin * this->pelec->charge->ngmc, {0.0, 0.0}); + std::vector> kin_g_space(PARAM.inp.nspin * this->chr->ngmc, {0.0, 0.0}); std::vector*> kin_g; for (int is = 0; is < PARAM.inp.nspin; is++) { - kin_g.push_back(kin_g_space.data() + is * this->pelec->charge->ngmc); - this->pw_rhod->real2recip(this->pelec->charge->kin_r_save[is], kin_g[is]); + kin_g.push_back(kin_g_space.data() + is * this->chr->ngmc); + this->pw_rhod->real2recip(this->chr->kin_r_save[is], kin_g[is]); } ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-TAU-DENSITY.restart", PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local, diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 9c95dd5b2d..1e7f23c02c 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -2,16 +2,31 @@ #define ESOLVER_FP_H #include "esolver.h" + +//! plane wave basis #include "module_basis/module_pw/pw_basis.h" + +//! symmetry analysis #include "module_cell/module_symmetry/symmetry.h" + +//! electronic states #include "module_elecstate/elecstate.h" + +//! charge extrapolation #include "module_elecstate/module_charge/charge_extra.h" + +//! solvation model #include "module_hamilt_general/module_surchem/surchem.h" + +//! local pseudopotential #include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h" + +//! structure factor related to plane wave basis #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include + //! The First-Principles (FP) Energy Solver Class /** * This class represents components that needed in @@ -22,7 +37,7 @@ namespace ModuleESolver { -class ESolver_FP : public ESolver +class ESolver_FP: public ESolver { public: //! Constructor @@ -49,39 +64,34 @@ class ESolver_FP : public ESolver //! ------------------------------------------------------------------------------ elecstate::ElecState* pelec = nullptr; ///< Electronic states - //! ------------------------------------------------------------------------------ + //! K points in Brillouin zone + K_Vectors kv; //! Electorn charge density Charge chr; - //! Structure factors that used with plane-wave basis set - Structure_Factor sf; - - //! K points in Brillouin zone - K_Vectors kv; - - //! Plane-wave basis set for charge density + //! pw_rho: Plane-wave basis set for charge density + //! pw_rhod: same as pw_rho for NCPP. Here 'd' stands for 'dense', + //! dense grid for for uspp, used for ultrasoft augmented charge density. + //! charge density and potential are defined on dense grids, + //! but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi> ModulePW::PW_Basis* pw_rho; + ModulePW::PW_Basis* pw_rhod; //! dense grid for USPP + ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class //! parallel for rho grid Parallel_Grid Pgrid; - //! pointer to local pseudopotential - pseudopot_cell_vl locpp; + //! Structure factors that used with plane-wave basis set + Structure_Factor sf; - /** - * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense' - * dense grid for for uspp, used for ultrasoft augmented charge density. - * charge density and potential are defined on dense grids, - * but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi> - */ - ModulePW::PW_Basis* pw_rhod; - ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class + //! local pseudopotentials + pseudopot_cell_vl locpp; - //! Charge extrapolation + //! charge extrapolation method Charge_Extra CE; - // solvent model + //! solvent model surchem solvent; }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 7dd356430a..75e6ff91bc 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -72,7 +72,7 @@ class ESolver_KS : public ESolver_FP //! PW for wave functions, only used in KSDFT, not in OFDFT ModulePW::PW_Basis_K* pw_wfc = nullptr; - //! Charge mixing method, only used in KDSFT, not in OFDFT + //! Charge mixing method Charge_Mixing* p_chgmix = nullptr; //! nonlocal pseudopotentials From 9f0f6892e0cd034a0001ebf3ae6a14daccc5b9cc Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 10:57:59 +0800 Subject: [PATCH 4/6] will delete pelec->charge in near future, we just directly use chr defined in esolver_fp --- source/module_esolver/esolver_fp.cpp | 34 +++++++++---------- source/module_esolver/esolver_ks.cpp | 20 +++++------ source/module_esolver/esolver_ks_lcao.cpp | 16 ++++----- .../module_esolver/esolver_ks_lcao_tddft.cpp | 8 ++--- source/module_esolver/esolver_ks_lcaopw.cpp | 4 +-- source/module_esolver/esolver_ks_pw.cpp | 14 ++++---- source/module_esolver/esolver_of.cpp | 34 +++++++++---------- .../module_esolver/esolver_of_interface.cpp | 4 +-- source/module_esolver/esolver_of_tool.cpp | 12 +++---- source/module_esolver/esolver_sdft_pw.cpp | 4 +-- source/module_esolver/lcao_before_scf.cpp | 4 +-- source/module_esolver/lcao_others.cpp | 8 ++--- 12 files changed, 81 insertions(+), 81 deletions(-) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index ada41292ca..c92c2ebd34 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -155,13 +155,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso double* data = nullptr; if (PARAM.inp.dm_to_rho) { - data = this->chr->rho[is]; - this->pw_rhod->real2recip(this->chr->rho[is], this->chr->rhog[is]); + data = this->chr.rho[is]; + this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]); } else { - data = this->chr->rho_save[is]; - this->pw_rhod->real2recip(this->chr->rho_save[is], this->chr->rhog_save[is]); + data = this->chr.rho_save[is]; + this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); } std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(Pgrid, @@ -178,7 +178,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso { fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(Pgrid, - this->chr->kin_r_save[is], + this->chr.kin_r_save[is], is, PARAM.inp.nspin, istep, @@ -219,7 +219,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso fn, istep, this->pw_rhod, - this->chr, + &this->chr, &(ucell), this->pelec->pot->get_fixed_v(), this->solvent); @@ -228,11 +228,11 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso // 6) write ELF if (PARAM.inp.out_elf[0] > 0) { - this->chr->cal_elf = true; + this->chr.cal_elf = true; Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->chr), this->pw_rhod, ucell.symm); + srho.begin(is, this->chr, this->pw_rhod, ucell.symm); } std::string out_dir =PARAM.globalv.global_out_dir; @@ -244,8 +244,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso out_dir, istep, PARAM.inp.nspin, - this->chr->rho, - this->chr->kin_r, + this->chr.rho, + this->chr.kin_r, this->pw_rhod, this->Pgrid, &(ucell), @@ -294,7 +294,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->CE.update_all_dis(ucell); this->CE.extrapolate_charge(&(this->Pgrid), ucell, - this->chr, + &this->chr, &(this->sf), GlobalV::ofs_running, GlobalV::ofs_warning); @@ -329,7 +329,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; ModuleIO::write_vdata_palgrid(this->Pgrid, - this->chr->rho[is], + this->chr.rho[is], is, PARAM.inp.nspin, istep, @@ -370,8 +370,8 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) { std::complex** rhog_tot - = (PARAM.inp.dm_to_rho) ? this->chr->rhog : this->chr->rhog_save; - double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr->rho : this->chr->rho_save; + = (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save; + double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save; for (int is = 0; is < PARAM.inp.nspin; is++) { this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]); @@ -388,12 +388,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& if (XC_Functional::get_ked_flag()) { - std::vector> kin_g_space(PARAM.inp.nspin * this->chr->ngmc, {0.0, 0.0}); + std::vector> kin_g_space(PARAM.inp.nspin * this->chr.ngmc, {0.0, 0.0}); std::vector*> kin_g; for (int is = 0; is < PARAM.inp.nspin; is++) { - kin_g.push_back(kin_g_space.data() + is * this->chr->ngmc); - this->pw_rhod->real2recip(this->chr->kin_r_save[is], kin_g[is]); + kin_g.push_back(kin_g_space.data() + is * this->chr.ngmc); + this->pw_rhod->real2recip(this->chr.kin_r_save[is], kin_g[is]); } ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-TAU-DENSITY.restart", PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local, diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 9cb971a798..84c8f2525b 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -368,7 +368,7 @@ void ESolver_KS::hamilt2density(UnitCell& ucell, const int istep, con // double drho = this->estate.caldr2(); // EState should be used after it is constructed. - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec); hsolver_error = 0.0; if (iter == 1 && PARAM.inp.calculation != "nscf") { @@ -390,7 +390,7 @@ void ESolver_KS::hamilt2density(UnitCell& ucell, const int istep, con this->hamilt2density_single(ucell, istep, iter, diag_ethr); - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec); hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, PARAM.inp.esolver_type, @@ -522,7 +522,7 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const in } // 1) save input rho - this->pelec->charge->save_rho_before_sum_band(); + this->chr.save_rho_before_sum_band(); } template @@ -552,9 +552,9 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i // compute magnetization, only for LSDA(spin==2) ucell.magnet.compute_magnetization(ucell.omega, - this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, + this->chr.nrxx, + this->chr.nxyz, + this->chr.rho, this->pelec->nelec_spin.data()); if (PARAM.globalv.ks_run) @@ -629,11 +629,11 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i } else { - p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing + p_chgmix->mix_rho(&this->chr); // update chr->rho by mixing } if (PARAM.inp.scf_thr_type == 2) { - pelec->charge->renormalize_rho(); // renormalize rho in R-space would + this->chr.renormalize_rho(); // renormalize rho in R-space would // induce a error in K-space } //----------charge mixing done----------- @@ -645,7 +645,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i // be careful! conv_esolver is bool, not double !! Maybe a bug 20250302 by mohan MPI_Bcast(&conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD); - MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD); + MPI_Bcast(this->chr.rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD); #endif // update potential @@ -677,7 +677,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i double dkin = 0.0; // for meta-GGA if (XC_Functional::get_ked_flag()) { - dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec); + dkin = p_chgmix->get_dkin(&this->chr, PARAM.inp.nelec); } diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index b9f865b593..a861a5a7d1 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -193,7 +193,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); // 8) inititlize the charge density - this->pelec->charge->allocate(PARAM.inp.nspin); + this->chr.allocate(PARAM.inp.nspin); this->pelec->omega = ucell.omega; // 9) initialize the potential @@ -469,7 +469,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) *this->pw_rho, *this->pw_rhod, this->locpp.vloc, - *this->pelec->charge, + this->chr, this->GG, this->GK, this->kv, @@ -497,7 +497,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) *this->pw_rho, *this->pw_rhod, this->locpp.vloc, - *this->pelec->charge, + this->chr, this->GG, this->GK, this->kv, @@ -628,7 +628,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const elecstate::cal_ux(ucell); //! update the potentials by using new electron charge density - this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); + this->pelec->pot->update_from_charge(&this->chr, &ucell); //! compute the correction energy for metals this->pelec->f_en.descf = this->pelec->cal_delta_escf(); @@ -665,7 +665,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const GlobalC::dftu.set_dmr(dynamic_cast*>(this->pelec)->get_DM()); } // Calculate U and J if Yukawa potential is used - GlobalC::dftu.cal_slater_UJ(ucell, this->pelec->charge->rho, this->pw_rho->nrxx); + GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx); } #ifdef __DEEPKS @@ -763,7 +763,7 @@ void ESolver_KS_LCAO::hamilt2density_single(UnitCell& ucell, int istep, Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); + srho.begin(is, this->chr, this->pw_rho, ucell.symm); } // 12) calculate delta energy @@ -783,7 +783,7 @@ void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const if (!conv_esolver) { elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); + this->pelec->pot->update_from_charge(&this->chr, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } else @@ -866,7 +866,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& { for (int is = 0; is < PARAM.inp.nspin; ++is) { - GlobalC::restart.save_disk("charge", is, this->pelec->charge->nrxx, this->pelec->charge->rho[is]); + GlobalC::restart.save_disk("charge", is, this->chr.nrxx, this->chr.rho[is]); } } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index f8b86fae26..664e666aac 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -154,7 +154,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(UnitCell& ucell, Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(pelec->charge), pw_rho, ucell.symm); + srho.begin(is, this->chr, pw_rho, ucell.symm); } } @@ -206,7 +206,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, if (!conv_esolver) { elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); + this->pelec->pot->update_from_charge(&this->chr, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } else @@ -365,8 +365,8 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep, std::stringstream ss_dipole; ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; ModuleIO::write_dipole(ucell, - pelec->charge->rho_save[is], - pelec->charge->rhopw, + this->chr.rho_save[is], + this->chr.rhopw, is, istep, ss_dipole.str()); diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index 29f6bfb94b..8439265bfd 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -163,7 +163,7 @@ namespace ModuleESolver Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); + srho.begin(is, this->chr, this->pw_rhod, ucell.symm); } // deband is calculated from "output" charge density calculated @@ -255,7 +255,7 @@ namespace ModuleESolver *this->pw_rho, *this->pw_rhod, this->locpp.vloc, - *this->pelec->charge, + this->chr, this->kv, this->pelec->wg #ifdef __EXX diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 0f6d52b6f0..5042f740e7 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -164,7 +164,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p } //! 4) inititlize the charge density. - this->pelec->charge->allocate(PARAM.inp.nspin); + this->chr.allocate(PARAM.inp.nspin); //! 5) set the cell volume variable in pelec this->pelec->omega = ucell.omega; @@ -266,7 +266,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); + srho.begin(is, this->chr, this->pw_rhod, ucell.symm); } //---------------------------------------------------------- @@ -470,7 +470,7 @@ void ESolver_KS_PW::hamilt2density_single(UnitCell& ucell, Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); + srho.begin(is, this->chr, this->pw_rhod, ucell.symm); } // deband is calculated from "output" charge density calculated @@ -488,7 +488,7 @@ void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, cons if (!conv_esolver) { elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); + this->pelec->pot->update_from_charge(&this->chr, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); #ifdef __MPI MPI_Bcast(&(this->pelec->f_en.descf), 1, MPI_DOUBLE, 0, BP_WORLD); @@ -605,7 +605,7 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const this->kv.wk, this->pw_big->bz, this->pw_big->nbz, - this->pelec->charge->ngmc, + this->chr.ngmc, &ucell, this->psi, this->pw_rhod, @@ -876,10 +876,10 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) // generate training data for ML-KEDF if (PARAM.inp.of_ml_gene_data == 1) { - this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); + this->pelec->pot->update_from_charge(this->chr, &ucell); ML_data ml_data; - ml_data.set_para(this->pelec->charge->nrxx, + ml_data.set_para(this->chr.nrxx, PARAM.inp.nelec, PARAM.inp.of_tf_weight, PARAM.inp.of_vw_weight, diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 0ddb31d507..1849bbb0f3 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -161,7 +161,7 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep) #ifdef __MLKEDF // for ML KEDF test - if (PARAM.inp.of_ml_local_test) this->ml_->localTest(pelec->charge->rho, this->pw_rho); + if (PARAM.inp.of_ml_local_test) this->ml_->localTest(this->chr.rho, this->pw_rho); #endif @@ -264,7 +264,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(pelec->charge), this->pw_rho, ucell.symm); + srho.begin(is, this->chr, this->pw_rho, ucell.symm); } for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -275,15 +275,15 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) { // Here we initialize rho to be uniform, // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements. - pelec->charge->rho[is][ibs] = this->nelec_[is] / this->pelec->omega; - this->pphi_[is][ibs] = sqrt(pelec->charge->rho[is][ibs]); + this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega; + this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); } } else { for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) { - this->pphi_[is][ibs] = sqrt(pelec->charge->rho[is][ibs]); + this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); } } } @@ -315,8 +315,8 @@ void ESolver_OF::update_potential(UnitCell& ucell) // (1) get dL/dphi elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(pelec->charge, &ucell); // Hartree + XC + external - this->kinetic_potential(pelec->charge->rho, + this->pelec->pot->update_from_charge(&this->chr, &ucell); // Hartree + XC + external + this->kinetic_potential(this->chr.rho, this->pphi_, this->pelec->pot->get_effective_v()); // (kinetic + Hartree + XC + external) * 2 * phi for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -405,7 +405,7 @@ void ESolver_OF::update_rho() { this->pphi_[is][ir] = this->pphi_[is][ir] * cos(this->theta_[is]) + this->pdirect_[is][ir] * sin(this->theta_[is]); - pelec->charge->rho[is][ir] = this->pphi_[is][ir] * this->pphi_[is][ir]; + this->chr.rho[is][ir] = this->pphi_[is][ir] * this->pphi_[is][ir]; } } // // ------------ turn on symmetry may cause instability in optimization ------------ @@ -414,10 +414,10 @@ void ESolver_OF::update_rho() // Symmetry_rho srho; // for (int is = 0; is < PARAM.inp.nspin; is++) // { - // srho.begin(is, *(pelec->charge), this->pw_rho, Pgrid, ucell.symm); + // srho.begin(is, *(this->chr), this->pw_rho, Pgrid, ucell.symm); // for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) // { - // this->pphi_[is][ibs] = sqrt(pelec->charge->rho[is][ibs]); + // this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]); // } // } // } @@ -494,7 +494,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso //------------------------------------------------------------------ if (PARAM.inp.out_elf[0] > 0) { - this->kinetic_energy_density(this->pelec->charge->rho, this->pphi_, this->pelec->charge->kin_r); + this->kinetic_energy_density(this->chr.rho, this->pphi_, this->chr.kin_r); } //------------------------------------------------------------------ @@ -506,7 +506,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso // should not be here? mohan note 2025-03-03 for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { - this->pelec->charge->rho_save[0][ir] = this->pelec->charge->rho[0][ir]; + this->chr.rho_save[0][ir] = this->chr.rho[0][ir]; } #ifdef __MLKEDF @@ -515,7 +515,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso //------------------------------------------------------------------ if (this->of_kinetic_ == "ml") { - this->tf_->get_energy(this->pelec->charge->rho); + this->tf_->get_energy(this->chr.rho); std::cout << "ML Term = " << this->ml_->ml_energy << " Ry, TF Term = " << this->tf_->tf_energy @@ -532,8 +532,8 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso //------------------------------------------------------------------ if (PARAM.inp.of_ml_gene_data) { - this->pelec->pot->update_from_charge(pelec->charge, &ucell); // Hartree + XC + external - this->kinetic_potential(pelec->charge->rho, this->pphi_, this->pelec->pot->get_effective_v()); // (kinetic + Hartree + XC + external) * 2 * phi + this->pelec->pot->update_from_charge(this->chr, &ucell); // Hartree + XC + external + this->kinetic_potential(this->chr.rho, this->pphi_, this->pelec->pot->get_effective_v()); // (kinetic + Hartree + XC + external) * 2 * phi const double* vr_eff = this->pelec->pot->get_effective_v(0); for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) @@ -547,7 +547,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso // ================= std::cout << "Generating Training data..." << std::endl; std::cout << "mu = " << this->pelec->eferm.get_efval(0) << std::endl; - this->ml_->generateTrainData(pelec->charge->rho, *(this->wt_), *(this->tf_), this->pw_rho, vr_eff); + this->ml_->generateTrainData(this->chr.rho, *(this->wt_), *(this->tf_), this->pw_rho, vr_eff); } #endif @@ -580,7 +580,7 @@ double ESolver_OF::cal_energy() for (int is = 0; is < PARAM.inp.nspin; ++is) { pseudopot_energy += this->inner_product(this->pelec->pot->get_fixed_v(), - pelec->charge->rho[is], + this->chr.rho[is], this->pw_rho->nrxx, this->dV_); } diff --git a/source/module_esolver/esolver_of_interface.cpp b/source/module_esolver/esolver_of_interface.cpp index 057b764118..5df1051af2 100644 --- a/source/module_esolver/esolver_of_interface.cpp +++ b/source/module_esolver/esolver_of_interface.cpp @@ -236,13 +236,13 @@ void ESolver_OF::kinetic_stress(ModuleBase::matrix& kinetic_stress_) if (this->of_kinetic_ == "wt") { - this->wt_->get_stress(pelec->charge->rho, this->pw_rho, PARAM.inp.of_vw_weight); + this->wt_->get_stress(this->chr.rho, this->pw_rho, PARAM.inp.of_vw_weight); kinetic_stress_ += this->wt_->stress; } if (this->of_kinetic_ == "lkt") { - this->lkt_->get_stress(pelec->charge->rho, this->pw_rho); + this->lkt_->get_stress(this->chr.rho, this->pw_rho); kinetic_stress_ += this->lkt_->stress; } if (this->of_kinetic_ == "ml") diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index c6f2725a1d..0fa78e97b3 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -20,7 +20,7 @@ void ESolver_OF::init_elecstate(UnitCell& ucell) if (this->pelec == nullptr) { this->pelec = new elecstate::ElecState((Charge*)(&chr), this->pw_rho, pw_big); - this->pelec->charge->allocate(PARAM.inp.nspin); + this->chr.allocate(PARAM.inp.nspin); } this->pelec->omega = ucell.omega; @@ -402,15 +402,15 @@ void ESolver_OF::print_info(const bool conv_esolver) // min/max(dE/dPhi)" << endl; } // ============ used to compare with PROFESS3.0 ================ - // double minDen = pelec->charge->rho[0][0]; - // double maxDen = pelec->charge->rho[0][0]; + // double minDen = this->chr.rho[0][0]; + // double maxDen = this->chr.rho[0][0]; // double minPot = this->pdEdphi_[0][0]; // double maxPot = this->pdEdphi_[0][0]; // for (int i = 0; i < this->pw_rho->nrxx; ++i) // { - // if (pelec->charge->rho[0][i] < minDen) minDen = - // pelec->charge->rho[0][i]; if (pelec->charge->rho[0][i] > maxDen) - // maxDen = pelec->charge->rho[0][i]; if (this->pdEdphi_[0][i] < minPot) + // if (this->chr.rho[0][i] < minDen) minDen = + // this->chr.rho[0][i]; if (this->chr.rho[0][i] > maxDen) + // maxDen = this->chr.rho[0][i]; if (this->pdEdphi_[0][i] < minPot) // minPot = this->pdEdphi_[0][i]; if (this->pdEdphi_[0][i] > maxPot) // maxPot = this->pdEdphi_[0][i]; // } diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index a733578820..e013372a11 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -191,7 +191,7 @@ void ESolver_SDFT_PW::hamilt2density_single(UnitCell& ucell, int iste Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); + srho.begin(is, this->chr, this->pw_rho, ucell.symm); } this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); } @@ -248,7 +248,7 @@ void ESolver_SDFT_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& this->pw_wfc, *this->kspw_psi, this->stowf, - this->pelec->charge, + &this->chr, &this->locpp, &this->ppcell, ucell); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index c6e020ca79..e991ac60f4 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -302,7 +302,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(this->Pgrid, - this->pelec->charge->rho[is], + this->chr.rho[is], is, PARAM.inp.nspin, istep, @@ -322,7 +322,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); + srho.begin(is, this->chr, this->pw_rho, ucell.symm); } this->p_hamilt->non_first_scf = istep; diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 8b211063cb..5426ce2294 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -279,7 +279,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) if (PARAM.globalv.gamma_only_local) { ISC.begin(this->GG, - this->pelec->charge->rho, + this->chr.rho, this->pelec->wg, this->pelec->eferm.get_all_ef(), this->pw_rhod->nrxx, @@ -307,8 +307,8 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) else { ISC.begin(this->GK, - this->pelec->charge->rho, - this->pelec->charge->rhog, + this->chr.rho, + this->chr.rhog, this->pelec->wg, this->pelec->eferm.get_all_ef(), this->pw_rhod, @@ -335,7 +335,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->kv, PARAM.inp.if_separate_k, &this->Pgrid, - this->pelec->charge->ngmc); + this->chr.ngmc); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); } From 27eae7cdcb82a4f5eb6f57725a36e54f3b26ac24 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 17:48:26 +0800 Subject: [PATCH 5/6] fix problem in ML_KEDF --- source/module_esolver/esolver_ks_pw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 5042f740e7..db5c5b613a 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -876,7 +876,7 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) // generate training data for ML-KEDF if (PARAM.inp.of_ml_gene_data == 1) { - this->pelec->pot->update_from_charge(this->chr, &ucell); + this->pelec->pot->update_from_charge(&this->chr, &ucell); ML_data ml_data; ml_data.set_para(this->chr.nrxx, From 82ba0b6f4d59c644d56adac697e0d6b7b4bd6d55 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 18:02:16 +0800 Subject: [PATCH 6/6] update this->chr --- source/module_esolver/esolver_of.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 1849bbb0f3..c79e04e2ea 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -532,7 +532,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso //------------------------------------------------------------------ if (PARAM.inp.of_ml_gene_data) { - this->pelec->pot->update_from_charge(this->chr, &ucell); // Hartree + XC + external + this->pelec->pot->update_from_charge(&this->chr, &ucell); // Hartree + XC + external this->kinetic_potential(this->chr.rho, this->pphi_, this->pelec->pot->get_effective_v()); // (kinetic + Hartree + XC + external) * 2 * phi const double* vr_eff = this->pelec->pot->get_effective_v(0);