diff --git a/source/Makefile.Objects b/source/Makefile.Objects index b64cc03f89..51c670dca8 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -259,6 +259,7 @@ OBJS_ESOLVER=esolver.o\ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ esolver_ks_lcao_tddft.o\ lcao_before_scf.o\ + lcao_after_scf.o\ esolver_gets.o\ lcao_others.o\ diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 395350f1d5..6618006588 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -17,6 +17,7 @@ if(ENABLE_LCAO) esolver_ks_lcao.cpp esolver_ks_lcao_tddft.cpp lcao_before_scf.cpp + lcao_after_scf.cpp esolver_gets.cpp lcao_others.cpp ) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 989e4cbfba..bd65cbe03d 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -127,15 +127,15 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp) } //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. -void ESolver_FP::after_scf(UnitCell& ucell, const int istep) +void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { ModuleBase::TITLE("ESolver_FP", "after_scf"); // 0) output convergence information - ModuleIO::output_convergence_after_scf(this->conv_esolver, this->pelec->f_en.etot); + ModuleIO::output_convergence_after_scf(conv_esolver, this->pelec->f_en.etot); // 1) write fermi energy - ModuleIO::output_efermi(this->conv_esolver, this->pelec->eferm.ef); + ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef); // 2) update delta rho for charge extrapolation CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); @@ -286,12 +286,12 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) return; } -void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter) +void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { //! output charge density if (PARAM.inp.out_chg[0] != -1) { - if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || this->conv_esolver) + 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; diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 355d6af5b5..9c95dd5b2d 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -39,10 +39,10 @@ class ESolver_FP : public ESolver virtual void before_scf(UnitCell& ucell, const int istep); //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void after_scf(UnitCell& ucell, const int istep); + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver); //! Something to do after hamilt2density function in each iter loop. - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter); + 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. diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 9ee009193c..50c5a1facf 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -437,7 +437,7 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); // 4) SCF iterations - this->conv_esolver = false; + bool conv_esolver = false; this->niter = this->maxniter; this->diag_ethr = PARAM.inp.pw_diag_thr; for (int iter = 1; iter <= this->maxniter; ++iter) @@ -449,10 +449,10 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) this->hamilt2density(ucell, istep, iter, diag_ethr); // 7) finish scf iterations - this->iter_finish(ucell, istep, iter); + this->iter_finish(ucell, istep, iter, conv_esolver); // 8) check convergence - if (this->conv_esolver || this->oscillate_esolver) + if (conv_esolver || this->oscillate_esolver) { this->niter = iter; if (this->oscillate_esolver) @@ -464,7 +464,7 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) } // end scf iterations // 9) after scf - this->after_scf(ucell, istep); + this->after_scf(ucell, istep, conv_esolver); ModuleBase::timer::tick(this->classname, "runner"); return; @@ -524,7 +524,7 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const in } template -void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter) +void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver) { if (PARAM.inp.out_bandgap) { @@ -578,30 +578,30 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i } #endif - this->conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged); + conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged); // add energy threshold for SCF convergence if (this->scf_ene_thr > 0.0) { // calculate energy of output charge density - this->update_pot(ucell, istep, iter); + this->update_pot(ucell, istep, iter, conv_esolver); this->pelec->cal_energies(2); // 2 means Kohn-Sham functional // now, etot_old is the energy of input density, while etot is the energy of output density this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; // output etot_delta GlobalV::ofs_running << " DeltaE_womix = " << this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV" << std::endl; - if (iter > 1 && this->conv_esolver == 1) // only check when density is converged + if (iter > 1 && conv_esolver == 1) // only check when density is converged { // update the convergence flag - this->conv_esolver + conv_esolver = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr); } } // If drho < hsolver_error in the first iter or drho < scf_thr, we // do not change rho. - if (drho < hsolver_error || this->conv_esolver || PARAM.inp.calculation == "nscf") + if (drho < hsolver_error || conv_esolver || PARAM.inp.calculation == "nscf") { if (drho < hsolver_error) { @@ -635,14 +635,16 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i #ifdef __MPI MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, BP_WORLD); - MPI_Bcast(&this->conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD); + + // 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); #endif // update potential // Hamilt should be used after it is constructed. // this->phamilt->update(conv_esolver); - this->update_pot(ucell, istep, iter); + this->update_pot(ucell, istep, iter, conv_esolver); // 1 means Harris-Foulkes functional // 2 means Kohn-Sham functional @@ -671,7 +673,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec); } - this->pelec->print_etot(ucell.magnet,this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr); + this->pelec->print_etot(ucell.magnet,conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr); // Json, need to be moved to somewhere else #ifdef __RAPIDJSON @@ -691,17 +693,17 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i std::cout << " SCF restart after this step!" << std::endl; } - ESolver_FP::iter_finish(ucell, istep, iter); + ESolver_FP::iter_finish(ucell, istep, iter, conv_esolver); } //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. template -void ESolver_KS::after_scf(UnitCell& ucell, const int istep) +void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { ModuleBase::TITLE("ESolver_KS", "after_scf"); // 1) call after_scf() of ESolver_FP - ESolver_FP::after_scf(ucell, istep); + ESolver_FP::after_scf(ucell, istep, conv_esolver); // 2) write eigenvalues if (istep % PARAM.inp.out_interval == 0) diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 00e1352c50..652c193d5c 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP virtual void iter_init(UnitCell& ucell, const int istep, const int iter); //! Something to do after hamilt2density function in each iter loop. - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; // calculate electron density from a specific Hamiltonian with ethr virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr); @@ -52,10 +52,10 @@ class ESolver_KS : public ESolver_FP void hamilt2density(UnitCell& ucell, const int istep, const int iter, const double ethr); //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void after_scf(UnitCell& ucell, const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; //! It should be replaced by a function in Hamilt Class - virtual void update_pot(UnitCell& ucell, const int istep, const int iter){}; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver){}; //! Hamiltonian hamilt::Hamilt* p_hamilt = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index ddefbef490..3f21c2e33b 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -22,14 +22,18 @@ #include "module_io/to_wannier90_lcao_in_pw.h" #include "module_io/write_HS.h" #include "module_io/write_dmr.h" -#include "module_io/write_eband_terms.hpp" #include "module_io/write_elecstat_pot.h" #include "module_io/write_istate_info.h" #include "module_io/write_proj_band_lcao.h" -#include "module_io/write_vxc.hpp" #include "module_io/write_wfc_nao.h" #include "module_parameter/parameter.h" + +//be careful of hpp, there may be multiple definitions of functions, 20250302, mohan +#include "module_io/write_eband_terms.hpp" +#include "module_io/write_vxc.hpp" +#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" + //--------------temporary---------------------------- #include "module_base/global_function.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" @@ -37,16 +41,11 @@ #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need DeePKS_init -#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" #include -#ifdef __EXX -#include "module_io/restart_exx_csr.h" -#include "module_ri/RPA_LRI.h" -#endif #ifdef __DEEPKS #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" @@ -59,9 +58,6 @@ #include "module_elecstate/elecstate_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hsolver/hsolver_lcao.h" -// function used by deepks -// #include "module_elecstate/cal_dm.h" -//--------------------------------------------------- // test RDMFT #include "module_rdmft/rdmft.h" @@ -776,11 +772,11 @@ void ESolver_KS_LCAO::hamilt2density_single(UnitCell& ucell, int istep, //! 1) print potential //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const int iter) +void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) { ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); - if (!this->conv_esolver) + if (!conv_esolver) { elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); @@ -801,7 +797,7 @@ void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const //! 4) output charge density and density matrix //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& iter) +void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); @@ -850,7 +846,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& } // call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(ucell, istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); // 1) mix density matrix if mixing_restart + mixing_dmr + not first // mixing_restart at every iter @@ -884,7 +880,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& this->scf_ene_thr, iter, istep, - this->conv_esolver) + conv_esolver) : this->exc->exx_iter_finish(this->kv, ucell, *this->p_hamilt, @@ -893,410 +889,18 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& this->scf_ene_thr, iter, istep, - this->conv_esolver); + conv_esolver); } } #endif // 6) use the converged occupation matrix for next MD/Relax SCF calculation - if (PARAM.inp.dft_plus_u && this->conv_esolver) + if (PARAM.inp.dft_plus_u && conv_esolver) { GlobalC::dftu.initialed_locale = true; } } -//------------------------------------------------------------------------------ -//! the 14th function of ESolver_KS_LCAO: after_scf -//! mohan add 2024-05-11 -//! 1) call after_scf() of ESolver_KS -//! 2) write density matrix for sparse matrix -//! 4) write density matrix -//! 5) write Exx matrix -//! 6) write Hamiltonian and Overlap matrix -//! 7) write wavefunctions -//! 11) write deepks information -//! 12) write rpa information -//! 13) write HR in npz format -//! 14) write dm in npz format -//! 15) write md related -//! 16) write spin constrian MW? -//! 17) delete grid -//! 18) write quasi-orbitals -//------------------------------------------------------------------------------ -template -void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); - - // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 - if (PARAM.inp.out_elf[0] > 0) - { - assert(this->psi != nullptr); - this->pelec->cal_tau(*(this->psi)); - } - - //! 2) call after_scf() of ESolver_KS - ESolver_KS::after_scf(ucell, istep); - - //! 3) write density matrix for sparse matrix - ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), - this->pv, - PARAM.inp.out_dm1, - false, - PARAM.inp.out_app_flag, - ucell.get_iat2iwt(), - &ucell.nat, - istep); - - //! 4) write density matrix - if (PARAM.inp.out_dm) - { - std::vector efermis(PARAM.inp.nspin == 2 ? 2 : 1); - for (int ispin = 0; ispin < efermis.size(); ispin++) - { - efermis[ispin] = this->pelec->eferm.get_efval(ispin); - } - const int precision = 3; - ModuleIO::write_dmk(dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(), - precision, - efermis, - &(ucell), - this->pv); - } - -#ifdef __EXX - //! 5) write Hexx matrix for NSCF (see `out_chg` in docs/advanced/input_files/input-main.md) - if (PARAM.inp.calculation != "nscf") - { - if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.out_chg[0] - && istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14 - { - const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) - { - ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exd->get_Hexxs()); - } - else - { - ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exc->get_Hexxs()); - } - } - } -#endif - - // 6) write Hamiltonian and Overlap matrix - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (PARAM.inp.out_mat_hs[0]) - { - this->p_hamilt->updateHk(ik); - } - bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core - // calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) - { - hamilt::MatrixBlock h_mat; - hamilt::MatrixBlock s_mat; - - this->p_hamilt->matrix(h_mat, s_mat); - - if (PARAM.inp.out_mat_hs[0]) - { - ModuleIO::save_mat(istep, - h_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "H", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - ModuleIO::save_mat(istep, - s_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "S", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - } - } - } - - // 7) write wavefunctions - if (elecstate::ElecStateLCAO::out_wfc_lcao && (istep % PARAM.inp.out_interval == 0)) - { - ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO::out_wfc_lcao, - this->psi[0], - this->pelec->ekb, - this->pelec->wg, - this->pelec->klist->kvec_c, - this->pv, - istep); - } - - //! 8) Write DeePKS information -#ifdef __DEEPKS - if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) - { - hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); - std::shared_ptr ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); - LCAO_Deepks_Interface LDI(ld_shared_ptr); - - LDI.out_deepks_labels(this->pelec->f_en.etot, - this->pelec->klist->get_nks(), - ucell.nat, - PARAM.globalv.nlocal, - this->pelec->ekb, - this->pelec->klist->kvec_d, - ucell, - orb_, - this->gd, - &(this->pv), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM(), - p_ham_deepks); - } -#endif - - //! 9) Perform RDMFT calculations - /******** test RDMFT *********/ - if (PARAM.inp.rdmft == true) // rdmft, added by jghan, 2024-10-17 - { - ModuleBase::matrix occ_number_ks(this->pelec->wg); - for (int ik = 0; ik < occ_number_ks.nr; ++ik) - { - for (int inb = 0; inb < occ_number_ks.nc; ++inb) - { - occ_number_ks(ik, inb) /= this->kv.wk[ik]; - } - } - this->rdmft_solver.update_elec(ucell, occ_number_ks, *(this->psi)); - - //! initialize the gradients of Etotal with respect to occupation numbers and wfc, - //! and set all elements to 0. - ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true); - psi::Psi dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis(), this->kv.ngk, true); - dE_dWfc.zero_out(); - - double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc); - } - /******** test RDMFT *********/ - -#ifdef __EXX - // 10) Write RPA information. - if (PARAM.inp.rpa) - { - // ModuleRPA::DFT_RPA_interface - // rpa_interface(GlobalC::exx_info.info_global); - RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); - rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), - MPI_COMM_WORLD, - ucell, - this->kv, - orb_); - rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); - rpa_lri_double.out_for_RPA(ucell, this->pv, *(this->psi), this->pelec); - } -#endif - - // 11) write HR in npz format. - if (PARAM.inp.out_hr_npz) - { - this->p_hamilt->updateHk(0); // first k point, up spin - hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); - std::string zipname = "output_HR0.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); - - if (PARAM.inp.nspin == 2) - { - this->p_hamilt->updateHk(this->kv.get_nks() / 2); // the other half of k points, down spin - hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); - zipname = "output_HR1.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); - } - } - - // 12) write density matrix in the 'npz' format. - if (PARAM.inp.out_dm_npz) - { - const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - std::string zipname = "output_DM0.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(1))); - - if (PARAM.inp.nspin == 2) - { - zipname = "output_DM1.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(2))); - } - } - - //! 13) Print out information every 'out_interval' steps. - if (PARAM.inp.calculation != "md" || istep % PARAM.inp.out_interval == 0) - { - //! Print out sparse matrix - ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2, - PARAM.inp.out_mat_dh, - PARAM.inp.out_mat_t, - PARAM.inp.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->pv, - this->GK, - two_center_bundle_, - orb_, - ucell, - this->gd, - this->kv, - this->p_hamilt); - - //! Perform Mulliken charge analysis - if (PARAM.inp.out_mul) - { - ModuleIO::cal_mag(&(this->pv), - this->p_hamilt, - this->kv, - this->pelec, - this->two_center_bundle_, - this->orb_, - ucell, - this->gd, - istep, - true); - } - } - - //! 14) Print out atomic magnetization only when 'spin_constraint' is on. - if (PARAM.inp.sc_mag_switch) - { - spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); - sc.cal_mi_lcao(istep); - sc.print_Mi(GlobalV::ofs_running); - sc.print_Mag_Force(GlobalV::ofs_running); - } - - //! 15) Clean up RA. - //! this should be last function and put it in the end, mohan request 2024-11-28 - if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress) - { - RA.delete_grid(); - } - - //! 16) Print out quasi-orbitals. - if (PARAM.inp.qo_switch) - { - toQO tqo(PARAM.inp.qo_basis, PARAM.inp.qo_strategy, PARAM.inp.qo_thr, PARAM.inp.qo_screening_coeff); - tqo.initialize(PARAM.globalv.global_out_dir, - PARAM.inp.pseudo_dir, - PARAM.inp.orbital_dir, - &ucell, - this->kv.kvec_d, - GlobalV::ofs_running, - GlobalV::MY_RANK, - GlobalV::NPROC); - tqo.calculate(); - } - - //! 17) Print out kinetic matrix. - if (PARAM.inp.out_mat_tk[0]) - { - hamilt::HS_Matrix_K hsk(&pv, true); - hamilt::HContainer hR(&pv); - hamilt::Operator* ekinetic - = new hamilt::EkineticNew>(&hsk, - this->kv.kvec_d, - &hR, - &ucell, - orb_.cutoffs(), - &this->gd, - two_center_bundle_.kinetic_orb.get()); - - const int nspin_k = (PARAM.inp.nspin == 2 ? 2 : 1); - for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik) - { - ekinetic->init(ik); - ModuleIO::save_mat(0, - hsk.get_hk(), - PARAM.globalv.nlocal, - false, - PARAM.inp.out_mat_tk[1], - 1, - PARAM.inp.out_app_flag, - "T", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - } - - // where is new? mohan ask 2024-11-28 - delete ekinetic; - } - - //! 18) Wannier 90 function, added by jingan in 2018.11.7 - if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) - { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wave function to Wannier90"); - if (PARAM.inp.wannier_method == 1) - { - toWannier90_LCAO_IN_PW myWannier(PARAM.inp.out_wannier_mmn, - PARAM.inp.out_wannier_amn, - PARAM.inp.out_wannier_unk, - PARAM.inp.out_wannier_eig, - PARAM.inp.out_wannier_wvfn_formatted, - PARAM.inp.nnkpfile, - PARAM.inp.wannier_spin); - myWannier.set_tpiba_omega(ucell.tpiba, ucell.omega); - myWannier.calculate(ucell, - this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->sf, - this->kv, - this->psi, - &(this->pv)); - } - else if (PARAM.inp.wannier_method == 2) - { - toWannier90_LCAO myWannier(PARAM.inp.out_wannier_mmn, - PARAM.inp.out_wannier_amn, - PARAM.inp.out_wannier_unk, - PARAM.inp.out_wannier_eig, - PARAM.inp.out_wannier_wvfn_formatted, - PARAM.inp.nnkpfile, - PARAM.inp.wannier_spin, - orb_); - - myWannier.calculate(ucell, this->gd, this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); - } - - //! 19) berry phase calculations, added by jingan - if (PARAM.inp.calculation == "nscf" && berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) - { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation"); - berryphase bp(&(this->pv)); - bp.lcao_init(ucell, this->gd, this->kv, this->GridT, orb_); - // additional step before calling - // macroscopic_polarization (why capitalize - // the function name?) - bp.Macroscopic_polarization(ucell, this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); - } - - ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); -} - template class ESolver_KS_LCAO; template class ESolver_KS_LCAO, double>; template class ESolver_KS_LCAO, std::complex>; diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 8a620291f7..4b4a3e530d 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -20,6 +20,7 @@ // added by jghan for rdmft calculation #include "module_rdmft/rdmft.h" + #include namespace LR @@ -53,11 +54,11 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; - virtual void after_scf(UnitCell& ucell, const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; virtual void others(UnitCell& ucell, const int istep) override; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 11792f6a2b..46d456c3f1 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -163,7 +163,11 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(UnitCell& ucell, } template -void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep, int& iter) +void ESolver_KS_LCAO_TDDFT::iter_finish( + UnitCell& ucell, + const int istep, + int& iter, + bool& conv_esolver) { // print occupation of each band if (iter == 1 && istep <= 2) @@ -189,14 +193,17 @@ void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep << std::endl; } - ESolver_KS_LCAO, double>::iter_finish(ucell, istep, iter); + ESolver_KS_LCAO, double>::iter_finish(ucell, istep, iter, conv_esolver); } template -void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const int iter) +void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, + const int istep, + const int iter, + const bool conv_esolver) { // Calculate new potential according to new Charge Density - if (!this->conv_esolver) + if (!conv_esolver) { elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); @@ -214,7 +221,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const int nlocal = PARAM.globalv.nlocal; // store wfc and Hk laststep - if (istep >= (PARAM.inp.init_wfc == "file" ? 0 : 1) && this->conv_esolver) + if (istep >= (PARAM.inp.init_wfc == "file" ? 0 : 1) && conv_esolver) { if (this->psi_laststep == nullptr) { @@ -305,7 +312,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, } // print "eigen value" for tddft - if (this->conv_esolver) + if (conv_esolver) { GlobalV::ofs_running << "---------------------------------------------------------------" "---------------------------------" @@ -331,10 +338,10 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, } template -void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) +void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { - ModuleBase::TITLE("ESolver_KS_LCAO_TDDFT", "after_scf"); - ModuleBase::timer::tick("ESolver_KS_LCAO_TDDFT", "after_scf"); + ModuleBase::TITLE("ESolver_LCAO_TDDFT", "after_scf"); + ModuleBase::timer::tick("ESolver_LCAO_TDDFT", "after_scf"); for (int is = 0; is < PARAM.inp.nspin; is++) { @@ -366,9 +373,9 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) orb_, this->RA); } - ESolver_KS_LCAO, double>::after_scf(ucell, istep); + ESolver_KS_LCAO, double>::after_scf(ucell, istep, conv_esolver); - ModuleBase::timer::tick("ESolver_KS_LCAO_TDDFT", "after_scf"); + ModuleBase::timer::tick("ESolver_LCAO_TDDFT", "after_scf"); } template diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index c5fdcbc497..9232590fbf 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -60,11 +60,11 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl protected: virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; - virtual void after_scf(UnitCell& ucell, const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; //! wave functions of last time step psi::Psi>* psi_laststep = nullptr; diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index 3f7296395b..29f6bfb94b 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -175,12 +175,12 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::iter_finish(UnitCell& ucell, const int istep, int& iter) + void ESolver_KS_LIP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { - ESolver_KS_PW::iter_finish(ucell, istep, iter); + ESolver_KS_PW::iter_finish(ucell, istep, iter, conv_esolver); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx && this->conv_esolver) + if (GlobalC::exx_info.info_global.cal_exx && conv_esolver) { // no separate_loop case if (!GlobalC::exx_info.info_global.separate_loop) @@ -198,7 +198,7 @@ namespace ModuleESolver iter = 0; std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; this->two_level_step++; - this->conv_esolver = false; + conv_esolver = false; } } // has separate_loop case @@ -206,7 +206,7 @@ namespace ModuleESolver else if (this->two_level_step == GlobalC::exx_info.info_global.hybrid_step || (iter == 1 && this->two_level_step != 0)) { - this->conv_esolver = true; + conv_esolver = true; } else { @@ -230,7 +230,7 @@ namespace ModuleESolver << (double)(t_end.tv_sec - t_start.tv_sec) + (double)(t_end.tv_usec - t_start.tv_usec) / 1000000.0 << std::defaultfloat << " (s)" << std::endl; - this->conv_esolver = false; + conv_esolver = false; } } #endif diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index 99f527081e..c73b74a38c 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -27,7 +27,7 @@ namespace ModuleESolver protected: virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; /// All the other interfaces except this one are the same as ESolver_KS_PW. virtual void hamilt2density_single(UnitCell& ucell, diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 74d3904a65..15c0bd740b 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -9,9 +9,7 @@ #include "module_hamilt_general/module_ewald/H_Ewald_pw.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" -//-----force------------------- #include "module_hamilt_pw/hamilt_pwdft/forces.h" -//-----stress------------------ #include "module_hamilt_pw/hamilt_pwdft/stress_pw.h" //--------------------------------------------------- #include "module_base/formatter.h" @@ -129,6 +127,7 @@ void ESolver_KS_PW::allocate_hamilt(const UnitCell& ucell) { this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell); } + template void ESolver_KS_PW::deallocate_hamilt() { @@ -282,21 +281,29 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) 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); } + //---------------------------------------------------------- //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho + //---------------------------------------------------------- elecstate::cal_ux(ucell); + //---------------------------------------------------------- //! calculate the total local pseudopotential in real space + //---------------------------------------------------------- this->pelec ->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm, (void*)this->pw_wfc); + //---------------------------------------------------------- //! output the initial charge density + //---------------------------------------------------------- if (PARAM.inp.out_chg[0] == 2) { for (int is = 0; is < PARAM.inp.nspin; is++) @@ -314,7 +321,9 @@ void ESolver_KS_PW::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++) @@ -334,19 +343,23 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) } } + //---------------------------------------------------------- //! 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->pelec->charge), this->pw_rhod, ucell.symm); } + //---------------------------------------------------------- // 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); @@ -555,9 +568,9 @@ void ESolver_KS_PW::hamilt2density_single(UnitCell& ucell, // Temporary, it should be rewritten with Hamilt class. template -void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, const int iter) +void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) { - if (!this->conv_esolver) + if (!conv_esolver) { elecstate::cal_ux(ucell); this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); @@ -573,15 +586,14 @@ 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) +void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { // 1) Call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(ucell, istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); // 2) Update USPP-related quantities - // D in uspp need vloc, thus needs update when veff updated - // calculate the effective coefficient matrix for non-local pseudopotential - // projectors + // 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) { @@ -589,20 +601,20 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); } - // 4) Print out electronic wavefunctions + // 3) Print out electronic wavefunctions in pw basis if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2) { - if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || this->conv_esolver) + if (iter % PARAM.inp.out_freq_elec == 0 || + iter == PARAM.inp.scf_nmax || + conv_esolver) { std::stringstream ssw; ssw << PARAM.globalv.global_out_dir << "WAVEFUNC"; - // mohan update 2011-02-21 // qianrui update 2020-10-17 ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc); - // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"write wave - // functions into file WAVEFUNC.dat"); } } + // 4) check if oscillate for delta_spin method if (PARAM.inp.sc_mag_switch) { @@ -621,7 +633,7 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int } template -void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) +void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { ModuleBase::TITLE("ESolver_KS_PW", "after_scf"); ModuleBase::timer::tick("ESolver_KS_PW", "after_scf"); @@ -633,9 +645,9 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) } // 2) call after_scf() of ESolver_KS - ESolver_KS::after_scf(ucell, istep); + ESolver_KS::after_scf(ucell, istep, conv_esolver); - // 3) output wavefunctions + // 3) output wavefunctions in pw basis if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2) { std::stringstream ssw; @@ -643,7 +655,8 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc); } - // 4) Transfer data from GPU to CPU + // 4) transfer data from GPU to CPU in pw basis + // a question: the wavefunctions have been output, then the data transfer occurs? mohan 20250302 if (this->device == base_device::GpuDevice) { castmem_2d_d2h_op()(this->psi[0].get_pointer() - this->psi[0].get_psi_bias(), @@ -651,7 +664,7 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) this->psi[0].size()); } - // 5) Calculate band-decomposed (partial) charge density + // 5) calculate band-decomposed (partial) charge density in pw basis const std::vector bands_to_print = PARAM.inp.bands_to_print; if (bands_to_print.size() > 0) { @@ -678,7 +691,7 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) PARAM.inp.if_separate_k); } - //! 6) Calculate Wannier functions + //! 6) calculate Wannier functions in PW basis if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wannier functions calculation"); @@ -734,6 +747,7 @@ template void ESolver_KS_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { Forces ff(ucell.nat); + if (this->__kspw_psi != nullptr && PARAM.inp.precision == "single") { delete reinterpret_cast, Device>*>(this->__kspw_psi); @@ -763,6 +777,7 @@ template void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) { Stress_PW ss(this->pelec); + if (this->__kspw_psi != nullptr && PARAM.inp.precision == "single") { delete reinterpret_cast, Device>*>(this->__kspw_psi); @@ -930,12 +945,31 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) { this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); - ML_data ml_data; - ml_data.set_para(this->pelec->charge->nrxx, PARAM.inp.nelec, PARAM.inp.of_tf_weight, PARAM.inp.of_vw_weight, - PARAM.inp.of_ml_chi_p, PARAM.inp.of_ml_chi_q, PARAM.inp.of_ml_chi_xi, PARAM.inp.of_ml_chi_pnl, PARAM.inp.of_ml_chi_qnl, - PARAM.inp.of_ml_nkernel, PARAM.inp.of_ml_kernel, PARAM.inp.of_ml_kernel_scaling, PARAM.inp.of_ml_yukawa_alpha, PARAM.inp.of_ml_kernel_file, ucell.omega, this->pw_rho); - ml_data.generateTrainData_KS(this->kspw_psi, this->pelec, this->pw_wfc, this->pw_rho, ucell, this->pelec->pot->get_effective_v(0)); - } + ML_data ml_data; + ml_data.set_para(this->pelec->charge->nrxx, + PARAM.inp.nelec, + PARAM.inp.of_tf_weight, + PARAM.inp.of_vw_weight, + PARAM.inp.of_ml_chi_p, + PARAM.inp.of_ml_chi_q, + PARAM.inp.of_ml_chi_xi, + PARAM.inp.of_ml_chi_pnl, + PARAM.inp.of_ml_chi_qnl, + PARAM.inp.of_ml_nkernel, + PARAM.inp.of_ml_kernel, + PARAM.inp.of_ml_kernel_scaling, + PARAM.inp.of_ml_yukawa_alpha, + PARAM.inp.of_ml_kernel_file, + ucell.omega, + this->pw_rho); + + ml_data.generateTrainData_KS(this->kspw_psi, + this->pelec, + this->pw_wfc, + this->pw_rho, + ucell, + this->pelec->pot->get_effective_v(0)); + } #endif } diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index 25cf73dc77..47ddb82e2a 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -36,11 +36,11 @@ class ESolver_KS_PW : public ESolver_KS virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; - virtual void after_scf(UnitCell& ucell, const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; virtual void others(UnitCell& ucell, const int istep) override; diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index fcf326543e..9160a735f1 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -164,6 +164,9 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep) if (PARAM.inp.of_ml_local_test) this->ml_->localTest(pelec->charge->rho, this->pw_rho); #endif + + bool conv_esolver = false; // this conv_esolver is added by mohan 20250302 + while (true) { // once we get a new rho and phi, update potential @@ -174,8 +177,9 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep) this->energy_last_ = this->energy_current_; this->energy_current_ = this->cal_energy(); + // check if the job is done - if (this->check_exit()) + if (this->check_exit(conv_esolver)) { break; } @@ -188,10 +192,10 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep) this->iter_++; - ESolver_FP::iter_finish(ucell, istep, this->iter_); + ESolver_FP::iter_finish(ucell, istep, this->iter_, conv_esolver); } - this->after_opt(istep, ucell); + this->after_opt(istep, ucell, conv_esolver); ModuleBase::timer::tick("ESolver_OF", "runner"); } @@ -434,9 +438,9 @@ void ESolver_OF::update_rho() * * @return exit or not */ -bool ESolver_OF::check_exit() +bool ESolver_OF::check_exit(bool& conv_esolver) { - this->conv_esolver = false; + conv_esolver = false; bool potConv = false; bool potHold = false; // if normdLdphi nearly remains unchanged bool energyConv = false; @@ -457,12 +461,12 @@ bool ESolver_OF::check_exit() energyConv = true; } - this->conv_esolver = (this->of_conv_ == "energy" && energyConv) || (this->of_conv_ == "potential" && potConv) + conv_esolver = (this->of_conv_ == "energy" && energyConv) || (this->of_conv_ == "potential" && potConv) || (this->of_conv_ == "both" && potConv && energyConv); - this->print_info(); + this->print_info(conv_esolver); - if (this->conv_esolver || this->iter_ >= this->max_iter_) + if (conv_esolver || this->iter_ >= this->max_iter_) { return true; } @@ -488,7 +492,7 @@ bool ESolver_OF::check_exit() * @param istep * @param ucell */ -void ESolver_OF::after_opt(const int istep, UnitCell& ucell) +void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_esolver) { ModuleBase::TITLE("ESolver_OF", "after_opt"); ModuleBase::timer::tick("ESolver_OF", "after_opt"); @@ -538,7 +542,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) } #endif // 2) call after_scf() of ESolver_FP - ESolver_FP::after_scf(ucell, istep); + ESolver_FP::after_scf(ucell, istep, conv_esolver); ModuleBase::timer::tick("ESolver_OF", "after_opt"); } diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index 8c69a477c5..f0aaf8ca7e 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -89,8 +89,8 @@ class ESolver_OF : public ESolver_FP void update_potential(UnitCell& ucell); void optimize(UnitCell& ucell); void update_rho(); - bool check_exit(); - void after_opt(const int istep, UnitCell& ucell); + bool check_exit(bool& conv_esolver); + void after_opt(const int istep, UnitCell& ucell, const bool conv_esolver); // ============================ tools =============================== // --------------------- initialize --------------------------------- @@ -110,7 +110,7 @@ class ESolver_OF : public ESolver_FP void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell); // --------------------- output the necessary information ----------- - void print_info(); + void print_info(const bool conv_esolver); // --------------------- interface to blas -------------------------- double inner_product(double* pa, double* pb, int length, double dV = 1) diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 5ac09580e4..c6f2725a1d 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -386,7 +386,7 @@ void ESolver_OF::test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& * @brief Print nessecary information to the screen, * and write the components of the total energy into running_log. */ -void ESolver_OF::print_info() +void ESolver_OF::print_info(const bool conv_esolver) { if (this->iter_ == 0) { @@ -431,8 +431,11 @@ void ESolver_OF::print_info() std::vector titles; std::vector energies_Ry; std::vector energies_eV; - if ((PARAM.inp.printe > 0 - && ((this->iter_ + 1) % PARAM.inp.printe == 0 || this->conv_esolver || this->iter_ == PARAM.inp.scf_nmax)) || PARAM.inp.init_chg == "file") + if ((PARAM.inp.printe > 0 && + ((this->iter_ + 1) % PARAM.inp.printe == 0 || + conv_esolver || + this->iter_ == PARAM.inp.scf_nmax)) || + PARAM.inp.init_chg == "file") { titles.push_back("E_Total"); energies_Ry.push_back(this->pelec->f_en.etot); diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index a4c20d37d8..3049ea68a7 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -125,20 +125,20 @@ void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) } template -void ESolver_SDFT_PW::iter_finish(UnitCell& ucell, const int istep, int& iter) +void ESolver_SDFT_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { // call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(ucell, istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); } template -void ESolver_SDFT_PW::after_scf(UnitCell& ucell, const int istep) +void ESolver_SDFT_PW::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { ModuleBase::TITLE("ESolver_SDFT_PW", "after_scf"); ModuleBase::timer::tick("ESolver_SDFT_PW", "after_scf"); // 1) call after_scf() of ESolver_KS_PW - ESolver_KS_PW::after_scf(ucell, istep); + ESolver_KS_PW::after_scf(ucell, istep, conv_esolver); ModuleBase::timer::tick("ESolver_SDFT_PW", "after_scf"); } diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 3287fc9ecb..4cb0fe0ddd 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -39,9 +39,9 @@ class ESolver_SDFT_PW : public ESolver_KS_PW virtual void others(UnitCell& ucell, const int istep) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; - virtual void after_scf(UnitCell& ucell, const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; virtual void after_all_runners(UnitCell& ucell) override; diff --git a/source/module_esolver/lcao_after_scf.cpp b/source/module_esolver/lcao_after_scf.cpp new file mode 100644 index 0000000000..d7dc57e1de --- /dev/null +++ b/source/module_esolver/lcao_after_scf.cpp @@ -0,0 +1,494 @@ +#include "esolver_ks_lcao.h" + +#include "module_base/formatter.h" +#include "module_base/global_variable.h" +#include "module_base/tool_title.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_hamilt_lcao/module_dftu/dftu.h" +#include "module_io/berryphase.h" +#include "module_io/cube_io.h" +#include "module_io/dos_nao.h" +#include "module_io/io_dmk.h" +#include "module_io/io_npz.h" +#include "module_io/nscf_band.h" +#include "module_io/output_dmk.h" +#include "module_io/output_log.h" +#include "module_io/output_mat_sparse.h" +#include "module_io/output_mulliken.h" +#include "module_io/output_sk.h" +#include "module_io/to_qo.h" +#include "module_io/to_wannier90_lcao.h" +#include "module_io/to_wannier90_lcao_in_pw.h" +#include "module_io/write_HS.h" +#include "module_io/write_dmr.h" +#include "module_io/write_elecstat_pot.h" +#include "module_io/write_istate_info.h" +#include "module_io/write_proj_band_lcao.h" +#include "module_io/write_wfc_nao.h" +#include "module_parameter/parameter.h" + + +//--------------temporary---------------------------- +#include "module_base/global_function.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_ux.h" +#include "module_elecstate/module_charge/symmetry_rho.h" +#include "module_elecstate/occupy.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need DeePKS_init +#include "module_hamilt_lcao/module_dftu/dftu.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/print_info.h" + + +//mohan add 20250302 +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h" + + +#include +#ifdef __EXX +#include "module_io/restart_exx_csr.h" +#include "module_ri/RPA_LRI.h" +#endif + +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#include "module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h" +#endif +//-----force& stress------------------- +#include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" + +//-----HSolver ElecState Hamilt-------- +#include "module_elecstate/elecstate_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_hsolver/hsolver_lcao.h" +// function used by deepks +// #include "module_elecstate/cal_dm.h" +//--------------------------------------------------- + +// test RDMFT +#include "module_rdmft/rdmft.h" + +#include + +namespace ModuleESolver +{ + +template +void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); + + //------------------------------------------------------------------ + // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 + //------------------------------------------------------------------ + if (PARAM.inp.out_elf[0] > 0) + { + assert(this->psi != nullptr); + this->pelec->cal_tau(*(this->psi)); + } + + //------------------------------------------------------------------ + //! 2) call after_scf() of ESolver_KS + //------------------------------------------------------------------ + ESolver_KS::after_scf(ucell, istep, conv_esolver); + + //------------------------------------------------------------------ + //! 3) write density matrix for sparse matrix + //------------------------------------------------------------------ + ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), + this->pv, + PARAM.inp.out_dm1, + false, + PARAM.inp.out_app_flag, + ucell.get_iat2iwt(), + &ucell.nat, + istep); + + //------------------------------------------------------------------ + //! 4) write density matrix + //------------------------------------------------------------------ + if (PARAM.inp.out_dm) + { + std::vector efermis(PARAM.inp.nspin == 2 ? 2 : 1); + for (int ispin = 0; ispin < efermis.size(); ispin++) + { + efermis[ispin] = this->pelec->eferm.get_efval(ispin); + } + const int precision = 3; + ModuleIO::write_dmk(dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(), + precision, + efermis, + &(ucell), + this->pv); + } + +#ifdef __EXX + //------------------------------------------------------------------ + //! 5) write Hexx matrix for NSCF (see `out_chg` in docs/advanced/input_files/input-main.md) + //------------------------------------------------------------------ + if (PARAM.inp.calculation != "nscf") + { + if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.out_chg[0] + && istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14 + { + const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); + if (GlobalC::exx_info.info_ri.real_number) + { + ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exd->get_Hexxs()); + } + else + { + ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exc->get_Hexxs()); + } + } + } +#endif + + //------------------------------------------------------------------ + // 6) write Hamiltonian and Overlap matrix + //------------------------------------------------------------------ + for (int ik = 0; ik < this->kv.get_nks(); ++ik) + { + if (PARAM.inp.out_mat_hs[0]) + { + this->p_hamilt->updateHk(ik); + } + bool bit = false; // LiuXh, 2017-03-21 + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) + { + hamilt::MatrixBlock h_mat; + hamilt::MatrixBlock s_mat; + + this->p_hamilt->matrix(h_mat, s_mat); + + if (PARAM.inp.out_mat_hs[0]) + { + ModuleIO::save_mat(istep, + h_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "H", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + ModuleIO::save_mat(istep, + s_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "S", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + } + } + } + + //------------------------------------------------------------------ + // 7) write electronic wavefunctions + //------------------------------------------------------------------ + if (elecstate::ElecStateLCAO::out_wfc_lcao && (istep % PARAM.inp.out_interval == 0)) + { + ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO::out_wfc_lcao, + this->psi[0], + this->pelec->ekb, + this->pelec->wg, + this->pelec->klist->kvec_c, + this->pv, + istep); + } + + //------------------------------------------------------------------ + //! 8) write DeePKS information + //------------------------------------------------------------------ +#ifdef __DEEPKS + if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) + { + hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); + std::shared_ptr ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); + LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); + + deepks_interface.out_deepks_labels(this->pelec->f_en.etot, + this->pelec->klist->get_nks(), + ucell.nat, + PARAM.globalv.nlocal, + this->pelec->ekb, + this->pelec->klist->kvec_d, + ucell, + orb_, + this->gd, + &(this->pv), + *(this->psi), + dynamic_cast*>(this->pelec)->get_DM(), + p_ham_deepks); + } +#endif + + + //------------------------------------------------------------------ + //! 9) Perform RDMFT calculations + // rdmft, added by jghan, 2024-10-17 + //------------------------------------------------------------------ + if (PARAM.inp.rdmft == true) + { + ModuleBase::matrix occ_num(this->pelec->wg); + for (int ik = 0; ik < occ_num.nr; ++ik) + { + for (int inb = 0; inb < occ_num.nc; ++inb) + { + occ_num(ik, inb) /= this->kv.wk[ik]; + } + } + this->rdmft_solver.update_elec(ucell, occ_num, *(this->psi)); + + //! initialize the gradients of Etotal with respect to occupation numbers and wfc, + //! and set all elements to 0. + //! dedocc = d E/d Occ_Num + ModuleBase::matrix dedocc(this->pelec->wg.nr, this->pelec->wg.nc, true); + + //! dedwfc = d E/d wfc + psi::Psi dedwfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis(), this->kv.ngk, true); + dedwfc.zero_out(); + + double etot_rdmft = this->rdmft_solver.run(dedocc, dedwfc); + } + + +#ifdef __EXX + //------------------------------------------------------------------ + // 10) Write RPA information. + //------------------------------------------------------------------ + if (PARAM.inp.rpa) + { + RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); + rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), + MPI_COMM_WORLD, + ucell, + this->kv, + orb_); + rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); + rpa_lri_double.out_for_RPA(ucell, this->pv, *(this->psi), this->pelec); + } +#endif + + + //------------------------------------------------------------------ + // 11) write HR in npz format. + //------------------------------------------------------------------ + if (PARAM.inp.out_hr_npz) + { + this->p_hamilt->updateHk(0); // first k point, up spin + hamilt::HamiltLCAO, double>* p_ham_lcao + = dynamic_cast, double>*>(this->p_hamilt); + std::string zipname = "output_HR0.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); + + if (PARAM.inp.nspin == 2) + { + this->p_hamilt->updateHk(this->kv.get_nks() / 2); // the other half of k points, down spin + hamilt::HamiltLCAO, double>* p_ham_lcao + = dynamic_cast, double>*>(this->p_hamilt); + zipname = "output_HR1.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); + } + } + + //------------------------------------------------------------------ + // 12) write density matrix in the 'npz' format. + //------------------------------------------------------------------ + if (PARAM.inp.out_dm_npz) + { + const elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + std::string zipname = "output_DM0.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(1))); + + if (PARAM.inp.nspin == 2) + { + zipname = "output_DM1.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(2))); + } + } + + //------------------------------------------------------------------ + //! 13) Print out information every 'out_interval' steps. + //------------------------------------------------------------------ + if (PARAM.inp.calculation != "md" || istep % PARAM.inp.out_interval == 0) + { + //! Print out sparse matrix + ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2, + PARAM.inp.out_mat_dh, + PARAM.inp.out_mat_t, + PARAM.inp.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->pv, + this->GK, + two_center_bundle_, + orb_, + ucell, + this->gd, + this->kv, + this->p_hamilt); + + //! Perform Mulliken charge analysis + if (PARAM.inp.out_mul) + { + ModuleIO::cal_mag(&(this->pv), + this->p_hamilt, + this->kv, + this->pelec, + this->two_center_bundle_, + this->orb_, + ucell, + this->gd, + istep, + true); + } + } + + //------------------------------------------------------------------ + //! 14) Print out atomic magnetization only when 'spin_constraint' is on. + //------------------------------------------------------------------ + if (PARAM.inp.sc_mag_switch) + { + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + sc.cal_mi_lcao(istep); + sc.print_Mi(GlobalV::ofs_running); + sc.print_Mag_Force(GlobalV::ofs_running); + } + + //------------------------------------------------------------------ + //! 15) Clean up RA. + //! this should be last function and put it in the end, mohan request 2024-11-28 + //------------------------------------------------------------------ + if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress) + { + RA.delete_grid(); + } + + //------------------------------------------------------------------ + //! 16) Print out quasi-orbitals. + //------------------------------------------------------------------ + if (PARAM.inp.qo_switch) + { + toQO tqo(PARAM.inp.qo_basis, PARAM.inp.qo_strategy, PARAM.inp.qo_thr, PARAM.inp.qo_screening_coeff); + tqo.initialize(PARAM.globalv.global_out_dir, + PARAM.inp.pseudo_dir, + PARAM.inp.orbital_dir, + &ucell, + this->kv.kvec_d, + GlobalV::ofs_running, + GlobalV::MY_RANK, + GlobalV::NPROC); + tqo.calculate(); + } + + //------------------------------------------------------------------ + //! 17) Print out kinetic matrix. + //------------------------------------------------------------------ + if (PARAM.inp.out_mat_tk[0]) + { + hamilt::HS_Matrix_K hsk(&pv, true); + hamilt::HContainer hR(&pv); + hamilt::Operator* ekinetic + = new hamilt::EkineticNew>(&hsk, + this->kv.kvec_d, + &hR, + &ucell, + orb_.cutoffs(), + &this->gd, + two_center_bundle_.kinetic_orb.get()); + + const int nspin_k = (PARAM.inp.nspin == 2 ? 2 : 1); + for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik) + { + ekinetic->init(ik); + ModuleIO::save_mat(0, + hsk.get_hk(), + PARAM.globalv.nlocal, + false, + PARAM.inp.out_mat_tk[1], + 1, + PARAM.inp.out_app_flag, + "T", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + } + + delete ekinetic; + } + + //------------------------------------------------------------------ + //! 18) wannier90 interface, added by jingan in 2018.11.7 + //------------------------------------------------------------------ + if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) + { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wave function to Wannier90"); + if (PARAM.inp.wannier_method == 1) + { + toWannier90_LCAO_IN_PW myWannier(PARAM.inp.out_wannier_mmn, + PARAM.inp.out_wannier_amn, + PARAM.inp.out_wannier_unk, + PARAM.inp.out_wannier_eig, + PARAM.inp.out_wannier_wvfn_formatted, + PARAM.inp.nnkpfile, + PARAM.inp.wannier_spin); + myWannier.set_tpiba_omega(ucell.tpiba, ucell.omega); + myWannier.calculate(ucell, + this->pelec->ekb, + this->pw_wfc, + this->pw_big, + this->sf, + this->kv, + this->psi, + &(this->pv)); + } + else if (PARAM.inp.wannier_method == 2) + { + toWannier90_LCAO myWannier(PARAM.inp.out_wannier_mmn, + PARAM.inp.out_wannier_amn, + PARAM.inp.out_wannier_unk, + PARAM.inp.out_wannier_eig, + PARAM.inp.out_wannier_wvfn_formatted, + PARAM.inp.nnkpfile, + PARAM.inp.wannier_spin, + orb_); + + myWannier.calculate(ucell, this->gd, this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); + } + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); + } + + //------------------------------------------------------------------ + //! 19) berry phase calculations, added by jingan + //------------------------------------------------------------------ + if (PARAM.inp.calculation == "nscf" && berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) + { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation"); + berryphase bp(&(this->pv)); + bp.lcao_init(ucell, this->gd, this->kv, this->GridT, orb_); + // additional step before calling + // macroscopic_polarization (why capitalize + // the function name?) + bp.Macroscopic_polarization(ucell, this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); + } + + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); +} + +template class ESolver_KS_LCAO; +template class ESolver_KS_LCAO, double>; +template class ESolver_KS_LCAO, std::complex>; +} // namespace ModuleESolver diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h index 660aa0ceab..1433024bf1 100644 --- a/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h @@ -25,16 +25,16 @@ class GintAtom // getter functions const Atom* get_atom() const { return atom_; }; - const int get_ia() const { return ia_; }; - const int get_iat() const { return iat_; }; + int get_ia() const { return ia_; }; + int get_iat() const { return iat_; }; const Vec3i& get_bgrid_idx() const { return biggrid_idx_; }; const Vec3i& get_unitcell_idx() const { return unitcell_idx_; }; const Vec3i& get_R() const { return unitcell_idx_; }; const Vec3d& get_tau_in_bgrid() const { return tau_in_biggrid_; }; const Numerical_Orbital* get_orb() const { return orb_; }; - const int get_nw() const { return atom_->nw; }; - const double get_rcut() const { return orb_->getRcut(); }; + int get_nw() const { return atom_->nw; }; + double get_rcut() const { return orb_->getRcut(); }; /** * @brief Get the wave function values of the atom at a meshgrid. @@ -115,4 +115,4 @@ class GintAtom }; -} // namespace ModuleGint \ No newline at end of file +} // namespace ModuleGint diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h index a39042e558..f24d1194b4 100644 --- a/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h +++ b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h @@ -17,14 +17,14 @@ class LocalCellInfo std::shared_ptr unitcell_info); // getter functions - const int get_startidx_bx() const { return startidx_bx_; }; - const int get_startidx_by() const { return startidx_by_; }; - const int get_startidx_bz() const { return startidx_bz_; }; - const int get_nbx() const { return nbx_; }; - const int get_nby() const { return nby_; }; - const int get_nbz() const { return nbz_; }; - const int get_bgrids_num() const { return nbxyz_; }; - const int get_mgrids_num() const { return nmxyz_; }; + int get_startidx_bx() const { return startidx_bx_; }; + int get_startidx_by() const { return startidx_by_; }; + int get_startidx_bz() const { return startidx_bz_; }; + int get_nbx() const { return nbx_; }; + int get_nby() const { return nby_; }; + int get_nbz() const { return nbz_; }; + int get_bgrids_num() const { return nbxyz_; }; + int get_mgrids_num() const { return nmxyz_; }; std::shared_ptr get_unitcell_info() const { return unitcell_info_; }; std::shared_ptr get_bgrid_info() const { return unitcell_info_->get_bgrid_info(); }; @@ -123,4 +123,4 @@ class LocalCellInfo }; -} // namespace ModuleGint \ No newline at end of file +} // namespace ModuleGint diff --git a/source/module_io/output_log.cpp b/source/module_io/output_log.cpp index 191dc413a2..7bf3a08d99 100644 --- a/source/module_io/output_log.cpp +++ b/source/module_io/output_log.cpp @@ -7,7 +7,7 @@ namespace ModuleIO { -void output_convergence_after_scf(bool& convergence, double& energy, std::ofstream& ofs_running) +void output_convergence_after_scf(const bool &convergence, double& energy, std::ofstream& ofs_running) { if (convergence) { @@ -43,7 +43,7 @@ void output_after_relax(bool conv_ion, bool conv_esolver, std::ofstream& ofs_run } } -void output_efermi(bool& convergence, double& efermi, std::ofstream& ofs_running) +void output_efermi(const bool &convergence, double& efermi, std::ofstream& ofs_running) { if (convergence && PARAM.inp.out_level != "m") { @@ -311,4 +311,4 @@ void write_head(std::ofstream& ofs_running, const int& istep, const int& iter, c << " ELEC=" << std::setw(4) << iter << "--------------------------------\n"; } -}// namespace ModuleIO \ No newline at end of file +}// namespace ModuleIO diff --git a/source/module_io/output_log.h b/source/module_io/output_log.h index 68c93ade9d..2bbc90cd94 100644 --- a/source/module_io/output_log.h +++ b/source/module_io/output_log.h @@ -14,7 +14,7 @@ namespace ModuleIO /// @param convergence if is convergence /// @param energy the total energy in Ry /// @param ofs_running the output stream -void output_convergence_after_scf(bool& convergence, double& energy, std::ofstream& ofs_running = GlobalV::ofs_running); +void output_convergence_after_scf(const bool&convergence, double& energy, std::ofstream& ofs_running = GlobalV::ofs_running); /// @brief output after relaxation /// @param conv_ion if is convergence for ions @@ -26,7 +26,7 @@ void output_after_relax(bool conv_ion, bool conv_esolver, std::ofstream& ofs_run /// @param convergence if is convergence /// @param efermi /// @param ofs_running the output stream -void output_efermi(bool& convergence, double& efermi, std::ofstream& ofs_running = GlobalV::ofs_running); +void output_efermi(const bool &convergence, double& efermi, std::ofstream& ofs_running = GlobalV::ofs_running); /// @brief calculate and output the vacuum level /// We first determine the vacuum direction, then get the vacuum position based on the minimum of charge density, @@ -75,4 +75,4 @@ void write_head(std::ofstream& ofs_running, const int& istep, const int& iter, c } // namespace ModuleIO -#endif \ No newline at end of file +#endif