From 285e9fb9d6dc7db97f2b8431327f38ff39097d54 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 09:35:13 +0800 Subject: [PATCH 01/10] 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 02/10] 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 03/10] 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 04/10] 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 05/10] 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 06/10] 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); From 508c542775bcef02d5c742ad6e7723a807419cea Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 18:21:54 +0800 Subject: [PATCH 07/10] fix some format issues --- source/module_elecstate/cal_ux.cpp | 103 +++++++++++++++----------- source/module_elecstate/magnetism.cpp | 37 ++++++--- source/module_elecstate/magnetism.h | 9 +-- 3 files changed, 93 insertions(+), 56 deletions(-) diff --git a/source/module_elecstate/cal_ux.cpp b/source/module_elecstate/cal_ux.cpp index 441aacceae..65e3b6bc72 100644 --- a/source/module_elecstate/cal_ux.cpp +++ b/source/module_elecstate/cal_ux.cpp @@ -4,62 +4,81 @@ namespace elecstate { void cal_ux(UnitCell& ucell) { + if (PARAM.inp.nspin != 4) { return; } - double amag, uxmod; + const double absolute_mag_thr = 1.0e-6; + + double amag = 0.0; + double uxmod = 0.0; + int starting_it = 0; int starting_ia = 0; - bool is_paraller; + bool is_paraller = false; + // do not sign feature in teh general case ucell.magnet.lsign_ = false; ModuleBase::GlobalFunc::ZEROS(ucell.magnet.ux_, 3); - for (int it = 0; it < ucell.ntype; it++) { - for (int ia = 0; ia < ucell.atoms[it].na; ia++) { + for (int it = 0; it < ucell.ntype; it++) + { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + // m_loc_: local magnetization vector for each atom amag = pow(ucell.atoms[it].m_loc_[ia].x, 2) + pow(ucell.atoms[it].m_loc_[ia].y, 2) + pow(ucell.atoms[it].m_loc_[ia].z, 2); - if (amag > 1e-6) { - ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x; - ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y; - ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z; - starting_it = it; - starting_ia = ia; - ucell.magnet.lsign_ = true; - break; - } - } - if (ucell.magnet.lsign_) { - break; - } - } + + if (amag > absolute_mag_thr) + { + ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x; + ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y; + ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z; + starting_it = it; + starting_ia = ia; + ucell.magnet.lsign_ = true; + break; + } + } + if (ucell.magnet.lsign_) + { + break; + } + } + // whether the initial magnetizations is parallel - for (int it = starting_it; it < ucell.ntype; it++) { - for (int ia = 0; ia < ucell.atoms[it].na; ia++) { - if (it > starting_it || ia > starting_ia) { - ucell.magnet.lsign_ - = ucell.magnet.lsign_ - && judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]); - } - } - } - if (ucell.magnet.lsign_) { - uxmod = pow(ucell.magnet.ux_[0], 2) + pow(ucell.magnet.ux_[1], 2) - + pow(ucell.magnet.ux_[2], 2); - if (uxmod < 1e-6) { - ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod"); - } - for (int i = 0; i < 3; i++) { - ucell.magnet.ux_[i] *= 1 / sqrt(uxmod); - } - // std::cout<<" Fixed quantization axis for GGA: " - //< starting_it || ia > starting_ia) + { + ucell.magnet.lsign_ + = ucell.magnet.lsign_ + && judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]); + } + } + } + + if (ucell.magnet.lsign_) + { + uxmod = pow(ucell.magnet.ux_[0], 2) + + pow(ucell.magnet.ux_[1], 2) + + pow(ucell.magnet.ux_[2], 2); + + if (uxmod < absolute_mag_thr) + { + ModuleBase::WARNING_QUIT("elecstate::cal_ux", "wrong uxmod"); + } + for (int i = 0; i < 3; i++) + { + ucell.magnet.ux_[i] *= 1 / sqrt(uxmod); + } + } + return; } bool judge_parallel(double a[3], ModuleBase::Vector3 b) { @@ -71,4 +90,4 @@ bool judge_parallel(double a[3], ModuleBase::Vector3 b) { jp = (fabs(cross) < 1e-6); return jp; } -} \ No newline at end of file +} diff --git a/source/module_elecstate/magnetism.cpp b/source/module_elecstate/magnetism.cpp index b41f5261ad..4cf77c3ff0 100644 --- a/source/module_elecstate/magnetism.cpp +++ b/source/module_elecstate/magnetism.cpp @@ -22,6 +22,8 @@ void Magnetism::compute_magnetization(const double& omega, const double* const * rho, double* nelec_spin) { + assert(nxyz>0); + if (PARAM.inp.nspin==2) { this->tot_magnetization = 0.00; @@ -57,36 +59,53 @@ void Magnetism::compute_magnetization(const double& omega, // noncolliear : else if(PARAM.inp.nspin==4) { - for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] = 0.00; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] = 0.00; + } + this->abs_magnetization = 0.00; for (int ir=0; irtot_magnetization_nc[i] += rho[i+1][ir]; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] += rho[i+1][ir]; + } this->abs_magnetization += std::abs(diff); } #ifdef __MPI Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3); Parallel_Reduce::reduce_pool(this->abs_magnetization); #endif - for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] *= omega/ nxyz; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] *= omega/ nxyz; + } + this->abs_magnetization *= omega/ nxyz; - GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<tot_magnetization_nc[0]<<'\t'<tot_magnetization_nc[1]<<'\t'<tot_magnetization_nc[2]<<'\n'; + + GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)" + <<'\t'<tot_magnetization_nc[0] + <<'\t'<tot_magnetization_nc[1] + <<'\t'<tot_magnetization_nc[2]<<'\n'; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization); } return; } + bool Magnetism::judge_parallel(double a[3], ModuleBase::Vector3 b) { bool jp=false; - double cross; - cross = pow((a[1]*b.z-a[2]*b.y),2) + pow((a[2]*b.x-a[0]*b.z),2) + pow((a[0]*b.y-a[1]*b.x),2); + double cross=0.0; + cross = pow((a[1]*b.z-a[2]*b.y),2) + + pow((a[2]*b.x-a[0]*b.z),2) + + pow((a[0]*b.y-a[1]*b.x),2); + jp = (fabs(cross)<1e-6); return jp; } diff --git a/source/module_elecstate/magnetism.h b/source/module_elecstate/magnetism.h index 408ba9d039..8a841c7006 100644 --- a/source/module_elecstate/magnetism.h +++ b/source/module_elecstate/magnetism.h @@ -13,7 +13,7 @@ class Magnetism Magnetism(); ~Magnetism(); - // notice : becast is done in unitcell + // notice : bcast (MPI operation) is done in unitcell double *start_magnetization=nullptr; // tot_magnetization : majority spin - minority spin (nelup - neldw). @@ -28,10 +28,9 @@ class Magnetism double* nelec_spin = nullptr); - ModuleBase::Vector3 *m_loc_=nullptr;//magnetization for each element along c-axis - double *angle1_=nullptr; //angle between c-axis and real spin std::vector - double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane - //void cal_ux(const int ntype); + ModuleBase::Vector3 *m_loc_=nullptr; //magnetization for each element along c-axis + double *angle1_=nullptr; //angle between c-axis and real spin std::vector + double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane double ux_[3]={0.0}; bool lsign_=false; From 5af7ded7025aaa321f3ce07c85940526a5933ea4 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 18:28:20 +0800 Subject: [PATCH 08/10] update cal_ux.cpp --- source/module_elecstate/cal_ux.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/module_elecstate/cal_ux.cpp b/source/module_elecstate/cal_ux.cpp index 65e3b6bc72..f8e4b1c82a 100644 --- a/source/module_elecstate/cal_ux.cpp +++ b/source/module_elecstate/cal_ux.cpp @@ -32,17 +32,23 @@ void cal_ux(UnitCell& ucell) { + pow(ucell.atoms[it].m_loc_[ia].y, 2) + pow(ucell.atoms[it].m_loc_[ia].z, 2); + // find the first atom (it,ia) whose magnetism is not zero + // compute ux if (amag > absolute_mag_thr) { ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x; ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y; ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z; + starting_it = it; starting_ia = ia; + ucell.magnet.lsign_ = true; break; } } + + // if any atom has magnetism, then break the for iteration if (ucell.magnet.lsign_) { break; @@ -63,6 +69,8 @@ void cal_ux(UnitCell& ucell) { } } + // if all of the atoms have the same parallel magnetism direction, + // then set the direction to a unit vector if (ucell.magnet.lsign_) { uxmod = pow(ucell.magnet.ux_[0], 2) @@ -73,6 +81,8 @@ void cal_ux(UnitCell& ucell) { { ModuleBase::WARNING_QUIT("elecstate::cal_ux", "wrong uxmod"); } + + // reset the magnetism for each direction for (int i = 0; i < 3; i++) { ucell.magnet.ux_[i] *= 1 / sqrt(uxmod); From ac6fb24b24ca93904b69deb442f279dce4955bd4 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Thu, 6 Mar 2025 18:36:23 +0800 Subject: [PATCH 09/10] update fp --- source/module_esolver/esolver_fp.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index c92c2ebd34..1f31453ebd 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -258,6 +258,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_FP", "before_scf"); + // if the cell has changed if (ucell.cell_parameter_updated) { // only G-vector and K-vector are changed due to the change of lattice @@ -266,6 +267,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rho->collect_local_pw(); this->pw_rho->collect_uniqgg(); + // if double grid used in USPP, update related quantities in dense grid if (PARAM.globalv.double_grid) { this->pw_rhod->initgrids(ucell.lat0, ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz); @@ -273,35 +275,40 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rhod->collect_uniqgg(); } + // reset local pseudopotentials this->locpp.init_vloc(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); this->pelec->omega = ucell.omega; + // perform symmetry analysis if (ModuleSymmetry::Symmetry::symm_flag == 1) { ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } + // reset k-points kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } + //---------------------------------------------------------- // charge extrapolation + //---------------------------------------------------------- if (ucell.ionic_position_updated) { this->CE.update_all_dis(ucell); - this->CE.extrapolate_charge(&(this->Pgrid), + this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, - &(this->sf), + &this->sf, GlobalV::ofs_running, GlobalV::ofs_warning); } //---------------------------------------------------------- - // about vdw, jiyy add vdwd3 and linpz add vdwd2 + //! calculate D2 or D3 vdW //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) @@ -309,19 +316,22 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pelec->f_en.evdw = vdw_solver->get_energy(); } - // calculate ewald energy + //---------------------------------------------------------- + //! calculate ewald energy + //---------------------------------------------------------- if (!PARAM.inp.test_skip_ewald) { this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac); } //---------------------------------------------------------- - //! cal_ux should be called before init_scf because - //! the direction of ux is used in noncoline_rho + //! set direction of magnetism, used in non-collinear case //---------------------------------------------------------- elecstate::cal_ux(ucell); + //---------------------------------------------------------- //! output the initial charge density + //---------------------------------------------------------- if (PARAM.inp.out_chg[0] == 2) { for (int is = 0; is < PARAM.inp.nspin; is++) @@ -339,7 +349,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) } } + //---------------------------------------------------------- //! output total local potential of the initial charge density + //---------------------------------------------------------- if (PARAM.inp.out_pot == 3) { for (int is = 0; is < PARAM.inp.nspin; is++) From 10e51bbfedc8e558a566dfef70921e91d0761bf5 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Fri, 7 Mar 2025 16:07:27 +0800 Subject: [PATCH 10/10] small updates --- source/module_esolver/esolver_ks_lcao.cpp | 4 +++- source/module_esolver/lcao_others.cpp | 21 ++++++++++++--------- source/module_esolver/pw_others.cpp | 11 ++++++++--- 3 files changed, 23 insertions(+), 13 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index b123a354fe..1e6cc8d827 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -160,7 +160,9 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa #ifdef __EXX // 5) initialize exx // PLEASE simplify the Exx_Global interface - if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" + if (PARAM.inp.calculation == "scf" + || PARAM.inp.calculation == "relax" + || PARAM.inp.calculation == "cell-relax" || PARAM.inp.calculation == "md") { if (GlobalC::exx_info.info_global.cal_exx) diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 5426ce2294..dc7b522c43 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -21,15 +21,16 @@ #include "module_base/formatter.h" #include "module_elecstate/elecstate_lcao.h" #include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_hamilt_general/module_ewald/H_Ewald_pw.h" -#include "module_hamilt_general/module_vdw/vdw.h" + #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" + #include "module_io/read_wfc_nao.h" #include "module_io/write_elecstat_pot.h" #include "module_io/write_wfc_nao.h" + #ifdef __EXX #include "module_io/restart_exx_csr.h" #endif @@ -263,6 +264,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->psi, this->pelec); } + //========================================================= // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho @@ -271,14 +273,15 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // pelec should be initialized before these calculations this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); - IState_Charge ISC(this->psi, &(this->pv)); + IState_Charge chr_i(this->psi, &(this->pv)); if (PARAM.globalv.gamma_only_local) { - ISC.begin(this->GG, + chr_i.begin(this->GG, this->chr.rho, this->pelec->wg, this->pelec->eferm.get_all_ef(), @@ -306,7 +309,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - ISC.begin(this->GK, + chr_i.begin(this->GK, this->chr.rho, this->chr.rhog, this->pelec->wg, @@ -342,10 +345,10 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) else if (cal_type == "get_wf") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); - IState_Envelope IEP(this->pelec); + IState_Envelope wavefunc(this->pelec); if (PARAM.globalv.gamma_only_local) { - IEP.begin(ucell, + wavefunc.begin(ucell, this->psi, this->pw_rhod, this->pw_wfc, @@ -367,7 +370,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - IEP.begin(ucell, + wavefunc.begin(ucell, this->psi, this->pw_rhod, this->pw_wfc, @@ -391,7 +394,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); } ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index c36ddae6e7..5461bed0a0 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -58,7 +58,8 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) const std::string cal_type = PARAM.inp.calculation; - if (cal_type == "test_memory") { + if (cal_type == "test_memory") + { Cal_Test::test_memory(ucell.nat, ucell.ntype, ucell.GGT, @@ -66,7 +67,9 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); - } else if (cal_type == "gen_bessel") { + } + else if (cal_type == "gen_bessel") + { Numerical_Descriptor nc; nc.output_descriptor(ucell, this->psi[0], @@ -75,7 +78,9 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) PARAM.inp.bessel_descriptor_tolerence, this->kv.get_nks()); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "GENERATE DESCRIPTOR FOR DEEPKS"); - } else { + } + else + { ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); }