From ceb2c46f1b1ba7f8314f21fc3c5c0e570c39c41d Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 27 Feb 2025 22:03:17 +0800 Subject: [PATCH 1/4] fix the R range of data-HR --- .../hamilt_lcaodft/spar_hsr.cpp | 8 ++++++-- .../hamilt_lcaodft/spar_hsr.h | 17 +++++++++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index e2ab76647bf..9481ab3ecc4 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -23,7 +23,7 @@ void sparse_format::cal_HSR(const UnitCell& ucell, ) { ModuleBase::TITLE("sparse_format", "cal_HSR"); - sparse_format::set_R_range(HS_Arrays.all_R_coor, grid); + // sparse_format::set_R_range(HS_Arrays.all_R_coor, grid); const int nspin = PARAM.inp.nspin; @@ -33,6 +33,8 @@ void sparse_format::cal_HSR(const UnitCell& ucell, = dynamic_cast, double>*>( p_ham); + HS_Arrays.all_R_coor = get_R_range(*(p_ham_lcao->getHR())); + if (TD_Velocity::tddft_velocity) { sparse_format::cal_HContainer_td( pv, @@ -57,7 +59,9 @@ void sparse_format::cal_HSR(const UnitCell& ucell, hamilt::HamiltLCAO, std::complex>* p_ham_lcao = dynamic_cast, - std::complex>*>(p_ham); + std::complex>*>(p_ham); + + HS_Arrays.all_R_coor = get_R_range(*(p_ham_lcao->getHR())); sparse_format::cal_HContainer_cd(pv, current_spin, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index bf803866788..fd62200a416 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -4,6 +4,23 @@ #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" namespace sparse_format { + template + std::set> get_R_range(const hamilt::HContainer& hR) + { + std::set> all_R_coor; + for (int iap = 0;iap < hR.size_atom_pairs(); ++iap) + { + const hamilt::AtomPair& atom_pair = hR.get_atom_pair(iap); + for (int iR = 0; iR < atom_pair.get_R_size(); ++iR) + { + const auto& r_index = atom_pair.get_R_index(iR); + Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); + all_R_coor.insert(dR); + } + } + return all_R_coor; + }; + using TAC = std::pair>; void cal_HSR(const UnitCell& ucell, const Parallel_Orbitals& pv, From daab36c32893fb3fc1aede3953ae7d333c8cd73e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 26 Mar 2025 17:17:39 +0800 Subject: [PATCH 2/4] feature: write Vxc(R) --- source/module_esolver/esolver_ks_lcao.cpp | 21 ++ .../operator_lcao/op_exx_lcao.cpp | 16 ++ .../operator_lcao/op_exx_lcao.h | 16 ++ .../operator_lcao/op_exx_lcao.hpp | 29 ++- source/module_io/write_vxc.hpp | 15 +- source/module_io/write_vxc_r.hpp | 202 ++++++++++++++++++ 6 files changed, 271 insertions(+), 28 deletions(-) create mode 100644 source/module_io/write_vxc_r.hpp diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 9549cc511e9..9c5326e0201 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -33,6 +33,7 @@ //be careful of hpp, there may be multiple definitions of functions, 20250302, mohan #include "module_io/write_eband_terms.hpp" #include "module_io/write_vxc.hpp" +#include "module_io/write_vxc_r.hpp" #include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" //--------------temporary---------------------------- @@ -509,6 +510,26 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr, this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr +#endif + ); + + ModuleIO::write_Vxc_R(PARAM.inp.nspin, + &this->pv, + ucell, + this->sf, + this->solvent, + *this->pw_rho, + *this->pw_rhod, + this->locpp.vloc, + this->chr, + this->GG, + this->GK, + this->kv, + orb_.cutoffs(), + this->gd +#ifdef __EXX + , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr, + this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr #endif ); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp index 56f9520b700..4e7b2b553fa 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp @@ -4,6 +4,22 @@ namespace hamilt { + RI::Cell_Nearest init_cell_nearest(const UnitCell& ucell, const std::array& Rs_period) + { + RI::Cell_Nearest cell_nearest; + std::map> atoms_pos; + for (int iat = 0; iat < ucell.nat; ++iat) { + atoms_pos[iat] = RI_Util::Vector3_to_array3( + ucell.atoms[ucell.iat2it[iat]] + .tau[ucell.iat2ia[iat]]); + } + const std::array, 3> latvec + = { RI_Util::Vector3_to_array3(ucell.a1), + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3) }; + cell_nearest.init(atoms_pos, latvec, Rs_period); + return cell_nearest; + } template<> void OperatorEXX>::add_loaded_Hexx(const int ik) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h index cc68a11b06d..b112a552cfe 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h @@ -73,6 +73,22 @@ class OperatorEXX> : public OperatorLCAO std::vector>> Hexxc_k_load; }; +using TAC = std::pair>; + +RI::Cell_Nearest init_cell_nearest(const UnitCell& ucell, const std::array& Rs_period); + +// allocate according to the read-in HexxR, used in nscf +template +void reallocate_hcontainer(const std::vector>>>& Hexxs, + HContainer* hR, + const RI::Cell_Nearest* const cell_nearest = nullptr); + +/// allocate according to BvK cells, used in scf +template +void reallocate_hcontainer(const int nat, HContainer* hR, + const std::array& Rs_period, + const RI::Cell_Nearest* const cell_nearest = nullptr); + } // namespace hamilt #endif // __EXX #include "op_exx_lcao.hpp" diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp index 5c4af43f406..8863b8b19c7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp @@ -12,10 +12,12 @@ namespace hamilt { using TAC = std::pair>; + // allocate according to the read-in HexxR, used in nscf template - inline void reallocate_hcontainer(const std::vector>>>& Hexxs, - HContainer* hR) + void reallocate_hcontainer(const std::vector>>>& Hexxs, + HContainer* hR, + const RI::Cell_Nearest* const cell_nearest) { auto* pv = hR->get_paraV(); bool need_allocate = false; @@ -27,7 +29,10 @@ namespace hamilt const int& iat1 = Htmp2.first.first; if (pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0) { - const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3(Htmp2.first.second); + const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3( + (cell_nearest ? + cell_nearest->get_cell_nearest_discrete(iat0, iat1, Htmp2.first.second) + : Htmp2.first.second)); BaseMatrix* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z); if (HlocR == nullptr) { // add R to HContainer @@ -40,11 +45,12 @@ namespace hamilt } if (need_allocate) { hR->allocate(nullptr, true); } } + /// allocate according to BvK cells, used in scf template - inline void reallocate_hcontainer(const int nat, HContainer* hR, + void reallocate_hcontainer(const int nat, HContainer* hR, const std::array& Rs_period, - const RI::Cell_Nearest* const cell_nearest = nullptr) + const RI::Cell_Nearest* const cell_nearest) { auto* pv = hR->get_paraV(); auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period); @@ -154,18 +160,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, const std::array Rs_period = { this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2] }; if (this->use_cell_nearest) { - // set cell_nearest - std::map> atoms_pos; - for (int iat = 0; iat < ucell.nat; ++iat) { - atoms_pos[iat] = RI_Util::Vector3_to_array3( - ucell.atoms[ucell.iat2it[iat]] - .tau[ucell.iat2ia[iat]]); - } - const std::array, 3> latvec - = { RI_Util::Vector3_to_array3(ucell.a1), - RI_Util::Vector3_to_array3(ucell.a2), - RI_Util::Vector3_to_array3(ucell.a3) }; - this->cell_nearest.init(atoms_pos, latvec, Rs_period); + this->cell_nearest = init_cell_nearest(ucell, Rs_period); reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest); } else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); } diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp index fa9da0ec629..201e8cbaae1 100644 --- a/source/module_io/write_vxc.hpp +++ b/source/module_io/write_vxc.hpp @@ -30,16 +30,6 @@ struct TGint> namespace ModuleIO { -inline void gint_vl(Gint_Gamma& gg, Gint_inout& io) -{ - gg.cal_vlocal(&io, false); -}; - -inline void gint_vl(Gint_k& gk, Gint_inout& io, ModuleBase::matrix& wg) -{ - gk.cal_gint(&io); -}; - inline void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d) { std::ofstream ofs; @@ -134,6 +124,8 @@ std::vector orbital_energy(const int ik, const int nbands, const std::ve return e; } +#ifndef SET_GINT_POINTER_H +#define SET_GINT_POINTER_H // mohan update 2024-04-01 template void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint::type*& gint); @@ -153,6 +145,7 @@ void set_gint_pointer>(Gint_Gamma& gint_gamma, { gint = &gint_k; } +#endif inline void write_orb_energy(const K_Vectors& kv, const int nspin0, const int nbands, @@ -223,7 +216,7 @@ void write_Vxc(const int nspin, // 2. allocate AO-matrix // R (the number of hR: 1 for nspin=1, 4; 2 for nspin=2) int nspin0 = (nspin == 2) ? 2 : 1; - std::vector> vxcs_R_ao(nspin0, hamilt::HContainer(pv)); + std::vector> vxcs_R_ao(nspin0, hamilt::HContainer(ucell, pv)); for (int is = 0; is < nspin0; ++is) { vxcs_R_ao[is].set_zero(); if (std::is_same::value) { vxcs_R_ao[is].fix_gamma(); } diff --git a/source/module_io/write_vxc_r.hpp b/source/module_io/write_vxc_r.hpp new file mode 100644 index 00000000000..ee180da40a1 --- /dev/null +++ b/source/module_io/write_vxc_r.hpp @@ -0,0 +1,202 @@ +#ifndef __WRITE_VXC_R_H_ +#define __WRITE_VXC_R_H_ +#include "module_parameter/parameter.h" +#include "module_base/scalapack_connector.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h" +#ifdef __EXX +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" +#include "module_ri/RI_2D_Comm.h" +#endif + +#ifndef TGINT_H +#define TGINT_H +template +struct TGint; + +template <> +struct TGint +{ + using type = Gint_Gamma; +}; + +template <> +struct TGint> +{ + using type = Gint_k; +}; +#endif + +namespace ModuleIO +{ + +#ifndef SET_GINT_POINTER_H +#define SET_GINT_POINTER_H +template +void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint::type*& gint); + +template <> +void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint::type*& gint) +{ + gint = &gint_gamma; +} + +template <> +void set_gint_pointer>(Gint_Gamma& gint_gamma, + Gint_k& gint_k, + typename TGint>::type*& gint) +{ + gint = &gint_k; +} +#endif + +template std::set> get_R_range(const hamilt::HContainer& hR) +{ + std::set> all_R_coor; + + return all_R_coor; +} + +template +std::map, std::map>> +cal_HR_sparse(const hamilt::HContainer& hR, + const int current_spin, + const double sparse_thr); + +template <> +std::map, std::map>> +cal_HR_sparse(const hamilt::HContainer& hR, + const int current_spin, + const double sparse_thr) +{ + std::map, std::map>> target; + sparse_format::cal_HContainer_d(*hR.get_paraV(), current_spin, sparse_thr, hR, target); + return target; +} +template <> +std::map, std::map>>> +cal_HR_sparse(const hamilt::HContainer>& hR, + const int current_spin, + const double sparse_thr) +{ + std::map, std::map>>> target; + sparse_format::cal_HContainer_cd(*hR.get_paraV(), current_spin, sparse_thr, hR, target); + return target; +} + +/// @brief write the Vxc matrix in KS orbital representation, usefull for GW calculation +/// including terms: local/semi-local XC, EXX, DFTU +template +void write_Vxc_R(const int nspin, + const Parallel_Orbitals* pv, + const UnitCell& ucell, + Structure_Factor& sf, + surchem& solvent, + const ModulePW::PW_Basis& rho_basis, + const ModulePW::PW_Basis& rhod_basis, + const ModuleBase::matrix& vloc, + const Charge& chg, + Gint_Gamma& gint_gamma, + Gint_k& gint_k, + const K_Vectors& kv, + const std::vector& orb_cutoff, + Grid_Driver& gd, +#ifdef __EXX + const std::vector>>>* const Hexxd, + const std::vector>>>>* const Hexxc, +#endif +const double sparse_thr=1e-10) +{ + ModuleBase::TITLE("ModuleIO", "write_Vxc_R"); + // 1. real-space xc potential + // ModuleBase::matrix vr_xc(nspin, chg.nrxx); + double etxc = 0.0; + double vtxc = 0.0; + // elecstate::PotXC* potxc(&rho_basis, &etxc, vtxc, nullptr); + // potxc.cal_v_eff(&chg, &ucell, vr_xc); + elecstate::Potential* potxc + = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc); + std::vector compnents_list = {"xc"}; + potxc->pot_register(compnents_list); + potxc->update_from_charge(&chg, &ucell); + + // 2. allocate H(R) + // (the number of hR: 1 for nspin=1, 4; 2 for nspin=2) + int nspin0 = (nspin == 2) ? 2 : 1; + std::vector> vxcs_R_ao(nspin0, hamilt::HContainer(ucell, pv)); // call move constructor +#ifdef __EXX + std::array Rs_period = { kv.nmp[0], kv.nmp[1], kv.nmp[2] }; + const auto cell_nearest = hamilt::init_cell_nearest(ucell, Rs_period); +#endif + for (int is = 0; is < nspin0; ++is) + { + if (std::is_same::value) { vxcs_R_ao[is].fix_gamma(); } +#ifdef __EXX + if (GlobalC::exx_info.info_global.cal_exx) + { + GlobalC::exx_info.info_ri.real_number ? + hamilt::reallocate_hcontainer(*Hexxd, &vxcs_R_ao[is], &cell_nearest) : + hamilt::reallocate_hcontainer(*Hexxc, &vxcs_R_ao[is], &cell_nearest); + } +#endif + } + + // 3. calculate the Vxc(R) + hamilt::HS_Matrix_K vxc_k_ao(pv, 1); // only hk is needed, sk is skipped + typename TGint::type* gint = nullptr; + set_gint_pointer(gint_gamma, gint_k, gint); + std::vector>*> vxcs_op_ao(nspin0); + for (int is = 0; is < nspin0; ++is) + { + vxcs_op_ao[is] = new hamilt::Veff>(gint, + &vxc_k_ao, kv.kvec_d, potxc, &vxcs_R_ao[is], &ucell, orb_cutoff, &gd, nspin); + vxcs_op_ao[is]->contributeHR(); +#ifdef __EXX + if (GlobalC::exx_info.info_global.cal_exx) + { + GlobalC::exx_info.info_ri.real_number ? + RI_2D_Comm::add_HexxR(is, GlobalC::exx_info.info_global.hybrid_alpha, *Hexxd, *pv, ucell.get_npol(), vxcs_R_ao[is], &cell_nearest) : + RI_2D_Comm::add_HexxR(is, GlobalC::exx_info.info_global.hybrid_alpha, *Hexxc, *pv, ucell.get_npol(), vxcs_R_ao[is], &cell_nearest); + } +#endif + } + + // test: fold Vxc(R) and check whether it is equal to Vxc(k) + // for (int ik = 0; ik < kv.get_nks(); ++ik) + // { + // vxc_k_ao.set_zero_hk(); + // dynamic_cast*>(vxcs_op_ao[kv.isk[ik]])->contributeHk(ik); + + // // output Vxc(k) (test) + // const TK* const hk = vxc_k_ao.get_hk(); + // std::cout << "ik=" << ik << ", Vxc(K): " << std::endl; + // for (int i = 0; i < pv->get_row_size(); i++) + // { + // for (int j = 0; j < pv->get_col_size(); j++) + // { + // std::cout << hk[j * pv->get_row_size() + i] << " "; + // } + // std::cout << std::endl; + // } + // } + + // 4. write Vxc(R) in csr format + for (int is = 0; is < nspin0; ++is) + { + std::set> all_R_coor = sparse_format::get_R_range(vxcs_R_ao[is]); + const std::string filename = "Vxc_R_spin" + std::to_string(is); + ModuleIO::save_sparse( + cal_HR_sparse(vxcs_R_ao[is], is, sparse_thr), + all_R_coor, + sparse_thr, + false, //binary + PARAM.globalv.global_out_dir + filename + ".csr", + *pv, + filename, + -1, + true); //all-reduce + } +} +} // namespace ModuleIO +#endif From 988db114c85eac27084eb8dfbfb6565412bfa53b Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 26 Mar 2025 18:02:51 +0800 Subject: [PATCH 3/4] add parameter out_mat_xc2 and update tests --- docs/advanced/input_files/input-main.md | 8 ++++ source/module_esolver/esolver_ks_lcao.cpp | 44 ++++++++++--------- source/module_io/read_input_item_output.cpp | 6 +++ source/module_io/test/read_input_ptest.cpp | 1 + source/module_parameter/input_parameter.h | 1 + tests/integrate/307_NO_GO_OXC/INPUT | 1 + tests/integrate/307_NO_GO_OXC/Vxc_R_spin0.ref | 7 +++ tests/integrate/307_NO_GO_OXC/result.ref | 1 + tests/integrate/tools/catch_properties.sh | 8 ++++ 9 files changed, 56 insertions(+), 21 deletions(-) create mode 100644 tests/integrate/307_NO_GO_OXC/Vxc_R_spin0.ref diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 2a5df1f162e..f8eec32d73d 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -153,6 +153,7 @@ - [out\_mat\_t](#out_mat_t) - [out\_mat\_dh](#out_mat_dh) - [out\_mat\_xc](#out_mat_xc) + - [out\_mat\_xc2](#out_mat_xc2) - [out\_eband\_terms](#out_eband_terms) - [out\_hr\_npz/out\_dm\_npz](#out_hr_npzout_dm_npz) - [dm\_to\_rho](#dm_to_rho) @@ -1799,6 +1800,13 @@ These variables are used to control the output of properties. The band (KS orbital) energy for each (k-point, spin, band) will be printed in the file `OUT.${suffix}/vxc_out.dat`. If EXX is calculated, the local and EXX part of band energy will also be printed in `OUT.${suffix}/vxc_local_out.dat`and `OUT.${suffix}/vxc_exx_out.dat`, respectively. All the `vxc*_out.dat` files contains 3 integers (nk, nspin, nband) followed by nk\*nspin\*nband lines of energy Hartree and eV. - **Default**: False +### out_mat_xc2 + +- **Type**: Boolean +- **Availability**: Numerical atomic orbital (NAO) basis +- **Description**: Whether to print the exchange-correlation matrices in **numerical orbital representation** (unit: Ry): $\braket{\phi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\phi_j}(\mathbf{R})$ in CSR format (the same format as [out_mat_hs2](../elec_properties/hs_matrix.md#out_mat_hs2)) in the directory `OUT.${suffix}`. (Note that currently DeePKS term is not included. ) The files are named `Vxc_R_spin$s`. +- **Default**: False + ### out_eband_terms - **Type**: Boolean diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 9c5326e0201..c2a90f83aeb 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -489,30 +489,32 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) if (PARAM.inp.out_mat_xc) { ModuleIO::write_Vxc(PARAM.inp.nspin, - PARAM.globalv.nlocal, - GlobalV::DRANK, - &this->pv, - *this->psi, - ucell, - this->sf, - this->solvent, - *this->pw_rho, - *this->pw_rhod, - this->locpp.vloc, - this->chr, - this->GG, - this->GK, - this->kv, - orb_.cutoffs(), - this->pelec->wg, - this->gd + PARAM.globalv.nlocal, + GlobalV::DRANK, + &this->pv, + *this->psi, + ucell, + this->sf, + this->solvent, + *this->pw_rho, + *this->pw_rhod, + this->locpp.vloc, + this->chr, + this->GG, + this->GK, + this->kv, + orb_.cutoffs(), + this->pelec->wg, + this->gd #ifdef __EXX - , - this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr, - this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr + , + this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr, + this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr #endif ); - + } + if (PARAM.inp.out_mat_xc2) + { ModuleIO::write_Vxc_R(PARAM.inp.nspin, &this->pv, ucell, diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index e2c9dfb0685..0614ca1d08e 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -298,6 +298,12 @@ void ReadInput::item_output() read_sync_bool(input.out_mat_xc); this->add_item(item); } + { + Input_Item item("out_mat_xc2"); + item.annotation = "output exchange-correlation matrix in NAO representation"; + read_sync_bool(input.out_mat_xc2); + this->add_item(item); + } { Input_Item item("out_eband_terms"); item.annotation = "output the band energy terms separately"; diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 90f985c307f..d551f4f5fc5 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -208,6 +208,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.out_mat_hs[1], 8); EXPECT_EQ(param.inp.out_mat_hs2, 0); EXPECT_FALSE(param.inp.out_mat_xc); + EXPECT_FALSE(param.inp.out_mat_xc2); EXPECT_FALSE(param.inp.out_eband_terms); EXPECT_EQ(param.inp.out_interval, 1); EXPECT_EQ(param.inp.out_app_flag, 0); diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 4e0e8537e41..77ee42a2ed2 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -376,6 +376,7 @@ struct Input_para bool out_mat_dh = false; bool out_mat_xc = false; ///< output exchange-correlation matrix in ///< KS-orbital representation. + bool out_mat_xc2 = false; ///< output exchange-correlation matrix Vxc(R) in NAO representation. bool out_eband_terms = false; ///< output the band energy terms separately bool out_hr_npz = false; ///< output exchange-correlation matrix in ///< KS-orbital representation. diff --git a/tests/integrate/307_NO_GO_OXC/INPUT b/tests/integrate/307_NO_GO_OXC/INPUT index be4873f986f..cf821c8d81b 100644 --- a/tests/integrate/307_NO_GO_OXC/INPUT +++ b/tests/integrate/307_NO_GO_OXC/INPUT @@ -23,6 +23,7 @@ smearing_method gauss smearing_sigma 0.002 out_mat_xc 1 +out_mat_xc2 1 out_eband_terms 1 #Parameters (5.Mixing) diff --git a/tests/integrate/307_NO_GO_OXC/Vxc_R_spin0.ref b/tests/integrate/307_NO_GO_OXC/Vxc_R_spin0.ref new file mode 100644 index 00000000000..68418006680 --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/Vxc_R_spin0.ref @@ -0,0 +1,7 @@ +STEP: 0 +Matrix Dimension of Vxc_R_spin0(R): 26 +Matrix number of Vxc_R_spin0(R): 1 +0 0 0 132 + -6.99425500e-01 1.54630627e-01 -2.02590600e-02 -9.98372482e-02 1.54630627e-01 -3.63859528e-01 -9.98178110e-02 -1.36032667e-01 -6.10221093e-01 -2.05221537e-01 4.84073458e-03 3.86425460e-02 -6.54685745e-02 1.71711614e-02 -6.10221093e-01 -2.05221537e-01 4.84073458e-03 3.86425460e-02 -6.54685745e-02 1.71711614e-02 -6.10221093e-01 -2.05221537e-01 4.84073458e-03 3.86425460e-02 -6.54685745e-02 1.71711614e-02 -2.05221537e-01 -4.39113944e-01 -1.04968231e-03 -6.55024474e-02 2.36444413e-02 -6.86923692e-02 -2.05221537e-01 -4.39113944e-01 -1.04968231e-03 -6.55024474e-02 2.36444413e-02 -6.86923692e-02 -2.05221537e-01 -4.39113944e-01 -1.04968231e-03 -6.55024474e-02 2.36444413e-02 -6.86923692e-02 -6.98160574e-01 -4.48366310e-03 4.84073458e-03 -1.04968231e-03 -6.99822165e-01 -1.71459508e-02 6.86343042e-02 -3.03794506e-04 4.84073458e-03 -1.04968231e-03 -6.99822165e-01 -1.71459508e-02 6.86343042e-02 -3.03794506e-04 -6.98160574e-01 -4.48366310e-03 4.84073458e-03 -1.04968231e-03 -6.99822165e-01 -1.71459508e-02 6.86343042e-02 -3.03794506e-04 -2.02590600e-02 -9.98178110e-02 -6.99426855e-01 1.54629055e-01 -9.98372482e-02 -1.36032667e-01 1.54629055e-01 -3.63863991e-01 3.86425460e-02 -6.55024474e-02 -1.71459508e-02 -6.10223512e-01 -2.05219163e-01 -4.84033739e-03 3.86425460e-02 -6.55024474e-02 -1.71459508e-02 -6.10223512e-01 -2.05219163e-01 -4.84033739e-03 3.86425460e-02 -6.55024474e-02 -1.71459508e-02 -6.10223512e-01 -2.05219163e-01 -4.84033739e-03 -6.54685745e-02 2.36444413e-02 6.86343042e-02 -2.05219163e-01 -4.39128101e-01 1.04971294e-03 -6.54685745e-02 2.36444413e-02 6.86343042e-02 -2.05219163e-01 -4.39128101e-01 1.04971294e-03 -6.54685745e-02 2.36444413e-02 6.86343042e-02 -2.05219163e-01 -4.39128101e-01 1.04971294e-03 -4.48366310e-03 -6.98168157e-01 1.71711614e-02 -6.86923692e-02 -3.03794506e-04 -4.84033739e-03 1.04971294e-03 -6.99823286e-01 1.71711614e-02 -6.86923692e-02 -3.03794506e-04 -4.84033739e-03 1.04971294e-03 -6.99823286e-01 -4.48366310e-03 -6.98168157e-01 1.71711614e-02 -6.86923692e-02 -3.03794506e-04 -4.84033739e-03 1.04971294e-03 -6.99823286e-01 + 0 1 13 14 0 1 13 14 2 5 12 15 18 25 3 6 10 16 19 23 4 7 9 17 20 22 2 5 12 15 18 25 3 6 10 16 19 23 4 7 9 17 20 22 8 21 4 7 9 17 20 22 3 6 10 16 19 23 11 24 2 5 12 15 18 25 0 1 13 14 0 1 13 14 2 5 12 15 18 25 3 6 10 16 19 23 4 7 9 17 20 22 2 5 12 15 18 25 3 6 10 16 19 23 4 7 9 17 20 22 8 21 4 7 9 17 20 22 3 6 10 16 19 23 11 24 2 5 12 15 18 25 + 0 4 8 14 20 26 32 38 44 46 52 58 60 66 70 74 80 86 92 98 104 110 112 118 124 126 132 diff --git a/tests/integrate/307_NO_GO_OXC/result.ref b/tests/integrate/307_NO_GO_OXC/result.ref index 80390d1cada..aeec3dc93ec 100644 --- a/tests/integrate/307_NO_GO_OXC/result.ref +++ b/tests/integrate/307_NO_GO_OXC/result.ref @@ -1,6 +1,7 @@ etotref -204.1062806122546 etotperatomref -102.0531403061 CompareVXC_pass 0 +CompareVXC_R_pass 0 CompareOrbXC_pass 0 CompareOrbKinetic_pass 0 CompareOrbVL_pass 0 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index b4447dc163c..c72afea9b34 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -42,6 +42,7 @@ has_cond=$(get_input_key_value "cal_cond" "INPUT") has_hs=$(get_input_key_value "out_mat_hs" "INPUT") has_hs2=$(get_input_key_value "out_mat_hs2" "INPUT") has_xc=$(get_input_key_value "out_mat_xc" "INPUT") +has_xc2=$(get_input_key_value "out_mat_xc2" "INPUT") has_eband_separate=$(get_input_key_value "out_eband_terms" "INPUT") has_r=$(get_input_key_value "out_mat_r" "INPUT") deepks_out_labels=$(get_input_key_value "deepks_out_labels" "INPUT") @@ -257,6 +258,13 @@ if ! test -z "$has_xc" && [ $has_xc == 1 ]; then echo "CompareOrbXC_pass $?" >>$1 fi +if ! test -z "$has_xc2" && [ $has_xc2 == 1 ]; then + xc2ref=Vxc_R_spin0.ref + xc2cal=OUT.autotest/Vxc_R_spin0.csr + python3 ../tools/CompareFile.py $xc2ref $xc2cal 8 + echo "CompareVXC_R_pass $?" >>$1 +fi + if ! test -z "$has_eband_separate" && [ $has_eband_separate == 1 ]; then ekref=kinetic_out.ref ekcal=OUT.autotest/kinetic_out.dat From 62db5442197b3504d3e9b5b0b62d60bd05b05a77 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 26 Mar 2025 18:26:47 +0800 Subject: [PATCH 4/4] fix compile --- source/module_io/write_vxc_r.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/source/module_io/write_vxc_r.hpp b/source/module_io/write_vxc_r.hpp index ee180da40a1..29866484151 100644 --- a/source/module_io/write_vxc_r.hpp +++ b/source/module_io/write_vxc_r.hpp @@ -5,6 +5,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h" +#include "module_io/write_HS_sparse.h" #ifdef __EXX #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_ri/RI_2D_Comm.h"