diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 5442ad010c..537ad1a4dd 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -238,6 +238,7 @@ OBJS_ELECSTAT=elecstate.o\ gatefield.o\ potential_new.o\ potential_types.o\ + pot_sep.o\ pot_local.o\ pot_local_paw.o\ H_Hartree_pw.o\ @@ -248,7 +249,7 @@ OBJS_ELECSTAT=elecstate.o\ cal_nelec_nband.o\ read_pseudo.o\ cal_wfc.o\ - pot_sep.o\ + setup_estate_pw.o\ OBJS_ELECSTAT_LCAO=elecstate_lcao.o\ elecstate_lcao_cal_tau.o\ @@ -714,6 +715,8 @@ OBJS_SRCPW=H_Ewald_pw.o\ charge_mixing_rho.o\ charge_mixing_uspp.o\ fp_energy.o\ + setup_pot.o\ + setup_pwrho.o\ forces.o\ forces_us.o\ forces_nl.o\ diff --git a/source/source_esolver/esolver.h b/source/source_esolver/esolver.h index 58238dcb1b..6716ea0c96 100644 --- a/source/source_esolver/esolver.h +++ b/source/source_esolver/esolver.h @@ -17,6 +17,9 @@ class ESolver virtual ~ESolver() { + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** } //! initialize the energy solver by using input parameters and cell modules diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 9d8326d2a3..f6aff301bd 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -1,12 +1,10 @@ #include "esolver_fp.h" -#include "source_base/global_variable.h" #include "source_estate/cal_ux.h" #include "source_estate/module_charge/symmetry_rho.h" #include "source_estate/read_pseudo.h" #include "source_hamilt/module_ewald/H_Ewald_pw.h" #include "source_hamilt/module_vdw/vdw.h" -#include "source_pw/module_pwdft/global.h" #include "source_io/cif_io.h" #include "source_io/cube_io.h" // use write_vdata_palgrid #include "source_io/json_output/init_info.h" @@ -15,9 +13,11 @@ #include "source_io/print_info.h" #include "source_io/rhog_io.h" #include "source_io/module_parameter/parameter.h" -#include "source_cell/k_vector_utils.h" #include "source_io/ctrl_output_fp.h" +#include "source_pw/module_pwdft/setup_pwrho.h" // mohan 20251005 +#include "source_hamilt/module_xc/xc_functional.h" // mohan 20251005 + namespace ModuleESolver { @@ -27,117 +27,36 @@ ESolver_FP::ESolver_FP() ESolver_FP::~ESolver_FP() { - if (pw_rho_flag == true) - { - delete this->pw_rho; - this->pw_rho_flag = false; - } - if (PARAM.globalv.double_grid) - { - delete pw_rhod; - } - delete this->pelec; + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + // mohan add 20251005 + pw::teardown_pwrho(this->pw_rho_flag, PARAM.globalv.double_grid, this->pw_rho, this->pw_rhod); + + delete this->pelec; } void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_FP", "before_all_runners"); - std::string fft_device = inp.device; - std::string fft_precison = inp.precision; - // LCAO basis doesn't support GPU acceleration on FFT currently - if(inp.basis_type == "lcao") - { - fft_device = "cpu"; - } - if ((inp.precision=="single") || (inp.precision=="mixing")) - { - fft_precison = "mixing"; - } - else if (inp.precision=="double") - { - fft_precison = "double"; - } - #if (not defined(__ENABLE_FLOAT_FFTW) and (defined(__CUDA) || defined(__RCOM))) - if (fft_device == "gpu") - { - fft_precison = "double"; - } - #endif - pw_rho = new ModulePW::PW_Basis_Big(fft_device, fft_precison); - pw_rho_flag = true; - if (PARAM.globalv.double_grid) - { - pw_rhod = new ModulePW::PW_Basis_Big(fft_device, fft_precison); - } - else - { - pw_rhod = pw_rho; - } - pw_big = static_cast(pw_rhod); - pw_big->setbxyz(inp.bx, inp.by, inp.bz); - sf.set(pw_rhod, inp.nbspline); - //! 1) read pseudopotentials + //! read pseudopotentials elecstate::read_pseudo(GlobalV::ofs_running, ucell); - //! 2) initialie the plane wave basis for rho -#ifdef __MPI - this->pw_rho->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); -#endif - if (this->classname == "ESolver_OF" || inp.of_ml_gene_data == 1) - { - this->pw_rho->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); - } + // setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials + pw::setup_pwrho(ucell, PARAM.globalv.double_grid, this->pw_rho_flag, + this->pw_rho, this->pw_rhod, this->pw_big, + this->classname, inp); - if (inp.nx * inp.ny * inp.nz == 0) - { - this->pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, 4.0 * inp.ecutwfc); - } - else - { - this->pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.nx, inp.ny, inp.nz); - } + // setup the structure factors + this->sf.set(this->pw_rhod, inp.nbspline); - this->pw_rho->initparameters(false, 4.0 * inp.ecutwfc); - this->pw_rho->fft_bundle.initfftmode(inp.fft_mode); - this->pw_rho->setuptransform(); - this->pw_rho->collect_local_pw(); - this->pw_rho->collect_uniqgg(); - - //! 3) initialize the double grid (for uspp) if necessary - if ( PARAM.globalv.double_grid) - { - ModulePW::PW_Basis_Sup* pw_rhod_sup = static_cast(pw_rhod); -#ifdef __MPI - this->pw_rhod->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); -#endif - if (this->classname == "ESolver_OF") - { - this->pw_rhod->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); - } - if (inp.ndx * inp.ndy * inp.ndz == 0) - { - this->pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ecutrho); - } - else - { - this->pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ndx, inp.ndy, inp.ndz); - } - this->pw_rhod->initparameters(false, inp.ecutrho); - this->pw_rhod->fft_bundle.initfftmode(inp.fft_mode); - pw_rhod_sup->setuptransform(this->pw_rho); - this->pw_rhod->collect_local_pw(); - this->pw_rhod->collect_uniqgg(); - } ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?"); - //! 4) print some information - ModuleIO::print_rhofft(this->pw_rhod, this->pw_rho, this->pw_big, GlobalV::ofs_running); - - //! 5) initialize the charge extrapolation method if necessary + //! initialize the charge extrapolation method if necessary this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap); return; @@ -148,16 +67,16 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso { ModuleBase::TITLE("ESolver_FP", "after_scf"); - // 1) output convergence information + //! Output convergence information ModuleIO::output_convergence_after_scf(conv_esolver, this->pelec->f_en.etot); - // 2) write fermi energy + //! Write Fermi energy ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef); - // 3) update delta_rho for charge extrapolation + //! Update delta_rho for charge extrapolation CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); - // 4) print out charge density, potential, elf, etc. + //! print out charge density, potential, elf, etc. ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, this->chr, this->solvent, this->Pgrid, istep); @@ -202,9 +121,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } - //---------------------------------------------------------- // charge extrapolation - //---------------------------------------------------------- if (ucell.ionic_position_updated) { this->CE.update_all_dis(ucell); @@ -216,33 +133,23 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) GlobalV::ofs_warning); } - //---------------------------------------------------------- //! calculate D2 or D3 vdW - //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_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); } - //---------------------------------------------------------- //! set direction of magnetism, used in non-collinear case - //---------------------------------------------------------- elecstate::cal_ux(ucell); - - - //---------------------------------------------------------- //! output the initial charge density - //---------------------------------------------------------- const int nspin = PARAM.inp.nspin; if (PARAM.inp.out_chg[0] == 2) { @@ -271,9 +178,7 @@ 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 < nspin; is++) @@ -308,7 +213,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { - //! output charge density + //! output charge density in G-space, or if available, kinetic energy density in G-space if (PARAM.inp.out_chg[0] != -1) { if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) @@ -352,10 +257,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& void ESolver_FP::after_all_runners(UnitCell& ucell) { + // print out the final total energy GlobalV::ofs_running << "\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; + } } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_fp.h b/source/source_esolver/esolver_fp.h index 7d4f34736e..b2bb8f065e 100644 --- a/source/source_esolver/esolver_fp.h +++ b/source/source_esolver/esolver_fp.h @@ -47,9 +47,7 @@ class ESolver_FP: public ESolver virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver); - //! ------------------------------------------------------------------------------ //! These pointers will be deleted in the free_pointers() function every ion step. - //! ------------------------------------------------------------------------------ elecstate::ElecState* pelec = nullptr; ///< Electronic states //! K points in Brillouin zone @@ -82,7 +80,7 @@ class ESolver_FP: public ESolver //! solvent model surchem solvent; - int pw_rho_flag = false; ///< flag for pw_rho, 0: not initialized, 1: initialized + bool pw_rho_flag = false; ///< flag for pw_rho, 0: not initialized, 1: initialized //! the start time of scf iteration #ifdef __MPI diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index ba919399e7..b477c89656 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -26,7 +26,6 @@ #include "source_io/json_output/output_info.h" - namespace ModuleESolver { @@ -39,7 +38,10 @@ ESolver_KS::ESolver_KS() template ESolver_KS::~ESolver_KS() { - delete this->psi; + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + delete this->psi; delete this->pw_wfc; delete this->p_hamilt; delete this->p_chgmix; diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index e7c212f829..b7ef11ed55 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -97,6 +97,9 @@ ESolver_KS_LCAO::ESolver_KS_LCAO() template ESolver_KS_LCAO::~ESolver_KS_LCAO() { + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** } template diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 857da73913..c68ff0955b 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -57,7 +57,10 @@ ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT() template ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() { - delete psi_laststep; + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + delete psi_laststep; if (Hk_laststep != nullptr) { for (int ik = 0; ik < this->kv.get_nks(); ++ik) diff --git a/source/source_esolver/esolver_ks_lcaopw.cpp b/source/source_esolver/esolver_ks_lcaopw.cpp index 1b0b7dae67..8597fa6847 100644 --- a/source/source_esolver/esolver_ks_lcaopw.cpp +++ b/source/source_esolver/esolver_ks_lcaopw.cpp @@ -51,6 +51,9 @@ namespace ModuleESolver template ESolver_KS_LIP::~ESolver_KS_LIP() { + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** delete this->psi_local; // delete Hamilt this->deallocate_hamilt(); diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 4bf5fa6a6c..6bb78b5c25 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -1,31 +1,22 @@ #include "esolver_ks_pw.h" -#include "source_base/global_variable.h" -#include "source_base/kernels/math_kernel_op.h" -#include "source_base/memory.h" #include "source_estate/cal_ux.h" #include "source_estate/elecstate_pw.h" -#include "source_estate/elecstate_pw_sdft.h" -#include "source_estate/elecstate_tools.h" #include "source_estate/module_charge/symmetry_rho.h" -#include "source_hamilt/module_ewald/H_Ewald_pw.h" -#include "source_hamilt/module_vdw/vdw.h" + #include "source_hsolver/diago_iter_assist.h" #include "source_hsolver/hsolver_pw.h" + #include "source_hsolver/kernels/hegvd_op.h" #include "source_io/module_parameter/parameter.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/onsite_projector.h" #include "source_lcao/module_dftu/dftu.h" #include "source_pw/module_pwdft/VSep_in_pw.h" -#include "source_pw/module_pwdft/forces.h" #include "source_pw/module_pwdft/hamilt_pw.h" -#include "source_pw/module_pwdft/stress_pw.h" -#include - -#include -#include +#include "source_pw/module_pwdft/forces.h" +#include "source_pw/module_pwdft/stress_pw.h" #ifdef __DSP #include "source_base/kernels/dsp/dsp_connector.h" @@ -33,6 +24,8 @@ #include +#include "source_pw/module_pwdft/setup_pot.h" // mohan add 20250929 +#include "source_estate/setup_estate_pw.h" // mohan add 20251005 #include "source_io/ctrl_output_pw.h" // mohan add 20250927 namespace ModuleESolver @@ -49,20 +42,12 @@ ESolver_KS_PW::ESolver_KS_PW() template ESolver_KS_PW::~ESolver_KS_PW() { - // delete Hamilt + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + // delete Hamilt this->deallocate_hamilt(); - if (this->vsep_cell != nullptr) - { - delete this->vsep_cell; - } - - if (this->pelec != nullptr) - { - delete reinterpret_cast*>(this->pelec); - this->pelec = nullptr; - } - if (PARAM.inp.device == "gpu" || PARAM.inp.precision == "single") { delete this->kspw_psi; @@ -95,84 +80,18 @@ void ESolver_KS_PW::deallocate_hamilt() template void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) { - // 1) call before_all_runners() of ESolver_KS + //! Call before_all_runners() of ESolver_KS ESolver_KS::before_all_runners(ucell, inp); - // 2) initialize ElecState, set pelec pointer - if (this->pelec == nullptr) - { - if (inp.esolver_type == "sdft") - { - //! SDFT only supports double precision currently - this->pelec = new elecstate::ElecStatePW_SDFT, Device>(this->pw_wfc, - &(this->chr), - &(this->kv), - &ucell, - &(this->ppcell), - this->pw_rhod, - this->pw_rho, - this->pw_big); - } - else - { - this->pelec = new elecstate::ElecStatePW(this->pw_wfc, - &(this->chr), - &(this->kv), - &ucell, - &this->ppcell, - this->pw_rhod, - this->pw_rho, - this->pw_big); - } - } - - //! set the cell volume variable in pelec - this->pelec->omega = ucell.omega; - - //! 3) inititlize the charge density. - this->chr.allocate(inp.nspin); - - // 3.5) initialize DFT-1/2 - if (PARAM.inp.dfthalf_type > 0) - { - this->vsep_cell = new VSep; - this->vsep_cell->init_vsep(*this->pw_rhod, ucell.sep_cell); - } - - - //! 4) initialize the potential. - if (this->pelec->pot == nullptr) - { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &ucell, - &this->locpp.vloc, - &(this->sf), - &(this->solvent), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc), - this->vsep_cell); - } - - //! 5) Initalize local pseudopotential - this->locpp.init_vloc(ucell, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - - //! 6) Initalize non-local pseudopotential - this->ppcell.init(ucell, &this->sf, this->pw_wfc); - this->ppcell.init_vnl(ucell, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); + //! setup and allocation for pelec, charge density, potentials, etc. + elecstate::setup_estate_pw(ucell, this->kv, this->sf, this->pelec, this->chr, + this->locpp, this->ppcell, this->vsep_cell, this->pw_wfc, this->pw_rho, + this->pw_rhod, this->pw_big, this->solvent, inp); - //! 7) Allocate and initialize psi + //! Allocate and initialize psi this->p_psi_init = new psi::PSIInit(inp.init_wfc, - inp.ks_solver, - inp.basis_type, - GlobalV::MY_RANK, - ucell, - this->sf, - this->kv, - this->ppcell, - *this->pw_wfc); + inp.ks_solver, inp.basis_type, GlobalV::MY_RANK, ucell, + this->sf, this->kv, this->ppcell, *this->pw_wfc); allocate_psi(this->psi, this->kv.get_nks(), this->kv.ngk, PARAM.globalv.nbands_l, this->pw_wfc->npwk_max); @@ -184,18 +103,7 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); - //! 8) setup occupations - if (inp.ocp) - { - elecstate::fixed_weights(inp.ocp_kb, - inp.nbands, - inp.nelec, - this->pelec->klist, - this->pelec->wg, - this->pelec->skip_weights); - } - - // 9) initialize exx pw + //! Initialize exx pw if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax" || inp.calculation == "md") { @@ -218,14 +126,10 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); ModuleBase::timer::tick("ESolver_KS_PW", "before_scf"); - //---------------------------------------------------------- - //! 1) Call before_scf() of ESolver_KS - //---------------------------------------------------------- + //! Call before_scf() of ESolver_KS ESolver_KS::before_scf(ucell, istep); - //---------------------------------------------------------- - //! 2) Init variables (cell changed) - //---------------------------------------------------------- + //! Init variables (once the cell has changed) if (ucell.cell_parameter_updated) { this->ppcell.rescale_vnl(ucell.omega); @@ -240,119 +144,29 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->p_psi_init->prepare_init(PARAM.inp.pw_seed); } - //---------------------------------------------------------- - //! 3) init Hamiltonian (cell changed) - //---------------------------------------------------------- + //! Init Hamiltonian (cell changed) //! Operators in HamiltPW should be reallocated once cell changed //! delete Hamilt if not first scf this->deallocate_hamilt(); - // allocate HamiltPW + //! Allocate HamiltPW this->allocate_hamilt(ucell); - //---------------------------------------------------------- - //! 4) DFT-1/2 calculations, sep potential need to generate before effective potential calculation - //---------------------------------------------------------- - if (PARAM.inp.dfthalf_type > 0) - { - this->vsep_cell->generate_vsep_r(this->pw_rhod[0], this->sf.strucFac, ucell.sep_cell); - } - - //---------------------------------------------------------- - //! 5) Renew local pseudopotential - //---------------------------------------------------------- - this->pelec - ->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm, (void*)this->pw_wfc); - - //---------------------------------------------------------- - //! 6) Symmetrize the charge density (rho) - //---------------------------------------------------------- - - //! Symmetry_rho should behind init_scf, because charge should be - //! initialized first. liuyu comment: Symmetry_rho should be - //! located between init_rho and v_of_rho? - Symmetry_rho srho; - for (int is = 0; is < PARAM.inp.nspin; is++) - { - srho.begin(is, this->chr, this->pw_rhod, ucell.symm); - } + //! Setup potentials (local, non-local, sc, +U, DFT-1/2) + pw::setup_pot(istep, ucell, this->kv, this->sf, this->pelec, this->Pgrid, + this->chr, this->locpp, this->ppcell, this->vsep_cell, + this->kspw_psi, this->p_hamilt, this->pw_wfc, this->pw_rhod, PARAM.inp); - //---------------------------------------------------------- - //! 7) Calculate the effective potential with rho - //---------------------------------------------------------- - //! liuyu move here 2023-10-09 - //! D in uspp need vloc, thus behind init_scf() - //! calculate the effective coefficient matrix - //! for non-local pseudopotential projectors - ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - - this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); - - //---------------------------------------------------------- - //! 8) Onsite projectors - //---------------------------------------------------------- - if (PARAM.inp.onsite_radius > 0) - { - auto* onsite_p = projectors::OnsiteProjector::get_instance(); - onsite_p->init(PARAM.inp.orbital_dir, - &ucell, - *(this->kspw_psi), - this->kv, - *(this->pw_wfc), - this->sf, - PARAM.inp.onsite_radius, - PARAM.globalv.nqx, - PARAM.globalv.dq, - this->pelec->wg, - this->pelec->ekb); - } - - //---------------------------------------------------------- - //! 9) Spin-constrained algorithms - //---------------------------------------------------------- - if (PARAM.inp.sc_mag_switch) - { - spinconstrain::SpinConstrain>& sc - = spinconstrain::SpinConstrain>::getScInstance(); - sc.init_sc(PARAM.inp.sc_thr, - PARAM.inp.nsc, - PARAM.inp.nsc_min, - PARAM.inp.alpha_trial, - PARAM.inp.sccut, - PARAM.inp.sc_drop_thr, - ucell, - nullptr, - PARAM.inp.nspin, - this->kv, - this->p_hamilt, - this->kspw_psi, - this->pelec, - this->pw_wfc); - } - - //---------------------------------------------------------- - //! 10) DFT+U algorithm - //---------------------------------------------------------- - if (PARAM.inp.dft_plus_u) - { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - dftu->init(ucell, nullptr, this->kv.get_nks()); - } - - //---------------------------------------------------------- - //! 10) Initialize wave functions - //---------------------------------------------------------- + //! Initialize wave functions if (!this->already_initpsi) { this->p_psi_init->initialize_psi(this->psi, this->kspw_psi, this->p_hamilt, GlobalV::ofs_running); this->already_initpsi = true; } - //---------------------------------------------------------- - //! 11) Exx calculations - //---------------------------------------------------------- - if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" - || PARAM.inp.calculation == "md") + //! Exx calculations + 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 && PARAM.inp.basis_type == "pw") { @@ -368,7 +182,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) template void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const int iter) { - // call iter_init() of ESolver_KS + // Call iter_init() of ESolver_KS ESolver_KS::iter_init(ucell, istep, iter); if (iter == 1) @@ -377,10 +191,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const this->p_chgmix->mixing_restart_step = PARAM.inp.scf_nmax + 1; } - //---------------------------------------------------------- - // for mixing restart - // the details should move somewhere else, 2025-04-13 - //---------------------------------------------------------- + // For mixing restart if (iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0) { this->p_chgmix->init_mixing(); @@ -419,16 +230,12 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const } } - //---------------------------------------------------------- // mohan move harris functional to here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) - //---------------------------------------------------------- this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell); - //---------------------------------------------------------- // update local occupations for DFT+U // should before lambda loop in DeltaSpin - //---------------------------------------------------------- if (PARAM.inp.dft_plus_u && (iter != 1 || istep != 0)) { auto* dftu = ModuleDFTU::DFTU::get_instance(); @@ -547,38 +354,28 @@ void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, cons template void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { - //---------------------------------------------------------- // Related to EXX - //---------------------------------------------------------- if (GlobalC::exx_info.info_global.cal_exx && !exx_helper.op_exx->first_iter) { this->pelec->set_exx(exx_helper.cal_exx_energy(kspw_psi)); } - //---------------------------------------------------------- - // deband is calculated from "output" charge density calculated - // in sum_band - // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))' - //---------------------------------------------------------- + // deband is calculated from "output" charge density this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); - //---------------------------------------------------------- - // 1) Call iter_finish() of ESolver_KS - //---------------------------------------------------------- + // Call iter_finish() of ESolver_KS ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); - //---------------------------------------------------------- - // 2) Update USPP-related quantities // D in USPP needs vloc, thus needs update when veff updated // calculate the effective coefficient matrix for non-local // pp projectors, liuyu 2023-10-24 - //---------------------------------------------------------- if (PARAM.globalv.use_uspp) { ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); } + // Related to EXX if (GlobalC::exx_info.info_global.cal_exx) { if (GlobalC::exx_info.info_global.separate_loop) @@ -621,17 +418,15 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int } } - //---------------------------------------------------------- - // 4) check if oscillate for delta_spin method - //---------------------------------------------------------- + // check if oscillate for delta_spin method if (PARAM.inp.sc_mag_switch) { spinconstrain::SpinConstrain>& sc = spinconstrain::SpinConstrain>::getScInstance(); if (!sc.higher_mag_prec) { - sc.higher_mag_prec - = this->p_chgmix->if_scf_oscillate(iter, this->drho, PARAM.inp.sc_os_ndim, PARAM.inp.scf_os_thr); + sc.higher_mag_prec = this->p_chgmix->if_scf_oscillate(iter, + this->drho, PARAM.inp.sc_os_ndim, PARAM.inp.scf_os_thr); if (sc.higher_mag_prec) { // if oscillate, increase the precision of magnetization and do mixing_restart in next iteration this->p_chgmix->mixing_restart_step = iter + 1; @@ -639,6 +434,7 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int } } + // the output quantities ModuleIO::ctrl_iter_pw(istep, iter, conv_esolver, this->psi, this->kv, this->pw_wfc, PARAM.inp); } @@ -649,25 +445,18 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const ModuleBase::TITLE("ESolver_KS_PW", "after_scf"); ModuleBase::timer::tick("ESolver_KS_PW", "after_scf"); - //------------------------------------------------------------------ - // 1) since ESolver_KS::psi is hidden by ESolver_KS_PW::psi, + // Since ESolver_KS::psi is hidden by ESolver_KS_PW::psi, // we need to copy the data from ESolver_KS::psi to ESolver_KS_PW::psi. - // This part needs to be removed when we have a better design. // sunliang 2025-04-10 - //------------------------------------------------------------------ if (PARAM.inp.out_elf[0] > 0) { this->ESolver_KS::psi = new psi::Psi(this->psi[0]); } - //------------------------------------------------------------------ - // 2) call after_scf() of ESolver_KS - //------------------------------------------------------------------ + // Call 'after_scf' of ESolver_KS ESolver_KS::after_scf(ucell, istep, conv_esolver); - //------------------------------------------------------------------ - // 3) transfer data from GPU to CPU in pw basis - //------------------------------------------------------------------ + // Transfer data from GPU to CPU in pw basis if (this->device == base_device::GpuDevice) { castmem_2d_d2h_op()(this->psi[0].get_pointer() - this->psi[0].get_psi_bias(), @@ -675,6 +464,7 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const this->psi[0].size()); } + // Output quantities ModuleIO::ctrl_scf_pw(istep, ucell, this->pelec, this->chr, this->kv, this->pw_wfc, this->pw_rho, this->pw_rhod, this->pw_big, this->psi, this->kspw_psi, this->__kspw_psi, this->ctx, this->Pgrid, PARAM.inp); @@ -740,9 +530,6 @@ void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& s template void ESolver_KS_PW::after_all_runners(UnitCell& ucell) { - //---------------------------------------------------------- - //! 1) after_all_runners in ESolver_KS - //---------------------------------------------------------- ESolver_KS::after_all_runners(ucell); ModuleIO::ctrl_runner_pw(ucell, this->pelec, this->pw_wfc, @@ -750,6 +537,8 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) this->kspw_psi, this->__kspw_psi, this->sf, this->ppcell, this->solvent, this->ctx, this->Pgrid, PARAM.inp); + elecstate::teardown_estate_pw(this->pelec, this->vsep_cell); + } template class ESolver_KS_PW, base_device::DEVICE_CPU>; diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 45f26babcb..6a53152030 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -25,7 +25,10 @@ ESolver_OF::ESolver_OF() ESolver_OF::~ESolver_OF() { - delete psi_; + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** + delete psi_; delete[] this->pphi_; for (int i = 0; i < PARAM.inp.nspin; ++i) diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 42fde3ebe5..170147ba06 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/source/source_esolver/esolver_sdft_pw.cpp @@ -27,6 +27,9 @@ ESolver_SDFT_PW::ESolver_SDFT_PW() template ESolver_SDFT_PW::~ESolver_SDFT_PW() { + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** } template diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index 651c602822..e22dbfc3b0 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -37,6 +37,7 @@ list(APPEND objects cal_nelec_nband.cpp read_pseudo.cpp cal_wfc.cpp + setup_estate_pw.cpp ) if(ENABLE_LCAO) diff --git a/source/source_estate/setup_estate_pw.cpp b/source/source_estate/setup_estate_pw.cpp new file mode 100644 index 0000000000..51bd14e8b1 --- /dev/null +++ b/source/source_estate/setup_estate_pw.cpp @@ -0,0 +1,194 @@ +#include "source_estate/setup_estate_pw.h" +#include "source_estate/elecstate_pw.h" // init of pelec +#include "source_estate/elecstate_pw_sdft.h" // init of pelec for sdft +#include "source_estate/elecstate_tools.h" // occupations + +template +void elecstate::setup_estate_pw(UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K* pw_wfc, // pw for wfc + ModulePW::PW_Basis* pw_rho, // pw for rho + ModulePW::PW_Basis* pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp) // input parameters +{ + ModuleBase::TITLE("elecstate", "setup_estate_pw"); + + //! Initialize ElecState, set pelec pointer + if (pelec == nullptr) + { + if (inp.esolver_type == "sdft") + { + //! SDFT only supports double precision currently + pelec = new elecstate::ElecStatePW_SDFT, Device>(pw_wfc, + &chr, &kv, &ucell, &ppcell, + pw_rhod, pw_rho, pw_big); + } + else + { + pelec = new elecstate::ElecStatePW(pw_wfc, + &chr, &kv, &ucell, &ppcell, + pw_rhod, pw_rho, pw_big); + } + } + + //! Set the cell volume variable in pelec + pelec->omega = ucell.omega; + + //! Inititlize the charge density. + chr.allocate(inp.nspin); + + //! Initialize DFT-1/2 + if (PARAM.inp.dfthalf_type > 0) + { + vsep_cell = new VSep; + vsep_cell->init_vsep(*pw_rhod, ucell.sep_cell); + } + + //! Initialize the potential. + if (pelec->pot == nullptr) + { + pelec->pot = new elecstate::Potential(pw_rhod, + pw_rho, &ucell, &locpp.vloc, &sf, + &solvent, &(pelec->f_en.etxc), &(pelec->f_en.vtxc), vsep_cell); + } + + //! Initalize local pseudopotential + locpp.init_vloc(ucell, pw_rhod); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); + + //! Initalize non-local pseudopotential + ppcell.init(ucell, &sf, pw_wfc); + ppcell.init_vnl(ucell, pw_rhod); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); + + //! Setup occupations + if (inp.ocp) + { + elecstate::fixed_weights(inp.ocp_kb, + inp.nbands, + inp.nelec, + pelec->klist, + pelec->wg, + pelec->skip_weights); + } + + return; +} + + +template +void elecstate::teardown_estate_pw(elecstate::ElecState* &pelec, VSep* &vsep_cell) +{ + ModuleBase::TITLE("elecstate", "teardown_estate_pw"); + + if (vsep_cell != nullptr) + { + delete vsep_cell; + } + + // mohan update 20251005 to increase the security level + if (pelec != nullptr) + { + auto* pw_elec = dynamic_cast*>(pelec); + if (pw_elec) + { + delete pw_elec; + pelec = nullptr; + } + else + { + ModuleBase::WARNING_QUIT("elecstate::teardown_estate_pw", "Invalid ElecState type"); + } + } +} + + +template void elecstate::setup_estate_pw, base_device::DEVICE_CPU>( + UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + ModulePW::PW_Basis *pw_rho, // pw for rho + ModulePW::PW_Basis *pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp); // input parameters + +template void elecstate::setup_estate_pw, base_device::DEVICE_CPU>( + UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + ModulePW::PW_Basis *pw_rho, // pw for rho + ModulePW::PW_Basis *pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp); // input parameters + + +template void elecstate::teardown_estate_pw, base_device::DEVICE_CPU>( + elecstate::ElecState* &pelec, VSep* &vsep_cell); + +template void elecstate::teardown_estate_pw, base_device::DEVICE_CPU>( + elecstate::ElecState* &pelec, VSep* &vsep_cell); + + +#if ((defined __CUDA) || (defined __ROCM)) + +template void elecstate::setup_estate_pw, base_device::DEVICE_GPU>( + UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + ModulePW::PW_Basis *pw_rho, // pw for rho + ModulePW::PW_Basis *pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp); // input parameters + +template void elecstate::setup_estate_pw, base_device::DEVICE_GPU>( + UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + ModulePW::PW_Basis *pw_rho, // pw for rho + ModulePW::PW_Basis *pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp); // input parameters + +template void elecstate::teardown_estate_pw, base_device::DEVICE_GPU>( + elecstate::ElecState* &pelec, VSep* &vsep_cell); + +template void elecstate::teardown_estate_pw, base_device::DEVICE_GPU>( + elecstate::ElecState* &pelec, VSep* &vsep_cell); + +#endif diff --git a/source/source_estate/setup_estate_pw.h b/source/source_estate/setup_estate_pw.h new file mode 100644 index 0000000000..44864ad588 --- /dev/null +++ b/source/source_estate/setup_estate_pw.h @@ -0,0 +1,37 @@ +#ifndef SETUP_ESTATE_PW_H +#define SETUP_ESTATE_PW_H + +#include "source_base/module_device/device.h" // use Device +#include "source_cell/unitcell.h" +#include "source_cell/klist.h" +#include "source_pw/module_pwdft/structure_factor.h" +#include "source_estate/elecstate.h" +#include "source_pw/module_pwdft/VL_in_pw.h" +#include "source_pw/module_pwdft/VSep_in_pw.h" + +namespace elecstate +{ + +template +void setup_estate_pw(UnitCell& ucell, // unitcell + K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState* &pelec, // pointer of electrons + Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* &vsep_cell, // U-1/2 method + ModulePW::PW_Basis_K* pw_wfc, // pw for wfc + ModulePW::PW_Basis* pw_rho, // pw for rho + ModulePW::PW_Basis* pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* pw_big, // pw for big grid + surchem &solvent, // solvent + const Input_para& inp); // input parameters + +template +void teardown_estate_pw(elecstate::ElecState* &pelec, VSep* &vsep_cell); + +} + + +#endif diff --git a/source/source_lcao/module_deltaspin/init_sc.cpp b/source/source_lcao/module_deltaspin/init_sc.cpp index fbba82a839..ab192fe133 100644 --- a/source/source_lcao/module_deltaspin/init_sc.cpp +++ b/source/source_lcao/module_deltaspin/init_sc.cpp @@ -11,7 +11,7 @@ void spinconstrain::SpinConstrain::init_sc(double sc_thr_in, const UnitCell& ucell, Parallel_Orbitals* ParaV_in, int nspin_in, - K_Vectors& kv_in, + const K_Vectors& kv_in, void* p_hamilt_in, void* psi_in, elecstate::ElecState* pelec_in, @@ -34,4 +34,4 @@ void spinconstrain::SpinConstrain::init_sc(double sc_thr_in, } template class spinconstrain::SpinConstrain>; -template class spinconstrain::SpinConstrain; \ No newline at end of file +template class spinconstrain::SpinConstrain; diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp index 26cf70cea5..6213f528f7 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.cpp +++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp @@ -487,7 +487,7 @@ double SpinConstrain::get_sc_drop_thr() } template -void SpinConstrain::set_solver_parameters(K_Vectors& kv_in, +void SpinConstrain::set_solver_parameters(const K_Vectors& kv_in, void* p_hamilt_in, void* psi_in, elecstate::ElecState* pelec_in) @@ -603,4 +603,4 @@ void SpinConstrain::print_Mag_Force(std::ofstream& ofs_running) template class SpinConstrain>; template class SpinConstrain; -} // namespace spinconstrain \ No newline at end of file +} // namespace spinconstrain diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index 8f510db57c..bcc4ed9ec1 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -36,7 +36,7 @@ class SpinConstrain const UnitCell& ucell, Parallel_Orbitals* ParaV_in, int nspin_in, - K_Vectors& kv_in, + const K_Vectors& kv_in, void* p_hamilt_in, void* psi_in, elecstate::ElecState* pelec_in, @@ -201,7 +201,7 @@ class SpinConstrain /// @brief set orbital parallel info void set_ParaV(Parallel_Orbitals* ParaV_in); /// @brief set parameters for solver - void set_solver_parameters(K_Vectors& kv_in, + void set_solver_parameters(const K_Vectors& kv_in, void* p_hamilt_in, void* psi_in, elecstate::ElecState* pelec_in); diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index 8f49801c48..fdfe37828a 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -12,6 +12,8 @@ list(APPEND objects operator_pw/op_exx_pw.cpp operator_pw/exx_pw_ace.cpp operator_pw/exx_pw_pot.cpp + setup_pot.cpp + setup_pwrho.cpp forces_nl.cpp forces_cc.cpp forces_scc.cpp diff --git a/source/source_pw/module_pwdft/forces_scc.cpp b/source/source_pw/module_pwdft/forces_scc.cpp index 9279a0bbea..1484166b4e 100644 --- a/source/source_pw/module_pwdft/forces_scc.cpp +++ b/source/source_pw/module_pwdft/forces_scc.cpp @@ -75,63 +75,65 @@ void Forces::cal_force_scc(ModuleBase::matrix& forcescc, int igg0 = 0; const int ig0 = rho_basis->ig_gge0; - if (rho_basis->gg_uniq[0] < 1.0e-8) { - igg0 = 1; -} + if (rho_basis->gg_uniq[0] < 1.0e-8) + { + igg0 = 1; + } double fact = 2.0; - for (int nt = 0; nt < ucell_in.ntype; nt++) { - // Here we compute the G.ne.0 term - const int mesh = ucell_in.atoms[nt].ncpp.msh; - this->deriv_drhoc_scc(numeric, - mesh, - ucell_in.atoms[nt].ncpp.r.data(), - ucell_in.atoms[nt].ncpp.rab.data(), - ucell_in.atoms[nt].ncpp.rho_at.data(), - rhocgnt.data(), - rho_basis, - ucell_in); - int iat = 0; - for (int it = 0; it < ucell_in.ntype; it++) { - for (int ia = 0; ia < ucell_in.atoms[it].na; ia++) { - if (nt == it) { - const ModuleBase::Vector3 pos - = ucell_in.atoms[it].tau[ia]; - double &force0 = forcescc(iat, 0), - &force1 = forcescc(iat, 1), - &force2 = forcescc(iat, 2); + for (int nt = 0; nt < ucell_in.ntype; nt++) + { + // Here we compute the G.ne.0 term + const int mesh = ucell_in.atoms[nt].ncpp.msh; + this->deriv_drhoc_scc(numeric, + mesh, + ucell_in.atoms[nt].ncpp.r.data(), + ucell_in.atoms[nt].ncpp.rab.data(), + ucell_in.atoms[nt].ncpp.rho_at.data(), + rhocgnt.data(), + rho_basis, + ucell_in); + int iat = 0; + for (int it = 0; it < ucell_in.ntype; it++) { + for (int ia = 0; ia < ucell_in.atoms[it].na; ia++) { + if (nt == it) { + const ModuleBase::Vector3 pos + = ucell_in.atoms[it].tau[ia]; + double &force0 = forcescc(iat, 0), + &force1 = forcescc(iat, 1), + &force2 = forcescc(iat, 2); #ifdef _OPENMP #pragma omp parallel for reduction(+ : force0) reduction(+ : force1) reduction(+ : force2) #endif - for (int ig = 0; ig < rho_basis->npw; ++ig) { - if (ig == ig0) { - continue; -} - const ModuleBase::Vector3 gv - = rho_basis->gcar[ig]; - const double rhocgntigg - = rhocgnt[rho_basis->ig2igg[ig]]; - const double arg = ModuleBase::TWO_PI * (gv * pos); - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - const std::complex cpm - = std::complex(sinp, cosp) * conj(psic[ig]); - - force0 += fact * rhocgntigg * ucell_in.tpiba - * gv.x * cpm.real(); - force1 += fact * rhocgntigg * ucell_in.tpiba - * gv.y * cpm.real(); - force2 += fact * rhocgntigg * ucell_in.tpiba - * gv.z * cpm.real(); - } - } - iat++; - } - } - } + for (int ig = 0; ig < rho_basis->npw; ++ig) { + if (ig == ig0) { + continue; + } + const ModuleBase::Vector3 gv + = rho_basis->gcar[ig]; + const double rhocgntigg + = rhocgnt[rho_basis->ig2igg[ig]]; + const double arg = ModuleBase::TWO_PI * (gv * pos); + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + const std::complex cpm + = std::complex(sinp, cosp) * conj(psic[ig]); + + force0 += fact * rhocgntigg * ucell_in.tpiba + * gv.x * cpm.real(); + force1 += fact * rhocgntigg * ucell_in.tpiba + * gv.y * cpm.real(); + force2 += fact * rhocgntigg * ucell_in.tpiba + * gv.z * cpm.real(); + } + } + iat++; + } + } + } - Parallel_Reduce::reduce_pool(forcescc.c, forcescc.nr * forcescc.nc); + Parallel_Reduce::reduce_pool(forcescc.c, forcescc.nr * forcescc.nc); ModuleBase::timer::tick("Forces", "cal_force_scc"); return; @@ -163,10 +165,12 @@ void Forces::deriv_drhoc_scc(const bool& numeric, /// /// G=0 term /// - if (rho_basis->gg_uniq[0] < 1.0e-8) { - drhocg[0] = 0.0; + if (rho_basis->gg_uniq[0] < 1.0e-8) + { + drhocg[0] = 0.0; igl0 = 1; - } else { + } else + { igl0 = 0; } @@ -178,17 +182,19 @@ void Forces::deriv_drhoc_scc(const bool& numeric, #ifdef _OPENMP #pragma omp parallel for #endif - for (int igl = igl0; igl < rho_basis->ngg; igl++) { - gx_arr[igl] = sqrt(rho_basis->gg_uniq[igl]) * ucell_in.tpiba; - } + for (int igl = igl0; igl < rho_basis->ngg; igl++) + { + gx_arr[igl] = sqrt(rho_basis->gg_uniq[igl]) * ucell_in.tpiba; + } double *r_d = nullptr; double *rhoc_d = nullptr; double *rab_d = nullptr; double *aux_d = nullptr; double *drhocg_d = nullptr; - if (this->device == base_device::GpuDevice) { - resmem_var_op()(r_d, mesh); + if (this->device == base_device::GpuDevice) + { + resmem_var_op()(r_d, mesh); resmem_var_op()(rhoc_d, mesh); resmem_var_op()(rab_d, mesh); @@ -204,14 +210,36 @@ void Forces::deriv_drhoc_scc(const bool& numeric, syncmem_var_h2d_op()(rhoc_d, rhoc, mesh); } - if(this->device == base_device::GpuDevice) { + if(this->device == base_device::GpuDevice) + { hamilt::cal_stress_drhoc_aux_op()( - r_d,rhoc_d,gx_arr_d+igl0,rab_d,drhocg_d+igl0,mesh,igl0,rho_basis->ngg-igl0,ucell_in.omega,2); + r_d, + rhoc_d, + gx_arr_d+igl0, + rab_d, + drhocg_d+igl0, + mesh, + igl0, + rho_basis->ngg-igl0, + ucell_in.omega, + 2); + syncmem_var_d2h_op()(drhocg+igl0, drhocg_d+igl0, rho_basis->ngg-igl0); - } else { + } + else + { hamilt::cal_stress_drhoc_aux_op()( - r,rhoc,gx_arr.data()+igl0,rab,drhocg+igl0,mesh,igl0,rho_basis->ngg-igl0,ucell_in.omega,2); + r, + rhoc, + gx_arr.data()+igl0, + rab, + drhocg+igl0, + mesh, + igl0, + rho_basis->ngg-igl0, + ucell_in.omega, + 2); } @@ -226,4 +254,4 @@ void Forces::deriv_drhoc_scc(const bool& numeric, template class Forces; #if ((defined __CUDA) || (defined __ROCM)) template class Forces; -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp new file mode 100644 index 0000000000..b65d25cdb4 --- /dev/null +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -0,0 +1,196 @@ +#include "source_pw/module_pwdft/setup_pot.h" + +#include "source_estate/module_charge/symmetry_rho.h" +#include "source_lcao/module_deltaspin/spin_constrain.h" +#include "source_pw/module_pwdft/onsite_projector.h" +#include "source_lcao/module_dftu/dftu.h" +#include "source_pw/module_pwdft/VSep_in_pw.h" + +template +void pw::setup_pot(const int istep, + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi* kspw_psi, // electronic wave functions + hamilt::Hamilt* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp) // input parameters +{ + ModuleBase::TITLE("pw", "setup_pot"); + //---------------------------------------------------------- + //! 1) Renew local pseudopotential + //---------------------------------------------------------- + pelec->init_scf(istep, + ucell, + para_grid, + sf.strucFac, + locpp.numeric, + ucell.symm, + (void*)pw_wfc); + + //---------------------------------------------------------- + //! 2) Symmetrize the charge density (rho) + //---------------------------------------------------------- + + //! Symmetry_rho should behind init_scf, because charge should be + //! initialized first. liuyu comment: Symmetry_rho should be + //! located between init_rho and v_of_rho? + Symmetry_rho srho; + for (int is = 0; is < inp.nspin; is++) + { + srho.begin(is, chr, pw_rhod, ucell.symm); + } + + //---------------------------------------------------------- + //! 3) Calculate the effective potential with rho + //---------------------------------------------------------- + //! liuyu move here 2023-10-09 + //! D in uspp need vloc, thus behind init_scf() + //! calculate the effective coefficient matrix + //! for non-local pseudopotential projectors + ModuleBase::matrix veff = pelec->pot->get_effective_v(); + + ppcell.cal_effective_D(veff, pw_rhod, ucell); + + //---------------------------------------------------------- + //! 4) Onsite projectors + //---------------------------------------------------------- + if (PARAM.inp.onsite_radius > 0) + { + auto* onsite_p = projectors::OnsiteProjector::get_instance(); + onsite_p->init(PARAM.inp.orbital_dir, + &ucell, + *(kspw_psi), + kv, + *(pw_wfc), + sf, + PARAM.inp.onsite_radius, + PARAM.globalv.nqx, + PARAM.globalv.dq, + pelec->wg, + pelec->ekb); + } + + //---------------------------------------------------------- + //! 5) Spin-constrained algorithms + //---------------------------------------------------------- + if (PARAM.inp.sc_mag_switch) + { + spinconstrain::SpinConstrain>& sc + = spinconstrain::SpinConstrain>::getScInstance(); + sc.init_sc(PARAM.inp.sc_thr, + PARAM.inp.nsc, + PARAM.inp.nsc_min, + PARAM.inp.alpha_trial, + PARAM.inp.sccut, + PARAM.inp.sc_drop_thr, + ucell, + nullptr, + PARAM.inp.nspin, + kv, + p_hamilt, + kspw_psi, + pelec, + pw_wfc); + } + + //---------------------------------------------------------- + //! 6) DFT+U algorithm + //---------------------------------------------------------- + if (PARAM.inp.dft_plus_u) + { + auto* dftu = ModuleDFTU::DFTU::get_instance(); + dftu->init(ucell, nullptr, kv.get_nks()); + } + + //---------------------------------------------------------- + //! 7) DFT-1/2 calculations, sep potential need to generate + // before effective potential calculation + //---------------------------------------------------------- + if (PARAM.inp.dfthalf_type > 0) + { + vsep_cell->generate_vsep_r(pw_rhod[0], sf.strucFac, ucell.sep_cell); + } + + return; +} + +template void pw::setup_pot, base_device::DEVICE_CPU>( + const int istep, // ionic step + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions + hamilt::Hamilt, base_device::DEVICE_CPU>* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp); // input parameters + + +template void pw::setup_pot, base_device::DEVICE_CPU>( + const int istep, // ionic step + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi, base_device::DEVICE_CPU>* kspw_psi, // electronic wave functions + hamilt::Hamilt, base_device::DEVICE_CPU>* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp); // input parameters + +#if ((defined __CUDA) || (defined __ROCM)) + +template void pw::setup_pot, base_device::DEVICE_GPU>( + const int istep, // ionic step + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions + hamilt::Hamilt, base_device::DEVICE_GPU>* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp); // input parameters + +template void pw::setup_pot, base_device::DEVICE_GPU>( + const int istep, // ionic step + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi, base_device::DEVICE_GPU>* kspw_psi, // electronic wave functions + hamilt::Hamilt, base_device::DEVICE_GPU>* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp); // input parameters + +#endif diff --git a/source/source_pw/module_pwdft/setup_pot.h b/source/source_pw/module_pwdft/setup_pot.h new file mode 100644 index 0000000000..b89e33c21e --- /dev/null +++ b/source/source_pw/module_pwdft/setup_pot.h @@ -0,0 +1,36 @@ +#ifndef SETUP_POT_H +#define SETUP_POT_H + +#include "source_base/module_device/device.h" // use Device +#include "source_cell/unitcell.h" +#include "source_cell/klist.h" +#include "source_pw/module_pwdft/structure_factor.h" +#include "source_estate/elecstate.h" +#include "source_pw/module_pwdft/VL_in_pw.h" +#include "source_hamilt/hamilt.h" + +namespace pw +{ + +template +void setup_pot(const int istep, + UnitCell& ucell, // unitcell + const K_Vectors &kv, // kpoints + Structure_Factor &sf, // structure factors + elecstate::ElecState *pelec, // pointer of electrons + const Parallel_Grid ¶_grid, // parallel of FFT grids + const Charge &chr, // charge density + pseudopot_cell_vl &locpp, // local pseudopotentials + pseudopot_cell_vnl &ppcell, // non-local pseudopotentials + VSep* vsep_cell, // U-1/2 method + psi::Psi* kspw_psi, // electronic wave functions + hamilt::Hamilt* p_hamilt, // hamiltonian + ModulePW::PW_Basis_K *pw_wfc, // pw for wfc + const ModulePW::PW_Basis *pw_rhod, // pw for rhod + const Input_para& inp); // input parameters + +} + + + +#endif diff --git a/source/source_pw/module_pwdft/setup_pwrho.cpp b/source/source_pw/module_pwdft/setup_pwrho.cpp new file mode 100644 index 0000000000..e08ee71d23 --- /dev/null +++ b/source/source_pw/module_pwdft/setup_pwrho.cpp @@ -0,0 +1,141 @@ +#include "source_pw/module_pwdft/setup_pwrho.h" +#include "source_io/print_info.h" // use print_rhofft +#include "source_base/parallel_comm.h" // use POOL_WORLD + +void pw::setup_pwrho( + UnitCell& ucell, // unitcell + const bool double_grid, // for USPP + bool &pw_rho_flag, // flag for allocation of pw_rho + ModulePW::PW_Basis* &pw_rho, // pw for rhod + ModulePW::PW_Basis* &pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* &pw_big, // pw for rhod + const std::string &classname, + const Input_para& inp) // input parameters * +{ + ModuleBase::TITLE("pw", "setup_pwrho"); + + std::string fft_device = inp.device; + std::string fft_precision = inp.precision; + + // LCAO basis doesn't support GPU acceleration on FFT currently + if(inp.basis_type == "lcao") + { + fft_device = "cpu"; + } + + // single, double, or mixing precision calculations + if ((inp.precision=="single") || (inp.precision=="mixing")) + { + fft_precision = "mixing"; + } + else if (inp.precision=="double") + { + fft_precision = "double"; + } + + // for GPU +#if (not defined(__ENABLE_FLOAT_FFTW) and (defined(__CUDA) || defined(__RCOM))) + if (fft_device == "gpu") + { + fft_precision = "double"; + } +#endif + + // initialize pw_rho + pw_rho = new ModulePW::PW_Basis_Big(fft_device, fft_precision); + pw_rho_flag = true; + + // initialize pw_rhod + if (double_grid) + { + pw_rhod = new ModulePW::PW_Basis_Big(fft_device, fft_precision); + } + else + { + pw_rhod = pw_rho; + } + + // initialize pw_big + pw_big = static_cast(pw_rhod); + pw_big->setbxyz(inp.bx, inp.by, inp.bz); + + //! initialie the plane wave basis for rho +#ifdef __MPI + pw_rho->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); +#endif + + //! for OFDFT calculations + if (classname == "ESolver_OF" || inp.of_ml_gene_data == 1) + { + pw_rho->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); + } + + //! initialize the FFT grid + if (inp.nx * inp.ny * inp.nz == 0) + { + pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, 4.0 * inp.ecutwfc); + } + else + { + pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.nx, inp.ny, inp.nz); + } + + pw_rho->initparameters(false, 4.0 * inp.ecutwfc); + pw_rho->fft_bundle.initfftmode(inp.fft_mode); + pw_rho->setuptransform(); + pw_rho->collect_local_pw(); + pw_rho->collect_uniqgg(); + + //! initialize the double grid (for uspp) if necessary + if (double_grid) + { + ModulePW::PW_Basis_Sup* pw_rhod_sup = static_cast(pw_rhod); +#ifdef __MPI + pw_rhod->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); +#endif + if (classname == "ESolver_OF") + { + pw_rhod->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); + } + if (inp.ndx * inp.ndy * inp.ndz == 0) + { + pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ecutrho); + } + else + { + pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ndx, inp.ndy, inp.ndz); + } + pw_rhod->initparameters(false, inp.ecutrho); + pw_rhod->fft_bundle.initfftmode(inp.fft_mode); + pw_rhod_sup->setuptransform(pw_rho); + pw_rhod->collect_local_pw(); + pw_rhod->collect_uniqgg(); + } + + ModuleIO::print_rhofft(pw_rhod, pw_rho, pw_big, GlobalV::ofs_running); + + return; +} + + +void pw::teardown_pwrho(bool &pw_rho_flag, + const bool double_grid, + ModulePW::PW_Basis* &pw_rho, // pw for rhod + ModulePW::PW_Basis* &pw_rhod) // pw for rhod +{ + if (pw_rho_flag == true) + { + delete pw_rho; + pw_rho = nullptr; + pw_rho_flag = false; + } + + if (double_grid == true) + { + delete pw_rhod; + pw_rhod = nullptr; + } + + return; +} + diff --git a/source/source_pw/module_pwdft/setup_pwrho.h b/source/source_pw/module_pwdft/setup_pwrho.h new file mode 100644 index 0000000000..820edb9d2d --- /dev/null +++ b/source/source_pw/module_pwdft/setup_pwrho.h @@ -0,0 +1,33 @@ +#ifndef SETUP_PWRHO_H +#define SETUP_PWRHO_H + +#include "source_base/module_device/device.h" // use Device +#include "source_cell/unitcell.h" // use UnitCell +#include "source_pw/module_pwdft/structure_factor.h" // use Structure_Factor +#include "source_basis/module_pw/pw_basis.h" // use PW_Basis +#include "source_io/module_parameter/input_parameter.h" // use Input_para + +namespace pw +{ + +void setup_pwrho( + UnitCell& ucell, // unitcell + const bool double_grid, // for USPP + bool &pw_rho_flag, // flag for allocation of pw_rho + ModulePW::PW_Basis* &pw_rho, // pw for rhod + ModulePW::PW_Basis* &pw_rhod, // pw for rhod + ModulePW::PW_Basis_Big* &pw_big, // pw for rhod + const std::string &classname, + const Input_para& inp); // input parameters * + + +void teardown_pwrho(bool &pw_rho_flag, + const bool double_grid, + ModulePW::PW_Basis* &pw_rho, // pw for rhod + ModulePW::PW_Basis* &pw_rhod); // pw for rhod + +} + + + +#endif diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index 027bc56039..a9ac9decb5 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -1,5 +1,5 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_source_pw_HAMILT_PWDFT_STRESS_PW_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_source_pw_HAMILT_PWDFT_STRESS_PW_H +#ifndef STRESS_PW_H +#define STRESS_PW_H #include "source_estate/elecstate.h" #include "source_pw/module_pwdft/VL_in_pw.h"