From 27af2b214a96526aaf597da9fe837db5c368828d Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sun, 19 Nov 2023 02:02:38 +0800 Subject: [PATCH] matrix-lapack, reuse DM, fix wf-sign-dependence --- source/module_beyonddft/AX/test/AX_test.cpp | 2 +- .../dm_trans/test/dm_trans_test.cpp | 2 +- source/module_beyonddft/esolver_lrtd_lcao.h | 6 +- source/module_beyonddft/esolver_lrtd_lcao.hpp | 17 ++--- source/module_beyonddft/hamilt_casida.hpp | 73 ++++++++++++++++-- source/module_beyonddft/hsolver_lrtd.cpp | 11 +++ source/module_beyonddft/hsolver_lrtd.h | 3 +- .../operator_casida/operator_lr_exx.h | 18 +++-- .../operator_casida/operator_lr_exx.hpp | 31 +++++++- .../operator_casida/operator_lr_hxc.h | 20 ++--- .../operator_casida/operator_lr_hxc.hpp | 75 ++++++++++--------- .../potentials/pot_hxc_lrtd.hpp | 2 +- source/module_beyonddft/utils/lr_util.h | 9 ++- .../utils/lr_util_algorithms.hpp | 35 ++++++++- source/module_hamilt_general/hamilt.h | 2 +- 15 files changed, 221 insertions(+), 85 deletions(-) diff --git a/source/module_beyonddft/AX/test/AX_test.cpp b/source/module_beyonddft/AX/test/AX_test.cpp index 36e647d25f..d117dd0841 100644 --- a/source/module_beyonddft/AX/test/AX_test.cpp +++ b/source/module_beyonddft/AX/test/AX_test.cpp @@ -2,7 +2,7 @@ #include "mpi.h" #include "../AX.h" -#include "module_beyonddft/utils/lr_util.h" +#include "module_beyonddft/utils/lr_util_algorithms.hpp" struct matsize { diff --git a/source/module_beyonddft/dm_trans/test/dm_trans_test.cpp b/source/module_beyonddft/dm_trans/test/dm_trans_test.cpp index 4678c91f54..7b6ee7f654 100644 --- a/source/module_beyonddft/dm_trans/test/dm_trans_test.cpp +++ b/source/module_beyonddft/dm_trans/test/dm_trans_test.cpp @@ -2,7 +2,7 @@ #include "mpi.h" #include "../dm_trans.h" #ifdef __MPI -#include "module_beyonddft/utils/lr_util.h" +#include "module_beyonddft/utils/lr_util_algorithms.hpp" #endif struct matsize { diff --git a/source/module_beyonddft/esolver_lrtd_lcao.h b/source/module_beyonddft/esolver_lrtd_lcao.h index 8d1ee9a0fc..73582cb74d 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.h +++ b/source/module_beyonddft/esolver_lrtd_lcao.h @@ -45,7 +45,6 @@ namespace ModuleESolver delete this->phsol; delete this->pot; delete this->psi_ks; - delete this->DM_trans; delete this->X; } @@ -76,10 +75,7 @@ namespace ModuleESolver //pelec in ESolver_FP // const psi::Psi* psi_ks = nullptr; psi::Psi* psi_ks = nullptr; - ModuleBase::matrix eig_ks; - /// transition density matrix in AO representation - elecstate::DensityMatrix* DM_trans = nullptr; - // energy of ground state is in pelec->ekb + ModuleBase::matrix eig_ks;///< energy of ground state /// @brief Excited state info. size: nstates * nks * (nocc(local) * nvirt (local)) psi::Psi* X; diff --git a/source/module_beyonddft/esolver_lrtd_lcao.hpp b/source/module_beyonddft/esolver_lrtd_lcao.hpp index b71bb203e6..d0625c687d 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.hpp +++ b/source/module_beyonddft/esolver_lrtd_lcao.hpp @@ -255,11 +255,10 @@ ModuleESolver::ESolver_LRTD::ESolver_LRTD(Input& inp, UnitCell& ucell) : // if EXX from scratch, init 2-center integral and calclate Cs, Vs #ifdef __EXX - bool exx_from_scratch = false; - if (exx_from_scratch && xc_kernel == "hf") + if (xc_kernel == "hf") { - int Lmax = GlobalC::exx_info.info_ri.abfs_Lmax; #ifndef USE_NEW_TWO_CENTER + int Lmax = GlobalC::exx_info.info_ri.abfs_Lmax; this->orb_con.set_orb_tables(GlobalV::ofs_running, GlobalC::UOT, GlobalC::ORB, @@ -271,10 +270,11 @@ ModuleESolver::ESolver_LRTD::ESolver_LRTD(Input& inp, UnitCell& ucell) : ucell.infoNL.Beta); #else two_center_bundle.reset(new TwoCenterBundle); - two_center_bundle->build(ucell.ntype, ucell.orbital_fn, ucell.infoNL.Beta, + two_center_bundle->build(ucell.ntype, ucell.orbital_fn, nullptr/*ucell.infoNL.Beta*/, GlobalV::deepks_setorb, &ucell.descriptor_file); GlobalC::UOT.two_center_bundle = std::move(two_center_bundle); #endif + std::cout << GlobalC::exx_info.info_ri.dm_threshold << std::endl; this->exx_lri = std::make_shared>(GlobalC::exx_info.info_ri); this->exx_lri->init(MPI_COMM_WORLD, this->kv); // using GlobalC::ORB this->exx_lri->cal_exx_ions(); @@ -388,18 +388,15 @@ void ModuleESolver::ESolver_LRTD::init_A(hamilt::HContainer* pHR_ pHR->allocate(true); } pHR->set_paraV(&this->paraMat_); - this->DM_trans = new elecstate::DensityMatrix(&this->kv, &this->paraMat_, this->nspin); - this->DM_trans->init_DMR(*pHR); - this->p_hamilt = new hamilt::HamiltCasidaLR(xc_kernel, this->nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, this->psi_ks, this->eig_ks, this->DM_trans, pHR, + this->p_hamilt = new hamilt::HamiltCasidaLR(xc_kernel, this->nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, this->psi_ks, this->eig_ks, pHR, #ifdef __EXX this->exx_lri.get(), #endif - this->gint, this->pot, this->kv, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); + this->gint, this->pot, this->kv, & this->paraX_, & this->paraC_, & this->paraMat_); // init HSolver - this->phsol = new hsolver::HSolverLR(); + this->phsol = new hsolver::HSolverLR(this->npairs); this->phsol->set_diagethr(0, 0, std::max(1e-13, lr_thr)); - // if (!pHR_in) delete pHR; } template diff --git a/source/module_beyonddft/hamilt_casida.hpp b/source/module_beyonddft/hamilt_casida.hpp index 921e089040..96e39c1cc5 100644 --- a/source/module_beyonddft/hamilt_casida.hpp +++ b/source/module_beyonddft/hamilt_casida.hpp @@ -20,7 +20,7 @@ namespace hamilt const UnitCell& ucell_in, const psi::Psi* psi_ks_in, const ModuleBase::matrix& eig_ks, - elecstate::DensityMatrix* DM_trans_in, + // elecstate::DensityMatrix* DM_trans_in, HContainer*& hR_in, #ifdef __EXX Exx_LRI* exx_lri_in, @@ -28,22 +28,28 @@ namespace hamilt TGint* gint_in, elecstate::PotHxcLR* pot_in, const K_Vectors& kv_in, - const std::vector p2d_in) + Parallel_2D* pX_in, + Parallel_2D* pc_in, + Parallel_Orbitals* pmat_in) : nocc(nocc), nvirt(nvirt), pX(pX_in) { ModuleBase::TITLE("HamiltCasidaLR", "HamiltCasidaLR"); this->classname = "HamiltCasidaLR"; assert(hR_in != nullptr); this->hR = new HContainer(std::move(*hR_in)); + this->DM_trans.resize(1); + this->DM_trans[0] = new elecstate::DensityMatrix(&kv_in, pmat_in, nspin); + this->DM_trans[0]->init_DMR(*this->hR); // add the diag operator (the first one) - this->ops = new OperatorLRDiag(eig_ks, p2d_in.at(0), kv_in.nks, nspin, nocc, nvirt); + this->ops = new OperatorLRDiag(eig_ks, pX_in, kv_in.nks, nspin, nocc, nvirt); //add Hxc operator - OperatorLRHxc* lr_hxc = new OperatorLRHxc(nspin, naos, nocc, nvirt, psi_ks_in, DM_trans_in, this->hR, gint_in, pot_in, kv_in.kvec_d, p2d_in); + OperatorLRHxc* lr_hxc = new OperatorLRHxc(nspin, naos, nocc, nvirt, psi_ks_in, + this->DM_trans, this->hR, gint_in, pot_in, kv_in, pX_in, pc_in, pmat_in); this->ops->add(lr_hxc); #ifdef __EXX if (xc_kernel == "hf") - { - //add Exx operator - Operator* lr_exx = new OperatorLREXX(nspin, naos, nocc, nvirt, ucell_in, psi_ks_in, DM_trans_in, exx_lri_in, kv_in, p2d_in); + { //add Exx operator + Operator* lr_exx = new OperatorLREXX(nspin, naos, nocc, nvirt, ucell_in, psi_ks_in, + this->DM_trans, exx_lri_in, kv_in, pX_in, pc_in, pmat_in); this->ops->add(lr_exx); } #endif @@ -55,10 +61,63 @@ namespace hamilt delete this->ops; } delete this->hR; + for (auto& d : this->DM_trans)delete d; }; HContainer* getHR() { return this->hR; } + + virtual std::vector matrix() override + { + ModuleBase::TITLE("HamiltCasidaLR", "matrix"); + int npairs = this->nocc * this->nvirt; + std::vector Amat_full(npairs * npairs, 0.0); + for (int lj = 0;lj < this->pX->get_col_size();++lj) + for (int lb = 0;lb < this->pX->get_row_size();++lb) + {//calculate A^{ai} for each bj + int b = pX->local2global_row(lb); + int j = pX->local2global_col(lj); + int bj = j * nvirt + b; + psi::Psi X_bj(1, 1, this->pX->get_local_size()); + X_bj.zero_out(); + X_bj(0, 0, lj * this->pX->get_row_size() + lb) = this->one(); + psi::Psi A_aibj(1, 1, this->pX->get_local_size()); + A_aibj.zero_out(); + Operator* node(this->ops); + while (node != nullptr) + { + node->act(X_bj, A_aibj, 1); + node = (Operator*)(node->next_op); + } + // reduce ai for a fixed bj + LR_Util::gather_2d_to_full(*this->pX, A_aibj.get_pointer(), Amat_full.data() + bj * npairs, false, this->nvirt, this->nocc); + } + // reduce all bjs + MPI_Allreduce(MPI_IN_PLACE, Amat_full.data(), npairs * npairs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + // output Amat + std::cout << "Amat_full:" << std::endl; + for (int i = 0;i < npairs;++i) + { + for (int j = 0;j < npairs;++j) + { + std::cout << Amat_full[i * npairs + j] << " "; + } + std::cout << std::endl; + } + return Amat_full; + } private: + int nocc; + int nvirt; + Parallel_2D* pX = nullptr; + T one(); HContainer* hR = nullptr; + /// transition density matrix in AO representation + /// Hxc only: size=1, calculate on the same address for each bands + /// Hxc+Exx: size=nbands, store the result of each bands for common use + std::vector*> DM_trans; }; + + template<> double HamiltCasidaLR::one() { return 1.0; } + template<> std::complex HamiltCasidaLR>::one() { return std::complex(1.0, 0.0); } + } \ No newline at end of file diff --git a/source/module_beyonddft/hsolver_lrtd.cpp b/source/module_beyonddft/hsolver_lrtd.cpp index fffada7a25..af408914f1 100644 --- a/source/module_beyonddft/hsolver_lrtd.cpp +++ b/source/module_beyonddft/hsolver_lrtd.cpp @@ -1,6 +1,7 @@ #include "hsolver_lrtd.h" #include "module_hsolver/diago_david.h" #include "module_hsolver/diago_cg.h" +#include "module_beyonddft/utils/lr_util.h" namespace hsolver { @@ -30,6 +31,16 @@ namespace hsolver this->pdiagh = new DiagoCG(precondition.data()); this->pdiagh->method = this->method; } + else if (this->method == "lapack") + { + std::vector Amat_full = pHamilt->matrix(); + eigenvalue.resize(npairs); + LR_Util::diag_lapack(npairs, Amat_full.data(), eigenvalue.data()); + std::cout << "eigenvalues:" << std::endl; + for (auto& e : eigenvalue)std::cout << e << " "; + std::cout << std::endl; + return; + } else throw std::runtime_error("HSolverLR::solve: method not implemented"); diff --git a/source/module_beyonddft/hsolver_lrtd.h b/source/module_beyonddft/hsolver_lrtd.h index 52e2a69263..40cbdf46e1 100644 --- a/source/module_beyonddft/hsolver_lrtd.h +++ b/source/module_beyonddft/hsolver_lrtd.h @@ -9,8 +9,9 @@ namespace hsolver { private: using Real = typename GetTypeReal::type; + const int npairs = 0; public: - HSolverLR() {}; + HSolverLR(const int npairs_in) :npairs(npairs_in) {}; virtual Real set_diagethr(const int istep, const int iter, const Real ethr) override { this->diag_ethr = ethr; diff --git a/source/module_beyonddft/operator_casida/operator_lr_exx.h b/source/module_beyonddft/operator_casida/operator_lr_exx.h index 1f5b8bdf32..8af8ebde65 100644 --- a/source/module_beyonddft/operator_casida/operator_lr_exx.h +++ b/source/module_beyonddft/operator_casida/operator_lr_exx.h @@ -21,14 +21,16 @@ namespace hamilt const int& nvirt, const UnitCell& ucell_in, const psi::Psi* psi_ks_in, - elecstate::DensityMatrix* DM_trans_in, + std::vector*>& DM_trans_in, // HContainer* hR_in, Exx_LRI* exx_lri_in, const K_Vectors& kv_in, - const std::vector p2d_in /*< 2d-block parallel info of {X, c, matrix}*/) + Parallel_2D* pX_in, + Parallel_2D* pc_in, + Parallel_Orbitals* pmat_in) : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), psi_ks(psi_ks_in), DM_trans(DM_trans_in), exx_lri(exx_lri_in), kv(kv_in), - pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)), ucell(ucell_in) + pX(pX_in), pc(pc_in), pmat(pmat_in), ucell(ucell_in) { ModuleBase::TITLE("OperatorLREXX", "OperatorLREXX"); this->nks = std::is_same::value ? 1 : this->kv.kvec_d.size(); @@ -38,14 +40,15 @@ namespace hamilt this->is_first_node = false; // reduce psi_ks for later use - this->psi_ks_full.resize(this->nsk, this->nbands_ks, this->naos); - LR_Util::gather_2d_to_full(*this->pc, this->psi_ks->get_pointer(), this->psi_ks_full.get_pointer(), false, this->naos, this->nbands_ks); + this->psi_ks_full.resize(this->nsk, this->psi_ks->get_nbands(), this->naos); + LR_Util::gather_2d_to_full(*this->pc, this->psi_ks->get_pointer(), this->psi_ks_full.get_pointer(), false, this->naos, this->psi_ks->get_nbands()); // get cells in BvK supercell const TC period = RI_Util::get_Born_vonKarmen_period(kv_in); this->BvK_cells = RI_Util::get_Born_von_Karmen_cells(period); this->allocate_Ds_onebase(); + this->exx_lri->Hexxs.resize(this->nspin); }; void init(const int ik_in) override {}; @@ -60,14 +63,13 @@ namespace hamilt int naos; int nocc; int nvirt; - int nbands_ks; const K_Vectors& kv; /// ground state wavefunction const psi::Psi* psi_ks = nullptr; psi::Psi psi_ks_full; /// transition density matrix - elecstate::DensityMatrix* DM_trans; + std::vector*>& DM_trans; /// density matrix of a certain (i, a, k), with full naos*naos size for each key /// D^{iak}_{\mu\nu}(k): 1/N_k * c^*_{ak,\mu} c_{ik,\nu} @@ -92,7 +94,7 @@ namespace hamilt ///parallel info Parallel_2D* pc = nullptr; Parallel_2D* pX = nullptr; - Parallel_2D* pmat = nullptr; + Parallel_Orbitals* pmat = nullptr; // allocate Ds_onebase diff --git a/source/module_beyonddft/operator_casida/operator_lr_exx.hpp b/source/module_beyonddft/operator_casida/operator_lr_exx.hpp index a7d937b37e..c2eeaeb988 100644 --- a/source/module_beyonddft/operator_casida/operator_lr_exx.hpp +++ b/source/module_beyonddft/operator_casida/operator_lr_exx.hpp @@ -32,7 +32,7 @@ namespace hamilt int iat2 = ucell.itia2iat(it2, ia2); auto& D2d = this->Ds_onebase[is][iat1][std::make_pair(iat2, cell)]; for (int iw1 = 0;iw1 < ucell.atoms[it1].nw;++iw1) - for (int iw2 = 0;iw1 < ucell.atoms[it2].nw;++iw2) + for (int iw2 = 0;iw2 < ucell.atoms[it2].nw;++iw2) D2d(iw1, iw2) = this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1)) * this->psi_ks_full(ik, iv, ucell.itiaiw2iwt(it2, ia2, iw2)); } } @@ -80,7 +80,8 @@ namespace hamilt // 1. set_Ds (once) // convert to vector for the interface of RI_2D_Comm::split_m2D_ktoR (interface will be unified to ct::Tensor) - std::vector> DMk_trans_vector = this->DM_trans->get_DMK_vector(); + std::cout << "ib=" << ib << std::endl; + std::vector> DMk_trans_vector = this->DM_trans[ib]->get_DMK_vector(); assert(DMk_trans_vector.size() == this->nsk); std::vector*> DMk_trans_pointer(this->nsk); for (int is = 0;is < this->nsk;++is) DMk_trans_pointer[is] = &DMk_trans_vector[is]; @@ -88,6 +89,30 @@ namespace hamilt std::vector>>> Ds_trans = RI_2D_Comm::split_m2D_ktoR(this->kv, DMk_trans_pointer, *this->pmat); + // output DM_trans + GlobalV::ofs_running << "DM_trans in OperatorLREXX::act" << std::endl; + for (int is = 0;is < this->nsk;++is) + { + GlobalV::ofs_running << "is = " << is << std::endl; + for (auto& mabR : Ds_trans[is + ]) + { + int ia = mabR.first; + GlobalV::ofs_running << "ia = " << ia << std::endl; + for (auto& mbR : mabR.second) + { + GlobalV::ofs_running << "ib = " << mbR.first.first << ", R=" << mbR.first.second[0] << + " " << mbR.first.second[1] << " " << mbR.first.second[2] << ", size=" << mbR.second.shape[0] << ", " << mbR.second.shape[1] << std::endl; + auto& t = mbR.second; + for (int i = 0;i < t.shape[0];++i) + { + for (int j = 0;j < t.shape[1];++j) + GlobalV::ofs_running << t(i, j) << " "; + GlobalV::ofs_running << std::endl; + } + } + } + } // 2. cal_Hs for (int is = 0;is < this->nsk;++is) { @@ -110,7 +135,7 @@ namespace hamilt for (int is = 0;is < this->nspin;++is) { this->cal_DM_onebase(this->pX->local2global_col(io), this->pX->local2global_row(iv), ik, is); //set Ds_onebase - psi_out_bfirst(ik, io * this->pX->get_row_size() + iv) = + psi_out_bfirst(ik, io * this->pX->get_row_size() + iv) += this->exx_lri->exx_lri.post_2D.cal_energy(this->Ds_onebase[is], this->exx_lri->Hexxs[is]); } } diff --git a/source/module_beyonddft/operator_casida/operator_lr_hxc.h b/source/module_beyonddft/operator_casida/operator_lr_hxc.h index c6d0b1915e..86ea71e2fa 100644 --- a/source/module_beyonddft/operator_casida/operator_lr_hxc.h +++ b/source/module_beyonddft/operator_casida/operator_lr_hxc.h @@ -28,18 +28,20 @@ namespace hamilt const int& nocc, const int& nvirt, const psi::Psi* psi_ks_in, - elecstate::DensityMatrix* DM_trans_in, + std::vector*>& DM_trans_in, HContainer* hR_in, typename TGint::type* gint_in, elecstate::PotHxcLR* pot_in, - const std::vector>& kvec_d_in, - const std::vector p2d_in /*< 2d-block parallel info of {X, c, matrix}*/) + const K_Vectors& kv_in, + Parallel_2D* pX_in, + Parallel_2D* pc_in, + Parallel_Orbitals* pmat_in) : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), - psi_ks(psi_ks_in), DM_trans(DM_trans_in), hR(hR_in), gint(gint_in), pot(pot_in), kvec_d(kvec_d_in), - pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) + psi_ks(psi_ks_in), DM_trans(DM_trans_in), hR(hR_in), gint(gint_in), pot(pot_in), kv(kv_in), + pX(pX_in), pc(pc_in), pmat(pmat_in) { ModuleBase::TITLE("OperatorLRHxc", "OperatorLRHxc"); - this->nks = std::is_same::value ? 1 : kvec_d.size(); + this->nks = std::is_same::value ? 1 : this->kv.kvec_d.size(); this->nsk = std::is_same::value ? nspin : nks; this->cal_type = calculation_type::lcao_gint; this->act_type = 2; @@ -58,12 +60,12 @@ namespace hamilt int naos; int nocc; int nvirt; - const std::vector>& kvec_d; + const K_Vectors& kv; /// ground state wavefunction const psi::Psi* psi_ks = nullptr; /// transition density matrix - elecstate::DensityMatrix* DM_trans; + std::vector*>& DM_trans; /// transition hamiltonian in AO representation hamilt::HContainer* hR = nullptr; @@ -71,7 +73,7 @@ namespace hamilt //parallel info Parallel_2D* pc = nullptr; Parallel_2D* pX = nullptr; - Parallel_2D* pmat = nullptr; + Parallel_Orbitals* pmat = nullptr; elecstate::PotHxcLR* pot = nullptr; diff --git a/source/module_beyonddft/operator_casida/operator_lr_hxc.hpp b/source/module_beyonddft/operator_casida/operator_lr_hxc.hpp index c3de8f08f1..0acd5bf7da 100644 --- a/source/module_beyonddft/operator_casida/operator_lr_hxc.hpp +++ b/source/module_beyonddft/operator_casida/operator_lr_hxc.hpp @@ -2,7 +2,7 @@ #include "operator_lr_hxc.h" #include #include "module_base/blas_connector.h" -#include "module_beyonddft/utils/lr_util.h" +#include "module_beyonddft/utils/lr_util_algorithms.hpp" // #include "module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" #include "module_beyonddft/dm_trans/dm_trans.h" @@ -15,14 +15,30 @@ namespace hamilt void OperatorLRHxc::act(const psi::Psi& psi_in, psi::Psi& psi_out, const int nbands) const { ModuleBase::TITLE("OperatorLRHxc", "act"); - assert(nbands <= psi_in.get_nbands()); + + /// initialize transion density matrix + if (this->next_op != nullptr) + { + int prev_size = this->DM_trans.size(); + std::cout << "prev_size= " << prev_size << std::endl; + std::cout << "nbands=" << nbands << std::endl; + if (prev_size > nbands)for (int ib = nbands;ib < prev_size;++ib)delete this->DM_trans[ib]; + this->DM_trans.resize(nbands); + for (int ib = prev_size;ib < nbands;++ib) + { + this->DM_trans[ib] = new elecstate::DensityMatrix(&this->kv, this->pmat, this->nspin); + this->DM_trans[ib]->init_DMR(*this->hR); + } + } + psi::Psi psi_in_bfirst = LR_Util::k1_to_bfirst_wrapper(psi_in, this->nsk, this->pX->get_local_size()); psi::Psi psi_out_bfirst = LR_Util::k1_to_bfirst_wrapper(psi_out, this->nsk, this->pX->get_local_size()); const int& lgd = gint->gridt->lgd; for (int ib = 0;ib < nbands;++ib) { + int ib_dm = (this->next_op == nullptr) ? 0 : ib; GlobalV::ofs_running << "ib=" << ib << std::endl; psi_in_bfirst.fix_b(ib); psi_out_bfirst.fix_b(ib); @@ -58,7 +74,7 @@ namespace hamilt if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos, *pmat); #else std::vector dm_trans_2d = cal_dm_trans_blas(psi_in_bfirst, psi_ks, nocc, nvirt); - if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos) + if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos); #endif // GlobalV::ofs_running << "1. Dm_trans: " << std::endl; // GlobalV::ofs_running << "local row:" << pmat->get_row_size() << " local col:" << pmat->get_col_size() << std::endl; @@ -70,19 +86,19 @@ namespace hamilt // } // GlobalV::ofs_running << std::endl; // tensor to vector, then set DMK - for (int isk = 0;isk < this->nsk;++isk)this->DM_trans->set_DMK_pointer(isk, dm_trans_2d[isk].data()); + for (int isk = 0;isk < this->nsk;++isk)this->DM_trans[ib_dm]->set_DMK_pointer(isk, dm_trans_2d[isk].data()); // use cal_DMR to get DMR form DMK by FT - this->DM_trans->cal_DMR(); //DM_trans->get_DMR_vector() is 2d-block parallized + this->DM_trans[ib_dm]->cal_DMR(); //DM_trans->get_DMR_vector() is 2d-block parallized GlobalV::ofs_running << "return cal_DMR (outside)" << std::endl; // 2d block to grid // new interface: transfer_DM2DtoGrid, set DMRGint for the next step Gint // this->gint->transfer_DM2DtoGrid(this->gint->get_DMRGint());//err? - this->gint->transfer_DM2DtoGrid(this->DM_trans->get_DMR_vector()); - // this->gint->set_DMRGint(this->DM_trans->get_DMR_vector()); + this->gint->transfer_DM2DtoGrid(this->DM_trans[ib_dm]->get_DMR_vector()); + // this->gint->set_DMRGint(this->DM_trans[ib_dm]->get_DMR_vector()); // GlobalV::ofs_running << "return transfer_DMR (outside)" << std::endl; // GlobalV::ofs_running << "2. DM(R) (2d): " << std::endl; - // for (auto& dr : this->DM_trans->get_DMR_vector()) + // for (auto& dr : this->DM_trans[ib_dm]->get_DMR_vector()) // { // for (int ia = 0;ia < GlobalC::ucell.nat;ia++) // for (int ja = 0;ja < GlobalC::ucell.nat;ja++) @@ -163,46 +179,37 @@ namespace hamilt GlobalV::ofs_running << "4.Vxc" << std::endl; std::vector v_hxc_2d(nsk, ct::Tensor(ct::DataTypeToEnum::value, ct::DeviceTypeToEnum::value, { pmat->get_col_size(), pmat->get_row_size() })); for (auto& v : v_hxc_2d) v.zero(); - // auto setter = [this](const int& iw1_all, const int& iw2_all, const double& v, double* out) { - // const int ir = this->pmat->global2local_row(iw1_all); - // const int ic = this->pmat->global2local_col(iw2_all); - // out[ic * this->pmat->nrow + ir] += v; - // }; // V(R) for each spin for (int is = 0;is < nspin;++is) { double* vr_hxc_is = &vr_hxc.c[is * this->pot->nrxx]; //v(r) at current spin - - double tot_vreff = 0.0; - for (int i = 0;i < this->pot->nrxx;++i)tot_vreff += vr_hxc_is[i]; - // std::cout << "total vreff in LR = " << tot_vreff << std::endl;// normal - Gint_inout inout_vlocal(vr_hxc_is, is, Gint_Tools::job_type::vlocal); this->gint->get_hRGint()->set_zero(); this->gint->cal_gint(&inout_vlocal); } + this->hR->set_zero(); // clear hR for each bands this->gint->transfer_pvpR(this->hR); - // GlobalV::ofs_running << "4.V(R):" << std::endl; - // for (int ia = 0;ia < GlobalC::ucell.nat;ia++) - // for (int ja = 0;ja < GlobalC::ucell.nat;ja++) - // { - // auto ap = this->hR->find_pair(ia, ja); - // GlobalV::ofs_running << "R-index size of atom pair(" << ia << ", " << ja << "): " << ap->get_R_size() << std::endl; - // for (int iR = 0;iR < ap->get_R_size();++iR) - // { - // GlobalV::ofs_running << "R(" << ap->get_R_index(iR)[0] << ", " << ap->get_R_index(iR)[1] << ", " << ap->get_R_index(iR)[2] << "): "; - // auto ptr = ap->get_HR_values(iR).get_pointer(); - // int size = ap->get_size(); - // for (int i = 0;i < size;++i)GlobalV::ofs_running << ptr[i] << " "; - // GlobalV::ofs_running << std::endl; - // } - // } + GlobalV::ofs_running << "4.V(R):" << std::endl; + for (int ia = 0;ia < GlobalC::ucell.nat;ia++) + for (int ja = 0;ja < GlobalC::ucell.nat;ja++) + { + auto ap = this->hR->find_pair(ia, ja); + GlobalV::ofs_running << "R-index size of atom pair(" << ia << ", " << ja << "): " << ap->get_R_size() << std::endl; + for (int iR = 0;iR < ap->get_R_size();++iR) + { + GlobalV::ofs_running << "R(" << ap->get_R_index(iR)[0] << ", " << ap->get_R_index(iR)[1] << ", " << ap->get_R_index(iR)[2] << "): "; + auto ptr = ap->get_HR_values(iR).get_pointer(); + int size = ap->get_size(); + for (int i = 0;i < size;++i)GlobalV::ofs_running << ptr[i] << " "; + GlobalV::ofs_running << std::endl; + } + } // V(R)->V(k) int nrow = ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER() ? this->pmat->get_row_size() : this->pmat->get_col_size(); for (int isk = 0;isk < this->nsk;++isk) { // this->gint->vl_grid_to_2D(this->gint->get_pvpR_grid(), *pmat, lgd, (is == 0), v_hxc_2d[is].c, setter); - hamilt::folding_HR(*this->hR, v_hxc_2d[isk].data(), this->kvec_d[isk], nrow, 1); // V(R) -> V(k) + hamilt::folding_HR(*this->hR, v_hxc_2d[isk].data(), this->kv.kvec_d[isk], nrow, 1); // V(R) -> V(k) } GlobalV::ofs_running << "4. V(k)" << std::endl; for (int isk = 0;isk < this->nsk;++isk) diff --git a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp index decb3a1104..4a71075514 100644 --- a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp +++ b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp @@ -43,7 +43,7 @@ namespace elecstate ModuleBase::timer::tick("PotHxcLR", "cal_v_eff"); const int nspin = v_eff.nr; v_eff += H_Hartree_pw::v_hartree(*ucell, const_cast(this->rho_basis_), v_eff.nr, rho); - if (xc_kernel == "rpa") return; + if (xc_kernel == "rpa" || xc_kernel == "hf") return; else if (XC_Functional::get_func_type() == 1)//LDA if (1 == nspin)// for LDA-spin0, just f*rho for (int ir = 0;ir < nrxx;++ir) diff --git a/source/module_beyonddft/utils/lr_util.h b/source/module_beyonddft/utils/lr_util.h index 7e6514b3af..d901e64744 100644 --- a/source/module_beyonddft/utils/lr_util.h +++ b/source/module_beyonddft/utils/lr_util.h @@ -104,6 +104,11 @@ namespace LR_Util template void gather_2d_to_full(const Parallel_2D& pv, const T* submat, T* fullmat, bool col_first, int global_nrow, int global_ncol); #endif + + ///=================diago-lapack==================== + /// @brief diagonalize a hermitian matrix + template + void diag_lapack(const int& n, T* mat, double* eig); } -#include "lr_util_algorithms.hpp" -#include "lr_util_physics.hpp" \ No newline at end of file +// #include "lr_util_algorithms.hpp" +// #include "lr_util_physics.hpp" \ No newline at end of file diff --git a/source/module_beyonddft/utils/lr_util_algorithms.hpp b/source/module_beyonddft/utils/lr_util_algorithms.hpp index e59f13c69b..cad4ecd4c5 100644 --- a/source/module_beyonddft/utils/lr_util_algorithms.hpp +++ b/source/module_beyonddft/utils/lr_util_algorithms.hpp @@ -4,6 +4,7 @@ #include #include "module_cell/unitcell.h" #include "module_base/constants.h" +#include "module_base/lapack_connector.h" #include "module_base/scalapack_connector.h" @@ -261,7 +262,7 @@ namespace LR_Util void setup_2d_division(Parallel_2D& pv, int nb, int gr, int gc, const MPI_Comm& comm_2D_in, const int& blacs_ctxt_in) { - ModuleBase::TITLE("ESolver_LRTD", "setup_2d_division"); + ModuleBase::TITLE("LR_Util", "setup_2d_division"); std::ofstream ofs(""); pv.set_block_size(nb); @@ -304,4 +305,34 @@ namespace LR_Util MPI_Allreduce(MPI_IN_PLACE, fullmat, global_nrow * global_ncol, get_mpi_datatype(), MPI_SUM, pv.comm_2D); }; #endif -} \ No newline at end of file + + template<> + void diag_lapack(const int& n, double* mat, double* eig) + { + ModuleBase::TITLE("LR_Util", "diag_lapack"); + int info = 0; + char jobz = 'V', uplo = 'U'; + double work_tmp; + constexpr int minus_one = -1; + dsyev_(&jobz, &uplo, &n, mat, &n, eig, &work_tmp, &minus_one, &info); // get best lwork + const int lwork = work_tmp; + double* work2 = new double[lwork]; + dsyev_(&jobz, &uplo, &n, mat, &n, eig, work2, &lwork, &info); + if (info) std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + delete[] work2; + } + template<> + void diag_lapack>(const int& n, std::complex* mat, double* eig) + { + ModuleBase::TITLE("LR_Util", "diag_lapack>"); + int lwork = 2 * n; + std::complex* work2 = new std::complex[lwork]; + double* rwork = new double[3 * n - 2]; + int info = 0; + char jobz = 'V', uplo = 'U'; + zheev_(&jobz, &uplo, &n, mat, &n, eig, work2, &lwork, rwork, &info); + if (info) std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + delete[] rwork; + delete[] work2; + } +} diff --git a/source/module_hamilt_general/hamilt.h b/source/module_hamilt_general/hamilt.h index 8ddd1e3e15..535e74da8a 100644 --- a/source/module_hamilt_general/hamilt.h +++ b/source/module_hamilt_general/hamilt.h @@ -38,7 +38,7 @@ class Hamilt /// core function: return H(k) and S(k) matrixs for direct solving eigenvalues. virtual void matrix(MatrixBlock> &hk_in, MatrixBlock> &sk_in){return;} virtual void matrix(MatrixBlock &hk_in, MatrixBlock &sk_in){return;} - + virtual std::vector matrix() { return std::vector(); } std::string classname = "none"; int non_first_scf=0;