diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index a861a5a7d1..b123a354fe 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -223,7 +223,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa ucell.nat, orb_.Alpha[0].getTotal_nchi() * ucell.nat, ld.lmaxd, - ld.inl_l, + ld.inl2l, *orb_.Alpha, ld.pdm); } diff --git a/source/module_esolver/lcao_after_scf.cpp b/source/module_esolver/lcao_after_scf.cpp index d8284ae7ab..73fb309a05 100644 --- a/source/module_esolver/lcao_after_scf.cpp +++ b/source/module_esolver/lcao_after_scf.cpp @@ -1,5 +1,4 @@ #include "esolver_ks_lcao.h" - #include "module_base/formatter.h" #include "module_base/global_variable.h" #include "module_base/tool_title.h" @@ -28,7 +27,6 @@ #include "module_io/write_wfc_nao.h" #include "module_parameter/parameter.h" - //--------------temporary---------------------------- #include "module_base/global_function.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" @@ -40,11 +38,9 @@ #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" - -//mohan add 20250302 +// mohan add 20250302 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h" - #include #ifdef __EXX #include "module_io/restart_exx_csr.h" @@ -95,7 +91,6 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const //------------------------------------------------------------------ ESolver_KS::after_scf(ucell, istep, conv_esolver); - //------------------------------------------------------------------ //! 3) write density matrix for sparse matrix in LCAO basis //------------------------------------------------------------------ @@ -133,13 +128,10 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const //------------------------------------------------------------------ if (PARAM.inp.calculation != "nscf") { - if (GlobalC::exx_info.info_global.cal_exx && - PARAM.inp.out_chg[0] && - istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14 + if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.out_chg[0] + && istep % PARAM.inp.out_interval == 0) // Peize Lin add if 2022.11.14 { - const std::string file_name_exx = PARAM.globalv.global_out_dir + - "HexxR" + - std::to_string(GlobalV::MY_RANK); + const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); if (GlobalC::exx_info.info_ri.real_number) { ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exd->get_Hexxs()); @@ -152,52 +144,52 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const } #endif - //------------------------------------------------------------------ - // 6) write Hamiltonian and Overlap matrix in LCAO basis - //------------------------------------------------------------------ - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (PARAM.inp.out_mat_hs[0]) - { - this->p_hamilt->updateHk(ik); - } - bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core - // calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) - { - hamilt::MatrixBlock h_mat; - hamilt::MatrixBlock s_mat; - - this->p_hamilt->matrix(h_mat, s_mat); - - if (PARAM.inp.out_mat_hs[0]) - { - ModuleIO::save_mat(istep, - h_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "H", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - ModuleIO::save_mat(istep, - s_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "S", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - } - } - } + //------------------------------------------------------------------ + // 6) write Hamiltonian and Overlap matrix in LCAO basis + //------------------------------------------------------------------ + for (int ik = 0; ik < this->kv.get_nks(); ++ik) + { + if (PARAM.inp.out_mat_hs[0]) + { + this->p_hamilt->updateHk(ik); + } + bool bit = false; // LiuXh, 2017-03-21 + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) + { + hamilt::MatrixBlock h_mat; + hamilt::MatrixBlock s_mat; + + this->p_hamilt->matrix(h_mat, s_mat); + + if (PARAM.inp.out_mat_hs[0]) + { + ModuleIO::save_mat(istep, + h_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "H", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + ModuleIO::save_mat(istep, + s_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "S", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + } + } + } //------------------------------------------------------------------ // 7) write electronic wavefunctions in LCAO basis @@ -224,23 +216,22 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); deepks_interface.out_deepks_labels(this->pelec->f_en.etot, - this->pelec->klist->get_nks(), - ucell.nat, - PARAM.globalv.nlocal, - this->pelec->ekb, - this->pelec->klist->kvec_d, - ucell, - orb_, - this->gd, - &(this->pv), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM(), - p_ham_deepks, - GlobalV::MY_RANK); + this->pelec->klist->get_nks(), + ucell.nat, + PARAM.globalv.nlocal, + this->pelec->ekb, + this->pelec->klist->kvec_d, + ucell, + orb_, + this->gd, + &(this->pv), + *(this->psi), + dynamic_cast*>(this->pelec)->get_DM(), + p_ham_deepks, + GlobalV::MY_RANK); } #endif - //------------------------------------------------------------------ //! 9) Perform RDMFT calculations // rdmft, added by jghan, 2024-10-17 @@ -269,7 +260,6 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const double etot_rdmft = this->rdmft_solver.run(dedocc, dedwfc); } - #ifdef __EXX //------------------------------------------------------------------ // 10) Write RPA information in LCAO basis @@ -287,7 +277,6 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const } #endif - //------------------------------------------------------------------ // 11) write HR in npz format in LCAO basis //------------------------------------------------------------------ @@ -309,59 +298,59 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const } } - //------------------------------------------------------------------ - // 12) write density matrix in the 'npz' format in LCAO basis - //------------------------------------------------------------------ - if (PARAM.inp.out_dm_npz) - { - const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - std::string zipname = "output_DM0.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(1))); - - if (PARAM.inp.nspin == 2) - { - zipname = "output_DM1.npz"; - ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(2))); - } - } - - //------------------------------------------------------------------ - //! 13) Print out information every 'out_interval' steps. - //------------------------------------------------------------------ - if (PARAM.inp.calculation != "md" || istep % PARAM.inp.out_interval == 0) - { - //! Print out sparse matrix - ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2, - PARAM.inp.out_mat_dh, - PARAM.inp.out_mat_t, - PARAM.inp.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->pv, - this->GK, - two_center_bundle_, - orb_, - ucell, - this->gd, - this->kv, - this->p_hamilt); - - //! Perform Mulliken charge analysis in LCAO basis - if (PARAM.inp.out_mul) - { - ModuleIO::cal_mag(&(this->pv), - this->p_hamilt, - this->kv, - this->pelec, - this->two_center_bundle_, - this->orb_, - ucell, - this->gd, - istep, - true); - } - } + //------------------------------------------------------------------ + // 12) write density matrix in the 'npz' format in LCAO basis + //------------------------------------------------------------------ + if (PARAM.inp.out_dm_npz) + { + const elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); + std::string zipname = "output_DM0.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(1))); + + if (PARAM.inp.nspin == 2) + { + zipname = "output_DM1.npz"; + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(2))); + } + } + + //------------------------------------------------------------------ + //! 13) Print out information every 'out_interval' steps. + //------------------------------------------------------------------ + if (PARAM.inp.calculation != "md" || istep % PARAM.inp.out_interval == 0) + { + //! Print out sparse matrix + ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2, + PARAM.inp.out_mat_dh, + PARAM.inp.out_mat_t, + PARAM.inp.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->pv, + this->GK, + two_center_bundle_, + orb_, + ucell, + this->gd, + this->kv, + this->p_hamilt); + + //! Perform Mulliken charge analysis in LCAO basis + if (PARAM.inp.out_mul) + { + ModuleIO::cal_mag(&(this->pv), + this->p_hamilt, + this->kv, + this->pelec, + this->two_center_bundle_, + this->orb_, + ucell, + this->gd, + istep, + true); + } + } //------------------------------------------------------------------ //! 14) Print out atomic magnetization in LCAO basis @@ -379,92 +368,90 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const //! 15) Print out kinetic matrix in LCAO basis //------------------------------------------------------------------ if (PARAM.inp.out_mat_tk[0]) - { - hamilt::HS_Matrix_K hsk(&pv, true); - hamilt::HContainer hR(&pv); - hamilt::Operator* ekinetic - = new hamilt::EkineticNew>(&hsk, - this->kv.kvec_d, - &hR, - &ucell, - orb_.cutoffs(), - &this->gd, - two_center_bundle_.kinetic_orb.get()); - - const int nspin_k = (PARAM.inp.nspin == 2 ? 2 : 1); - for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik) - { - ekinetic->init(ik); - ModuleIO::save_mat(0, - hsk.get_hk(), - PARAM.globalv.nlocal, - false, - PARAM.inp.out_mat_tk[1], - 1, - PARAM.inp.out_app_flag, - "T", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); - } - - delete ekinetic; - } - - //------------------------------------------------------------------ - //! 16) wannier90 interface in LCAO basis + { + hamilt::HS_Matrix_K hsk(&pv, true); + hamilt::HContainer hR(&pv); + hamilt::Operator* ekinetic + = new hamilt::EkineticNew>(&hsk, + this->kv.kvec_d, + &hR, + &ucell, + orb_.cutoffs(), + &this->gd, + two_center_bundle_.kinetic_orb.get()); + + const int nspin_k = (PARAM.inp.nspin == 2 ? 2 : 1); + for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik) + { + ekinetic->init(ik); + ModuleIO::save_mat(0, + hsk.get_hk(), + PARAM.globalv.nlocal, + false, + PARAM.inp.out_mat_tk[1], + 1, + PARAM.inp.out_app_flag, + "T", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); + } + + delete ekinetic; + } + + //------------------------------------------------------------------ + //! 16) wannier90 interface in LCAO basis // added by jingan in 2018.11.7 //------------------------------------------------------------------ if (PARAM.inp.calculation == "nscf" && PARAM.inp.towannier90) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wave function to Wannier90"); if (PARAM.inp.wannier_method == 1) - { - toWannier90_LCAO_IN_PW 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.set_tpiba_omega(ucell.tpiba, ucell.omega); - wan.calculate(ucell, - this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->sf, - this->kv, - this->psi, - &(this->pv)); - } - else if (PARAM.inp.wannier_method == 2) - { - toWannier90_LCAO 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, - orb_); - - wan.calculate(ucell, this->gd, this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); - } + { + toWannier90_LCAO_IN_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.set_tpiba_omega(ucell.tpiba, ucell.omega); + wan.calculate(ucell, + this->pelec->ekb, + this->pw_wfc, + this->pw_big, + this->sf, + this->kv, + this->psi, + &(this->pv)); + } + else if (PARAM.inp.wannier_method == 2) + { + toWannier90_LCAO 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, + orb_); + + wan.calculate(ucell, this->gd, this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); + } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); } //------------------------------------------------------------------ - //! 17) berry phase calculations in LCAO basis + //! 17) berry phase calculations in LCAO basis // added by jingan //------------------------------------------------------------------ - if (PARAM.inp.calculation == "nscf" && - berryphase::berry_phase_flag && - ModuleSymmetry::Symmetry::symm_flag != 1) + if (PARAM.inp.calculation == "nscf" && berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation"); berryphase bp(&(this->pv)); bp.lcao_init(ucell, this->gd, this->kv, this->GridT, orb_); - // additional step before calling macroscopic_polarization + // additional step before calling macroscopic_polarization bp.Macroscopic_polarization(ucell, this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); } @@ -472,19 +459,19 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const //------------------------------------------------------------------ //! 18) calculate quasi-orbitals in LCAO basis //------------------------------------------------------------------ - if (PARAM.inp.qo_switch) - { - toQO tqo(PARAM.inp.qo_basis, PARAM.inp.qo_strategy, PARAM.inp.qo_thr, PARAM.inp.qo_screening_coeff); - tqo.initialize(PARAM.globalv.global_out_dir, - PARAM.inp.pseudo_dir, - PARAM.inp.orbital_dir, - &ucell, - this->kv.kvec_d, - GlobalV::ofs_running, - GlobalV::MY_RANK, - GlobalV::NPROC); - tqo.calculate(); - } + if (PARAM.inp.qo_switch) + { + toQO tqo(PARAM.inp.qo_basis, PARAM.inp.qo_strategy, PARAM.inp.qo_thr, PARAM.inp.qo_screening_coeff); + tqo.initialize(PARAM.globalv.global_out_dir, + PARAM.inp.pseudo_dir, + PARAM.inp.orbital_dir, + &ucell, + this->kv.kvec_d, + GlobalV::ofs_running, + GlobalV::MY_RANK, + GlobalV::NPROC); + tqo.calculate(); + } //------------------------------------------------------------------ //! 19) Clean up RA, which is used to serach for adjacent atoms diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index e5db652a05..ef5b86ef26 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -216,6 +216,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, DM_in, ld_in); this->getOperator()->add(deepks); + this->V_delta_R = dynamic_cast>*>(deepks)->get_V_delta_R(); } #endif @@ -340,6 +341,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, DM_in, ld_in); this->getOperator()->add(deepks); + this->V_delta_R = dynamic_cast>*>(deepks)->get_V_delta_R(); } #endif // TDDFT_velocity_gague diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index a52dcda24d..74a9f08d4a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -104,6 +104,13 @@ class HamiltLCAO : public Hamilt { return this->sR; } +#ifdef __DEEPKS + /// get V_delta_R pointer of *this->V_delta_R, which is a HContainer and contains V_delta(R) + HContainer*& get_V_delta_R() + { + return this->V_delta_R; + } +#endif /// refresh the status of HR void refresh() override; @@ -131,6 +138,10 @@ class HamiltLCAO : public Hamilt HContainer* hR = nullptr; HContainer* sR = nullptr; +#ifdef __DEEPKS + HContainer* V_delta_R = nullptr; +#endif + HS_Matrix_K* hsk = nullptr; // special case for NSPIN=2 , data of HR should be separated into two parts // save them in this->hRS2; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp index b7d3b24c1a..115033645b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp @@ -169,7 +169,7 @@ void hamilt::DeePKS>::contributeHR() DeePKS_domain::cal_pdm(this->ld->init_pdm, inlmax, this->ld->lmaxd, - this->ld->inl_l, + this->ld->inl2l, this->ld->inl_index, this->DM, this->ld->phialpha, @@ -182,7 +182,7 @@ void hamilt::DeePKS>::contributeHR() std::vector descriptor; DeePKS_domain::cal_descriptor(this->ucell->nat, inlmax, - this->ld->inl_l, + this->ld->inl2l, this->ld->pdm, descriptor, this->ld->des_per_atom); @@ -193,7 +193,7 @@ void hamilt::DeePKS>::contributeHR() this->ld->nmaxd, inlmax, this->ld->des_per_atom, - this->ld->inl_l, + this->ld->inl2l, descriptor, this->ld->gedm, this->ld->E_delta, @@ -204,7 +204,7 @@ void hamilt::DeePKS>::contributeHR() DeePKS_domain::cal_edelta_gedm(this->ucell->nat, inlmax, this->ld->des_per_atom, - this->ld->inl_l, + this->ld->inl2l, descriptor, this->ld->pdm, this->ld->model_deepks, @@ -488,13 +488,6 @@ void hamilt::DeePKS>::cal_HR_IJR(const double* hr_i } } -template -inline void get_h_delta_k(int ik, TK*& h_delta_k, LCAO_Deepks* ld_in) -{ - h_delta_k = ld_in->V_delta[ik].data(); - return; -} - // contributeHk() template void hamilt::DeePKS>::contributeHk(int ik) @@ -502,8 +495,7 @@ void hamilt::DeePKS>::contributeHk(int ik) ModuleBase::TITLE("DeePKS", "contributeHk"); ModuleBase::timer::tick("DeePKS", "contributeHk"); - TK* h_delta_k = nullptr; - get_h_delta_k(ik, h_delta_k, this->ld); + TK* h_delta_k = this->ld->V_delta[ik].data(); // set SK to zero and then calculate SK for each k vector ModuleBase::GlobalFunc::ZEROS(h_delta_k, this->hsk->get_size()); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h index 8f2b787f28..ccd1cc6f9b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h @@ -58,6 +58,11 @@ class DeePKS> : public OperatorLCAO * @param ik: the index of k-point */ virtual void contributeHk(int ik) override; + + HContainer* get_V_delta_R() const + { + return this->V_delta_R; + } #endif private: diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp index b792705d9e..1340f04e0c 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp @@ -20,7 +20,6 @@ template LCAO_Deepks::LCAO_Deepks() { inl_index = new ModuleBase::IntArray[1]; - inl_l = nullptr; gedm = nullptr; this->phialpha.resize(1); } @@ -30,7 +29,6 @@ template LCAO_Deepks::~LCAO_Deepks() { delete[] inl_index; - delete[] inl_l; //=======1. to use deepks, pdm is required========== pdm.clear(); @@ -106,7 +104,7 @@ void LCAO_Deepks::init(const LCAO_Orbitals& orb, // init pdm for (int inl = 0; inl < this->inlmax; inl++) { - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; pdm_size += nm * nm; this->pdm[inl] = torch::zeros({nm, nm}, torch::kFloat64); } @@ -120,9 +118,9 @@ void LCAO_Deepks::init(const LCAO_Orbitals& orb, pdm_size = pdm_size * pdm_size; this->des_per_atom = pdm_size; ofs << " Equivariant version, size of pdm matrices : " << pdm_size << std::endl; - for (int inl = 0; inl < this->inlmax; inl++) + for (int iat = 0; iat < nat; iat++) { - this->pdm[inl] = torch::zeros({pdm_size}, torch::kFloat64); + this->pdm[iat] = torch::zeros({pdm_size}, torch::kFloat64); } } @@ -142,9 +140,7 @@ void LCAO_Deepks::init_index(const int ntype, { delete[] this->inl_index; this->inl_index = new ModuleBase::IntArray[ntype]; - delete[] this->inl_l; - this->inl_l = new int[this->inlmax]; - ModuleBase::GlobalFunc::ZEROS(this->inl_l, this->inlmax); + this->inl2l.resize(this->inlmax, 0); int inl = 0; int alpha = 0; @@ -162,7 +158,7 @@ void LCAO_Deepks::init_index(const int ntype, for (int n = 0; n < orb.Alpha[0].getNchi(l); n++) { this->inl_index[it](ia, l, n) = inl; - this->inl_l[inl] = l; + this->inl2l[inl] = l; inl++; } } diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h index ac40e05897..e3058a15a7 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h @@ -73,7 +73,7 @@ class LCAO_Deepks int inlmax = 0; // tot. number {i,n,l} - atom, n, l int n_descriptor; // natoms * des_per_atom, size of descriptor(projector) basis set int des_per_atom; // \sum_L{Nchi(L)*(2L+1)} - int* inl_l; // inl_l[inl_index] = l of descriptor with inl_index + std::vector inl2l; // inl2l[inl] = l of descriptor with inl_index ModuleBase::IntArray* inl_index; // caoyu add 2021-05-07 bool init_pdm = false; // for DeePKS NSCF calculation, set init_pdm to skip the calculation of pdm in SCF iteration diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp index f84359bb09..a1af889b57 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -6,6 +6,7 @@ #include "module_base/tool_title.h" #include "module_elecstate/cal_dm.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_hamilt_lcao/module_hcontainer/output_hcontainer.h" #include "module_parameter/parameter.h" template @@ -43,7 +44,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, const int nmaxd = ld->nmaxd; const int des_per_atom = ld->des_per_atom; - const int* inl_l = ld->inl_l; + const std::vector inl2l = ld->inl2l; const ModuleBase::IntArray* inl_index = ld->inl_index; const std::vector*> phialpha = ld->phialpha; @@ -61,16 +62,16 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // this part is for integrated test of deepks // so it is printed no matter even if deepks_out_labels is not used DeePKS_domain::cal_pdm< - TK>(init_pdm, inlmax, lmaxd, inl_l, inl_index, dm, phialpha, ucell, orb, GridD, *ParaV, pdm); + TK>(init_pdm, inlmax, lmaxd, inl2l, inl_index, dm, phialpha, ucell, orb, GridD, *ParaV, pdm); - DeePKS_domain::check_pdm(inlmax, inl_l, pdm); // print out the projected dm for NSCF calculaiton + DeePKS_domain::check_pdm(inlmax, inl2l, pdm); // print out the projected dm for NSCF calculaiton std::vector descriptor; - DeePKS_domain::cal_descriptor(nat, inlmax, inl_l, pdm, descriptor, + DeePKS_domain::cal_descriptor(nat, inlmax, inl2l, pdm, descriptor, des_per_atom); // final descriptor DeePKS_domain::check_descriptor(inlmax, des_per_atom, - inl_l, + inl2l, ucell, PARAM.globalv.global_out_dir, descriptor, @@ -81,7 +82,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, LCAO_deepks_io::save_npy_d(nat, des_per_atom, inlmax, - inl_l, + inl2l, PARAM.inp.deepks_equiv, descriptor, PARAM.globalv.global_out_dir, @@ -99,7 +100,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, nmaxd, inlmax, des_per_atom, - inl_l, + inl2l, descriptor, ld->gedm, E_delta, @@ -110,7 +111,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, DeePKS_domain::cal_edelta_gedm(nat, inlmax, des_per_atom, - inl_l, + inl2l, descriptor, pdm, ld->model_deepks, @@ -130,7 +131,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, std::vector gevdm; if (PARAM.inp.deepks_scf) { - DeePKS_domain::cal_gevdm(nat, inlmax, inl_l, pdm, gevdm); + DeePKS_domain::cal_gevdm(nat, inlmax, inl2l, pdm, gevdm); } // Energy Part @@ -162,7 +163,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dm_vec, ucell, orb, *ParaV, GridD, gdmx); torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, inlmax, des_per_atom, inl_l, gevdm, gdmx, gvx, rank); + DeePKS_domain::cal_gvx(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmx, gvx, rank); const std::string file_gradvx = PARAM.globalv.global_out_dir + "deepks_gradvx.npy"; LCAO_deepks_io::save_tensor2npy(file_gradvx, gvx, rank); @@ -186,7 +187,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dm_vec, ucell, orb, *ParaV, GridD, gdmepsl); torch::Tensor gvepsl; - DeePKS_domain::cal_gvepsl(ucell.nat, inlmax, des_per_atom, inl_l, gevdm, gdmepsl, gvepsl, rank); + DeePKS_domain::cal_gvepsl(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmepsl, gvepsl, rank); const std::string file_gvepsl = PARAM.globalv.global_out_dir + "deepks_gvepsl.npy"; LCAO_deepks_io::save_tensor2npy(file_gvepsl, gvepsl, rank); @@ -252,7 +253,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, inlmax, nat, nks, - inl_l, + inl2l, kvec_d, phialpha, gevdm, @@ -278,17 +279,44 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } // end deepks_scf == 0 } // end bandgap label - // H(R) matrix part, not realized now - if (true) // should be modified later! + // H(R) matrix part, for HR, base will not be calculated since they are HContainer objects + if (PARAM.inp.deepks_v_delta < 0) { - const std::string file_hr = PARAM.globalv.global_out_dir + "deepks_hr.npy"; - const hamilt::HContainer& hR = *(p_ham->getHR()); + // set the output + const double sparse_threshold = 1e-10; + const int precision = 8; + const std::string file_hrtot = PARAM.globalv.global_out_dir + "deepks_hrtot.csr"; + hamilt::HContainer* hR_tot = (p_ham->getHR()); - // How to save H(R)? + if (rank == 0) + { + std::ofstream ofs_hr(file_hrtot, std::ios::out); + ofs_hr << "Matrix Dimension of H(R): " << hR_tot->get_nbasis() << std::endl; + ofs_hr << "Matrix number of H(R): " << hR_tot->size_R_loop() << std::endl; + hamilt::Output_HContainer out_hr(hR_tot, ofs_hr, sparse_threshold, precision); + out_hr.write(); + ofs_hr.close(); + } + + if (PARAM.inp.deepks_scf) + { + const std::string file_vdeltar = PARAM.globalv.global_out_dir + "deepks_hrdelta.csr"; + hamilt::HContainer* h_deltaR = p_ham->get_V_delta_R(); + + if (rank == 0) + { + std::ofstream ofs_hr(file_vdeltar, std::ios::out); + ofs_hr << "Matrix Dimension of H_delta(R): " << h_deltaR->get_nbasis() << std::endl; + ofs_hr << "Matrix number of H_delta(R): " << h_deltaR->size_R_loop() << std::endl; + hamilt::Output_HContainer out_hr(h_deltaR, ofs_hr, sparse_threshold, precision); + out_hr.write(); + ofs_hr.close(); + } + } } // H(k) matrix part - if (PARAM.inp.deepks_v_delta) + if (PARAM.inp.deepks_v_delta > 0) { std::vector h_tot(nks); std::vector> h_mat(nks, std::vector(ParaV->nloc)); @@ -338,7 +366,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, inlmax, nat, nks, - inl_l, + inl2l, kvec_d, phialpha, gevdm, @@ -387,7 +415,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, { LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, ParaV->nrow, dm->get_DMK_vector()); - DeePKS_domain::check_gedm(inlmax, inl_l, ld->gedm); + DeePKS_domain::check_gedm(inlmax, inl2l, ld->gedm); std::ofstream ofs("E_delta_bands.dat"); ofs << std::setprecision(10) << e_delta_band; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp index d69f9da9d4..c58672b180 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp @@ -77,7 +77,7 @@ void LCAO_deepks_io::load_npy_gedm(const int nat, void LCAO_deepks_io::save_npy_d(const int nat, const int des_per_atom, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const bool deepks_equiv, const std::vector& descriptor, const std::string& out_dir, @@ -98,7 +98,7 @@ void LCAO_deepks_io::save_npy_d(const int nat, for (int inl = 0; inl < inlmax; ++inl) { auto accessor = descriptor[inl].accessor(); - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; for (int im = 0; im < nm; im++) { npy_des.push_back(accessor[im]); diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h index 8ab0725a7f..4c536a4999 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h @@ -42,7 +42,7 @@ void load_npy_gedm(const int nat, const int des_per_atom, double** gedm, double& void save_npy_d(const int nat, const int des_per_atom, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const bool deepks_equiv, const std::vector& descriptor, const std::string& out_dir, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp b/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp index 7a3c18752b..7be7f58fc2 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp @@ -12,7 +12,7 @@ // Dimension is different for each inl, so there's a vector of tensors void DeePKS_domain::cal_gevdm(const int nat, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const std::vector& pdm, std::vector& gevdm) { @@ -26,7 +26,7 @@ void DeePKS_domain::cal_gevdm(const int nat, for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; // repeat each block for nm times in an additional dimension torch::Tensor tmp_x = pdm[inl].reshape({nm, nm}).unsqueeze(0).repeat({nm, 1, 1}); // torch::Tensor tmp_y = std::get<0>(torch::symeig(tmp_x, true)); @@ -127,7 +127,7 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, const int nmaxd, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& descriptor, double** gedm, double& E_delta, @@ -139,7 +139,7 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, LCAO_deepks_io::save_npy_d(nat, des_per_atom, inlmax, - inl_l, + inl2l, PARAM.inp.deepks_equiv, descriptor, PARAM.globalv.global_out_dir, @@ -169,7 +169,7 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, void DeePKS_domain::cal_edelta_gedm(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& descriptor, const std::vector& pdm, torch::jit::script::Module& model_deepks, @@ -201,7 +201,7 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, // gedm_tensor(Hartree) to gedm(Ry) for (int inl = 0; inl < inlmax; ++inl) { - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; auto accessor = gedm_tensor[inl].accessor(); for (int m1 = 0; m1 < nm; ++m1) { @@ -216,13 +216,13 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, return; } -void DeePKS_domain::check_gedm(const int inlmax, const int* inl_l, double** gedm) +void DeePKS_domain::check_gedm(const int inlmax, const std::vector& inl2l, double** gedm) { std::ofstream ofs("gedm.dat"); for (int inl = 0; inl < inlmax; inl++) { - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; for (int m1 = 0; m1 < nm; ++m1) { for (int m2 = 0; m2 < nm; ++m2) diff --git a/source/module_hamilt_lcao/module_deepks/deepks_basic.h b/source/module_hamilt_lcao/module_deepks/deepks_basic.h index 194d0ae014..bd8ff9ff08 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_basic.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_basic.h @@ -30,7 +30,7 @@ void load_model(const std::string& model_file, torch::jit::script::Module& model // calculate gevdm void cal_gevdm(const int nat, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const std::vector& pdm, std::vector& gevdm); @@ -38,19 +38,19 @@ void cal_gevdm(const int nat, void cal_edelta_gedm(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& descriptor, const std::vector& pdm, torch::jit::script::Module& model_deepks, double** gedm, double& E_delta); -void check_gedm(const int inlmax, const int* inl_l, double** gedm); +void check_gedm(const int inlmax, const std::vector& inl2l, double** gedm); void cal_edelta_gedm_equiv(const int nat, const int lmaxd, const int nmaxd, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& descriptor, double** gedm, double& E_delta, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_descriptor.cpp b/source/module_hamilt_lcao/module_deepks/deepks_descriptor.cpp index f29c1f0cf4..5bac498e6e 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_descriptor.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_descriptor.cpp @@ -37,7 +37,7 @@ void DeePKS_domain::cal_descriptor_equiv(const int nat, // calculates descriptors from projected density matrices void DeePKS_domain::cal_descriptor(const int nat, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const std::vector& pdm, std::vector& descriptor, const int des_per_atom = -1) @@ -53,7 +53,7 @@ void DeePKS_domain::cal_descriptor(const int nat, for (int inl = 0; inl < inlmax; ++inl) { - const int nm = 2 * inl_l[inl] + 1; + const int nm = 2 * inl2l[inl] + 1; pdm[inl].requires_grad_(true); descriptor.push_back(torch::ones({nm}, torch::requires_grad(true))); } @@ -74,7 +74,7 @@ void DeePKS_domain::cal_descriptor(const int nat, void DeePKS_domain::check_descriptor(const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const UnitCell& ucell, const std::string& out_dir, const std::vector& descriptor, @@ -104,7 +104,7 @@ void DeePKS_domain::check_descriptor(const int inlmax, int id = 0; for (int inl = 0; inl < inlmax / ucell.nat; inl++) { - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; const int ind = iat * inlmax / ucell.nat + inl; auto accessor = descriptor[ind].accessor(); for (int im = 0; im < nm; im++) diff --git a/source/module_hamilt_lcao/module_deepks/deepks_descriptor.h b/source/module_hamilt_lcao/module_deepks/deepks_descriptor.h index 3f365b9643..f4664c7f25 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_descriptor.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_descriptor.h @@ -30,14 +30,14 @@ namespace DeePKS_domain /// which are eigenvalues of pdm in blocks of I_n_l void cal_descriptor(const int nat, const int inlmax, - const int* inl_l, + const std::vector& inl2l, const std::vector& pdm, std::vector& descriptor, const int des_per_atom); /// print descriptors based on LCAO basis void check_descriptor(const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const UnitCell& ucell, const std::string& out_dir, const std::vector& descriptor, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_fpre.cpp b/source/module_hamilt_lcao/module_deepks/deepks_fpre.cpp index f6db839180..2616fb5d1c 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_fpre.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_fpre.cpp @@ -221,7 +221,7 @@ void DeePKS_domain::check_gdmx(const torch::Tensor& gdmx) void DeePKS_domain::cal_gvx(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& gevdm, const torch::Tensor& gdmx, torch::Tensor& gvx, @@ -249,7 +249,7 @@ void DeePKS_domain::cal_gvx(const int nat, for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; std::vector mmv; for (int m1 = 0; m1 < nm; ++m1) { diff --git a/source/module_hamilt_lcao/module_deepks/deepks_fpre.h b/source/module_hamilt_lcao/module_deepks/deepks_fpre.h index fa6db95b85..fcc10d13a9 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_fpre.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_fpre.h @@ -59,7 +59,7 @@ void check_gdmx(const torch::Tensor& gdmx); void cal_gvx(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& gevdm, const torch::Tensor& gdmx, torch::Tensor& gvx, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_hmat.cpp b/source/module_hamilt_lcao/module_deepks/deepks_hmat.cpp index 54d143687c..90936e88cc 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_hmat.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_hmat.cpp @@ -1,71 +1,70 @@ #ifdef __DEEPKS -#include "module_parameter/parameter.h" #include "LCAO_deepks.h" #include "module_base/parallel_reduce.h" +#include "module_parameter/parameter.h" template -void DeePKS_domain::collect_h_mat( - const Parallel_Orbitals &pv, - const std::vector>& h_in, - std::vector &h_out, - const int nlocal, - const int nks) +void DeePKS_domain::collect_h_mat(const Parallel_Orbitals& pv, + const std::vector>& h_in, + std::vector& h_out, + const int nlocal, + const int nks) { ModuleBase::TITLE("DeePKS_domain", "collect_h_tot"); - //construct the total H matrix - for (int k=0; k lineH(nlocal-i, TK(0.0)); + std::vector lineH(nlocal - i, TK(0.0)); ir = pv.global2local_row(i); - if (ir>=0) + if (ir >= 0) { // data collection - for (int j=i; j=0) + if (ic >= 0) { - int iic=0; + int iic = 0; if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) { - iic=ir+ic*pv.nrow; + iic = ir + ic * pv.nrow; } else { - iic=ir*pv.ncol+ic; + iic = ir * pv.ncol + ic; } - lineH[j-i] = h_in[k][iic]; + lineH[j - i] = h_in[k][iic]; } } } else { - //do nothing + // do nothing } - Parallel_Reduce::reduce_all(lineH.data(),nlocal-i); + Parallel_Reduce::reduce_all(lineH.data(), nlocal - i); - for (int j=i; j -void DeePKS_domain::check_h_mat( - const std::vector &H, - const std::string &h_file, - const int nlocal, - const int nks) +void DeePKS_domain::check_h_mat(const std::vector& H, const std::string& h_file, const int nlocal, const int nks) { std::ofstream ofs(h_file.c_str()); ofs << std::setprecision(10); - for (int k=0; k( - const Parallel_Orbitals &pv, - const std::vector>& h_in, - std::vector &h_out, - const int nlocal, - const int nks); +template void DeePKS_domain::collect_h_mat(const Parallel_Orbitals& pv, + const std::vector>& h_in, + std::vector& h_out, + const int nlocal, + const int nks); template void DeePKS_domain::collect_h_mat, ModuleBase::ComplexMatrix>( - const Parallel_Orbitals &pv, - const std::vector>>& h_in, - std::vector &h_out, - const int nlocal, - const int nks); + const Parallel_Orbitals& pv, + const std::vector>>& h_in, + std::vector& h_out, + const int nlocal, + const int nks); -template void DeePKS_domain::check_h_mat( - const std::vector &H, - const std::string &h_file, - const int nlocal, - const int nks); +template void DeePKS_domain::check_h_mat(const std::vector& H, + const std::string& h_file, + const int nlocal, + const int nks); -template void DeePKS_domain::check_h_mat( - const std::vector &H, - const std::string &h_file, - const int nlocal, - const int nks); +template void DeePKS_domain::check_h_mat(const std::vector& H, + const std::string& h_file, + const int nlocal, + const int nks); #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_hmat.h b/source/module_hamilt_lcao/module_deepks/deepks_hmat.h index 9668d5ede9..98d365e70a 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_hmat.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_hmat.h @@ -1,10 +1,10 @@ -#ifndef DEEPKS_HMAT_H -#define DEEPKS_HMAT_H +#ifndef DEEPKS_HMAT_H +#define DEEPKS_HMAT_H #ifdef __DEEPKS -#include "module_base/matrix.h" #include "module_base/complexmatrix.h" +#include "module_base/matrix.h" #include "module_base/timer.h" #include "module_basis/module_ao/parallel_orbitals.h" @@ -14,23 +14,18 @@ namespace DeePKS_domain { - //Collect data in h_in to matrix h_out. Note that left lower trianger in h_out is filled - template - void collect_h_mat( - const Parallel_Orbitals &pv, - const std::vector>& h_in, - std::vector &h_out, - const int nlocal, - const int nks); +// Collect data in h_in to matrix h_out. Note that left lower trianger in h_out is filled +template +void collect_h_mat(const Parallel_Orbitals& pv, + const std::vector>& h_in, + std::vector& h_out, + const int nlocal, + const int nks); - // write h_mat to file h_file for checking // not used in the code now - template - void check_h_mat( - const std::vector &H, - const std::string &h_file, - const int nlocal, - const int nks); -} +// write h_mat to file h_file for checking // not used in the code now +template +void check_h_mat(const std::vector& H, const std::string& h_file, const int nlocal, const int nks); +} // namespace DeePKS_domain #endif #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbpre.cpp b/source/module_hamilt_lcao/module_deepks/deepks_orbpre.cpp index 16e8ac7c54..d679177b5b 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_orbpre.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbpre.cpp @@ -22,7 +22,7 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, @@ -303,7 +303,7 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; std::vector mmv; for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d @@ -345,7 +345,7 @@ template void DeePKS_domain::cal_orbital_precalc( const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, @@ -362,7 +362,7 @@ template void DeePKS_domain::cal_orbital_precalc, ModuleBas const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbpre.h b/source/module_hamilt_lcao/module_deepks/deepks_orbpre.h index e044f9d54f..cf982ec879 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_orbpre.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbpre.h @@ -31,7 +31,7 @@ void cal_orbital_precalc(const std::vector& dm_hl, const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_pdm.cpp b/source/module_hamilt_lcao/module_deepks/deepks_pdm.cpp index 56eceb9045..fc7978a4db 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_pdm.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_pdm.cpp @@ -30,7 +30,7 @@ void DeePKS_domain::read_pdm(bool read_pdm_file, const int nat, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const Numerical_Orbital& alpha, std::vector& pdm) { @@ -47,7 +47,7 @@ void DeePKS_domain::read_pdm(bool read_pdm_file, { for (int inl = 0; inl < inlmax; inl++) { - int nm = inl_l[inl] * 2 + 1; + int nm = 2 * inl2l[inl] + 1; auto accessor = pdm[inl].accessor(); for (int m1 = 0; m1 < nm; m1++) { @@ -91,7 +91,7 @@ template void DeePKS_domain::cal_pdm(bool& init_pdm, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const ModuleBase::IntArray* inl_index, const elecstate::DensityMatrix* dm, const std::vector*> phialpha, @@ -116,7 +116,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, { for (int inl = 0; inl < inlmax; inl++) { - int nm = inl_l[inl] * 2 + 1; + int nm = 2 * inl2l[inl] + 1; pdm[inl] = torch::zeros({nm, nm}, torch::kFloat64); } } @@ -283,6 +283,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, dRy = dR2.y - dR1.y; dRz = dR2.z - dR1.z; } + // dm_R auto* tmp = dm->get_DMR_vector()[is]->find_matrix(ibt1, ibt2, dRx, dRy, dRz); if (tmp == nullptr) { @@ -325,7 +326,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, } // ad2 if (!PARAM.inp.deepks_equiv) { - int ib = 0, index = 0, inc = 1; + int index = 0, inc = 1; for (int L0 = 0; L0 <= orb.Alpha[0].getLmax(); ++L0) { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) @@ -346,7 +347,6 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, index++; } } - ib += nm; } } } @@ -380,7 +380,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, #ifdef __MPI for (int inl = 0; inl < inlmax; inl++) { - int pdm_size = (2 * inl_l[inl] + 1) * (2 * inl_l[inl] + 1); + int pdm_size = (2 * inl2l[inl] + 1) * (2 * inl2l[inl] + 1); Parallel_Reduce::reduce_all(pdm[inl].data_ptr(), pdm_size); } #endif @@ -388,7 +388,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, return; } -void DeePKS_domain::check_pdm(const int inlmax, const int* inl_l, const std::vector& pdm) +void DeePKS_domain::check_pdm(const int inlmax, const std::vector& inl2l, const std::vector& pdm) { const std::string file_projdm = PARAM.globalv.global_out_dir + "pdm.dat"; std::ofstream ofs(file_projdm.c_str()); @@ -396,7 +396,7 @@ void DeePKS_domain::check_pdm(const int inlmax, const int* inl_l, const std::vec ofs << std::setprecision(10); for (int inl = 0; inl < inlmax; inl++) { - const int nm = 2 * inl_l[inl] + 1; + const int nm = 2 * inl2l[inl] + 1; auto accessor = pdm[inl].accessor(); for (int m1 = 0; m1 < nm; m1++) { @@ -412,7 +412,7 @@ void DeePKS_domain::check_pdm(const int inlmax, const int* inl_l, const std::vec template void DeePKS_domain::cal_pdm(bool& init_pdm, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const ModuleBase::IntArray* inl_index, const elecstate::DensityMatrix* dm, const std::vector*> phialpha, @@ -426,7 +426,7 @@ template void DeePKS_domain::cal_pdm>( bool& init_pdm, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const ModuleBase::IntArray* inl_index, const elecstate::DensityMatrix, double>* dm, const std::vector*> phialpha, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_pdm.h b/source/module_hamilt_lcao/module_deepks/deepks_pdm.h index 44063a407d..ec9c4b37ac 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_pdm.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_pdm.h @@ -35,7 +35,7 @@ void read_pdm(bool read_pdm_file, const int nat, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const Numerical_Orbital& alpha, std::vector& pdm); @@ -48,7 +48,7 @@ template void cal_pdm(bool& init_pdm, const int inlmax, const int lmaxd, - const int* inl_l, + const std::vector& inl2l, const ModuleBase::IntArray* inl_index, const elecstate::DensityMatrix* dm, const std::vector*> phialpha, @@ -58,7 +58,7 @@ void cal_pdm(bool& init_pdm, const Parallel_Orbitals& pv, std::vector& pdm); -void check_pdm(const int inlmax, const int* inl_l, const std::vector& pdm); +void check_pdm(const int inlmax, const std::vector& inl2l, const std::vector& pdm); } // namespace DeePKS_domain #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_spre.cpp b/source/module_hamilt_lcao/module_deepks/deepks_spre.cpp index be81262fa1..dfc83cad1e 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_spre.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_spre.cpp @@ -216,7 +216,7 @@ void DeePKS_domain::check_gdmepsl(const torch::Tensor& gdmepsl) void DeePKS_domain::cal_gvepsl(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& gevdm, const torch::Tensor& gdmepsl, torch::Tensor& gvepsl, @@ -240,7 +240,7 @@ void DeePKS_domain::cal_gvepsl(const int nat, for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; std::vector mmv; for (int m1 = 0; m1 < nm; ++m1) { diff --git a/source/module_hamilt_lcao/module_deepks/deepks_spre.h b/source/module_hamilt_lcao/module_deepks/deepks_spre.h index 396c3d79cb..674c6bbf3b 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_spre.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_spre.h @@ -52,7 +52,7 @@ void check_gdmepsl(const torch::Tensor& gdmepsl); void cal_gvepsl(const int nat, const int inlmax, const int des_per_atom, - const int* inl_l, + const std::vector& inl2l, const std::vector& gevdm, const torch::Tensor& gdmepsl, torch::Tensor& gvepsl, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_vdpre.cpp b/source/module_hamilt_lcao/module_deepks/deepks_vdpre.cpp index 45358a225e..96da03df66 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_vdpre.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_vdpre.cpp @@ -27,7 +27,7 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, @@ -176,7 +176,7 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl_l[inl] + 1; + int nm = 2 * inl2l[inl] + 1; std::vector mmv; for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d @@ -505,7 +505,7 @@ template void DeePKS_domain::cal_v_delta_precalc(const int nlocal, const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, @@ -521,7 +521,7 @@ template void DeePKS_domain::cal_v_delta_precalc>( const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_vdpre.h b/source/module_hamilt_lcao/module_deepks/deepks_vdpre.h index 41abcf8b0c..00676e582e 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_vdpre.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_vdpre.h @@ -39,7 +39,7 @@ void cal_v_delta_precalc(const int nlocal, const int inlmax, const int nat, const int nks, - const int* inl_l, + const std::vector& inl2l, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, diff --git a/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt b/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt index 692ea81f6c..6dc611cbd3 100644 --- a/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_deepks/test/CMakeLists.txt @@ -21,6 +21,7 @@ add_executable( ../../../module_cell/read_pp_blps.cpp ../../../module_hamilt_pw/hamilt_pwdft/soc.cpp ../../../module_io/output.cpp + ../../../module_io/sparse_matrix.cpp ../../../module_elecstate/read_pseudo.cpp ../../../module_elecstate/cal_wfc.cpp ../../../module_elecstate/read_orb.cpp @@ -33,6 +34,7 @@ add_executable( ../../../module_hamilt_lcao/module_hcontainer/func_transfer.cpp ../../../module_hamilt_lcao/module_hcontainer/func_folding.cpp ../../../module_hamilt_lcao/module_hcontainer/transfer.cpp + ../../../module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp ../../../module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp ../../../module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp ../../../module_hamilt_general/operator.cpp diff --git a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp index 86247bc4eb..6bcb75d7bb 100644 --- a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp +++ b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp @@ -140,7 +140,7 @@ void test_deepks::check_pdm() DeePKS_domain::cal_pdm(this->ld.init_pdm, this->ld.inlmax, this->ld.lmaxd, - this->ld.inl_l, + this->ld.inl2l, this->ld.inl_index, p_elec_DM, this->ld.phialpha, @@ -149,7 +149,7 @@ void test_deepks::check_pdm() Test_Deepks::GridD, ParaO, this->ld.pdm); - DeePKS_domain::check_pdm(this->ld.inlmax, this->ld.inl_l, this->ld.pdm); + DeePKS_domain::check_pdm(this->ld.inlmax, this->ld.inl2l, this->ld.pdm); this->compare_with_ref("pdm.dat", "pdm_ref.dat"); } @@ -229,11 +229,11 @@ void test_deepks::check_descriptor(std::vector& descriptor) { DeePKS_domain::cal_descriptor(ucell.nat, this->ld.inlmax, - this->ld.inl_l, + this->ld.inl2l, this->ld.pdm, descriptor, this->ld.des_per_atom); - DeePKS_domain::check_descriptor(this->ld.inlmax, this->ld.des_per_atom, this->ld.inl_l, ucell, "./", descriptor, 0); + DeePKS_domain::check_descriptor(this->ld.inlmax, this->ld.des_per_atom, this->ld.inl2l, ucell, "./", descriptor, 0); this->compare_with_ref("deepks_desc.dat", "descriptor_ref.dat"); } @@ -241,9 +241,9 @@ template void test_deepks::check_gvx(torch::Tensor& gdmx) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl_l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, this->ld.inl_l, gevdm, gdmx, gvx, 0); + DeePKS_domain::cal_gvx(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, this->ld.inl2l, gevdm, gdmx, gvx, 0); DeePKS_domain::check_gvx(gvx, 0); for (int ia = 0; ia < ucell.nat; ia++) @@ -274,12 +274,12 @@ template void test_deepks::check_gvepsl(torch::Tensor& gdmepsl) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl_l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); torch::Tensor gvepsl; DeePKS_domain::cal_gvepsl(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, - this->ld.inl_l, + this->ld.inl2l, gevdm, gdmepsl, gvepsl, @@ -310,7 +310,7 @@ void test_deepks::check_edelta(std::vector& descriptor) this->ld.nmaxd, this->ld.inlmax, this->ld.des_per_atom, - this->ld.inl_l, + this->ld.inl2l, descriptor, this->ld.gedm, this->ld.E_delta, @@ -321,7 +321,7 @@ void test_deepks::check_edelta(std::vector& descriptor) DeePKS_domain::cal_edelta_gedm(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, - this->ld.inl_l, + this->ld.inl2l, descriptor, this->ld.pdm, this->ld.model_deepks, @@ -334,7 +334,7 @@ void test_deepks::check_edelta(std::vector& descriptor) ofs.close(); this->compare_with_ref("E_delta.dat", "E_delta_ref.dat"); - DeePKS_domain::check_gedm(this->ld.inlmax, this->ld.inl_l, this->ld.gedm); + DeePKS_domain::check_gedm(this->ld.inlmax, this->ld.inl2l, this->ld.gedm); this->compare_with_ref("gedm.dat", "gedm_ref.dat"); } diff --git a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp index 04fab988a7..c8bb58d95e 100644 --- a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp @@ -119,9 +119,7 @@ void Output_HContainer::write_single_R(int rx, int ry, int rz) this->_hcontainer->unfix_R(); } -// explicit instantiation of template class with double type template class Output_HContainer; -// to do: explicit instantiation of template class with std::complex type -// template class Output_HContainer>; +template class Output_HContainer>; } // namespace hamilt \ No newline at end of file