diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 885a61a6df..682eb2c3d6 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -478,7 +478,6 @@ void ESolver_KS_LCAO::after_all_runners() this->GG, this->GK, this->LM, - this->LOC, this->kv, this->pelec->wg, GlobalC::GridD); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 69d1844bd9..5f640d8535 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -175,7 +175,6 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // effective potential term Operator* veff = new Veff>(GG_in, - loc_in, LM_in, this->kv->kvec_d, pot_in, @@ -248,7 +247,6 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // Veff term this->getOperator() = new Veff>(GK_in, - loc_in, LM_in, kv->kvec_d, pot_in, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h index 2587fb54cc..f8c3aaa49d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h @@ -1,13 +1,12 @@ #ifndef VEFFLCAO_H #define VEFFLCAO_H #include "module_base/timer.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" #include "module_elecstate/potentials/potential_new.h" -#include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "operator_lcao.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_cell/unitcell.h" namespace hamilt { @@ -25,8 +24,8 @@ class Veff : public T /// @brief Effective potential class, used for calculating Hamiltonian with grid integration tools /// If user want to separate the contribution of V_{eff} into V_{H} and V_{XC} and V_{local pseudopotential} and so on, /// the user can separate the Potential class into different parts, and construct different Veff class for each part. -/// @tparam TK -/// @tparam TR +/// @tparam TK +/// @tparam TR template class Veff> : public OperatorLCAO { @@ -34,23 +33,17 @@ class Veff> : public OperatorLCAO /** * @brief Construct a new Veff object for multi-kpoint calculation * @param GK_in: the pointer of Gint_k object, used for grid integration - */ + */ Veff>(Gint_k* GK_in, - Local_Orbital_Charge* loc_in, - LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - elecstate::Potential* pot_in, - hamilt::HContainer* hR_in, - std::vector* hK_in, - const UnitCell* ucell_in, - Grid_Driver* GridD_in, - const Parallel_Orbitals* paraV) - : GK(GK_in), - loc(loc_in), - pot(pot_in), - ucell(ucell_in), - gd(GridD_in), - OperatorLCAO(LM_in, kvec_d_in, hR_in, hK_in) + LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + std::vector* hK_in, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV) + : GK(GK_in), pot(pot_in), ucell(ucell_in), gd(GridD_in), OperatorLCAO(LM_in, kvec_d_in, hR_in, hK_in) { this->cal_type = calculation_type::lcao_gint; @@ -60,20 +53,17 @@ class Veff> : public OperatorLCAO /** * @brief Construct a new Veff object for Gamma-only calculation * @param GG_in: the pointer of Gint_Gamma object, used for grid integration - */ + */ Veff>(Gint_Gamma* GG_in, - Local_Orbital_Charge* loc_in, - LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - elecstate::Potential* pot_in, - hamilt::HContainer* hR_in, - std::vector* hK_in, - const UnitCell* ucell_in, - Grid_Driver* GridD_in, - const Parallel_Orbitals* paraV - ) - : GG(GG_in), loc(loc_in), pot(pot_in), - OperatorLCAO(LM_in, kvec_d_in, hR_in, hK_in) + LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + std::vector* hK_in, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV) + : GG(GG_in), pot(pot_in), OperatorLCAO(LM_in, kvec_d_in, hR_in, hK_in) { this->cal_type = calculation_type::lcao_gint; this->initialize_HR(ucell_in, GridD_in, paraV); @@ -86,13 +76,14 @@ class Veff> : public OperatorLCAO /** * @brief contributeHR() is used to calculate the HR matrix * - * the contribution of V_{eff} is calculated by the contribution of V_{H} and V_{XC} and V_{local pseudopotential} and so on. - * grid integration is used to calculate the contribution Hamiltonian of effective potential + * the contribution of V_{eff} is calculated by the contribution of V_{H} and V_{XC} and V_{local pseudopotential} + * and so on. grid integration is used to calculate the contribution Hamiltonian of effective potential */ virtual void contributeHR() override; - - const UnitCell* ucell; - Grid_Driver* gd; + + const UnitCell* ucell; + Grid_Driver* gd; + private: // used for k-dependent grid integration. Gint_k* GK = nullptr; @@ -101,7 +92,6 @@ class Veff> : public OperatorLCAO Gint_Gamma* GG = nullptr; // Charge calculating method in LCAO base and contained grid base calculation: DM_R, DM, pvpR_reduced - Local_Orbital_Charge* loc = nullptr; elecstate::Potential* pot = nullptr; @@ -114,7 +104,6 @@ class Veff> : public OperatorLCAO * the size of HR will be fixed after initialization */ void initialize_HR(const UnitCell* ucell_in, Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); - }; } // namespace hamilt diff --git a/source/module_io/write_Vxc.hpp b/source/module_io/write_Vxc.hpp index d8b395e773..0992d3fca9 100644 --- a/source/module_io/write_Vxc.hpp +++ b/source/module_io/write_Vxc.hpp @@ -1,21 +1,24 @@ #pragma once -#include "module_psi/psi.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h" -#include "module_base/scalapack_connector.h" #include "module_base/parallel_reduce.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_psi/psi.h" #ifndef TGINT_H #define TGINT_H -template struct TGint; +template +struct TGint; template <> -struct TGint { +struct TGint +{ using type = Gint_Gamma; }; template <> -struct TGint> { +struct TGint> +{ using type = Gint_k; }; #endif @@ -23,340 +26,415 @@ struct TGint> { namespace ModuleIO { - inline void gint_vl( - Gint_Gamma& gg, - Gint_inout& io, - LCAO_Matrix& lm) - { - gg.cal_vlocal(&io, false); - }; +inline void gint_vl(Gint_Gamma& gg, Gint_inout& io, LCAO_Matrix& lm) +{ + gg.cal_vlocal(&io, false); +}; - inline void gint_vl( - Gint_k& gk, - Gint_inout& io, - LCAO_Matrix& lm, - ModuleBase::matrix& wg) - { - gk.cal_gint(&io); - }; +inline void gint_vl(Gint_k& gk, Gint_inout& io, LCAO_Matrix& lm, ModuleBase::matrix& wg) +{ + gk.cal_gint(&io); +}; - void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d) - { - std::ofstream ofs; +void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d) +{ + std::ofstream ofs; #ifdef __MPI - p2d.set(nbands, nbands, pv.nb, pv.comm_2D, pv.blacs_ctxt); + p2d.set(nbands, nbands, pv.nb, pv.comm_2D, pv.blacs_ctxt); #else - p2d.set_serial(nbands, nbands); + p2d.set_serial(nbands, nbands); #endif - } +} - std::vector> cVc( - std::complex* V, - std::complex* c, - int nbasis, - int nbands, - const Parallel_Orbitals& pv, - const Parallel_2D& p2d) - { - std::vector> Vc(pv.nloc_wfc, 0.0); - char transa = 'N'; - char transb = 'N'; - const std::complex alpha(1.0, 0.0); - const std::complex beta(0.0, 0.0); +std::vector> cVc(std::complex* V, + std::complex* c, + int nbasis, + int nbands, + const Parallel_Orbitals& pv, + const Parallel_2D& p2d) +{ + std::vector> Vc(pv.nloc_wfc, 0.0); + char transa = 'N'; + char transb = 'N'; + const std::complex alpha(1.0, 0.0); + const std::complex beta(0.0, 0.0); #ifdef __MPI - const int i1 = 1; - pzgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, - &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); + const int i1 = 1; + pzgemm_(&transa, + &transb, + &nbasis, + &nbands, + &nbasis, + &alpha, + V, + &i1, + &i1, + pv.desc, + c, + &i1, + &i1, + pv.desc_wfc, + &beta, + Vc.data(), + &i1, + &i1, + pv.desc_wfc); #else - zgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, - c, &nbasis, &beta, Vc.data(), &nbasis); + zgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); #endif - std::vector> cVc(p2d.nloc, 0.0); - transa = 'C'; + std::vector> cVc(p2d.nloc, 0.0); + transa = 'C'; #ifdef __MPI - pzgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, - pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); + pzgemm_(&transa, + &transb, + &nbands, + &nbands, + &nbasis, + &alpha, + c, + &i1, + &i1, + pv.desc_wfc, + Vc.data(), + &i1, + &i1, + pv.desc_wfc, + &beta, + cVc.data(), + &i1, + &i1, + p2d.desc); #else - zgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, - Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); + zgemm_(&transa, + &transb, + &nbands, + &nbands, + &nbasis, + &alpha, + c, + &nbasis, + Vc.data(), + &nbasis, + &beta, + cVc.data(), + &nbasis); #endif - return cVc; - } + return cVc; +} - std::vector cVc( - double* V, - double* c, - int nbasis, - int nbands, - const Parallel_Orbitals& pv, - const Parallel_2D& p2d) - { - std::vector Vc(pv.nloc_wfc, 0.0); - char transa = 'N'; - char transb = 'N'; - const double alpha = 1.0; - const double beta = 0.0; +std::vector cVc(double* V, + double* c, + int nbasis, + int nbands, + const Parallel_Orbitals& pv, + const Parallel_2D& p2d) +{ + std::vector Vc(pv.nloc_wfc, 0.0); + char transa = 'N'; + char transb = 'N'; + const double alpha = 1.0; + const double beta = 0.0; #ifdef __MPI - const int i1 = 1; - pdgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, - &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); + const int i1 = 1; + pdgemm_(&transa, + &transb, + &nbasis, + &nbands, + &nbasis, + &alpha, + V, + &i1, + &i1, + pv.desc, + c, + &i1, + &i1, + pv.desc_wfc, + &beta, + Vc.data(), + &i1, + &i1, + pv.desc_wfc); #else - dgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, - c, &nbasis, &beta, Vc.data(), &nbasis); + dgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); #endif - std::vector cVc(p2d.nloc, 0.0); - transa = 'T'; + std::vector cVc(p2d.nloc, 0.0); + transa = 'T'; #ifdef __MPI - pdgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, - pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); + pdgemm_(&transa, + &transb, + &nbands, + &nbands, + &nbasis, + &alpha, + c, + &i1, + &i1, + pv.desc_wfc, + Vc.data(), + &i1, + &i1, + pv.desc_wfc, + &beta, + cVc.data(), + &i1, + &i1, + p2d.desc); #else - dgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, - &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); + dgemm_(&transa, + &transb, + &nbands, + &nbands, + &nbasis, + &alpha, + c, + &nbasis, + Vc.data(), + &nbasis, + &beta, + cVc.data(), + &nbasis); #endif - return cVc; - } + return cVc; +} - inline double get_real(const std::complex& c) { return c.real(); } +inline double get_real(const std::complex& c) +{ + return c.real(); +} - inline double get_real(const double& d) { return d; } +inline double get_real(const double& d) +{ + return d; +} - template - double all_band_energy( - const int ik, - const std::vector& mat_mo, - const Parallel_2D& p2d, - const ModuleBase::matrix& wg) +template +double all_band_energy(const int ik, const std::vector& mat_mo, const Parallel_2D& p2d, const ModuleBase::matrix& wg) +{ + double e = 0.0; + for (int i = 0; i < p2d.get_row_size(); ++i) { - double e = 0.0; - for (int i = 0;i < p2d.get_row_size();++i) - { - for (int j = 0;j < p2d.get_col_size();++j) - { - if (p2d.local2global_row(i) == p2d.local2global_col(j)) - { - e += get_real(mat_mo[j * p2d.get_row_size() + i]) * wg(ik, p2d.local2global_row(i)); - } - } - } - Parallel_Reduce::reduce_all(e); - return e; + for (int j = 0; j < p2d.get_col_size(); ++j) + { + if (p2d.local2global_row(i) == p2d.local2global_col(j)) + { + e += get_real(mat_mo[j * p2d.get_row_size() + i]) * wg(ik, p2d.local2global_row(i)); + } + } } + Parallel_Reduce::reduce_all(e); + return e; +} - template - std::vector orbital_energy( - const int ik, - const int nbands, - const std::vector& mat_mo, - const Parallel_2D& p2d) - { +template +std::vector orbital_energy(const int ik, const int nbands, const std::vector& mat_mo, const Parallel_2D& p2d) +{ #ifdef __DEBUG - assert(nbands >= 0); + assert(nbands >= 0); #endif - std::vector e(nbands, 0.0); - for (int i = 0;i < nbands;++i) + std::vector e(nbands, 0.0); + for (int i = 0; i < nbands; ++i) + { + if (p2d.in_this_processor(i, i)) { - if (p2d.in_this_processor(i, i)) - { - const int index = p2d.global2local_col(i) * p2d.get_row_size() + p2d.global2local_row(i); - e[i] = get_real(mat_mo[index]); - } + const int index = p2d.global2local_col(i) * p2d.get_row_size() + p2d.global2local_row(i); + e[i] = get_real(mat_mo[index]); } - Parallel_Reduce::reduce_all(e.data(), nbands); - return e; } + Parallel_Reduce::reduce_all(e.data(), nbands); + return e; +} - // mohan update 2024-04-01 - template - void set_gint_pointer( - Gint_Gamma &gint_gamma, - Gint_k &gint_k, - typename TGint::type*& gint); +// mohan update 2024-04-01 +template +void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint::type*& gint); - // mohan update 2024-04-01 - template <> - void set_gint_pointer( - Gint_Gamma &gint_gamma, - Gint_k &gint_k, - typename TGint::type*& gint) - { - gint = &gint_gamma; - } +// mohan update 2024-04-01 +template <> +void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint::type*& gint) +{ + gint = &gint_gamma; +} - // mohan update 2024-04-01 - template <> - void set_gint_pointer>( - Gint_Gamma &gint_gamma, - Gint_k &gint_k, - typename TGint>::type*& gint) - { - gint = &gint_k; - } +// mohan update 2024-04-01 +template <> +void set_gint_pointer>(Gint_Gamma& gint_gamma, + Gint_k& gint_k, + typename TGint>::type*& gint) +{ + gint = &gint_k; +} - /// @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( - int nspin, - int nbasis, - int drank, - const psi::Psi& psi, - const UnitCell& ucell, - Structure_Factor& sf, - const ModulePW::PW_Basis& rho_basis, - const ModulePW::PW_Basis& rhod_basis, - const ModuleBase::matrix& vloc, - const Charge& chg, - Gint_Gamma &gint_gamma, // mohan add 2024-04-01 - Gint_k &gint_k, // mohan add 2024-04-01 - LCAO_Matrix& lm, - Local_Orbital_Charge& loc, - const K_Vectors& kv, - const ModuleBase::matrix& wg, - Grid_Driver& gd) - { - ModuleBase::TITLE("ModuleIO", "write_Vxc"); - const Parallel_Orbitals* pv = lm.ParaV; - int nbands = wg.nc; - // 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, &etxc, &vtxc); - std::vector compnents_list = { "xc" }; - potxc->pot_register(compnents_list); - potxc->update_from_charge(&chg, &ucell); +/// @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(int nspin, + int nbasis, + int drank, + const psi::Psi& psi, + const UnitCell& ucell, + Structure_Factor& sf, + const ModulePW::PW_Basis& rho_basis, + const ModulePW::PW_Basis& rhod_basis, + const ModuleBase::matrix& vloc, + const Charge& chg, + Gint_Gamma& gint_gamma, // mohan add 2024-04-01 + Gint_k& gint_k, // mohan add 2024-04-01 + LCAO_Matrix& lm, + const K_Vectors& kv, + const ModuleBase::matrix& wg, + Grid_Driver& gd) +{ + ModuleBase::TITLE("ModuleIO", "write_Vxc"); + const Parallel_Orbitals* pv = lm.ParaV; + int nbands = wg.nc; + // 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, &etxc, &vtxc); + std::vector compnents_list = {"xc"}; + potxc->pot_register(compnents_list); + potxc->update_from_charge(&chg, &ucell); - // 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)); - for (int is = 0;is < nspin0;++is) vxcs_R_ao[is].set_zero(); - // k (size for each k-point) - std::vector vxc_k_ao(pv->nloc); + // 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)); + for (int is = 0; is < nspin0; ++is) + vxcs_R_ao[is].set_zero(); + // k (size for each k-point) + std::vector vxc_k_ao(pv->nloc); - // 3. allocate operators and contribute HR - // op (corresponding to hR) - typename TGint::type* gint = nullptr; + // 3. allocate operators and contribute HR + // op (corresponding to hR) + typename TGint::type* gint = nullptr; - set_gint_pointer(gint_gamma, gint_k, gint); + 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, - &loc, - &lm, - kv.kvec_d, - potxc, - &vxcs_R_ao[is], - &vxc_k_ao, - &ucell, - &gd, - pv); + std::vector>*> vxcs_op_ao(nspin0); + for (int is = 0; is < nspin0; ++is) + { + vxcs_op_ao[is] = new hamilt::Veff>(gint, + &lm, + kv.kvec_d, + potxc, + &vxcs_R_ao[is], + &vxc_k_ao, + &ucell, + &gd, + pv); - vxcs_op_ao[is]->contributeHR(); - } - std::vector> e_orb_locxc; // orbital energy (local XC) - std::vector> e_orb_tot; // orbital energy (total) + vxcs_op_ao[is]->contributeHR(); + } + std::vector> e_orb_locxc; // orbital energy (local XC) + std::vector> e_orb_tot; // orbital energy (total) #ifdef __EXX - hamilt::OperatorEXX> vexx_op_ao(&lm, nullptr, &vxc_k_ao, kv); - std::vector vexxonly_k_ao(pv->nloc); - hamilt::OperatorEXX> vexxonly_op_ao(&lm, nullptr, &vexxonly_k_ao, kv); - std::vector> e_orb_exx; // orbital energy (EXX) + hamilt::OperatorEXX> vexx_op_ao(&lm, nullptr, &vxc_k_ao, kv); + std::vector vexxonly_k_ao(pv->nloc); + hamilt::OperatorEXX> vexxonly_op_ao(&lm, nullptr, &vexxonly_k_ao, kv); + std::vector> e_orb_exx; // orbital energy (EXX) #endif - hamilt::OperatorDFTU> vdftu_op_ao(&lm, kv.kvec_d, nullptr, &vxc_k_ao, kv.isk); + hamilt::OperatorDFTU> vdftu_op_ao(&lm, kv.kvec_d, nullptr, &vxc_k_ao, kv.isk); - //4. calculate and write the MO-matrix Exc - Parallel_2D p2d; - set_para2d_MO(*pv, nbands, p2d); + // 4. calculate and write the MO-matrix Exc + Parallel_2D p2d; + set_para2d_MO(*pv, nbands, p2d); - // ======test======= - // double total_energy = 0.0; - // double exx_energy = 0.0; - // ======test======= - for (int ik = 0;ik < kv.get_nks();++ik) - { - ModuleBase::GlobalFunc::ZEROS(vxc_k_ao.data(), pv->nloc); - int is = kv.isk[ik]; - dynamic_cast*>(vxcs_op_ao[is])->contributeHk(ik); - const std::vector& vlocxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + // ======test======= + // double total_energy = 0.0; + // double exx_energy = 0.0; + // ======test======= + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + ModuleBase::GlobalFunc::ZEROS(vxc_k_ao.data(), pv->nloc); + int is = kv.isk[ik]; + dynamic_cast*>(vxcs_op_ao[is])->contributeHk(ik); + const std::vector& vlocxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) - { - e_orb_locxc.emplace_back(orbital_energy(ik, nbands, vlocxc_k_mo, p2d)); - ModuleBase::GlobalFunc::ZEROS(vexxonly_k_ao.data(), pv->nloc); - vexx_op_ao.contributeHk(ik); - vexxonly_op_ao.contributeHk(ik); - std::vector vexx_k_mo = cVc(vexxonly_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); - e_orb_exx.emplace_back(orbital_energy(ik, nbands, vexx_k_mo, p2d)); - } - // ======test======= - // exx_energy += all_band_energy(ik, vexx_k_mo, p2d, wg); - // ======test======= -#endif - if (GlobalV::dft_plus_u) - { - vdftu_op_ao.contributeHk(ik); - } - const std::vector& vxc_tot_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); - e_orb_tot.emplace_back(orbital_energy(ik, nbands, vxc_tot_k_mo, p2d)); - // write - ModuleIO::save_mat(-1, vxc_tot_k_mo.data(), nbands, - false/*binary*/, GlobalV::out_ndigits, true/*triangle*/, false/*append*/, - "Vxc", "k-" + std::to_string(ik), p2d, drank); - // ======test======= - // total_energy += all_band_energy(ik, vxc_tot_k_mo, p2d, wg); - // ======test======= + if (GlobalC::exx_info.info_global.cal_exx) + { + e_orb_locxc.emplace_back(orbital_energy(ik, nbands, vlocxc_k_mo, p2d)); + ModuleBase::GlobalFunc::ZEROS(vexxonly_k_ao.data(), pv->nloc); + vexx_op_ao.contributeHk(ik); + vexxonly_op_ao.contributeHk(ik); + std::vector vexx_k_mo = cVc(vexxonly_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + e_orb_exx.emplace_back(orbital_energy(ik, nbands, vexx_k_mo, p2d)); } // ======test======= - // total_energy -= 0.5 * exx_energy; - // std::cout << "total energy: " << total_energy << std::endl; - // std::cout << "etxc: " << etxc << std::endl; - // std::cout << "vtxc_cal: " << total_energy - 0.5 * exx_energy << std::endl; - // std::cout << "vtxc_ref: " << vtxc << std::endl; - // std::cout << "exx_energy: " << 0.5 * exx_energy << std::endl; + // exx_energy += all_band_energy(ik, vexx_k_mo, p2d, wg); // ======test======= - delete potxc; - for (int is = 0;is < nspin0;++is) - { - delete vxcs_op_ao[is]; +#endif + if (GlobalV::dft_plus_u) + { + vdftu_op_ao.contributeHk(ik); } - // write the orbital energy for xc and exx in LibRPA format - auto write_orb_energy = [&kv, &nspin0, &nbands](const std::vector>& e_orb, - const std::string& label, const bool app = false) - { - assert(e_orb.size() == kv.get_nks()); - const int nk = kv.get_nks() / nspin0; - std::ofstream ofs; - ofs.open(GlobalV::global_out_dir + "vxc_" + (label == "" ? "out" : label + "_out"), - app ? std::ios::app : std::ios::out); - ofs << nk << "\n" << nspin0 << "\n" << nbands << "\n"; - ofs << std::scientific << std::setprecision(16); - for (int ik = 0;ik < nk;++ik) - { - for (int is = 0;is < nspin0;++is) - { - for (auto e : e_orb[is * nk + ik]) - { // Hartree and eV - ofs << e / 2. << "\t" << e * ModuleBase::Ry_to_eV << "\n"; - } - } - } - }; - if (GlobalV::MY_RANK == 0) + const std::vector& vxc_tot_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + e_orb_tot.emplace_back(orbital_energy(ik, nbands, vxc_tot_k_mo, p2d)); + // write + ModuleIO::save_mat(-1, + vxc_tot_k_mo.data(), + nbands, + false /*binary*/, + GlobalV::out_ndigits, + true /*triangle*/, + false /*append*/, + "Vxc", + "k-" + std::to_string(ik), + p2d, + drank); + // ======test======= + // total_energy += all_band_energy(ik, vxc_tot_k_mo, p2d, wg); + // ======test======= + } + // ======test======= + // total_energy -= 0.5 * exx_energy; + // std::cout << "total energy: " << total_energy << std::endl; + // std::cout << "etxc: " << etxc << std::endl; + // std::cout << "vtxc_cal: " << total_energy - 0.5 * exx_energy << std::endl; + // std::cout << "vtxc_ref: " << vtxc << std::endl; + // std::cout << "exx_energy: " << 0.5 * exx_energy << std::endl; + // ======test======= + delete potxc; + for (int is = 0; is < nspin0; ++is) + { + delete vxcs_op_ao[is]; + } + // write the orbital energy for xc and exx in LibRPA format + auto write_orb_energy = [&kv, &nspin0, &nbands](const std::vector>& e_orb, + const std::string& label, + const bool app = false) { + assert(e_orb.size() == kv.get_nks()); + const int nk = kv.get_nks() / nspin0; + std::ofstream ofs; + ofs.open(GlobalV::global_out_dir + "vxc_" + (label == "" ? "out" : label + "_out"), + app ? std::ios::app : std::ios::out); + ofs << nk << "\n" << nspin0 << "\n" << nbands << "\n"; + ofs << std::scientific << std::setprecision(16); + for (int ik = 0; ik < nk; ++ik) { - write_orb_energy(e_orb_tot, ""); -#ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) + for (int is = 0; is < nspin0; ++is) { - write_orb_energy(e_orb_locxc, "local"); - write_orb_energy(e_orb_exx, "exx"); + for (auto e: e_orb[is * nk + ik]) + { // Hartree and eV + ofs << e / 2. << "\t" << e * ModuleBase::Ry_to_eV << "\n"; + } } -#endif } + }; + if (GlobalV::MY_RANK == 0) + { + write_orb_energy(e_orb_tot, ""); +#ifdef __EXX + if (GlobalC::exx_info.info_global.cal_exx) + { + write_orb_energy(e_orb_locxc, "local"); + write_orb_energy(e_orb_exx, "exx"); + } +#endif } } +} // namespace ModuleIO