diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 3b95d05b35..351c8c340a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -245,10 +245,8 @@ OBJS_ESOLVER=esolver.o\ esolver_of.o\ esolver_of_tool.o\ esolver_of_interface.o\ - pw_fun.o\ pw_init_after_vc.o\ pw_init_globalc.o\ - pw_nscf.o\ pw_others.o\ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ @@ -259,7 +257,6 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ set_matrix_grid.o\ lcao_before_scf.o\ lcao_gets.o\ - lcao_nscf.o\ lcao_others.o\ lcao_init_after_vc.o\ lcao_fun.o\ diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 451c4025b1..ca7242248e 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -62,7 +62,7 @@ void Driver::driver_run() { { Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM); } - else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax") + else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") { Relax_Driver rl_driver; rl_driver.relax_driver(p_esolver); @@ -70,7 +70,6 @@ void Driver::driver_run() { else { //! supported "other" functions: - //! nscf(PW,LCAO), //! get_pchg(LCAO), //! test_memory(PW,LCAO), //! test_neighbour(LCAO), diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 5a517b23e0..f5672f4d99 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -1,7 +1,5 @@ #include "elecstate_lcao.h" -#include - #include "cal_dm.h" #include "module_base/timer.h" #include "module_elecstate/module_dm/cal_dm_psi.h" @@ -11,6 +9,8 @@ #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_parameter/parameter.h" +#include + namespace elecstate { @@ -21,34 +21,31 @@ void ElecStateLCAO>::psiToRho(const psi::Psicalculate_weights(); - - // the calculations of dm, and dm -> rho are, technically, two separate - // functionalities, as we cannot rule out the possibility that we may have a - // dm from other sources, such as read from file. However, since we are not - // separating them now, I opt to add a flag to control how dm is obtained as - // of now - if (!PARAM.inp.dm_to_rho) - { - this->calEBand(); - - ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); - - // this part for calculating DMK in 2d-block format, not used for charge - // now - // psi::Psi> dm_k_2d(); - - if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "scalapack_gvx" || PARAM.inp.ks_solver == "lapack" - || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp" - || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15 - { - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), - this->wg, - psi, - *(this->DM)); - this->DM->cal_DMR(); - } - } + // // the calculations of dm, and dm -> rho are, technically, two separate + // // functionalities, as we cannot rule out the possibility that we may have a + // // dm from other sources, such as read from file. However, since we are not + // // separating them now, I opt to add a flag to control how dm is obtained as + // // of now + // if (!PARAM.inp.dm_to_rho) + // { + // ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); + + // // this part for calculating DMK in 2d-block format, not used for charge + // // now + // // psi::Psi> dm_k_2d(); + + // if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == + // "scalapack_gvx" || PARAM.inp.ks_solver == "lapack" + // || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp" + // || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15 + // { + // elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), + // this->wg, + // psi, + // *(this->DM)); + // this->DM->cal_DMR(); + // } + // } for (int is = 0; is < PARAM.inp.nspin; is++) { @@ -83,23 +80,6 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); - this->calculate_weights(); - this->calEBand(); - - if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "scalapack_gvx" || PARAM.inp.ks_solver == "lapack" - || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp" || PARAM.inp.ks_solver == "cg_in_lcao") - { - ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); - - // get DMK in 2d-block format - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), - this->wg, - psi, - *(this->DM)); - this->DM->cal_DMR(); - ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); - } - for (int is = 0; is < PARAM.inp.nspin; is++) { ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], diff --git a/source/module_elecstate/elecstate_lcao.h b/source/module_elecstate/elecstate_lcao.h index e0fcf6affe..4beeb017f0 100644 --- a/source/module_elecstate/elecstate_lcao.h +++ b/source/module_elecstate/elecstate_lcao.h @@ -74,6 +74,8 @@ class ElecStateLCAO : public ElecState void dmToRho(std::vector pexsi_DM, std::vector pexsi_EDM); #endif + DensityMatrix* DM = nullptr; + protected: // calculate electronic charge density on grid points or density matrix in real space // the consequence charge density rho saved into rho_out, preparing for charge mixing. @@ -85,7 +87,6 @@ class ElecStateLCAO : public ElecState Gint_Gamma* gint_gamma = nullptr; // mohan add 2024-04-01 Gint_k* gint_k = nullptr; // mohan add 2024-04-01 - DensityMatrix* DM = nullptr; }; template diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index d73163da77..0070f1a193 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -89,9 +89,6 @@ void ElecStatePW::psiToRho(const psi::Psi& psi) ModuleBase::timer::tick("ElecStatePW", "psiToRho"); this->init_rho_data(); - this->calculate_weights(); - - this->calEBand(); for(int is=0; is::psiToRho(const psi::Psi& psi) if (GlobalV::MY_STOGROUP == 0) { - this->calEBand(); - for (int is = 0; is < nspin; is++) { setmem_var_op()(this->ctx, this->rho[is], 0, this->charge->nrxx); diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 54f6d602d4..491bb95268 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -10,10 +10,8 @@ list(APPEND objects esolver_of.cpp esolver_of_interface.cpp esolver_of_tool.cpp - pw_fun.cpp pw_init_after_vc.cpp pw_init_globalc.cpp - pw_nscf.cpp pw_others.cpp ) if(ENABLE_LCAO) @@ -26,7 +24,6 @@ if(ENABLE_LCAO) dftu_cal_occup_m.cpp lcao_before_scf.cpp lcao_gets.cpp - lcao_nscf.cpp lcao_others.cpp lcao_init_after_vc.cpp lcao_fun.cpp diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 02aa386dd8..cc0d56fc75 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -52,9 +52,6 @@ namespace ModuleESolver //! Electorn charge density Charge chr; - //! Non-Self-Consistant Filed (NSCF) calculations - virtual void nscf(){}; - //! Structure factors that used with plane-wave basis set Structure_Factor sf; diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 7cf729958a..d4560d2a7c 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -402,7 +402,7 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); hsolver_error = 0.0; - if (iter == 1) + if (iter == 1 && PARAM.inp.calculation != "nscf") { hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, PARAM.inp.esolver_type, diag_ethr, PARAM.inp.nelec); @@ -472,13 +472,10 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); + // 4) SCF iterations this->conv_esolver = false; this->niter = this->maxniter; - - // 4) SCF iterations this->diag_ethr = PARAM.inp.pw_diag_thr; - - std::cout << " * * * * * *\n << Start SCF iteration." << std::endl; for (int iter = 1; iter <= this->maxniter; ++iter) { // 6) initialization of SCF iterations @@ -500,7 +497,6 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) break; } } // end scf iterations - std::cout << " >> Leave SCF iteration.\n * * * * * *" << std::endl; // 15) after scf ModuleBase::timer::tick(this->classname, "after_scf"); @@ -632,7 +628,7 @@ void ESolver_KS::iter_finish(const int istep, int& iter) // If drho < hsolver_error in the first iter or drho < scf_thr, we // do not change rho. - if (drho < hsolver_error || this->conv_esolver) + if (drho < hsolver_error || this->conv_esolver || PARAM.inp.calculation == "nscf") { if (drho < hsolver_error) { diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 88d127b619..f18264f7f7 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -49,9 +49,6 @@ class ESolver_KS : public ESolver_FP // calculate electron density from a specific Hamiltonian with ethr virtual void hamilt2density_single(const int istep, const int iter, const double ethr); - // calculate electron states from a specific Hamiltonian - virtual void hamilt2estates(const double ethr) {}; - // calculate electron density from a specific Hamiltonian void hamilt2density(const int istep, const int iter, const double ethr); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 0249f7e3ae..78c424bfd0 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -1,7 +1,10 @@ #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_io/berryphase.h" #include "module_io/cube_io.h" #include "module_io/dos_nao.h" #include "module_io/nscf_band.h" @@ -10,6 +13,8 @@ #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_eband_terms.hpp" #include "module_io/write_elecstat_pot.h" @@ -47,7 +52,7 @@ #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" +// #include "module_elecstate/cal_dm.h" //--------------------------------------------------- #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" @@ -586,6 +591,14 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // and the ncalculate the charge density on grid. this->pelec->skip_weights = true; + this->pelec->calculate_weights(); + if (!PARAM.inp.dm_to_rho) + { + auto _pelec = dynamic_cast*>(this->pelec); + _pelec->calEBand(); + elecstate::cal_dm_psi(_pelec->DM->get_paraV_pointer(), _pelec->wg, *this->psi, *(_pelec->DM)); + _pelec->DM->cal_DMR(); + } this->pelec->psiToRho(*this->psi); this->pelec->skip_weights = false; @@ -714,9 +727,9 @@ void ESolver_KS_LCAO::hamilt2density_single(int istep, int iter, double // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; - + bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false); + hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge); // 5) what's the exd used for? #ifdef __EXX @@ -1192,6 +1205,53 @@ void ESolver_KS_LCAO::after_scf(const int istep) delete ekinetic; } + + // add 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 + .calculate(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(this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); + } + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); + } + + // add 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(this->kv, + this->GridT, + orb_); // additional step before calling + // macroscopic_polarization (why capitalize + // the function name?) + bp.Macroscopic_polarization(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"); + } } //------------------------------------------------------------------------------ diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index b5f7c6093b..75eced339a 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -39,8 +39,6 @@ class ESolver_KS_LCAO : public ESolver_KS { void after_all_runners() override; - void nscf() override; - void get_S(); void cal_mag(const int istep, const bool print = false); diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 0a7a88402b..2221e291ef 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -165,8 +165,9 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite this->pelec->f_en.demet = 0.0; if (this->psi != nullptr) { + bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverLCAO> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, false); + hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, skip_charge); } } diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index fc40a76eac..0eb22b4ebf 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -126,6 +126,7 @@ namespace ModuleESolver hsolver::DiagoIterAssist::SCF_ITER = iter; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; + bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; // It is not a good choice to overload another solve function here, this will spoil the concept of // multiple inheritance and polymorphism. But for now, we just do it in this way. @@ -138,7 +139,7 @@ namespace ModuleESolver } hsolver::HSolverLIP hsolver_lip_obj(this->pw_wfc); - hsolver_lip_obj.solve(this->p_hamilt, this->kspw_psi[0], this->pelec, psig.lock().get()[0], false); + hsolver_lip_obj.solve(this->p_hamilt, this->kspw_psi[0], this->pelec, psig.lock().get()[0], skip_charge); // add exx #ifdef __EXX diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 51eb2c35a0..8af664799d 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -1,14 +1,5 @@ #include "esolver_ks_pw.h" -#include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/elecond.h" -#include "module_io/input_conv.h" -#include "module_io/nscf_band.h" -#include "module_io/output_log.h" -#include "module_io/write_dos_pw.h" -#include "module_io/write_istate_info.h" -#include "module_io/write_wfc_pw.h" - #include //--------------temporary---------------------------- @@ -22,10 +13,13 @@ //-----stress------------------ #include "module_hamilt_pw/hamilt_pwdft/stress_pw.h" //--------------------------------------------------- +#include "module_base/formatter.h" +#include "module_base/global_variable.h" #include "module_base/memory.h" #include "module_base/module_device/device.h" #include "module_elecstate/elecstate_pw.h" #include "module_hamilt_general/module_vdw/vdw.h" +#include "module_hamilt_pw/hamilt_pwdft/elecond.h" #include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h" #include "module_hsolver/diago_iter_assist.h" #include "module_hsolver/hsolver_pw.h" @@ -34,17 +28,22 @@ #include "module_io/berryphase.h" #include "module_io/cube_io.h" #include "module_io/get_pchg_pw.h" +#include "module_io/input_conv.h" +#include "module_io/nscf_band.h" #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" +#include "module_io/output_log.h" #include "module_io/to_wannier90_pw.h" #include "module_io/winput.h" +#include "module_io/write_dos_pw.h" #include "module_io/write_elecstat_pot.h" +#include "module_io/write_istate_info.h" +#include "module_io/write_wfc_pw.h" #include "module_io/write_wfc_r.h" #include "module_parameter/parameter.h" #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" #endif -#include "module_base/formatter.h" #include #include @@ -113,6 +112,21 @@ ESolver_KS_PW::~ESolver_KS_PW() delete this->p_wf_init; } +template +void ESolver_KS_PW::allocate_hamilt() +{ + this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv); +} +template +void ESolver_KS_PW::deallocate_hamilt() +{ + if (this->p_hamilt != nullptr) + { + delete reinterpret_cast*>(this->p_hamilt); + this->p_hamilt = nullptr; + } +} + template void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) { @@ -339,7 +353,11 @@ void ESolver_KS_PW::hamilt2density_single(const int istep, const int hsolver::DiagoIterAssist::SCF_ITER = iter; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; + if (PARAM.inp.calculation != "nscf") + { + hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; + } + bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, &this->wf, @@ -361,7 +379,7 @@ void ESolver_KS_PW::hamilt2density_single(const int istep, const int this->pelec->ekb.c, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, - false); + skip_charge); this->init_psi = true; @@ -526,6 +544,31 @@ void ESolver_KS_PW::after_scf(const int istep) PARAM.globalv.global_out_dir, PARAM.inp.if_separate_k); } + + //! 6) calculate Wannier functions + if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) + { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wannier functions calculation"); + toWannier90_PW wan(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); + + wan.calculate(this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi); + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation"); + } + + //! 7) calculate Berry phase polarization + 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 polarization"); + berryphase bp; + bp.Macroscopic_polarization(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 polarization"); + } } template diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index d27b1198c5..f9476bafb8 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -33,10 +33,6 @@ class ESolver_KS_PW : public ESolver_KS virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; - virtual void hamilt2estates(const double ethr) override; - - virtual void nscf() override; - void after_all_runners() override; protected: diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index a0e93c325c..8ee6098fe2 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -174,7 +174,7 @@ void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, doub this->pelec->f_en.demet = 0.0; // choose if psi should be diag in subspace // be careful that istep start from 0 and iter start from 1 - if (istep == 0 && iter == 1) + if (istep == 0 && iter == 1 || PARAM.inp.calculation == "nscf") { hsolver::DiagoIterAssist::need_subspace = false; } @@ -183,8 +183,8 @@ void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, doub hsolver::DiagoIterAssist::need_subspace = true; } + bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; // hsolver only exists in this function @@ -206,7 +206,15 @@ void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, doub hsolver::DiagoIterAssist::need_subspace, this->init_psi); - hsolver_pw_sdft_obj.solve(this->p_hamilt, this->kspw_psi[0], this->psi[0], this->pelec, this->pw_wfc, this->stowf, istep, iter, false); + hsolver_pw_sdft_obj.solve(this->p_hamilt, + this->kspw_psi[0], + this->psi[0], + this->pelec, + this->pw_wfc, + this->stowf, + istep, + iter, + skip_charge); this->init_psi = true; // set_diagethr need it @@ -353,39 +361,8 @@ void ESolver_SDFT_PW::others(const int istep) { ModuleBase::TITLE("ESolver_SDFT_PW", "others"); - if (PARAM.inp.calculation == "nscf") - { - this->nscf(); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_SDFT_PW::others", "CALCULATION type not supported"); - } - - return; -} - -template -void ESolver_SDFT_PW::nscf() -{ - ModuleBase::TITLE("ESolver_SDFT_PW", "nscf"); - ModuleBase::timer::tick("ESolver_SDFT_PW", "nscf"); - - const int istep = 0; - - const int iter = 1; - - const double diag_thr = std::max(std::min(1e-5, 0.1 * PARAM.inp.scf_thr / std::max(1.0, PARAM.inp.nelec)), 1e-12); - - std::cout << " DIGA_THR : " << diag_thr << std::endl; - - this->before_scf(istep); - - this->hamilt2density_single(istep, iter, diag_thr); - - this->pelec->cal_energies(2); + ModuleBase::WARNING_QUIT("ESolver_SDFT_PW::others", "CALCULATION type not supported"); - ModuleBase::timer::tick("ESolver_SDFT_PW", "nscf"); return; } diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 9b00f33d78..63337942fc 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -37,8 +37,6 @@ class ESolver_SDFT_PW : public ESolver_KS_PW virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; - virtual void nscf() override; - virtual void others(const int istep) override; virtual void iter_finish(const int istep, int& iter) override; diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 2062f5afda..72de1c87de 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -286,6 +286,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) this->read_mat_npz(zipname, *(dm->get_DMR_pointer(2))); } + this->pelec->calculate_weights(); this->pelec->psiToRho(*this->psi); int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1; diff --git a/source/module_esolver/lcao_nscf.cpp b/source/module_esolver/lcao_nscf.cpp deleted file mode 100644 index 4bb4624e6a..0000000000 --- a/source/module_esolver/lcao_nscf.cpp +++ /dev/null @@ -1,245 +0,0 @@ -#include "module_elecstate/module_charge/symmetry_rho.h" -#include "module_esolver/esolver_ks_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" -#include "module_hamilt_lcao/module_dftu/dftu.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -// -#include "module_base/timer.h" -#include "module_cell/module_neighbor/sltk_atom_arrange.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_io/berryphase.h" -#include "module_io/cube_io.h" -#include "module_io/get_pchg_lcao.h" -#include "module_io/get_wf_lcao.h" -#include "module_io/to_wannier90_lcao.h" -#include "module_io/to_wannier90_lcao_in_pw.h" -#include "module_io/write_HS_R.h" -#include "module_io/write_elecstat_pot.h" -#include "module_parameter/parameter.h" -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif -#include "module_base/formatter.h" -#include "module_elecstate/elecstate_lcao.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_hamilt_general/module_ewald/H_Ewald_pw.h" -#include "module_hamilt_general/module_vdw/vdw.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" -#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" -#include "module_io/read_wfc_nao.h" -#include "module_io/write_elecstat_pot.h" -#include "module_io/write_wfc_nao.h" -#ifdef __EXX -#include "module_io/restart_exx_csr.h" -#endif - -namespace ModuleESolver -{ - -template -void ESolver_KS_LCAO::nscf() { - ModuleBase::TITLE("ESolver_KS_LCAO", "nscf"); - - std::cout << " NON-SELF CONSISTENT CALCULATIONS" << std::endl; - - time_t time_start = std::time(nullptr); - - // mohan add 2021-02-09 - // in ions, istep starts from 1, - // then when the istep is a variable of scf or nscf, - // istep becomes istep-1, this should be fixed in future - int istep = 0; - hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, true); - - time_t time_finish = std::time(nullptr); - ModuleBase::GlobalFunc::OUT_TIME("cal_bands", time_start, time_finish); - - GlobalV::ofs_running << " end of band structure calculation " << std::endl; - GlobalV::ofs_running << " band eigenvalue in this processor (eV) :" << std::endl; - - const int nspin = PARAM.inp.nspin; - const int nbands = PARAM.inp.nbands; - - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (nspin == 2) - { - if (ik == 0) - { - GlobalV::ofs_running << " spin up :" << std::endl; - } - if (ik == (this->kv.get_nks() / 2)) - { - GlobalV::ofs_running << " spin down :" << std::endl; - } - } - - GlobalV::ofs_running << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() << "): " << this->kv.kvec_c[ik].x - << " " << this->kv.kvec_c[ik].y << " " << this->kv.kvec_c[ik].z << std::endl; - - for (int ib = 0; ib < nbands; ++ib) - { - GlobalV::ofs_running << " spin" << this->kv.isk[ik] + 1 << "final_state " << ib + 1 << " " - << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " - << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; - } - GlobalV::ofs_running << std::endl; - } - if (PARAM.inp.out_bandgap) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "band gap calculation"); - if (!PARAM.globalv.two_fermi) { - this->pelec->cal_bandgap(); - GlobalV::ofs_running << " E_bandgap " << this->pelec->bandgap * ModuleBase::Ry_to_eV << " eV" << std::endl; - } - else - { - this->pelec->cal_bandgap_updw(); - GlobalV::ofs_running << " E_bandgap_up " << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" - << std::endl; - GlobalV::ofs_running << " E_bandgap_dw " << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" - << std::endl; - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "band gap calculation"); - } - - // add by jingan in 2018.11.7 - if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) - { -#ifdef __LCAO - 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.calculate(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(this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); -#endif - } - - // add by jingan - if (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(this->kv, - this->GridT, - orb_); // additional step before calling - // macroscopic_polarization (why capitalize - // the function name?) - bp.Macroscopic_polarization(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"); - } - - // below is for DeePKS NSCF calculation -#ifdef __DEEPKS - if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "DeepKS output"); - const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - this->dpks_cal_projected_DM(dm); - GlobalC::ld.cal_descriptor(GlobalC::ucell.nat); // final descriptor - GlobalC::ld.cal_gedm(GlobalC::ucell.nat); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "DeepKS output"); - } -#endif - - this->create_Output_Mat_Sparse(0).write(); - - // mulliken charge analysis - if (PARAM.inp.out_mul) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Mulliken charge analysis"); - elecstate::ElecStateLCAO* pelec_lcao - = dynamic_cast*>(this->pelec); - this->pelec->calculate_weights(); - this->pelec->calEBand(); - elecstate::cal_dm_psi(&(this->pv), pelec_lcao->wg, *(this->psi), *(pelec_lcao->get_DM())); - this->cal_mag(istep, true); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Mulliken charge analysis"); - } - - /// write potential - if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3) - { - for (int is = 0; is < PARAM.inp.nspin; is++) - { - std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_POT.cube"; - - ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - 0, - fn, - 0.0, // efermi - &(GlobalC::ucell), - 3, // precision - 0); // out_fermi - } - } - else if (PARAM.inp.out_pot == 2) - { - std::string fn = PARAM.globalv.global_out_dir + "/ElecStaticPot.cube"; - ModuleIO::write_elecstat_pot( -#ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, -#endif - fn, - 0, // istep - this->pw_rhod, - this->pelec->charge, - &(GlobalC::ucell), - this->pelec->pot->get_fixed_v()); - } - - // write wfc - if (PARAM.inp.out_wfc_lcao) - { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing wave function"); - ModuleIO::write_wfc_nao(PARAM.inp.out_wfc_lcao, - *this->psi, - this->pelec->ekb, - this->pelec->wg, - this->pelec->klist->kvec_c, - this->pv, - istep); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing wave function"); - } -} - -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_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 5b07c91f16..418d4db8f9 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -85,11 +85,7 @@ void ESolver_KS_LCAO::others(const int istep) // pelec should be initialized before these calculations this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm); // self consistent calculations for electronic ground state - if (PARAM.inp.calculation == "nscf") - { - this->nscf(); - } - else if (cal_type == "get_pchg") + if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); IState_Charge ISC(this->psi, &(this->pv)); diff --git a/source/module_esolver/pw_fun.cpp b/source/module_esolver/pw_fun.cpp deleted file mode 100644 index 81ca1d00be..0000000000 --- a/source/module_esolver/pw_fun.cpp +++ /dev/null @@ -1,107 +0,0 @@ -#include "esolver_ks_pw.h" - -#include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/elecond.h" -#include "module_io/input_conv.h" -#include "module_io/nscf_band.h" -#include "module_io/output_log.h" -#include "module_io/write_dos_pw.h" -#include "module_io/write_istate_info.h" -#include "module_io/write_wfc_pw.h" - -#include - -//--------------temporary---------------------------- -#include "module_elecstate/module_charge/symmetry_rho.h" -#include "module_elecstate/occupy.h" -#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/memory.h" -#include "module_base/module_device/device.h" -#include "module_elecstate/elecstate_pw.h" -#include "module_hamilt_general/module_vdw/vdw.h" -#include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h" -#include "module_hsolver/diago_iter_assist.h" -#include "module_hsolver/hsolver_pw.h" -#include "module_hsolver/kernels/dngvd_op.h" -#include "module_hsolver/kernels/math_kernel_op.h" -#include "module_io/berryphase.h" -#include "module_io/numerical_basis.h" -#include "module_io/numerical_descriptor.h" -#include "module_io/to_wannier90_pw.h" -#include "module_io/winput.h" -#include "module_io/write_elecstat_pot.h" -#include "module_io/write_wfc_r.h" -#include "module_parameter/parameter.h" -#ifdef USE_PAW -#include "module_cell/module_paw/paw_cell.h" -#endif -#include -#include -#include "module_base/formatter.h" - -namespace ModuleESolver { - -template -void ESolver_KS_PW::allocate_hamilt() -{ - this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv); -} -template -void ESolver_KS_PW::deallocate_hamilt() -{ - if (this->p_hamilt != nullptr) - { - delete reinterpret_cast*>(this->p_hamilt); - this->p_hamilt = nullptr; - } -} - - -template -void ESolver_KS_PW::hamilt2estates(const double ethr) -{ - hsolver::DiagoIterAssist::need_subspace = false; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - - hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, - &this->wf, - - PARAM.inp.calculation, - PARAM.inp.basis_type, - PARAM.inp.ks_solver, - PARAM.inp.use_paw, - PARAM.globalv.use_uspp, - PARAM.inp.nspin, - - hsolver::DiagoIterAssist::SCF_ITER, - hsolver::DiagoIterAssist::PW_DIAG_NMAX, - hsolver::DiagoIterAssist::PW_DIAG_THR, - - hsolver::DiagoIterAssist::need_subspace, - this->init_psi); - - hsolver_pw_obj.solve(this->p_hamilt, - this->kspw_psi[0], - this->pelec, - this->pelec->ekb.c, - GlobalV::RANK_IN_POOL, - GlobalV::NPROC_IN_POOL, - true); - - this->init_psi = true; -} - -template class ESolver_KS_PW, base_device::DEVICE_CPU>; -template class ESolver_KS_PW, base_device::DEVICE_CPU>; -#if ((defined __CUDA) || (defined __ROCM)) -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -#endif -} // namespace ModuleESolver diff --git a/source/module_esolver/pw_nscf.cpp b/source/module_esolver/pw_nscf.cpp deleted file mode 100644 index 09217ff02f..0000000000 --- a/source/module_esolver/pw_nscf.cpp +++ /dev/null @@ -1,206 +0,0 @@ -#include "esolver_ks_pw.h" -#include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/elecond.h" -#include "module_io/cube_io.h" -#include "module_io/input_conv.h" -#include "module_io/nscf_band.h" -#include "module_io/output_log.h" -#include "module_io/write_dos_pw.h" -#include "module_io/write_istate_info.h" -#include "module_io/write_wfc_pw.h" - -#include - -//--------------temporary---------------------------- -#include "module_elecstate/module_charge/symmetry_rho.h" -#include "module_elecstate/occupy.h" -#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/memory.h" -#include "module_base/module_device/device.h" -#include "module_elecstate/elecstate_pw.h" -#include "module_hamilt_general/module_vdw/vdw.h" -#include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h" -#include "module_hsolver/diago_iter_assist.h" -#include "module_hsolver/hsolver_pw.h" -#include "module_hsolver/kernels/dngvd_op.h" -#include "module_hsolver/kernels/math_kernel_op.h" -#include "module_io/berryphase.h" -#include "module_io/numerical_basis.h" -#include "module_io/numerical_descriptor.h" -#include "module_io/to_wannier90_pw.h" -#include "module_io/winput.h" -#include "module_io/write_elecstat_pot.h" -#include "module_io/write_wfc_r.h" -#include "module_parameter/parameter.h" -#ifdef USE_PAW -#include "module_cell/module_paw/paw_cell.h" -#endif -#include -#include -#include "module_base/formatter.h" - -namespace ModuleESolver { - - -template -void ESolver_KS_PW::nscf() { - ModuleBase::TITLE("ESolver_KS_PW", "nscf"); - ModuleBase::timer::tick("ESolver_KS_PW", "nscf"); - - //! 1) before_scf - const int istep_tmp = 0; - this->before_scf(istep_tmp); - - //! 2) setup the parameters for diagonalization - double diag_ethr = PARAM.inp.pw_diag_thr; - if (diag_ethr - 1e-2 > -1e-5) { - diag_ethr - = std::max(1e-13, - 0.1 * std::min(1e-2, PARAM.inp.scf_thr / PARAM.inp.nelec)); - } - GlobalV::ofs_running << " PW_DIAG_THR = " << diag_ethr << std::endl; - - this->hamilt2estates(diag_ethr); - - //! 3) calculate weights/Fermi energies - this->pelec->calculate_weights(); - - GlobalV::ofs_running << "\n End of Band Structure Calculation \n" - << std::endl; - - //! 4) print out band energies and weights - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing band energies"); - const int nspin = PARAM.inp.nspin; - const int nbands = PARAM.inp.nbands; - for (int ik = 0; ik < this->kv.get_nks(); ik++) { - if (nspin == 2) { - if (ik == 0) { - GlobalV::ofs_running << " spin up :" << std::endl; - } - if (ik == (this->kv.get_nks() / 2)) { - GlobalV::ofs_running << " spin down :" << std::endl; - } - } - - GlobalV::ofs_running - << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() - << "): " << this->kv.kvec_c[ik].x << " " << this->kv.kvec_c[ik].y - << " " << this->kv.kvec_c[ik].z << std::endl; - - for (int ib = 0; ib < nbands; ib++) { - GlobalV::ofs_running - << " spin" << this->kv.isk[ik] + 1 << "_final_band " << ib + 1 - << " " << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " - << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; - } - GlobalV::ofs_running << std::endl; - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing band energies"); - - //! 5) print out band gaps - if (PARAM.inp.out_bandgap) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing band gaps"); - if (!PARAM.globalv.two_fermi) { - this->pelec->cal_bandgap(); - GlobalV::ofs_running << " E_bandgap " - << this->pelec->bandgap * ModuleBase::Ry_to_eV - << " eV" << std::endl; - } else { - this->pelec->cal_bandgap_updw(); - GlobalV::ofs_running - << " E_bandgap_up " - << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" - << std::endl; - GlobalV::ofs_running - << " E_bandgap_dw " - << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" - << std::endl; - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing band gaps"); - } - - //! 6) calculate Wannier functions - if (PARAM.inp.towannier90) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wannier functions calculation"); - toWannier90_PW wan(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); - - wan.calculate(this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->kv, - this->psi); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation"); - } - - //! 7) calculate Berry phase polarization - if (berryphase::berry_phase_flag - && ModuleSymmetry::Symmetry::symm_flag != 1) { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase polarization"); - berryphase bp; - bp.Macroscopic_polarization(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 polarization"); - } - - /// write potential - if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3) - { - for (int is = 0; is < PARAM.inp.nspin; is++) - { - std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_POT.cube"; - - ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - 0, - fn, - 0.0, // efermi - &(GlobalC::ucell), - 3, // precision - 0); // out_fermi - } - } - else if (PARAM.inp.out_pot == 2) - { - std::string fn = PARAM.globalv.global_out_dir + "/ElecStaticPot.cube"; - ModuleIO::write_elecstat_pot( -#ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, -#endif - fn, - 0, // istep - this->pw_rhod, - this->pelec->charge, - &(GlobalC::ucell), - this->pelec->pot->get_fixed_v()); - } - - ModuleBase::timer::tick("ESolver_KS_PW", "nscf"); - return; -} - -template class ESolver_KS_PW, base_device::DEVICE_CPU>; -template class ESolver_KS_PW, base_device::DEVICE_CPU>; -#if ((defined __CUDA) || (defined __ROCM)) -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -#endif -} // namespace ModuleESolver diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index 23439501c1..5db594c225 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -67,10 +67,7 @@ void ESolver_KS_PW::others(const int istep) { PARAM.inp.bessel_descriptor_rcut, PARAM.inp.bessel_descriptor_tolerence, this->kv.get_nks()); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, - "GENERATE DESCRIPTOR FOR DEEPKS"); - } else if (cal_type == "nscf") { - this->nscf(); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "GENERATE DESCRIPTOR FOR DEEPKS"); } else { ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); diff --git a/source/module_hsolver/hsolver.cpp b/source/module_hsolver/hsolver.cpp index 1d0edb0e56..5b9a6bd40e 100644 --- a/source/module_hsolver/hsolver.cpp +++ b/source/module_hsolver/hsolver.cpp @@ -21,7 +21,14 @@ double set_diagethr_ks(const std::string basis_type, if (basis_type == "pw" && esolver_type == "ksdft") { // It is too complex now and should be modified. - if (iter == 1) + if (calculation_in == "nscf") + { + if (res_diag_ethr - 1e-2 > -1e-5) + { + res_diag_ethr = std::max(1e-13, 0.1 * std::min(1e-2, PARAM.inp.scf_thr / PARAM.inp.nelec)); + } + } + else if (iter == 1) { if (std::abs(res_diag_ethr - 1.0e-2) < 1.0e-6) { @@ -94,7 +101,11 @@ double set_diagethr_sdft(const std::string basis_type, if (basis_type == "pw" && esolver_type == "sdft") { - if (iter == 1) + if (calculation_in == "nscf") + { + res_diag_ethr = std::max(std::min(1e-5, 0.1 * PARAM.inp.scf_thr / std::max(1.0, PARAM.inp.nelec)), 1e-12); + } + else if (iter == 1) { if (istep == 0) { diff --git a/source/module_hsolver/hsolver_lcao.cpp b/source/module_hsolver/hsolver_lcao.cpp index 741398faad..2f9a1b4313 100644 --- a/source/module_hsolver/hsolver_lcao.cpp +++ b/source/module_hsolver/hsolver_lcao.cpp @@ -22,12 +22,14 @@ #ifdef __PEXSI #include "diago_pexsi.h" -#include "module_elecstate/elecstate_lcao.h" #endif #include "module_base/global_variable.h" #include "module_base/memory.h" #include "module_base/timer.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/module_dm/density_matrix.h" #include "module_hsolver/parallel_k2d.h" #include "module_parameter/parameter.h" @@ -73,6 +75,15 @@ void HSolverLCAO::solve(hamilt::Hamilt* pHamilt, "This method and KPAR setting is not supported for lcao basis in ABACUS!"); } + pes->calculate_weights(); + if (!PARAM.inp.dm_to_rho) + { + auto _pes = dynamic_cast*>(pes); + _pes->calEBand(); + elecstate::cal_dm_psi(_pes->DM->get_paraV_pointer(), _pes->wg, psi, *(_pes->DM)); + _pes->DM->cal_DMR(); + } + if (!skip_charge) { // used in scf calculation diff --git a/source/module_hsolver/hsolver_lcaopw.cpp b/source/module_hsolver/hsolver_lcaopw.cpp index f5dc48fe67..be509b5629 100644 --- a/source/module_hsolver/hsolver_lcaopw.cpp +++ b/source/module_hsolver/hsolver_lcaopw.cpp @@ -272,6 +272,8 @@ void HSolverLIP::solve(hamilt::Hamilt* pHamilt, // ESolver_KS_PW::p_hamilt eigenvalues.data(), pes->ekb.nr * pes->ekb.nc); + reinterpret_cast*>(pes)->calculate_weights(); + reinterpret_cast*>(pes)->calEBand(); if (skip_charge) { ModuleBase::timer::tick("HSolverLIP", "solve"); diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 5b384ea3c5..e9d64e8455 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -327,6 +327,8 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // pes->ekb.nr * pes->ekb.nc this->wfc_basis->nks * psi.get_nbands()); + reinterpret_cast*>(pes)->calculate_weights(); + reinterpret_cast*>(pes)->calEBand(); if (skip_charge) { ModuleBase::timer::tick("HSolverPW", "solve"); diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp index f802162f0f..8698324f79 100644 --- a/source/module_hsolver/hsolver_pw_sdft.cpp +++ b/source/module_hsolver/hsolver_pw_sdft.cpp @@ -80,6 +80,12 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt* pHamilt, // prepare sqrt{f(\hat{H})}|\chi> to calculate density, force and stress stoiter.calHsqrtchi(stowf); + + elecstate::ElecStatePW* pes_pw = static_cast*>(pes); + if (GlobalV::MY_STOGROUP == 0) + { + pes_pw->calEBand(); + } if (skip_charge) { ModuleBase::timer::tick("HSolverPW_SDFT", "solve"); @@ -87,7 +93,6 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt* pHamilt, } //(5) calculate new charge density // calculate KS rho. - elecstate::ElecStatePW* pes_pw = static_cast*>(pes); pes_pw->init_rho_data(); if (nbands > 0) { diff --git a/source/module_io/read_input_item_elec_stru.cpp b/source/module_io/read_input_item_elec_stru.cpp index e039b7f043..6b6b9944e7 100644 --- a/source/module_io/read_input_item_elec_stru.cpp +++ b/source/module_io/read_input_item_elec_stru.cpp @@ -492,6 +492,12 @@ void ReadInput::item_elec_stru() Input_Item item("scf_nmax"); item.annotation = "number of electron iterations"; read_sync_int(input.scf_nmax); + item.reset_value = [](const Input_Item& item, Parameter& para) { + if (para.input.calculation == "nscf") + { + para.input.scf_nmax = 1; + } + }; this->add_item(item); } { diff --git a/source/module_io/read_input_item_system.cpp b/source/module_io/read_input_item_system.cpp index af0456aa17..fb698ae877 100644 --- a/source/module_io/read_input_item_system.cpp +++ b/source/module_io/read_input_item_system.cpp @@ -549,6 +549,10 @@ void ReadInput::item_system() { ModuleBase::WARNING_QUIT("ReadInput", "dm_to_rho is not available for parallel calculations"); } + if (para.input.dm_to_rho && para.inp.gamma_only) + { + ModuleBase::WARNING_QUIT("ReadInput", "dm_to_rho is not available for gamma_only calculations"); + } if (para.input.dm_to_rho) { #ifndef __USECNPY diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp index faec37cbc3..ed3688ded5 100644 --- a/source/module_io/test_serial/read_input_item_test.cpp +++ b/source/module_io/test_serial/read_input_item_test.cpp @@ -934,6 +934,14 @@ TEST_F(InputTest, Item_test2) output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + param.input.dm_to_rho = true; + param.input.gamma_only = true; + GlobalV::NPROC = 1; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(0), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + #ifndef __USECNPY param.input.dm_to_rho = true; GlobalV::NPROC = 1; diff --git a/tests/integrate/101_PW_15_f_pseudopots/log b/tests/integrate/101_PW_15_f_pseudopots/log deleted file mode 100644 index e69de29bb2..0000000000