diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 4f135f8724..2c726c05e0 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -241,14 +241,25 @@ OBJS_ESOLVER=esolver.o\ esolver_of_tool.o\ esolver_of_interface.o\ print_funcs.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\ - esolver_ks_lcao_elec.o\ esolver_ks_lcao_tddft.o\ dpks_cal_e_delta_band.o\ dftu_cal_occup_m.o\ io_npz.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\ + cal_edm_tddft.o\ OBJS_GINT=gint.o\ gint_gamma_env.o\ diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 932d501851..8ddf0949c3 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -11,16 +11,27 @@ list(APPEND objects esolver_of_interface.cpp esolver_of_tool.cpp print_funcs.cpp + pw_fun.cpp + pw_init_after_vc.cpp + pw_init_globalc.cpp + pw_nscf.cpp + pw_others.cpp ) if(ENABLE_LCAO) list(APPEND objects esolver_ks_lcao.cpp - esolver_ks_lcao_elec.cpp esolver_ks_lcao_tddft.cpp dpks_cal_e_delta_band.cpp io_npz.cpp set_matrix_grid.cpp 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 + cal_edm_tddft.cpp ) endif() diff --git a/source/module_esolver/cal_edm_tddft.cpp b/source/module_esolver/cal_edm_tddft.cpp new file mode 100644 index 0000000000..e3fe4d0b2c --- /dev/null +++ b/source/module_esolver/cal_edm_tddft.cpp @@ -0,0 +1,291 @@ +#include "esolver_ks_lcao_tddft.h" + +#include "module_io/cal_r_overlap_R.h" +#include "module_io/dipole_io.h" +#include "module_io/rho_io.h" +#include "module_io/td_current_io.h" +#include "module_io/write_HS.h" +#include "module_io/write_HS_R.h" +#include "module_io/write_wfc_nao.h" + +//--------------temporary---------------------------- +#include "module_base/blas_connector.h" +#include "module_base/global_function.h" +#include "module_base/scalapack_connector.h" +#include "module_base/lapack_connector.h" +#include "module_elecstate/module_charge/symmetry_rho.h" +#include "module_elecstate/occupy.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" +#include "module_hamilt_lcao/module_tddft/td_velocity.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/print_info.h" + +//-----HSolver ElecState Hamilt-------- +#include "module_elecstate/elecstate_lcao.h" +#include "module_elecstate/elecstate_lcao_tddft.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_hsolver/hsolver_lcao.h" +#include "module_parameter/parameter.h" +#include "module_psi/psi.h" + +//-----force& stress------------------- +#include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" + +//--------------------------------------------------- + +namespace ModuleESolver +{ + +// use the original formula (Hamiltonian matrix) to calculate energy density +// matrix +void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() +{ + // mohan add 2024-03-27 + const int nlocal = GlobalV::NLOCAL; + assert(nlocal >= 0); + + dynamic_cast>*>(this->pelec) + ->get_DM() + ->EDMK.resize(kv.get_nks()); + for (int ik = 0; ik < kv.get_nks(); ++ik) { + std::complex* tmp_dmk + = dynamic_cast>*>(this->pelec)->get_DM()->get_DMK_pointer(ik); + + ModuleBase::ComplexMatrix& tmp_edmk + = dynamic_cast>*>(this->pelec)->get_DM()->EDMK[ik]; + + const Parallel_Orbitals* tmp_pv + = dynamic_cast>*>(this->pelec)->get_DM()->get_paraV_pointer(); + +#ifdef __MPI + + // mohan add 2024-03-27 + //! be careful, the type of nloc is 'long' + //! whether the long type is safe, needs more discussion + const long nloc = this->pv.nloc; + const int ncol = this->pv.ncol; + const int nrow = this->pv.nrow; + + tmp_edmk.create(ncol, nrow); + complex* Htmp = new complex[nloc]; + complex* Sinv = new complex[nloc]; + complex* tmp1 = new complex[nloc]; + complex* tmp2 = new complex[nloc]; + complex* tmp3 = new complex[nloc]; + complex* tmp4 = new complex[nloc]; + + ModuleBase::GlobalFunc::ZEROS(Htmp, nloc); + ModuleBase::GlobalFunc::ZEROS(Sinv, nloc); + ModuleBase::GlobalFunc::ZEROS(tmp1, nloc); + ModuleBase::GlobalFunc::ZEROS(tmp2, nloc); + ModuleBase::GlobalFunc::ZEROS(tmp3, nloc); + ModuleBase::GlobalFunc::ZEROS(tmp4, nloc); + + const int inc = 1; + + hamilt::MatrixBlock> h_mat; + hamilt::MatrixBlock> s_mat; + + p_hamilt->matrix(h_mat, s_mat); + zcopy_(&nloc, h_mat.p, &inc, Htmp, &inc); + zcopy_(&nloc, s_mat.p, &inc, Sinv, &inc); + + vector ipiv(nloc, 0); + int info = 0; + const int one_int = 1; + + pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, this->pv.desc, ipiv.data(), &info); + + int lwork = -1; + int liwork = -1; + + // if lwork == -1, then the size of work is (at least) of length 1. + std::vector> work(1, 0); + + // if liwork = -1, then the size of iwork is (at least) of length 1. + std::vector iwork(1, 0); + + pzgetri_(&nlocal, + Sinv, + &one_int, + &one_int, + this->pv.desc, + ipiv.data(), + work.data(), + &lwork, + iwork.data(), + &liwork, + &info); + + lwork = work[0].real(); + work.resize(lwork, 0); + liwork = iwork[0]; + iwork.resize(liwork, 0); + + pzgetri_(&nlocal, + Sinv, + &one_int, + &one_int, + this->pv.desc, + ipiv.data(), + work.data(), + &lwork, + iwork.data(), + &liwork, + &info); + + const char N_char = 'N'; + const char T_char = 'T'; + const complex one_float = {1.0, 0.0}; + const complex zero_float = {0.0, 0.0}; + const complex half_float = {0.5, 0.0}; + + pzgemm_(&N_char, + &N_char, + &nlocal, + &nlocal, + &nlocal, + &one_float, + tmp_dmk, + &one_int, + &one_int, + this->pv.desc, + Htmp, + &one_int, + &one_int, + this->pv.desc, + &zero_float, + tmp1, + &one_int, + &one_int, + this->pv.desc); + + pzgemm_(&N_char, + &N_char, + &nlocal, + &nlocal, + &nlocal, + &one_float, + tmp1, + &one_int, + &one_int, + this->pv.desc, + Sinv, + &one_int, + &one_int, + this->pv.desc, + &zero_float, + tmp2, + &one_int, + &one_int, + this->pv.desc); + + pzgemm_(&N_char, + &N_char, + &nlocal, + &nlocal, + &nlocal, + &one_float, + Sinv, + &one_int, + &one_int, + this->pv.desc, + Htmp, + &one_int, + &one_int, + this->pv.desc, + &zero_float, + tmp3, + &one_int, + &one_int, + this->pv.desc); + + pzgemm_(&N_char, + &N_char, + &nlocal, + &nlocal, + &nlocal, + &one_float, + tmp3, + &one_int, + &one_int, + this->pv.desc, + tmp_dmk, + &one_int, + &one_int, + this->pv.desc, + &zero_float, + tmp4, + &one_int, + &one_int, + this->pv.desc); + + pzgeadd_(&N_char, + &nlocal, + &nlocal, + &half_float, + tmp2, + &one_int, + &one_int, + this->pv.desc, + &half_float, + tmp4, + &one_int, + &one_int, + this->pv.desc); + + zcopy_(&nloc, tmp4, &inc, tmp_edmk.c, &inc); + + delete[] Htmp; + delete[] Sinv; + delete[] tmp1; + delete[] tmp2; + delete[] tmp3; + delete[] tmp4; +#else + // for serial version + tmp_edmk.create(this->pv.ncol, this->pv.nrow); + ModuleBase::ComplexMatrix Sinv(nlocal, nlocal); + ModuleBase::ComplexMatrix Htmp(nlocal, nlocal); + + hamilt::MatrixBlock> h_mat; + hamilt::MatrixBlock> s_mat; + + p_hamilt->matrix(h_mat, s_mat); + // cout<<"hmat "<* work = new std::complex[lwork]; + ModuleBase::GlobalFunc::ZEROS(work, lwork); + + int IPIV[nlocal]; + + LapackConnector::zgetrf(nlocal, nlocal, Sinv, nlocal, IPIV, &INFO); + LapackConnector::zgetri(nlocal, Sinv, nlocal, IPIV, work, lwork, &INFO); + // I just use ModuleBase::ComplexMatrix temporarily, and will change it + // to complex* + ModuleBase::ComplexMatrix tmp_dmk_base(nlocal, nlocal); + for (int i = 0; i < nlocal; i++) + { + for (int j = 0; j < nlocal; j++) + { + tmp_dmk_base(i, j) = tmp_dmk[i * nlocal + j]; + } + } + tmp_edmk = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); + delete[] work; +#endif + } + return; +} +} // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index adbbaa1b4e..fda4aefd18 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -98,7 +98,7 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 1) calculate overlap matrix S or initialize //! 2) init ElecState //! 3) init LCAO basis -//! 4) redundant ParaV and LM pointers +//! 4) redundant pv and LM pointers //! 5) initialize Density Matrix //! 6) initialize Hamilt in LCAO //! 7) initialize exx @@ -161,7 +161,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 3) init LCAO basis // reading the localized orbitals/projectors // construct the interpolation tables. - LCAO_domain::init_basis_lcao(this->ParaV, + LCAO_domain::init_basis_lcao(this->pv, inp.onsite_radius, inp.lcao_ecut, inp.lcao_dk, @@ -175,7 +175,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // DensityMatrix is allocated here, DMK is also initialized here // DMR is not initialized here, it will be constructed in each before_scf dynamic_cast*>(this->pelec) - ->init_DM(&this->kv, &(this->ParaV), GlobalV::NSPIN); + ->init_DM(&this->kv, &(this->pv), GlobalV::NSPIN); // this function should be removed outside of the function if (GlobalV::CALCULATION == "get_S") @@ -187,7 +187,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 6) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - LCAO_domain::divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, ParaV, this->kv.get_nks()); + LCAO_domain::divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, pv, this->kv.get_nks()); #ifdef __EXX // 7) initialize exx @@ -208,7 +208,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 8) initialize DFT+U if (GlobalV::dft_plus_u) { - GlobalC::dftu.init(ucell, &this->ParaV, this->kv.get_nks()); + GlobalC::dftu.init(ucell, &this->pv, this->kv.get_nks()); } // 9) initialize ppcell @@ -217,7 +217,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 10) initialize the HSolver if (this->phsol == nullptr) { - this->phsol = new hsolver::HSolverLCAO(&(this->ParaV)); + this->phsol = new hsolver::HSolverLCAO(&(this->pv)); this->phsol->method = GlobalV::KS_SOLVER; } @@ -279,54 +279,6 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell return; } -//------------------------------------------------------------------------------ -//! the 4th function of ESolver_KS_LCAO: init_after_vc -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ -template -void ESolver_KS_LCAO::init_after_vc(const Input_para& inp, UnitCell& ucell) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); - - ESolver_KS::init_after_vc(inp, ucell); - if (inp.mdp.md_prec_level == 2) - { - delete this->pelec; - this->pelec = new elecstate::ElecStateLCAO( - &(this->chr), - &(this->kv), - this->kv.get_nks(), - &(this->GG), // mohan add 2024-04-01 - &(this->GK), // mohan add 2024-04-01 - this->pw_rho, - this->pw_big); - - dynamic_cast*>(this->pelec) - ->init_DM(&this->kv, &this->ParaV, GlobalV::NSPIN); - - GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); - - this->pelec->charge->allocate(GlobalV::NSPIN); - this->pelec->omega = GlobalC::ucell.omega; - - // Initialize the potential. - if (this->pelec->pot == nullptr) - { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &GlobalC::ucell, - &(GlobalC::ppcell.vloc), - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); - } - } - - ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); - return; -} - //------------------------------------------------------------------------------ //! the 5th function of ESolver_KS_LCAO: cal_energy //! mohan add 2024-05-11 @@ -355,7 +307,7 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) GlobalV::CAL_STRESS, GlobalV::TEST_FORCE, GlobalV::TEST_STRESS, - this->ParaV, + this->pv, this->pelec, this->psi, this->GG, // mohan add 2024-04-01 @@ -477,13 +429,13 @@ void ESolver_KS_LCAO::after_all_runners() if (PARAM.inp.out_proj_band) // Projeced band structure added by jiyy-2022-4-20 { - ModuleIO::write_proj_band_lcao(this->psi, this->ParaV, this->pelec, this->kv, GlobalC::ucell, this->p_hamilt); + ModuleIO::write_proj_band_lcao(this->psi, this->pv, this->pelec, this->kv, GlobalC::ucell, this->p_hamilt); } if (PARAM.inp.out_dos) { ModuleIO::out_dos_nao(this->psi, - this->ParaV, + this->pv, this->pelec->ekb, this->pelec->wg, PARAM.inp.dos_edelta_ev, @@ -502,7 +454,7 @@ void ESolver_KS_LCAO::after_all_runners() ModuleIO::write_Vxc(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, - &this->ParaV, + &this->pv, *this->psi, GlobalC::ucell, this->sf, @@ -527,7 +479,7 @@ void ESolver_KS_LCAO::after_all_runners() ModuleIO::write_eband_terms(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, - &this->ParaV, + &this->pv, *this->psi, GlobalC::ucell, this->sf, @@ -783,11 +735,11 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) #ifdef __EXX if (GlobalC::exx_info.info_ri.real_number) { - this->exd->exx_hamilt2density(*this->pelec, this->ParaV, iter); + this->exd->exx_hamilt2density(*this->pelec, this->pv, iter); } else { - this->exc->exx_hamilt2density(*this->pelec, this->ParaV, iter); + this->exc->exx_hamilt2density(*this->pelec, this->pv, iter); } #endif @@ -894,7 +846,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), - this->ParaV, + this->pv, GlobalV::DRANK); ModuleIO::save_mat(istep, s_mat.p, @@ -905,13 +857,13 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), - this->ParaV, + this->pv, GlobalV::DRANK); } #ifdef __DEEPKS if(GlobalV::deepks_out_labels && GlobalV::deepks_v_delta) { - DeePKS_domain::save_h_mat(h_mat.p, this->ParaV.nloc); + DeePKS_domain::save_h_mat(h_mat.p, this->pv.nloc); } #endif } @@ -927,7 +879,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) this->pelec->ekb, this->pelec->wg, this->pelec->klist->kvec_c, - this->ParaV, + this->pv, istep); } @@ -997,7 +949,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) { ////////// for Add_Hexx_Type::k /* - hamilt::HS_Matrix_K Hexxk_save(&this->ParaV, 1); + hamilt::HS_Matrix_K Hexxk_save(&this->pv, 1); for (int ik = 0; ik < this->kv.get_nks(); ++ik) { Hexxk_save.set_zero_hk(); @@ -1009,7 +961,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) GlobalC::restart.save_disk("Hexx", ik, - this->ParaV.get_local_size(), + this->pv.get_local_size(), Hexxk_save.get_hk()); }*/ ////////// for Add_Hexx_Type:R @@ -1102,7 +1054,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) // 2) write density matrix for sparse matrix ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), - this->ParaV, + this->pv, PARAM.inp.out_dm1, false, GlobalV::out_app_flag, @@ -1134,7 +1086,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) precision, efermis, &(GlobalC::ucell), - this->ParaV); + this->pv); } #ifdef __EXX @@ -1180,7 +1132,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) GlobalC::ucell, GlobalC::ORB, GlobalC::GridD, - &(this->ParaV), + &(this->pv), *(this->psi), dynamic_cast*>(this->pelec)->get_DM(), GlobalV::deepks_v_delta); @@ -1200,7 +1152,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) MPI_COMM_WORLD, this->kv); rpa_lri_double.init(MPI_COMM_WORLD, this->kv); - rpa_lri_double.out_for_RPA(this->ParaV, *(this->psi), this->pelec); + rpa_lri_double.out_for_RPA(this->pv, *(this->psi), this->pelec); } #endif @@ -1322,89 +1274,6 @@ bool ESolver_KS_LCAO::do_after_converge(int& iter) return true; } -//------------------------------------------------------------------------------ -//! the 16th function of ESolver_KS_LCAO: create_Output_DM -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ - -//------------------------------------------------------------------------------ -//! the 17th function of ESolver_KS_LCAO: create_Output_DM1 -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ - -//------------------------------------------------------------------------------ -//! the 18th function of ESolver_KS_LCAO: create_Output_Mat_Sparse -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ -template -ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) -{ - return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, - hsolver::HSolverLCAO::out_mat_dh, - hsolver::HSolverLCAO::out_mat_t, - PARAM.inp.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->ParaV, - this->GK, // mohan add 2024-04-01 - two_center_bundle_, - GlobalC::GridD, // mohan add 2024-04-06 - this->kv, - this->p_hamilt); -} - -//------------------------------------------------------------------------------ -//! the 19th function of ESolver_KS_LCAO: md_skip_out -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ -template -bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, int interval) -{ - if (calculation == "md") - { - if (istep % interval != 0) - { - return true; - } - } - return false; -} - -template -void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) -{ - auto cell_index = CellIndex(GlobalC::ucell.get_atomLabels(), - GlobalC::ucell.get_atomCounts(), - GlobalC::ucell.get_lnchiCounts(), - GlobalV::NSPIN); - auto out_sk = ModuleIO::Output_Sk(this->p_hamilt, - &(this->ParaV), - GlobalV::NSPIN, - this->kv.get_nks()); - auto out_dmk = ModuleIO::Output_DMK(dynamic_cast*>(this->pelec)->get_DM(), - &(this->ParaV), - GlobalV::NSPIN, - this->kv.get_nks()); - auto mulp = ModuleIO::Output_Mulliken(&(out_sk), - &(out_dmk), - &(this->ParaV), - &cell_index, - this->kv.isk, - GlobalV::NSPIN); - auto atom_chg = mulp.get_atom_chg(); - /// used in updating mag info in STRU file - GlobalC::ucell.atom_mulliken = mulp.get_atom_mulliken(atom_chg); - if (print && GlobalV::MY_RANK == 0) - { - /// write the Orbital file - cell_index.write_orb_info(GlobalV::global_out_dir); - /// write mulliken.txt - mulp.write(istep, GlobalV::global_out_dir); - /// write atomic mag info in running log file - mulp.print_atom_mag(atom_chg, GlobalV::ofs_running); - } -} - //------------------------------------------------------------------------------ //! the 20th,21th,22th functions of ESolver_KS_LCAO //! mohan add 2024-05-11 diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 6c87f2d4b3..7e0791bb1e 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -68,7 +68,7 @@ class ESolver_KS_LCAO : public ESolver_KS { Record_adj RA; // 2d block-cyclic distribution info - Parallel_Orbitals ParaV; + Parallel_Orbitals pv; // used for k-dependent grid integration. Gint_k GK; diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp deleted file mode 100644 index b699efb457..0000000000 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ /dev/null @@ -1,751 +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/get_pchg.h" -#include "module_io/get_wf.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_parameter/parameter.h" -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif -#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/rho_io.h" -#include "module_io/write_pot.h" -#include "module_io/write_wfc_nao.h" -#include "module_io/read_wfc_nao.h" -#include "module_base/formatter.h" -#ifdef __EXX -#include "module_io/restart_exx_csr.h" -#endif - -namespace ModuleESolver -{ - -template -void ESolver_KS_LCAO::beforesolver(const int istep) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); - - // 1. prepare HS matrices, prepare grid integral - this->set_matrix_grid(this->RA); - - // 2. density matrix extrapolation - - // set the augmented orbitals index. - // after ParaO and GridT, - // this information is used to calculate - // the force. - - // init psi - if (this->psi == nullptr) - { - int nsk = 0; - int ncol = 0; - if (GlobalV::GAMMA_ONLY_LOCAL) - { - nsk = GlobalV::NSPIN; - ncol = this->ParaV.ncol_bands; - if (GlobalV::KS_SOLVER == "genelpa" - || GlobalV::KS_SOLVER == "lapack" - || GlobalV::KS_SOLVER == "pexsi" - || GlobalV::KS_SOLVER == "cusolver" - || GlobalV::KS_SOLVER == "cusolvermp") { - ncol = this->ParaV.ncol; - } - } - else - { - nsk = this->kv.get_nks(); -#ifdef __MPI - ncol = this->ParaV.ncol_bands; -#else - ncol = GlobalV::NBANDS; -#endif - } - this->psi = new psi::Psi(nsk, ncol, this->ParaV.nrow, nullptr); - } - - // init wfc from file - if(istep == 0 && PARAM.inp.init_wfc == "file") - { - if (! ModuleIO::read_wfc_nao(GlobalV::global_readin_dir, this->ParaV, *(this->psi), this->pelec)) - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", - "read wfc nao failed"); - } - } - - // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, *this->pw_rho, *this->pw_big); - - // init Hamiltonian - if (this->p_hamilt != nullptr) - { - delete this->p_hamilt; - this->p_hamilt = nullptr; - } - if (this->p_hamilt == nullptr) - { - elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); - this->p_hamilt = new hamilt::HamiltLCAO( - GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, - GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), - &this->ParaV, - this->pelec->pot, - this->kv, - two_center_bundle_, - DM -#ifdef __EXX - , GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step - , GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr - , GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs -#endif - ); - } - -#ifdef __DEEPKS - // for each ionic step, the overlap must be rebuilt - // since it depends on ionic positions - if (GlobalV::deepks_setorb) - { - const Parallel_Orbitals* pv = &this->ParaV; - // build and save at beginning - GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, - GlobalC::ucell, - GlobalC::ORB, - GlobalC::GridD, - *(two_center_bundle_.overlap_orb_alpha)); - - if (PARAM.inp.deepks_out_unittest) - { - GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD); - } - } -#endif - if (PARAM.inp.sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); - sc.init_sc(GlobalV::sc_thr, - PARAM.inp.nsc, - PARAM.inp.nsc_min, - PARAM.inp.alpha_trial, - PARAM.inp.sccut, - PARAM.inp.sc_mag_switch, - GlobalC::ucell, - PARAM.inp.sc_file, - GlobalV::NPOL, - &(this->ParaV), - GlobalV::NSPIN, - this->kv, - GlobalV::KS_SOLVER, - this->phsol, - this->p_hamilt, - this->psi, - this->pelec); - } - //========================================================= - // cal_ux should be called before init_scf because - // the direction of ux is used in noncoline_rho - //========================================================= - if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) - { - GlobalC::ucell.cal_ux(); - } - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); -} - -template -void ESolver_KS_LCAO::before_scf(const int istep) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); - - if (GlobalC::ucell.cell_parameter_updated) - { - this->init_after_vc(PARAM.inp, GlobalC::ucell); - } - if (GlobalC::ucell.ionic_position_updated) - { - this->CE.update_all_dis(GlobalC::ucell); - this->CE.extrapolate_charge( -#ifdef __MPI - &(GlobalC::Pgrid), -#endif - GlobalC::ucell, - this->pelec->charge, - &(this->sf)); - } - - //---------------------------------------------------------- - // about vdw, jiyy add vdwd3 and linpz add vdwd2 - //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); - if (vdw_solver != nullptr) - { - this->pelec->f_en.evdw = vdw_solver->get_energy(); - } - - this->beforesolver(istep); - - // Peize Lin add 2016-12-03 -#ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf - if (GlobalC::exx_info.info_ri.real_number) - { - this->exd->exx_beforescf(this->kv, *this->p_chgmix); - } - else - { - this->exc->exx_beforescf(this->kv, *this->p_chgmix); - } -#endif // __EXX - - this->pelec->init_scf(istep, this->sf.strucFac); - if (PARAM.inp.out_chg == 2) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - std::stringstream ss; - ss << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; - ModuleIO::write_rho( -#ifdef __MPI - this->pw_big->bz, // bz first, then nbz - this->pw_big->nbz, - this->pw_rho->nplane, - this->pw_rho->startz_current, -#endif - this->pelec->charge->rho[is], - is, - GlobalV::NSPIN, - 0, - ss.str(), - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pelec->eferm.ef, - &(GlobalC::ucell), - 11); - } - } - - ModuleIO::write_pot(GlobalV::out_pot, - GlobalV::NSPIN, - GlobalV::global_out_dir, -#ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, - this->pw_rho->nplane, - this->pw_rho->startz_current, -#endif - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pelec->pot->get_effective_v()); - - // initalize DMR - // DMR should be same size with Hamiltonian(R) - dynamic_cast*>(this->pelec) - ->get_DM() - ->init_DMR(*(dynamic_cast*>(this->p_hamilt)->getHR())); - // two cases are considered: - // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK - // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros - if(istep > 0) - { - dynamic_cast*>(this->pelec) - ->get_DM() - ->cal_DMR(); - } - - if (PARAM.inp.dm_to_rho) - { - std::string zipname = "output_DM0.npz"; - elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - this->read_mat_npz(zipname, *(dm->get_DMR_pointer(1))); - if (GlobalV::NSPIN == 2) - { - zipname = "output_DM1.npz"; - this->read_mat_npz(zipname, *(dm->get_DMR_pointer(2))); - } - - this->pelec->psiToRho(*this->psi); - - this->create_Output_Rho(0, istep).write(); - if (GlobalV::NSPIN == 2) - { - this->create_Output_Rho(1, istep).write(); - } - - return; - } - - // the electron charge density should be symmetrized, - // here is the initialization - Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); - } - - // 1. calculate ewald energy. - // mohan update 2021-02-25 - if (!PARAM.inp.test_skip_ewald) - { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho, this->sf.strucFac); - } - - this->p_hamilt->non_first_scf = istep; - - ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); - return; -} - -template -void ESolver_KS_LCAO::others(const int istep) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "others"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); - - const std::string cal_type = GlobalV::CALCULATION; - - if (cal_type == "get_S") { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing the overlap matrix"); - this->get_S(); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing the overlap matrix"); - - ModuleBase::QUIT(); - - // return; // use 'return' will cause segmentation fault. by mohan - // 2024-06-09 - } else if (cal_type == "test_memory") { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing memory"); - Cal_Test::test_memory(this->pw_rho, - this->pw_wfc, - this->p_chgmix->get_mixing_mode(), - this->p_chgmix->get_mixing_ndim()); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing memory"); - return; - } - else if (cal_type == "test_neighbour") - { - // test_search_neighbor(); - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing neighbour"); - if (GlobalV::SEARCH_RADIUS < 0) { - std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS - << std::endl; - std::cout << " please make sure search_radius > 0" << std::endl; - } - - atom_arrange::search(PARAM.inp.search_pbc, - GlobalV::ofs_running, - GlobalC::GridD, - GlobalC::ucell, - GlobalV::SEARCH_RADIUS, - GlobalV::test_atom_input, - true); - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing neighbour"); - return; - } - - this->beforesolver(istep); - // pelec should be initialized before these calculations - this->pelec->init_scf(istep, this->sf.strucFac); - // self consistent calculations for electronic ground state - if (GlobalV::CALCULATION == "nscf") - { - this->nscf(); - } else if (cal_type == "get_pchg") { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); - IState_Charge ISC(this->psi, &(this->ParaV)); - if (GlobalV::GAMMA_ONLY_LOCAL) { - ISC.begin(this->GG, - this->pelec->charge->rho, - this->pelec->wg, - this->pelec->eferm.get_all_ef(), - this->pw_rho->nrxx, - this->pw_rho->nplane, - this->pw_rho->startz_current, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pw_big->bz, - this->pw_big->nbz, - GlobalV::GAMMA_ONLY_LOCAL, - PARAM.inp.nbands_istate, - PARAM.inp.bands_to_print, - GlobalV::NBANDS, - GlobalV::nelec, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir, - GlobalV::MY_RANK, - GlobalV::ofs_warning, - &GlobalC::ucell, - &GlobalC::GridD, - this->kv); - } else { - ISC.begin(this->GK, - this->pelec->charge->rho, - this->pelec->charge->rhog, - this->pelec->wg, - this->pelec->eferm.get_all_ef(), - this->pw_rho, - this->pw_rho->nrxx, - this->pw_rho->nplane, - this->pw_rho->startz_current, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pw_big->bz, - this->pw_big->nbz, - GlobalV::GAMMA_ONLY_LOCAL, - PARAM.inp.nbands_istate, - PARAM.inp.bands_to_print, - GlobalV::NBANDS, - GlobalV::nelec, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir, - GlobalV::MY_RANK, - GlobalV::ofs_warning, - &GlobalC::ucell, - &GlobalC::GridD, - this->kv, - PARAM.inp.if_separate_k, - &GlobalC::Pgrid, - this->pelec->charge->ngmc); - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); - } else if (cal_type == "get_wf") { - std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); - IState_Envelope IEP(this->pelec); - if (GlobalV::GAMMA_ONLY_LOCAL) - { - IEP.begin(this->psi, - this->pw_rho, - this->pw_wfc, - this->pw_big, - this->ParaV, - this->GG, - PARAM.inp.out_wfc_pw, - this->wf.out_wfc_r, - this->kv, - GlobalV::nelec, - PARAM.inp.nbands_istate, - PARAM.inp.bands_to_print, - GlobalV::NBANDS, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir); - } - else - { - IEP.begin(this->psi, - this->pw_rho, - this->pw_wfc, - this->pw_big, - this->ParaV, - this->GK, - PARAM.inp.out_wfc_pw, - this->wf.out_wfc_r, - this->kv, - GlobalV::nelec, - PARAM.inp.nbands_istate, - PARAM.inp.bands_to_print, - GlobalV::NBANDS, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir); - } - std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting wave function"); - } else { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", - "CALCULATION type not supported"); - } - - ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); - return; -} - -template <> -void ESolver_KS_LCAO::get_S(void) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", "not implemented for"); -} - -template <> -void ESolver_KS_LCAO, double>::get_S(void) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); - // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); - - atom_arrange::search(PARAM.inp.search_pbc, - GlobalV::ofs_running, - GlobalC::GridD, - GlobalC::ucell, - GlobalV::SEARCH_RADIUS, - GlobalV::test_atom_input); - - this->RA.for_2d(this->ParaV, GlobalV::GAMMA_ONLY_LOCAL); - - if (this->p_hamilt == nullptr) { - this->p_hamilt = new hamilt::HamiltLCAO, double>( - &this->ParaV, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast, double>*>( - this->p_hamilt->ops) - ->contributeHR(); - } - - // mohan add 2024-06-09 - const std::string fn = GlobalV::global_out_dir + "SR.csr"; - - std::cout << " The file is saved in " << fn << std::endl; - - ModuleIO::output_SR(ParaV, GlobalC::GridD, this->p_hamilt, fn); - - return; -} - -template <> -void ESolver_KS_LCAO, std::complex>::get_S(void) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); - // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); - - atom_arrange::search(PARAM.inp.search_pbc, - GlobalV::ofs_running, - GlobalC::GridD, - GlobalC::ucell, - GlobalV::SEARCH_RADIUS, - GlobalV::test_atom_input); - - this->RA.for_2d(this->ParaV, GlobalV::GAMMA_ONLY_LOCAL); - if (this->p_hamilt == nullptr) { - this->p_hamilt = new hamilt::HamiltLCAO, - std::complex>( - &this->ParaV, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast< - hamilt::OperatorLCAO, std::complex>*>( - this->p_hamilt->ops) - ->contributeHR(); - } - - // mohan add 2024-06-09 - const std::string fn = GlobalV::global_out_dir + "SR.csr"; - - std::cout << " The file is saved in " << fn << std::endl; - - ModuleIO::output_SR(ParaV, GlobalC::GridD, this->p_hamilt, fn); - - return; -} - -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; - if (this->phsol != nullptr) - { - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER, true); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); - } - - 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 = GlobalV::NSPIN; - const int nbands = GlobalV::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 (!GlobalV::TWO_EFERMI) { - 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 (GlobalV::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->ParaV)); - } - 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); - - myWannier.calculate(this->pelec->ekb, this->kv, *(this->psi), &(this->ParaV)); - } - 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->ParaV)); - bp.lcao_init(this->kv, - this->GridT); // 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 (GlobalV::deepks_out_labels || GlobalV::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->ParaV), 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 - this->create_Output_Potential(0).write(); - - // 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->ParaV, - 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/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index ce73a4c46f..48f48829c1 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -84,7 +84,8 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& } // 4) read the local orbitals and construct the interpolation tables. - LCAO_domain::init_basis_lcao(this->ParaV, + // initialize the pv + LCAO_domain::init_basis_lcao(this->pv, inp.onsite_radius, inp.lcao_ecut, inp.lcao_dk, @@ -94,16 +95,16 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& two_center_bundle_); // 5) allocate H and S matrices according to computational resources - LCAO_domain::divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, ParaV, kv.get_nks()); + LCAO_domain::divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, this->pv, kv.get_nks()); // 6) initialize Density Matrix dynamic_cast>*>(this->pelec) - ->init_DM(&kv, &this->ParaV, GlobalV::NSPIN); + ->init_DM(&kv, &this->pv, GlobalV::NSPIN); // 7) initialize Hsolver if (this->phsol == nullptr) { - this->phsol = new hsolver::HSolverLCAO>(&this->ParaV); + this->phsol = new hsolver::HSolverLCAO>(&this->pv); this->phsol->method = GlobalV::KS_SOLVER; } @@ -136,7 +137,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons GlobalV::NBANDS, GlobalV::NLOCAL, this->p_hamilt, - this->ParaV, + this->pv, this->psi, this->psi_laststep, this->Hk_laststep, @@ -155,7 +156,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons GlobalV::NBANDS, GlobalV::NLOCAL, this->p_hamilt, - this->ParaV, + this->pv, this->psi, this->psi_laststep, this->Hk_laststep, @@ -267,7 +268,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), - this->ParaV, + this->pv, GlobalV::DRANK); ModuleIO::save_mat(istep, @@ -279,7 +280,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), - this->ParaV, + this->pv, GlobalV::DRANK); } } @@ -294,7 +295,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) this->pelec->ekb, this->pelec->wg, this->pelec->klist->kvec_c, - this->ParaV, + this->pv, istep); } @@ -313,9 +314,9 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) this->pelec->cal_converged(); } - const int nloc = this->ParaV.nloc; - const int ncol_nbands = this->ParaV.ncol_bands; - const int nrow = this->ParaV.nrow; + const int nloc = this->pv.nloc; + const int ncol_nbands = this->pv.ncol_bands; + const int nrow = this->pv.nrow; const int nbands = GlobalV::NBANDS; const int nlocal = GlobalV::NLOCAL; @@ -420,7 +421,9 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) } if (TD_Velocity::out_current == true) { - elecstate::DensityMatrix, double>* tmp_DM = dynamic_cast>*>(this->pelec)->get_DM(); + elecstate::DensityMatrix, double>* tmp_DM = + dynamic_cast>*>(this->pelec)->get_DM(); + ModuleIO::write_current(istep, this->psi, pelec, @@ -432,255 +435,4 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) ESolver_KS_LCAO, double>::after_scf(istep); } -// use the original formula (Hamiltonian matrix) to calculate energy density -// matrix -void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() -{ - // mohan add 2024-03-27 - const int nlocal = GlobalV::NLOCAL; - assert(nlocal >= 0); - - dynamic_cast>*>(this->pelec) - ->get_DM() - ->EDMK.resize(kv.get_nks()); - for (int ik = 0; ik < kv.get_nks(); ++ik) { - std::complex* tmp_dmk - = dynamic_cast>*>(this->pelec)->get_DM()->get_DMK_pointer(ik); - - ModuleBase::ComplexMatrix& tmp_edmk - = dynamic_cast>*>(this->pelec)->get_DM()->EDMK[ik]; - - const Parallel_Orbitals* tmp_pv - = dynamic_cast>*>(this->pelec)->get_DM()->get_paraV_pointer(); - -#ifdef __MPI - - // mohan add 2024-03-27 - //! be careful, the type of nloc is 'long' - //! whether the long type is safe, needs more discussion - const long nloc = this->ParaV.nloc; - const int ncol = this->ParaV.ncol; - const int nrow = this->ParaV.nrow; - - tmp_edmk.create(ncol, nrow); - complex* Htmp = new complex[nloc]; - complex* Sinv = new complex[nloc]; - complex* tmp1 = new complex[nloc]; - complex* tmp2 = new complex[nloc]; - complex* tmp3 = new complex[nloc]; - complex* tmp4 = new complex[nloc]; - - ModuleBase::GlobalFunc::ZEROS(Htmp, nloc); - ModuleBase::GlobalFunc::ZEROS(Sinv, nloc); - ModuleBase::GlobalFunc::ZEROS(tmp1, nloc); - ModuleBase::GlobalFunc::ZEROS(tmp2, nloc); - ModuleBase::GlobalFunc::ZEROS(tmp3, nloc); - ModuleBase::GlobalFunc::ZEROS(tmp4, nloc); - - const int inc = 1; - - hamilt::MatrixBlock> h_mat; - hamilt::MatrixBlock> s_mat; - - p_hamilt->matrix(h_mat, s_mat); - zcopy_(&nloc, h_mat.p, &inc, Htmp, &inc); - zcopy_(&nloc, s_mat.p, &inc, Sinv, &inc); - - vector ipiv(nloc, 0); - int info = 0; - const int one_int = 1; - - pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, this->ParaV.desc, ipiv.data(), &info); - - int lwork = -1; - int liwork = -1; - - // if lwork == -1, then the size of work is (at least) of length 1. - std::vector> work(1, 0); - - // if liwork = -1, then the size of iwork is (at least) of length 1. - std::vector iwork(1, 0); - - pzgetri_(&nlocal, - Sinv, - &one_int, - &one_int, - this->ParaV.desc, - ipiv.data(), - work.data(), - &lwork, - iwork.data(), - &liwork, - &info); - - lwork = work[0].real(); - work.resize(lwork, 0); - liwork = iwork[0]; - iwork.resize(liwork, 0); - - pzgetri_(&nlocal, - Sinv, - &one_int, - &one_int, - this->ParaV.desc, - ipiv.data(), - work.data(), - &lwork, - iwork.data(), - &liwork, - &info); - - const char N_char = 'N'; - const char T_char = 'T'; - const complex one_float = {1.0, 0.0}; - const complex zero_float = {0.0, 0.0}; - const complex half_float = {0.5, 0.0}; - - pzgemm_(&N_char, - &N_char, - &nlocal, - &nlocal, - &nlocal, - &one_float, - tmp_dmk, - &one_int, - &one_int, - this->ParaV.desc, - Htmp, - &one_int, - &one_int, - this->ParaV.desc, - &zero_float, - tmp1, - &one_int, - &one_int, - this->ParaV.desc); - - pzgemm_(&N_char, - &N_char, - &nlocal, - &nlocal, - &nlocal, - &one_float, - tmp1, - &one_int, - &one_int, - this->ParaV.desc, - Sinv, - &one_int, - &one_int, - this->ParaV.desc, - &zero_float, - tmp2, - &one_int, - &one_int, - this->ParaV.desc); - - pzgemm_(&N_char, - &N_char, - &nlocal, - &nlocal, - &nlocal, - &one_float, - Sinv, - &one_int, - &one_int, - this->ParaV.desc, - Htmp, - &one_int, - &one_int, - this->ParaV.desc, - &zero_float, - tmp3, - &one_int, - &one_int, - this->ParaV.desc); - - pzgemm_(&N_char, - &N_char, - &nlocal, - &nlocal, - &nlocal, - &one_float, - tmp3, - &one_int, - &one_int, - this->ParaV.desc, - tmp_dmk, - &one_int, - &one_int, - this->ParaV.desc, - &zero_float, - tmp4, - &one_int, - &one_int, - this->ParaV.desc); - - pzgeadd_(&N_char, - &nlocal, - &nlocal, - &half_float, - tmp2, - &one_int, - &one_int, - this->ParaV.desc, - &half_float, - tmp4, - &one_int, - &one_int, - this->ParaV.desc); - - zcopy_(&nloc, tmp4, &inc, tmp_edmk.c, &inc); - - delete[] Htmp; - delete[] Sinv; - delete[] tmp1; - delete[] tmp2; - delete[] tmp3; - delete[] tmp4; -#else - // for serial version - tmp_edmk.create(this->ParaV.ncol, this->ParaV.nrow); - ModuleBase::ComplexMatrix Sinv(nlocal, nlocal); - ModuleBase::ComplexMatrix Htmp(nlocal, nlocal); - - hamilt::MatrixBlock> h_mat; - hamilt::MatrixBlock> s_mat; - - p_hamilt->matrix(h_mat, s_mat); - // cout<<"hmat "<* work = new std::complex[lwork]; - ModuleBase::GlobalFunc::ZEROS(work, lwork); - - int IPIV[nlocal]; - - LapackConnector::zgetrf(nlocal, nlocal, Sinv, nlocal, IPIV, &INFO); - LapackConnector::zgetri(nlocal, Sinv, nlocal, IPIV, work, lwork, &INFO); - // I just use ModuleBase::ComplexMatrix temporarily, and will change it - // to complex* - ModuleBase::ComplexMatrix tmp_dmk_base(nlocal, nlocal); - for (int i = 0; i < nlocal; i++) - { - for (int j = 0; j < nlocal; j++) - { - tmp_dmk_base(i, j) = tmp_dmk[i * nlocal + j]; - } - } - tmp_edmk = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); - delete[] work; -#endif - } - return; -} } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 707515baf7..6380c6c229 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -99,60 +99,7 @@ ESolver_KS_PW::~ESolver_KS_PW() } template -void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell) -{ - // GlobalC is a historically left-over namespace, it is used to store global - // classes, including: pseudopot_cell_vnl: pseudopotential in cell, V - // non-local UnitCell: cell information with atomic properties Grid_Driver: - // Parallel_Grid: - // Parallel_Kpoints: - // Restart: - // Exx_Info: - // Exx_Lip: - - // GlobalC would be refactored out in the future. If there is better idea - // about how to organize information stored in classes above, please feel - // free to discuss with issue or pull request. - - if (this->psi != nullptr) - { - delete this->psi; - } - - //! init pseudopotential - ppcell.init(ucell.ntype, &this->sf, this->pw_wfc); - - //! initalize local pseudopotential - ppcell.init_vloc(ppcell.vloc, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - - //! Initalize non-local pseudopotential - ppcell.init_vnl(ucell, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); - - //! Allocate psi - this->p_wf_init->allocate_psi(this->psi, - this->kv.get_nkstot(), - this->kv.get_nks(), - this->kv.ngk.data(), - this->pw_wfc->npwk_max, - &(this->sf)); - - this->kspw_psi = GlobalV::device_flag == "gpu" || GlobalV::precision_flag == "single" - ? new psi::Psi(this->psi[0]) - : reinterpret_cast*>(this->psi); - - if (GlobalV::precision_flag == "single") - { - ModuleBase::Memory::record("Psi_single", sizeof(T) * this->psi[0].size()); - } - - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); -} - -template -void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) -{ +void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) { // 1) call before_all_runners() of ESolver_KS ESolver_KS::before_all_runners(inp, ucell); @@ -218,163 +165,7 @@ void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCel this->pelec->fixed_weights(PARAM.inp.ocp_kb, GlobalV::NBANDS, GlobalV::nelec); } } -template -void ESolver_KS_PW::allocate_hsolver() -{ - this->phsol = new hsolver::HSolverPW(this->pw_wfc, &this->wf); -} -template -void ESolver_KS_PW::deallocate_hsolver() -{ - delete reinterpret_cast*>(this->phsol); - this->phsol = nullptr; -} -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() -{ - delete reinterpret_cast*>(this->p_hamilt); - this->p_hamilt = nullptr; -} -template -void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& ucell) -{ - ModuleBase::TITLE("ESolver_KS_PW", "init_after_vc"); - ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); - - ESolver_KS::init_after_vc(inp, ucell); - - if (inp.mdp.md_prec_level == 2) - { - this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz); - - this->pw_wfc->initparameters(false, inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); - -#ifdef __MPI - if (PARAM.inp.pw_seed > 0) - { - MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - } - // qianrui add 2021-8-13 to make different kpar parameters can get the - // same results -#endif - - this->pw_wfc->setuptransform(); - - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - this->kv.ngk[ik] = this->pw_wfc->npwk[ik]; - } - - this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); - delete this->phsol; - this->allocate_hsolver(); - - delete this->pelec; - this->pelec = new elecstate::ElecStatePW(this->pw_wfc, - &(this->chr), - (K_Vectors*)(&(this->kv)), - &ucell, - &GlobalC::ppcell, - this->pw_rhod, - this->pw_rho, - this->pw_big); - - this->pelec->charge->allocate(GlobalV::NSPIN); - - //! setup cell volume - this->pelec->omega = ucell.omega; - - delete this->pelec->pot; - - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &ucell, - &GlobalC::ppcell.vloc, - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); - - // temporary - this->Init_GlobalC(inp, ucell, GlobalC::ppcell); - } - else - { - GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); - - this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); - - this->pw_wfc->initparameters(false, PARAM.inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); - - this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); - - this->p_wf_init->make_table(this->kv.get_nks(), &this->sf); - } - -#ifdef USE_PAW - if (GlobalV::use_paw) - { - GlobalC::paw_cell.set_libpaw_ecut(PARAM.inp.ecutwfc / 2.0, - PARAM.inp.ecutwfc / 2.0); // in Hartree - GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx, - this->pw_wfc->ny, - this->pw_wfc->nz, - this->pw_wfc->nx, - this->pw_wfc->ny, - this->pw_wfc->nz, - this->pw_wfc->startz, - this->pw_wfc->numz); - -#ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) - { - GlobalC::paw_cell.prepare_paw(); - } -#else - GlobalC::paw_cell.prepare_paw(); -#endif - GlobalC::paw_cell.set_sij(); - - std::vector> rhoijp; - std::vector> rhoijselect; - std::vector nrhoijsel; -#ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) - { - GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - - for (int iat = 0; iat < ucell.nat; iat++) - { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); - } - } -#else - GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - - for (int iat = 0; iat < ucell.nat; iat++) - { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); - } -#endif - } -#endif - - ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); -} template void ESolver_KS_PW::before_scf(const int istep) @@ -503,42 +294,6 @@ void ESolver_KS_PW::before_scf(const int istep) } } -template -void ESolver_KS_PW::others(const int istep) -{ - ModuleBase::TITLE("ESolver_KS_PW", "others"); - - const std::string cal_type = GlobalV::CALCULATION; - - if (cal_type == "test_memory") - { - Cal_Test::test_memory(this->pw_rho, - this->pw_wfc, - this->p_chgmix->get_mixing_mode(), - this->p_chgmix->get_mixing_ndim()); - } - else if (cal_type == "gen_bessel") - { - Numerical_Descriptor nc; - nc.output_descriptor(this->psi[0], - PARAM.inp.bessel_descriptor_lmax, - 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(); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); - } - - return; -} - template void ESolver_KS_PW::iter_init(const int istep, const int iter) { @@ -1090,145 +845,6 @@ void ESolver_KS_PW::after_all_runners() } } -template -void ESolver_KS_PW::hamilt2estates(const double ethr) -{ - if (this->phsol != nullptr) - { - hsolver::DiagoIterAssist::need_subspace = false; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - this->phsol->solve(this->p_hamilt, - this->kspw_psi[0], - this->pelec, - PARAM.inp.ks_solver, - PARAM.inp.calculation, - PARAM.inp.basis_type, - PARAM.inp.use_paw, - GlobalV::use_uspp, - GlobalV::RANK_IN_POOL, - GlobalV::NPROC_IN_POOL, - - hsolver::DiagoIterAssist::SCF_ITER, - hsolver::DiagoIterAssist::need_subspace, - hsolver::DiagoIterAssist::PW_DIAG_NMAX, - hsolver::DiagoIterAssist::PW_DIAG_THR, - - true); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); - } -} - -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 = GlobalV::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 / GlobalV::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 = GlobalV::NSPIN; - const int nbands = GlobalV::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 (!GlobalV::TWO_EFERMI) - { - 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 - this->create_Output_Potential(0).write(); - - 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)) diff --git a/source/module_esolver/io_npz.cpp b/source/module_esolver/io_npz.cpp index 6edb2bf4a5..2747f9d286 100644 --- a/source/module_esolver/io_npz.cpp +++ b/source/module_esolver/io_npz.cpp @@ -33,7 +33,7 @@ void ESolver_KS_LCAO::read_mat_npz(std::string& zipname, hamilt::HContai { ModuleBase::TITLE("LCAO_Hamilt","read_mat_npz"); - const Parallel_Orbitals* paraV = &(this->ParaV); + const Parallel_Orbitals* paraV = &(this->pv); #ifdef __USECNPY diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp new file mode 100644 index 0000000000..5b545ff419 --- /dev/null +++ b/source/module_esolver/lcao_before_scf.cpp @@ -0,0 +1,324 @@ +#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/get_pchg.h" +#include "module_io/get_wf.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_parameter/parameter.h" +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif +#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/rho_io.h" +#include "module_io/write_pot.h" +#include "module_io/write_wfc_nao.h" +#include "module_io/read_wfc_nao.h" +#include "module_base/formatter.h" +#ifdef __EXX +#include "module_io/restart_exx_csr.h" +#endif + +namespace ModuleESolver +{ + +template +void ESolver_KS_LCAO::beforesolver(const int istep) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); + + // 1. prepare HS matrices, prepare grid integral + this->set_matrix_grid(this->RA); + + // 2. density matrix extrapolation + + // set the augmented orbitals index. + // after ParaO and GridT, + // this information is used to calculate + // the force. + + // init psi + if (this->psi == nullptr) + { + int nsk = 0; + int ncol = 0; + if (GlobalV::GAMMA_ONLY_LOCAL) + { + nsk = GlobalV::NSPIN; + ncol = this->pv.ncol_bands; + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "lapack" + || GlobalV::KS_SOLVER == "pexsi" + || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp") { + ncol = this->pv.ncol; + } + } + else + { + nsk = this->kv.get_nks(); +#ifdef __MPI + ncol = this->pv.ncol_bands; +#else + ncol = GlobalV::NBANDS; +#endif + } + this->psi = new psi::Psi(nsk, ncol, this->pv.nrow, nullptr); + } + + // init wfc from file + if(istep == 0 && PARAM.inp.init_wfc == "file") + { + if (! ModuleIO::read_wfc_nao(GlobalV::global_readin_dir, this->pv, *(this->psi), this->pelec)) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", + "read wfc nao failed"); + } + } + + // prepare grid in Gint + LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, *this->pw_rho, *this->pw_big); + + // init Hamiltonian + if (this->p_hamilt != nullptr) + { + delete this->p_hamilt; + this->p_hamilt = nullptr; + } + if (this->p_hamilt == nullptr) + { + elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); + this->p_hamilt = new hamilt::HamiltLCAO( + GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, + GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), + &this->pv, + this->pelec->pot, + this->kv, + two_center_bundle_, + DM +#ifdef __EXX + , GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step + , GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr + , GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs +#endif + ); + } + +#ifdef __DEEPKS + // for each ionic step, the overlap must be rebuilt + // since it depends on ionic positions + if (GlobalV::deepks_setorb) + { + const Parallel_Orbitals* pv = &this->pv; + // build and save at beginning + GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD, + *(two_center_bundle_.overlap_orb_alpha)); + + if (PARAM.inp.deepks_out_unittest) + { + GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD); + } + } +#endif + if (PARAM.inp.sc_mag_switch) + { + SpinConstrain& sc = SpinConstrain::getScInstance(); + sc.init_sc(GlobalV::sc_thr, + PARAM.inp.nsc, + PARAM.inp.nsc_min, + PARAM.inp.alpha_trial, + PARAM.inp.sccut, + PARAM.inp.sc_mag_switch, + GlobalC::ucell, + PARAM.inp.sc_file, + GlobalV::NPOL, + &(this->pv), + GlobalV::NSPIN, + this->kv, + GlobalV::KS_SOLVER, + this->phsol, + this->p_hamilt, + this->psi, + this->pelec); + } + //========================================================= + // cal_ux should be called before init_scf because + // the direction of ux is used in noncoline_rho + //========================================================= + if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) + { + GlobalC::ucell.cal_ux(); + } + ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); +} + +template +void ESolver_KS_LCAO::before_scf(const int istep) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); + + if (GlobalC::ucell.cell_parameter_updated) + { + this->init_after_vc(PARAM.inp, GlobalC::ucell); + } + if (GlobalC::ucell.ionic_position_updated) + { + this->CE.update_all_dis(GlobalC::ucell); + this->CE.extrapolate_charge( +#ifdef __MPI + &(GlobalC::Pgrid), +#endif + GlobalC::ucell, + this->pelec->charge, + &(this->sf)); + } + + //---------------------------------------------------------- + // about vdw, jiyy add vdwd3 and linpz add vdwd2 + //---------------------------------------------------------- + auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); + if (vdw_solver != nullptr) + { + this->pelec->f_en.evdw = vdw_solver->get_energy(); + } + + this->beforesolver(istep); + + // Peize Lin add 2016-12-03 +#ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf + if (GlobalC::exx_info.info_ri.real_number) + { + this->exd->exx_beforescf(this->kv, *this->p_chgmix); + } + else + { + this->exc->exx_beforescf(this->kv, *this->p_chgmix); + } +#endif // __EXX + + this->pelec->init_scf(istep, this->sf.strucFac); + if (PARAM.inp.out_chg == 2) + { + for (int is = 0; is < GlobalV::NSPIN; is++) + { + std::stringstream ss; + ss << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; + ModuleIO::write_rho( +#ifdef __MPI + this->pw_big->bz, // bz first, then nbz + this->pw_big->nbz, + this->pw_rho->nplane, + this->pw_rho->startz_current, +#endif + this->pelec->charge->rho[is], + is, + GlobalV::NSPIN, + 0, + ss.str(), + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pelec->eferm.ef, + &(GlobalC::ucell), + 11); + } + } + + ModuleIO::write_pot(GlobalV::out_pot, + GlobalV::NSPIN, + GlobalV::global_out_dir, +#ifdef __MPI + this->pw_big->bz, + this->pw_big->nbz, + this->pw_rho->nplane, + this->pw_rho->startz_current, +#endif + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pelec->pot->get_effective_v()); + + // initalize DMR + // DMR should be same size with Hamiltonian(R) + dynamic_cast*>(this->pelec) + ->get_DM() + ->init_DMR(*(dynamic_cast*>(this->p_hamilt)->getHR())); + // two cases are considered: + // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK + // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros + if(istep > 0) + { + dynamic_cast*>(this->pelec) + ->get_DM() + ->cal_DMR(); + } + + if (PARAM.inp.dm_to_rho) + { + std::string zipname = "output_DM0.npz"; + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + this->read_mat_npz(zipname, *(dm->get_DMR_pointer(1))); + if (GlobalV::NSPIN == 2) + { + zipname = "output_DM1.npz"; + this->read_mat_npz(zipname, *(dm->get_DMR_pointer(2))); + } + + this->pelec->psiToRho(*this->psi); + + this->create_Output_Rho(0, istep).write(); + if (GlobalV::NSPIN == 2) + { + this->create_Output_Rho(1, istep).write(); + } + + return; + } + + // the electron charge density should be symmetrized, + // here is the initialization + Symmetry_rho srho; + for (int is = 0; is < GlobalV::NSPIN; is++) + { + srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + } + + // 1. calculate ewald energy. + // mohan update 2021-02-25 + if (!PARAM.inp.test_skip_ewald) + { + this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho, this->sf.strucFac); + } + + this->p_hamilt->non_first_scf = istep; + + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); + return; +} + + +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_fun.cpp b/source/module_esolver/lcao_fun.cpp new file mode 100644 index 0000000000..a11f0b02a1 --- /dev/null +++ b/source/module_esolver/lcao_fun.cpp @@ -0,0 +1,128 @@ +#include "esolver_ks_lcao.h" + +#include "module_base/global_variable.h" +#include "module_base/tool_title.h" +#include "module_io/dos_nao.h" +#include "module_io/nscf_band.h" +#include "module_io/output_dmk.h" +#include "module_io/output_log.h" +#include "module_io/output_mulliken.h" +#include "module_io/output_sk.h" +#include "module_io/to_qo.h" +#include "module_io/write_HS.h" +#include "module_io/write_istate_info.h" +#include "module_io/write_proj_band_lcao.h" +#include "module_parameter/parameter.h" + +//--------------temporary---------------------------- +#include + +#include "module_base/global_function.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/module_charge/symmetry_rho.h" +#include "module_elecstate/occupy.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag +#include "module_hamilt_lcao/module_dftu/dftu.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/print_info.h" + +//-----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" +//--------------------------------------------------- + +#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_io/io_dmk.h" +#include "module_io/write_dmr.h" +#include "module_io/write_wfc_nao.h" + +namespace ModuleESolver +{ + +//------------------------------------------------------------------------------ +//! the 18th function of ESolver_KS_LCAO: create_Output_Mat_Sparse +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ +template +ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) +{ + return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, + hsolver::HSolverLCAO::out_mat_dh, + hsolver::HSolverLCAO::out_mat_t, + PARAM.inp.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->pv, + this->GK, // mohan add 2024-04-01 + two_center_bundle_, + GlobalC::GridD, // mohan add 2024-04-06 + this->kv, + this->p_hamilt); +} + +//------------------------------------------------------------------------------ +//! the 19th function of ESolver_KS_LCAO: md_skip_out +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ +template +bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, int interval) +{ + if (calculation == "md") + { + if (istep % interval != 0) + { + return true; + } + } + return false; +} + +template +void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) +{ + auto cell_index = CellIndex(GlobalC::ucell.get_atomLabels(), + GlobalC::ucell.get_atomCounts(), + GlobalC::ucell.get_lnchiCounts(), + GlobalV::NSPIN); + auto out_sk = ModuleIO::Output_Sk(this->p_hamilt, + &(this->pv), + GlobalV::NSPIN, + this->kv.get_nks()); + auto out_dmk = ModuleIO::Output_DMK(dynamic_cast*>(this->pelec)->get_DM(), + &(this->pv), + GlobalV::NSPIN, + this->kv.get_nks()); + auto mulp = ModuleIO::Output_Mulliken(&(out_sk), + &(out_dmk), + &(this->pv), + &cell_index, + this->kv.isk, + GlobalV::NSPIN); + auto atom_chg = mulp.get_atom_chg(); + /// used in updating mag info in STRU file + GlobalC::ucell.atom_mulliken = mulp.get_atom_mulliken(atom_chg); + if (print && GlobalV::MY_RANK == 0) + { + /// write the Orbital file + cell_index.write_orb_info(GlobalV::global_out_dir); + /// write mulliken.txt + mulp.write(istep, GlobalV::global_out_dir); + /// write atomic mag info in running log file + mulp.print_atom_mag(atom_chg, GlobalV::ofs_running); + } +} + +//------------------------------------------------------------------------------ +//! the 20th,21th,22th functions of ESolver_KS_LCAO +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ +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_gets.cpp b/source/module_esolver/lcao_gets.cpp new file mode 100644 index 0000000000..8686a26192 --- /dev/null +++ b/source/module_esolver/lcao_gets.cpp @@ -0,0 +1,126 @@ +#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/write_HS_R.h" +#include "module_parameter/parameter.h" +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif +#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/rho_io.h" +#include "module_io/write_pot.h" +#include "module_io/write_wfc_nao.h" +#include "module_io/read_wfc_nao.h" +#include "module_base/formatter.h" +#ifdef __EXX +#include "module_io/restart_exx_csr.h" +#endif + +namespace ModuleESolver +{ + +template <> +void ESolver_KS_LCAO::get_S(void) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", "not implemented for"); +} + +template <> +void ESolver_KS_LCAO, double>::get_S(void) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); + // (1) Find adjacent atoms for each atom. + GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); + + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + GlobalC::GridD, + GlobalC::ucell, + GlobalV::SEARCH_RADIUS, + GlobalV::test_atom_input); + + this->RA.for_2d(this->pv, GlobalV::GAMMA_ONLY_LOCAL); + + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, double>( + &this->pv, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast, double>*>( + this->p_hamilt->ops) + ->contributeHR(); + } + + // mohan add 2024-06-09 + const std::string fn = GlobalV::global_out_dir + "SR.csr"; + + std::cout << " The file is saved in " << fn << std::endl; + + ModuleIO::output_SR(pv, GlobalC::GridD, this->p_hamilt, fn); + + return; +} + +template <> +void ESolver_KS_LCAO, std::complex>::get_S(void) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); + // (1) Find adjacent atoms for each atom. + GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); + + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + GlobalC::GridD, + GlobalC::ucell, + GlobalV::SEARCH_RADIUS, + GlobalV::test_atom_input); + + this->RA.for_2d(this->pv, GlobalV::GAMMA_ONLY_LOCAL); + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, + std::complex>( + &this->pv, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast< + hamilt::OperatorLCAO, std::complex>*>( + this->p_hamilt->ops) + ->contributeHR(); + } + + // mohan add 2024-06-09 + const std::string fn = GlobalV::global_out_dir + "SR.csr"; + + std::cout << " The file is saved in " << fn << std::endl; + + ModuleIO::output_SR(pv, GlobalC::GridD, this->p_hamilt, fn); + + return; +} + +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_init_after_vc.cpp b/source/module_esolver/lcao_init_after_vc.cpp new file mode 100644 index 0000000000..4f7aa982da --- /dev/null +++ b/source/module_esolver/lcao_init_after_vc.cpp @@ -0,0 +1,79 @@ +#include "esolver_ks_lcao.h" + +#include "module_base/global_variable.h" +#include "module_base/tool_title.h" +#include "module_parameter/parameter.h" + +//--------------temporary---------------------------- +#include + +#include "module_base/global_function.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/module_charge/symmetry_rho.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/print_info.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" +//--------------------------------------------------- + + +namespace ModuleESolver +{ + +//------------------------------------------------------------------------------ +//! the 4th function of ESolver_KS_LCAO: init_after_vc +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ +template +void ESolver_KS_LCAO::init_after_vc(const Input_para& inp, UnitCell& ucell) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); + + ESolver_KS::init_after_vc(inp, ucell); + if (inp.mdp.md_prec_level == 2) + { + delete this->pelec; + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), + &(this->kv), + this->kv.get_nks(), + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 + this->pw_rho, + this->pw_big); + + dynamic_cast*>(this->pelec) + ->init_DM(&this->kv, &this->pv, GlobalV::NSPIN); + + GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); + + this->pelec->charge->allocate(GlobalV::NSPIN); + this->pelec->omega = GlobalC::ucell.omega; + + // Initialize the potential. + if (this->pelec->pot == nullptr) + { + this->pelec->pot = new elecstate::Potential(this->pw_rhod, + this->pw_rho, + &GlobalC::ucell, + &(GlobalC::ppcell.vloc), + &(this->sf), + &(this->pelec->f_en.etxc), + &(this->pelec->f_en.vtxc)); + } + } + + ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); + return; +} + + +template class ESolver_KS_LCAO; +template class ESolver_KS_LCAO, double>; +template class ESolver_KS_LCAO, std::complex>; +} diff --git a/source/module_esolver/lcao_nscf.cpp b/source/module_esolver/lcao_nscf.cpp new file mode 100644 index 0000000000..f04a7af130 --- /dev/null +++ b/source/module_esolver/lcao_nscf.cpp @@ -0,0 +1,216 @@ +#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/get_pchg.h" +#include "module_io/get_wf.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_parameter/parameter.h" +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif +#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/rho_io.h" +#include "module_io/write_pot.h" +#include "module_io/write_wfc_nao.h" +#include "module_io/read_wfc_nao.h" +#include "module_base/formatter.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; + if (this->phsol != nullptr) + { + this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER, true); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); + } + + 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 = GlobalV::NSPIN; + const int nbands = GlobalV::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 (!GlobalV::TWO_EFERMI) { + 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 (GlobalV::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); + + 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); // 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 (GlobalV::deepks_out_labels || GlobalV::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 + this->create_Output_Potential(0).write(); + + // 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 new file mode 100644 index 0000000000..8e167507b7 --- /dev/null +++ b/source/module_esolver/lcao_others.cpp @@ -0,0 +1,210 @@ +#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/get_pchg.h" +#include "module_io/get_wf.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_parameter/parameter.h" +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif +#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/rho_io.h" +#include "module_io/write_pot.h" +#include "module_io/write_wfc_nao.h" +#include "module_io/read_wfc_nao.h" +#include "module_base/formatter.h" +#ifdef __EXX +#include "module_io/restart_exx_csr.h" +#endif + +namespace ModuleESolver +{ + +template +void ESolver_KS_LCAO::others(const int istep) +{ + ModuleBase::TITLE("ESolver_KS_LCAO", "others"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); + + const std::string cal_type = GlobalV::CALCULATION; + + if (cal_type == "get_S") { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing the overlap matrix"); + this->get_S(); + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing the overlap matrix"); + + ModuleBase::QUIT(); + + // return; // use 'return' will cause segmentation fault. by mohan + // 2024-06-09 + } else if (cal_type == "test_memory") { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing memory"); + Cal_Test::test_memory(this->pw_rho, + this->pw_wfc, + this->p_chgmix->get_mixing_mode(), + this->p_chgmix->get_mixing_ndim()); + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing memory"); + return; + } + else if (cal_type == "test_neighbour") + { + // test_search_neighbor(); + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing neighbour"); + if (GlobalV::SEARCH_RADIUS < 0) { + std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS + << std::endl; + std::cout << " please make sure search_radius > 0" << std::endl; + } + + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + GlobalC::GridD, + GlobalC::ucell, + GlobalV::SEARCH_RADIUS, + GlobalV::test_atom_input, + true); + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing neighbour"); + return; + } + + this->beforesolver(istep); + // pelec should be initialized before these calculations + this->pelec->init_scf(istep, this->sf.strucFac); + // self consistent calculations for electronic ground state + if (GlobalV::CALCULATION == "nscf") + { + this->nscf(); + } else if (cal_type == "get_pchg") { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); + IState_Charge ISC(this->psi, &(this->pv)); + if (GlobalV::GAMMA_ONLY_LOCAL) { + ISC.begin(this->GG, + this->pelec->charge->rho, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + PARAM.inp.nbands_istate, + PARAM.inp.bands_to_print, + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv); + } else { + ISC.begin(this->GK, + this->pelec->charge->rho, + this->pelec->charge->rhog, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho, + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + PARAM.inp.nbands_istate, + PARAM.inp.bands_to_print, + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv, + PARAM.inp.if_separate_k, + &GlobalC::Pgrid, + this->pelec->charge->ngmc); + } + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); + } else if (cal_type == "get_wf") { + std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); + IState_Envelope IEP(this->pelec); + if (GlobalV::GAMMA_ONLY_LOCAL) + { + IEP.begin(this->psi, + this->pw_rho, + this->pw_wfc, + this->pw_big, + this->pv, + this->GG, + PARAM.inp.out_wfc_pw, + this->wf.out_wfc_r, + this->kv, + GlobalV::nelec, + PARAM.inp.nbands_istate, + PARAM.inp.bands_to_print, + GlobalV::NBANDS, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir); + } + else + { + IEP.begin(this->psi, + this->pw_rho, + this->pw_wfc, + this->pw_big, + this->pv, + this->GK, + PARAM.inp.out_wfc_pw, + this->wf.out_wfc_r, + this->kv, + GlobalV::nelec, + PARAM.inp.nbands_istate, + PARAM.inp.bands_to_print, + GlobalV::NBANDS, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir); + } + std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting wave function"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", + "CALCULATION type not supported"); + } + + ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); + return; +} + +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/pw_fun.cpp b/source/module_esolver/pw_fun.cpp new file mode 100644 index 0000000000..284cc45713 --- /dev/null +++ b/source/module_esolver/pw_fun.cpp @@ -0,0 +1,110 @@ +#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/rho_io.h" +#include "module_io/to_wannier90_pw.h" +#include "module_io/winput.h" +#include "module_io/write_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_hsolver() +{ + this->phsol = new hsolver::HSolverPW(this->pw_wfc, &this->wf); +} +template +void ESolver_KS_PW::deallocate_hsolver() +{ + delete reinterpret_cast*>(this->phsol); + this->phsol = nullptr; +} +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() +{ + delete reinterpret_cast*>(this->p_hamilt); + this->p_hamilt = nullptr; +} + + +template +void ESolver_KS_PW::hamilt2estates(const double ethr) { + if (this->phsol != nullptr) { + hsolver::DiagoIterAssist::need_subspace = false; + hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; + this->phsol->solve(this->p_hamilt, + this->kspw_psi[0], + this->pelec, + PARAM.inp.ks_solver, + PARAM.inp.calculation, + PARAM.inp.basis_type, + PARAM.inp.use_paw, + GlobalV::use_uspp, + GlobalV::RANK_IN_POOL, + GlobalV::NPROC_IN_POOL, + + hsolver::DiagoIterAssist::SCF_ITER, + hsolver::DiagoIterAssist::need_subspace, + hsolver::DiagoIterAssist::PW_DIAG_NMAX, + hsolver::DiagoIterAssist::PW_DIAG_THR, + + true); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", + "HSolver has not been initialed!"); + } +} + +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_init_after_vc.cpp b/source/module_esolver/pw_init_after_vc.cpp new file mode 100644 index 0000000000..3e8359bfef --- /dev/null +++ b/source/module_esolver/pw_init_after_vc.cpp @@ -0,0 +1,208 @@ +#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/rho_io.h" +#include "module_io/to_wannier90_pw.h" +#include "module_io/winput.h" +#include "module_io/write_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::init_after_vc(const Input_para& inp, UnitCell& ucell) { + ModuleBase::TITLE("ESolver_KS_PW", "init_after_vc"); + ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); + + ESolver_KS::init_after_vc(inp, ucell); + + if (inp.mdp.md_prec_level == 2) { + this->pw_wfc->initgrids(ucell.lat0, + ucell.latvec, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz); + + this->pw_wfc->initparameters(false, + inp.ecutwfc, + this->kv.get_nks(), + this->kv.kvec_d.data()); + +#ifdef __MPI + if (PARAM.inp.pw_seed > 0) { + MPI_Allreduce(MPI_IN_PLACE, + &this->pw_wfc->ggecut, + 1, + MPI_DOUBLE, + MPI_MAX, + MPI_COMM_WORLD); + } + // qianrui add 2021-8-13 to make different kpar parameters can get the + // same results +#endif + + this->pw_wfc->setuptransform(); + + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { + this->kv.ngk[ik] = this->pw_wfc->npwk[ik]; + } + + this->pw_wfc->collect_local_pw(inp.erf_ecut, + inp.erf_height, + inp.erf_sigma); + + delete this->phsol; + this->allocate_hsolver(); + + delete this->pelec; + this->pelec + = new elecstate::ElecStatePW(this->pw_wfc, + &(this->chr), + (K_Vectors*)(&(this->kv)), + &ucell, + &GlobalC::ppcell, + this->pw_rhod, + this->pw_rho, + this->pw_big); + + this->pelec->charge->allocate(GlobalV::NSPIN); + + //! setup cell volume + this->pelec->omega = ucell.omega; + + delete this->pelec->pot; + + this->pelec->pot = new elecstate::Potential(this->pw_rhod, + this->pw_rho, + &ucell, + &GlobalC::ppcell.vloc, + &(this->sf), + &(this->pelec->f_en.etxc), + &(this->pelec->f_en.vtxc)); + + // temporary + this->Init_GlobalC(inp, ucell, GlobalC::ppcell); + } else { + GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, + "NON-LOCAL POTENTIAL"); + + this->pw_wfc->initgrids(ucell.lat0, + ucell.latvec, + this->pw_wfc->nx, + this->pw_wfc->ny, + this->pw_wfc->nz); + + this->pw_wfc->initparameters(false, + PARAM.inp.ecutwfc, + this->kv.get_nks(), + this->kv.kvec_d.data()); + + this->pw_wfc->collect_local_pw(inp.erf_ecut, + inp.erf_height, + inp.erf_sigma); + + this->p_wf_init->make_table(this->kv.get_nks(), &this->sf); + } + +#ifdef USE_PAW + if (GlobalV::use_paw) { + GlobalC::paw_cell.set_libpaw_ecut(PARAM.inp.ecutwfc / 2.0, + PARAM.inp.ecutwfc / 2.0); // in Hartree + GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx, + this->pw_wfc->ny, + this->pw_wfc->nz, + this->pw_wfc->nx, + this->pw_wfc->ny, + this->pw_wfc->nz, + this->pw_wfc->startz, + this->pw_wfc->numz); + +#ifdef __MPI + if (GlobalV::RANK_IN_POOL == 0) { + GlobalC::paw_cell.prepare_paw(); + } +#else + GlobalC::paw_cell.prepare_paw(); +#endif + GlobalC::paw_cell.set_sij(); + + std::vector> rhoijp; + std::vector> rhoijselect; + std::vector nrhoijsel; +#ifdef __MPI + if (GlobalV::RANK_IN_POOL == 0) { + GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); + + for (int iat = 0; iat < ucell.nat; iat++) { + GlobalC::paw_cell.set_rhoij(iat, + nrhoijsel[iat], + rhoijselect[iat].size(), + rhoijselect[iat].data(), + rhoijp[iat].data()); + } + } +#else + GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); + + for (int iat = 0; iat < ucell.nat; iat++) { + GlobalC::paw_cell.set_rhoij(iat, + nrhoijsel[iat], + rhoijselect[iat].size(), + rhoijselect[iat].data(), + rhoijp[iat].data()); + } +#endif + } +#endif + + ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); +} + +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_init_globalc.cpp b/source/module_esolver/pw_init_globalc.cpp new file mode 100644 index 0000000000..0c35e5dbe0 --- /dev/null +++ b/source/module_esolver/pw_init_globalc.cpp @@ -0,0 +1,113 @@ +#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/rho_io.h" +#include "module_io/to_wannier90_pw.h" +#include "module_io/winput.h" +#include "module_io/write_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::Init_GlobalC(const Input_para& inp, + UnitCell& ucell, + pseudopot_cell_vnl& ppcell) { + // GlobalC is a historically left-over namespace, it is used to store global + // classes, including: pseudopot_cell_vnl: pseudopotential in cell, V + // non-local UnitCell: cell information with atomic properties Grid_Driver: + // Parallel_Grid: + // Parallel_Kpoints: + // Restart: + // Exx_Info: + // Exx_Lip: + + // GlobalC would be refactored out in the future. If there is better idea + // about how to organize information stored in classes above, please feel + // free to discuss with issue or pull request. + + if (this->psi != nullptr) { + delete this->psi; + } + + //! init pseudopotential + ppcell.init(ucell.ntype, &this->sf, this->pw_wfc); + + //! initalize local pseudopotential + ppcell.init_vloc(ppcell.vloc, this->pw_rhod); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); + + //! Initalize non-local pseudopotential + ppcell.init_vnl(ucell, this->pw_rhod); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); + + //! Allocate psi + this->p_wf_init->allocate_psi(this->psi, + this->kv.get_nkstot(), + this->kv.get_nks(), + this->kv.ngk.data(), + this->pw_wfc->npwk_max, + &(this->sf)); + + this->kspw_psi + = GlobalV::device_flag == "gpu" || GlobalV::precision_flag == "single" + ? new psi::Psi(this->psi[0]) + : reinterpret_cast*>(this->psi); + + if (GlobalV::precision_flag == "single") { + ModuleBase::Memory::record("Psi_single", + sizeof(T) * this->psi[0].size()); + } + + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); +} + + + +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 new file mode 100644 index 0000000000..88c29cc04c --- /dev/null +++ b/source/module_esolver/pw_nscf.cpp @@ -0,0 +1,175 @@ +#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/rho_io.h" +#include "module_io/to_wannier90_pw.h" +#include "module_io/winput.h" +#include "module_io/write_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 = GlobalV::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 / GlobalV::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 = GlobalV::NSPIN; + const int nbands = GlobalV::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 (!GlobalV::TWO_EFERMI) { + 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 + this->create_Output_Potential(0).write(); + + 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 new file mode 100644 index 0000000000..ab0352fe3b --- /dev/null +++ b/source/module_esolver/pw_others.cpp @@ -0,0 +1,89 @@ +#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/rho_io.h" +#include "module_io/to_wannier90_pw.h" +#include "module_io/winput.h" +#include "module_io/write_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::others(const int istep) { + ModuleBase::TITLE("ESolver_KS_PW", "others"); + + const std::string cal_type = GlobalV::CALCULATION; + + if (cal_type == "test_memory") { + Cal_Test::test_memory(this->pw_rho, + this->pw_wfc, + this->p_chgmix->get_mixing_mode(), + this->p_chgmix->get_mixing_ndim()); + } else if (cal_type == "gen_bessel") { + Numerical_Descriptor nc; + nc.output_descriptor(this->psi[0], + PARAM.inp.bessel_descriptor_lmax, + 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(); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", + "CALCULATION type not supported"); + } + + 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/set_matrix_grid.cpp b/source/module_esolver/set_matrix_grid.cpp index ddc56a6b9b..207e9e7198 100644 --- a/source/module_esolver/set_matrix_grid.cpp +++ b/source/module_esolver/set_matrix_grid.cpp @@ -82,16 +82,17 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) dpsi_u.shrink_to_fit(); d2psi_u.clear(); d2psi_u.shrink_to_fit(); + // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - Parallel_Orbitals* pv = &this->ParaV; - ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL); + ra.for_2d(this->pv, GlobalV::GAMMA_ONLY_LOCAL); + if (!GlobalV::GAMMA_ONLY_LOCAL) { // need to first calculae lgd. // using GridT.init. - this->GridT.cal_nnrg(pv); + this->GridT.cal_nnrg(&this->pv); } ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid"); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index dcb8d3d9e3..4a6398aee2 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -140,8 +140,8 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol // setup_wd_division is not need to be covered in #ifdef __MPI, see its implementation LR_Util::setup_2d_division(this->paraMat_, 1, this->nbasis, this->nbasis); - this->paraMat_.atom_begin_row = std::move(ks_sol.ParaV.atom_begin_row); - this->paraMat_.atom_begin_col = std::move(ks_sol.ParaV.atom_begin_col); + this->paraMat_.atom_begin_row = std::move(ks_sol.pv.atom_begin_row); + this->paraMat_.atom_begin_col = std::move(ks_sol.pv.atom_begin_col); this->paraMat_.iat2iwt_ = ucell.get_iat2iwt(); LR_Util::setup_2d_division(this->paraC_, 1, this->nbasis, this->nbands @@ -165,7 +165,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol const int start_band = this->nocc_max - this->nocc; for (int ik = 0;ik < this->kv.get_nks();++ik) { - Cpxgemr2d(this->nbasis, this->nbands, &(*ks_sol.psi)(ik, 0, 0), 1, start_band + 1, ks_sol.ParaV.desc_wfc, + Cpxgemr2d(this->nbasis, this->nbands, &(*ks_sol.psi)(ik, 0, 0), 1, start_band + 1, ks_sol.pv.desc_wfc, &(*this->psi_ks)(ik, 0, 0), 1, 1, this->paraC_.desc, this->paraC_.blacs_ctxt); for (int ib = 0;ib < this->nbands;++ib) { this->eig_ks(ik, ib) = ks_sol.pelec->ekb(ik, start_band + ib); } }