From fd653dba42daf4e08c893df2ff834f479af0cdd9 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 7 Sep 2023 12:39:01 +0800 Subject: [PATCH] Procedure: build gamma_only AX --- source/module_beyonddft/esolver_lrtd_lcao.h | 11 ++- source/module_beyonddft/esolver_lrtd_lcao.hpp | 10 +++ source/module_beyonddft/hamilt_casida.hpp | 2 +- .../operator_casida/operatorA_hxc.h | 45 +++++----- .../operator_casida/operatorA_hxc.hpp | 85 ++++++++++++------- .../potentials/pot_hxc_lrtd.hpp | 6 +- 6 files changed, 98 insertions(+), 61 deletions(-) diff --git a/source/module_beyonddft/esolver_lrtd_lcao.h b/source/module_beyonddft/esolver_lrtd_lcao.h index 69b8ac4ad5..37a4a43266 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.h +++ b/source/module_beyonddft/esolver_lrtd_lcao.h @@ -16,7 +16,8 @@ #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_gint/grid_technique.h" -#include "module_elecstate/potentials/potential_new.h" +#include "module_beyonddft/potentials/pot_hxc_lrtd.hpp" +#include "module_beyonddft/hamilt_casida.hpp" namespace ModuleESolver { @@ -52,12 +53,12 @@ namespace ModuleESolver const UnitCell* p_ucell = nullptr; const Input* p_input = nullptr; - hamilt::Hamilt* p_hamilt = nullptr; + hamilt::HamiltCasidaLR* p_hamilt = nullptr; hsolver::HSolver* phsol = nullptr; // not to use ElecState because 2-particle state is quite different from 1-particle state. // implement a independent one (ExcitedState) to pack physical properties if needed. // put the components of ElecState here: - elecstate::PotBase* pot = nullptr; + elecstate::PotHxcLR* pot = nullptr; // ground state info //pelec in ESolver_FP @@ -67,7 +68,7 @@ namespace ModuleESolver /// @brief Excited state info. size: nstates*nks*nocc*nvirt std::vector> X; - //psi::Psi* AX = nullptr; + std::vector> AX; size_t nocc; size_t nvirt; @@ -96,6 +97,8 @@ namespace ModuleESolver Parallel_2D paraX_; /// @brief variables for parallel distribution of matrix in AO representation Parallel_2D paraMat_; + + /// move to hsolver::updatePsiK void init_X(); }; diff --git a/source/module_beyonddft/esolver_lrtd_lcao.hpp b/source/module_beyonddft/esolver_lrtd_lcao.hpp index 530664743c..6d70a6de6c 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.hpp +++ b/source/module_beyonddft/esolver_lrtd_lcao.hpp @@ -86,6 +86,16 @@ ModuleESolver::ESolver_LRTD::ESolver_LRTD(ModuleESolver::ESolver this->p_hamilt = new hamilt::HamiltCasidaLR(this->nks, this->nbasis, this->nocc, this->nvirt, &this->gint_k, this->pot, &this->gridt, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); + // init HSolver + // this->phsol = new hsolver::HSolver + + // try hPsi + int istate = 0; + int iks = 0; + using hpsi_info = typename hamilt::Operator::hpsi_info; + hpsi_info dav_info(&this->X[istate], psi::Range(1, iks, 0, this->nocc - 1), this->AX[istate].get_pointer()); + // this->p_hamilt->ops->hPsi(dav_info); + } template diff --git a/source/module_beyonddft/hamilt_casida.hpp b/source/module_beyonddft/hamilt_casida.hpp index e6eec83436..48fbeca9aa 100644 --- a/source/module_beyonddft/hamilt_casida.hpp +++ b/source/module_beyonddft/hamilt_casida.hpp @@ -13,7 +13,7 @@ namespace hamilt const int& nocc, const int& nvirt, TGint* gint_in, - elecstate::PotBase* pot_in, + elecstate::PotHxcLR* pot_in, const Grid_Technique* gt, const std::vector p2d_in) { diff --git a/source/module_beyonddft/operator_casida/operatorA_hxc.h b/source/module_beyonddft/operator_casida/operatorA_hxc.h index 300c68313f..73e2444bf2 100644 --- a/source/module_beyonddft/operator_casida/operatorA_hxc.h +++ b/source/module_beyonddft/operator_casida/operatorA_hxc.h @@ -1,16 +1,16 @@ +#pragma once #include "module_hamilt_general/operator.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/grid_technique.h" #include "module_hamilt_lcao/module_gint/gint_k.h" - +#include "module_beyonddft/potentials/pot_hxc_lrtd.hpp" namespace hamilt { - /// @brief kernel for Hxc operator - /// @tparam T - template - class OperatorA_Hxc : public Operator + /// @brief Hxc part of A operator + template + class OperatorA_Hxc : public Operator { public: OperatorA_Hxc(const int& nsk, @@ -18,12 +18,12 @@ namespace hamilt const int& nocc, const int& nvirt, Gint_Gamma* gg_in, - elecstate::PotBase* pot_in, + elecstate::PotHxcLR* pot_in, const Grid_Technique* gt_in/* grid parallel info*/, const std::vector p2d_in /*< 2d-block parallel info of {X, c, matrix}*/) : nsk(nsk), naos(naos), nocc(nocc), nvirt(nvirt), gg(gg_in), pot(pot_in), gridt(gt_in), - px(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) + pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) { this->cal_type = calculation_type::lcao_gint; this->is_first_node = true; @@ -33,12 +33,12 @@ namespace hamilt const int& nocc, const int& nvirt, Gint_k* gk_in, - elecstate::PotBase* pot_in, + elecstate::PotHxcLR* pot_in, const Grid_Technique* gt_in/* grid parallel info*/, const std::vector p2d_in /*< 2d-block parallel info of {X, c, matrix}*/) : nsk(nsk), naos(naos), nocc(nocc), nvirt(nvirt), gk(gk_in), pot(pot_in), gridt(gt_in), - px(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) + pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) { this->cal_type = calculation_type::lcao_gint; this->is_first_node = true; @@ -46,8 +46,14 @@ namespace hamilt void init(const int ik_in) override {}; //move hpsi and act to base class ! // void act() const override; // call gint - void act() const; - + virtual void act(const int nbands, + const int nbasis, + const int npol, + const FPTYPE* tmpsi_in, + FPTYPE* tmhpsi, + const int ngk_ik = 0)const override {}; + //tmp, for only one state + virtual psi::Psi act(const psi::Psi& psi_in) const override; private: //global sizes int nsk; //nspin*nkpoints @@ -55,21 +61,20 @@ namespace hamilt int nocc; int nvirt; + /// ground state wavefunction + const psi::Psi* psi_ks = nullptr; //parallel info - const Parallel_2D* pc = nullptr; - const Parallel_2D* px = nullptr; - const Parallel_2D* pmat = nullptr; + Parallel_2D* pc = nullptr; + Parallel_2D* pX = nullptr; + Parallel_2D* pmat = nullptr; const Grid_Technique* gridt = nullptr; - elecstate::PotBase* pot = nullptr; + elecstate::PotHxcLR* pot = nullptr; Gint_Gamma* gg = nullptr; Gint_k* gk = nullptr; - - /// \f[ \tilde{\rho}(r)=\sum_{\mu_j, \mu_b}\tilde{\rho}_{\mu_j,\mu_b}\phi_{\mu_b}(r)\phi_{\mu_j}(r) \f] void cal_rho_trans(); }; - - -} \ No newline at end of file +} +#include "module_beyonddft/operator_casida/operatorA_hxc.hpp" \ No newline at end of file diff --git a/source/module_beyonddft/operator_casida/operatorA_hxc.hpp b/source/module_beyonddft/operator_casida/operatorA_hxc.hpp index c2b709d293..d696ca7dfb 100644 --- a/source/module_beyonddft/operator_casida/operatorA_hxc.hpp +++ b/source/module_beyonddft/operator_casida/operatorA_hxc.hpp @@ -1,56 +1,75 @@ +#pragma once #include "operatorA_hxc.h" #include #include "module_base/blas_connector.h" -#include "utils/lr_util.h" - +#include "module_beyonddft/utils/lr_util.h" +#include "module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h" +#include "module_beyonddft/density/dm_trans.h" +#include "module_beyonddft/AX/AX.h" namespace hamilt { // for double - template - void OperatorA_Hxc::act() const + template + psi::Psi OperatorA_Hxc::act(const psi::Psi& psi_in) const { - auto block2grid = [](std::vector block, double*** grid)->void {}; + // gamma-only now // 1. transition density matrix - std::vector dm_trans_2d = cal_dm_trans_blas(*px, *pc); - double*** dm_trans_grid = LR_Util::new_p3(nspin, naos_local_grid, naos_local_grid); - block2grid(dm_trans_2d, dm_trans_grid); + // multi-k needs k-to-R FT +#ifdef __MPI + std::vector dm_trans_2d = cal_dm_trans_pblas(psi_in, *pX, *psi_ks, *pc, naos, nocc, nvirt, *pmat); +#else + std::vector dm_trans_2d = cal_dm_trans_blas(*pX, *pc); +#endif + double*** dm_trans_grid; + LR_Util::new_p3(dm_trans_grid, nsk, gridt->lgd, gridt->lgd); + //2d block to grid + DMgamma_2dtoGrid dm2g; +#ifdef __MPI + dm2g.setAlltoallvParameter(pmat->comm_2D, naos, pmat->blacs_ctxt, pmat->nb, gridt->lgd, gridt->trace_lo); +#endif + dm2g.cal_dk_gamma_from_2D(LR_Util::ten2mat_double(dm_trans_2d), dm_trans_grid, nsk, naos, gridt->lgd, GlobalV::ofs_running); // 2. transition electron density - double** rho_trans = LR_Util::new_p2(nspin, this->gint_g->nbxx); // is nbxx local grid num ? + double** rho_trans; + LR_Util::new_p2(rho_trans, nsk, this->pot->nrxx); Gint_inout inout_rho(dm_trans_grid, rho_trans, Gint_Tools::job_type::rho); - this->gint_g->cal_gint(&inout_rho); + this->gg->cal_gint(&inout_rho); // 3. v_hxc = f_hxc * rho_trans - this->pot->update_from_charge(rho_trans, GlobalC::ucell); + ModuleBase::matrix vr_hxc(nsk, this->pot->nrxx); //grid + this->pot->cal_v_eff(rho_trans, &GlobalC::ucell, vr_hxc); // 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r) - // loop for nspin, or use current spin (how?) - // results are stored in gint_g->pvpR_grid(gamma_only) + // loop for nsk, or use current spin (how?) + // results are stored in gg->pvpR_grid(gamma_only) // or gint_k->pvpR_reduced(multi_k) - - std::vector v_hxc_local(nspin); // 2D local matrix) - for (int is = 0;is < this->nspin;++is) + std::vector v_hxc_2d(nsk); + this->gg->init_pvpR_grid(gridt->lgd); + 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; + }; + for (int is = 0;is < this->nsk;++is) { - const double* v1_hxc_grid = this->pot->get_effective_v(is); - Gint_intout inout_vlocal(v1_hxc_grid, is, Gint_Tools::job_type::vlocal); - // this->gint_g->cal_gint(&inout); - bool new_e_iteration = false; // what is this? - this->gint_g->cal_vlocal(&inout_vlocal, new_e_iteration); - - // grid-to-2d needs refactor ! - v_hxc_2d[is].create(naos_local_row, naos_local_col); - //LR_Util::grid2block(this->px, this->pc, this->gint_g->pvpR_grid, v_hxc_local.c); + ModuleBase::GlobalFunc::ZEROS(this->gg->get_pvpR_grid(), gridt->lgd * gridt->lgd); + double* vr_hxc_is = &vr_hxc.c[is * this->pot->nrxx]; //current spin + Gint_inout inout_vlocal(vr_hxc_is, Gint_Tools::job_type::vlocal); + this->gg->cal_gint(&inout_vlocal); + v_hxc_2d[is].create(pmat->get_col_size(), pmat->get_row_size()); + this->gg->vl_grid_to_2D(this->gg->get_pvpR_grid(), *pmat, this->gridt->lgd, (is == 0), v_hxc_2d[is].c, setter); } + this->gg->delete_pvpR_grid(); + // clear useless matrices + LR_Util::delete_p3(dm_trans_grid, nsk, gridt->lgd); + LR_Util::delete_p2(rho_trans, nsk); // 5. [AX]^{Hxc}_{ai}=\sum_{\mu,\nu}c^*_{a,\mu,}V^{Hxc}_{\mu,\nu}c_{\nu,i} - // use 2 pzgemms - this->cal_AX_cVc(v_hxc_2d, this->px); - // result is in which "psi" ? X？ - - // final clear - LR_Util::delete_p3(dm_trans_grid); - LR_Util::delete_p2(rho_trans); - return; +#ifdef __MPI + return cal_AX_pblas(LR_Util::mat2ten_double(v_hxc_2d), *this->pmat, *this->psi_ks, *this->pc, naos, nocc, nvirt, *this->pX); +#else + return cal_AX_blas(LR_Util::mat2ten_double(v_hxc_2d), *this->psi_ks, nocc); +#endif } } \ No newline at end of file diff --git a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp index a7662f065f..99bf8f3f4e 100644 --- a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp +++ b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp @@ -28,8 +28,8 @@ namespace elecstate { for (auto k : this->kernel_hxc) delete k; } - - void cal_v_eff(const Charge* chg/*excited state*/, const UnitCell* ucell, ModuleBase::matrix& v_eff) override + void cal_v_eff(const Charge* chg/*excited state*/, const UnitCell* ucell, ModuleBase::matrix& v_eff) override {}; + void cal_v_eff(double** rho, const UnitCell* ucell, ModuleBase::matrix& v_eff) { ModuleBase::TITLE("PotHxcLR", "cal_v_eff"); ModuleBase::timer::tick("PotHxcLR", "cal_v_eff"); @@ -38,7 +38,7 @@ namespace elecstate if (XC_Functional::get_func_type() == 1) if (1 == nspin)// for LDA-spin0, just f*rho for (int ir = 0;ir < nrxx;++ir) - v_eff(0, ir) += (this->kernel_hxc[0]->get_kernel()["Hartree"](0, ir) + this->kernel_hxc[1]->get_kernel()["LDA-v2rho2"](0, ir)) * chg->rho[0][ir]; + v_eff(0, ir) += (this->kernel_hxc[0]->get_kernel()["Hartree"](0, ir) + this->kernel_hxc[1]->get_kernel()["LDA-v2rho2"](0, ir)) * rho[0][ir]; else //remain for spin2, 4 throw std::domain_error("GlobalV::NSPIN =" + std::to_string(GlobalV::NSPIN) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));