From f1bfb8bc9573aaf5b3baf736cf37f3225d84ba73 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 15 Sep 2023 01:04:41 +0800 Subject: [PATCH] param xc_kernel & fix gint-call and Kernels --- docs/advanced/input_files/input-main.md | 6 +++ source/module_beyonddft/esolver_lrtd_lcao.h | 10 ++-- source/module_beyonddft/esolver_lrtd_lcao.hpp | 39 +++++++-------- source/module_beyonddft/hamilt_casida.hpp | 3 +- .../operator_casida/operatorA_hxc.h | 10 ++-- .../operator_casida/operatorA_hxc.hpp | 27 +++++----- .../module_beyonddft/potentials/kernel_base.h | 3 +- .../potentials/kernel_hartree.hpp | 14 +++--- .../module_beyonddft/potentials/kernel_xc.hpp | 49 ++++++++++++------- .../potentials/pot_hxc_lrtd.hpp | 2 +- source/module_hamilt_general/operator.h | 1 + source/module_io/input.cpp | 7 ++- source/module_io/input.h | 1 + source/module_io/test/input_test.cpp | 1 + source/module_io/test/input_test_para.cpp | 1 + source/module_io/test/write_input_test.cpp | 1 + source/module_io/write_input.cpp | 1 + 17 files changed, 100 insertions(+), 76 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 294b127c20..45007f2095 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3146,4 +3146,10 @@ These parameters are used to solve the excited states using. e.g. lr-tddft - **Description**: The number of 2-particle states to be solved - **Default**: 0 +### xc_kernel + +- **Type**: String +- **Description**: The exchange-correlation kernel used in the calculation. Currently, only `LDA` is supported. +- **Default**: LDA + [back to top](#full-list-of-input-keywords) \ No newline at end of file diff --git a/source/module_beyonddft/esolver_lrtd_lcao.h b/source/module_beyonddft/esolver_lrtd_lcao.h index 37a4a43266..de60eebc0f 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.h +++ b/source/module_beyonddft/esolver_lrtd_lcao.h @@ -34,6 +34,8 @@ namespace ModuleESolver delete this->phsol; delete this->pot; delete this->psi_ks; + delete this->X; + delete this->AX; } ///input: input, call, basis(LCAO), psi(ground state), elecstate @@ -66,9 +68,9 @@ namespace ModuleESolver ModuleBase::matrix eig_ks; // energy of ground state is in pelec->ekb - /// @brief Excited state info. size: nstates*nks*nocc*nvirt - std::vector> X; - std::vector> AX; + /// @brief Excited state info. size: nstates * nks * (nocc(local) * nvirt (local)) + psi::Psi* X; + psi::Psi* AX; size_t nocc; size_t nvirt; @@ -85,7 +87,7 @@ namespace ModuleESolver // adj info Record_adj ra; // grid parallel info (no need for 2d-block distribution)? - Grid_Technique gridt; + // Grid_Technique gridt; // grid integration method(will be moved to OperatorKernelHxc) ModulePW::PW_Basis_Big bigpw; Gint_Gamma gint_g; diff --git a/source/module_beyonddft/esolver_lrtd_lcao.hpp b/source/module_beyonddft/esolver_lrtd_lcao.hpp index e7ad3442df..4bb4a3bc37 100644 --- a/source/module_beyonddft/esolver_lrtd_lcao.hpp +++ b/source/module_beyonddft/esolver_lrtd_lcao.hpp @@ -41,11 +41,9 @@ ModuleESolver::ESolver_LRTD::ESolver_LRTD(ModuleESolver::ESolver this->gint_g = std::move(ks_sol.UHM.GG); else this->gint_k = std::move(ks_sol.UHM.GK); - // can GlobalC be moved? try - - // grid parallel info - this->gridt = std::move(ks_sol.GridT); + // xc kernel + XC_Functional::set_xc_type(inp.xc_kernel); //init potential and calculate kernels using ground state charge if (this->pot == nullptr) { @@ -75,22 +73,22 @@ ModuleESolver::ESolver_LRTD::ESolver_LRTD(ModuleESolver::ESolver //init Hamiltonian if (typeid(FPTYPE) == typeid(double)) this->p_hamilt = new hamilt::HamiltCasidaLR(this->nks, this->nbasis, this->nocc, this->nvirt, this->psi_ks, - &this->gint_g, this->pot, &this->gridt, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); + &this->gint_g, this->pot, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); else if (typeid(FPTYPE) == typeid(std::complex)) this->p_hamilt = new hamilt::HamiltCasidaLR(this->nks, this->nbasis, this->nocc, this->nvirt, this->psi_ks, - &this->gint_k, this->pot, &this->gridt, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); + &this->gint_k, this->pot, std::vector({ &this->paraX_, &this->paraC_, &this->paraMat_ })); // init HSolver // this->phsol = new hsolver::HSolver /// =========just for test============== - // 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->paraX_.get_row_size() - 1), this->AX[istate].get_pointer()); - this->AX[0] = this->p_hamilt->ops->act(this->X[0]); - + // try act + for (int istate = 0;istate < nstates;++istate) + { + this->X->fix_b(istate); + this->AX->fix_b(istate); + this->p_hamilt->ops->act(*this->X, *this->AX); + } } template @@ -112,15 +110,12 @@ void ModuleESolver::ESolver_LRTD::init_X() int nsk = (nspin == 4) ? this->nks : this->nks * this->nspin; LR_Util::setup_2d_division(this->paraX_, 1, this->nvirt, this->nocc);//nvirt - row, nocc - col - for (int i = 0; i < this->nstates; i++) - { - this->X.emplace_back(nsk, this->paraX_.get_col_size(), this->paraX_.get_row_size()); - X[i].zero_out(); - } + this->X = new psi::Psi(nsk, this->nstates, this->paraX_.get_local_size(), nullptr, false); // band(state)-first + this->X->zero_out(); LR_Util::setup_2d_division(this->paraMat_, 1, this->nbasis, this->nbasis, this->paraX_.comm_2D, this->paraX_.blacs_ctxt); LR_Util::setup_2d_division(this->paraC_, 1, this->nbasis, this->nocc + this->nvirt, this->paraX_.comm_2D, this->paraX_.blacs_ctxt); - this->AX = this->X; - + // this->AX = new psi::Psi(*this->X); + this->AX = new psi::Psi(nsk, this->nstates, this->paraX_.get_local_size(), nullptr, false); // set the initial guess of X // if (E_{lumo}-E_{homo-1} < E_{lumo+1}-E{homo}), mode = 0, else 1(smaller first) bool ix_mode = 0; //default @@ -143,9 +138,11 @@ void ModuleESolver::ESolver_LRTD::init_X() // use unit vectors as the initial guess for (int i = 0; i < this->nstates; i++) { + this->X->fix_b(i); int row_global = std::get<0>(ix2iciv[i]); int col_global = std::get<1>(ix2iciv[i]); if (this->paraX_.in_this_processor(row_global, col_global)) - X[i](this->paraX_.global2local_row(row_global), this->paraX_.global2local_col(col_global)) = static_cast(1.0); + for (int isk = 0;isk < nsk;++isk) + (*X)(isk, this->paraX_.global2local_row(row_global) * this->paraX_.get_col_size() + this->paraX_.global2local_col(col_global)) = static_cast(1.0); } } diff --git a/source/module_beyonddft/hamilt_casida.hpp b/source/module_beyonddft/hamilt_casida.hpp index 32ff49b916..9c1137aced 100644 --- a/source/module_beyonddft/hamilt_casida.hpp +++ b/source/module_beyonddft/hamilt_casida.hpp @@ -15,7 +15,6 @@ namespace hamilt const psi::Psi* psi_ks_in, TGint* gint_in, elecstate::PotHxcLR* pot_in, - const Grid_Technique* gt, const std::vector p2d_in) { ModuleBase::TITLE("HamiltCasidaLR", "HamiltCasidaLR"); @@ -23,7 +22,7 @@ namespace hamilt // ops and opsd in base class may be unified in the future? //add Hxc operator (the first one) - this->ops = new OperatorA_Hxc(nsk, naos, nocc, nvirt, psi_ks_in, gint_in, pot_in, gt, p2d_in); + this->ops = new OperatorA_Hxc(nsk, naos, nocc, nvirt, psi_ks_in, gint_in, pot_in, p2d_in); //add Exx operator (remaining) // Operator* A_Exx = new OperatorA_Exx; // this->opsd->add(A_Exx); diff --git a/source/module_beyonddft/operator_casida/operatorA_hxc.h b/source/module_beyonddft/operator_casida/operatorA_hxc.h index 4919252e52..b25f325f50 100644 --- a/source/module_beyonddft/operator_casida/operatorA_hxc.h +++ b/source/module_beyonddft/operator_casida/operatorA_hxc.h @@ -20,10 +20,9 @@ namespace hamilt const psi::Psi* psi_ks_in, Gint_Gamma* gg_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), - psi_ks(psi_ks_in), gg(gg_in), pot(pot_in), gridt(gt_in), + psi_ks(psi_ks_in), gg(gg_in), pot(pot_in), pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) { ModuleBase::TITLE("OperatorA_Hxc", "OperatorA_Hxc(gamma)"); @@ -37,10 +36,9 @@ namespace hamilt const psi::Psi* psi_ks_in, Gint_k* gk_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), - psi_ks(psi_ks_in), gk(gk_in), pot(pot_in), gridt(gt_in), + psi_ks(psi_ks_in), gk(gk_in), pot(pot_in), pX(p2d_in.at(0)), pc(p2d_in.at(1)), pmat(p2d_in.at(2)) { ModuleBase::TITLE("OperatorA_Hxc", "OperatorA_Hxc(k)"); @@ -56,7 +54,8 @@ namespace hamilt 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; + // virtual psi::Psi act(const psi::Psi& psi_in) const override; + virtual void act(const psi::Psi& psi_in, psi::Psi& psi_out) const override; private: //global sizes int nsk; //nspin*nkpoints @@ -70,7 +69,6 @@ namespace hamilt Parallel_2D* pc = nullptr; Parallel_2D* pX = nullptr; Parallel_2D* pmat = nullptr; - const Grid_Technique* gridt = nullptr; elecstate::PotHxcLR* pot = nullptr; diff --git a/source/module_beyonddft/operator_casida/operatorA_hxc.hpp b/source/module_beyonddft/operator_casida/operatorA_hxc.hpp index dc9da1ba72..93587a36c5 100644 --- a/source/module_beyonddft/operator_casida/operatorA_hxc.hpp +++ b/source/module_beyonddft/operator_casida/operatorA_hxc.hpp @@ -10,8 +10,11 @@ namespace hamilt { // for double template - psi::Psi OperatorA_Hxc::act(const psi::Psi& psi_in) const + // psi::Psi OperatorA_Hxc::act(const psi::Psi& psi_in) const + void OperatorA_Hxc::act(const psi::Psi& psi_in, psi::Psi& psi_out) const { + ModuleBase::TITLE("OperatorA_Hxc", "act"); + const int& lgd = gg->gridt->lgd; // gamma-only now // 1. transition density matrix // multi-k needs k-to-R FT @@ -21,13 +24,13 @@ namespace hamilt std::vector dm_trans_2d = cal_dm_trans_blas(*pX, *pc, nocc, nvirt); #endif double*** dm_trans_grid; - LR_Util::new_p3(dm_trans_grid, nsk, gridt->lgd, gridt->lgd); + LR_Util::new_p3(dm_trans_grid, nsk, lgd, 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); + dm2g.setAlltoallvParameter(pmat->comm_2D, naos, pmat->blacs_ctxt, pmat->nb, lgd, gg->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); + dm2g.cal_dk_gamma_from_2D(LR_Util::ten2mat_double(dm_trans_2d), dm_trans_grid, nsk, naos, lgd, GlobalV::ofs_running); // 2. transition electron density double** rho_trans; @@ -44,7 +47,7 @@ namespace hamilt // results are stored in gg->pvpR_grid(gamma_only) // or gint_k->pvpR_reduced(multi_k) std::vector v_hxc_2d(nsk); - this->gg->init_pvpR_grid(gridt->lgd); + this->gg->init_pvpR_grid(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); @@ -52,28 +55,24 @@ namespace hamilt }; for (int is = 0;is < this->nsk;++is) { - ModuleBase::GlobalFunc::ZEROS(this->gg->get_pvpR_grid(), gridt->lgd * gridt->lgd); + ModuleBase::GlobalFunc::ZEROS(this->gg->get_pvpR_grid(), lgd * 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->vl_grid_to_2D(this->gg->get_pvpR_grid(), *pmat, 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_p3(dm_trans_grid, nsk, 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} #ifdef __MPI - psi::Psi AX(1, nsk, this->pX->get_local_size(), nullptr, false); - cal_AX_pblas(LR_Util::mat2ten_double(v_hxc_2d), *this->pmat, *this->psi_ks, *this->pc, naos, nocc, nvirt, *this->pX, AX); - return AX; + cal_AX_pblas(LR_Util::mat2ten_double(v_hxc_2d), *this->pmat, *this->psi_ks, *this->pc, naos, nocc, nvirt, *this->pX, psi_out); #else - psi::Psi AX(1, nsk, nocc * nvirt, nullptr, false); - cal_AX_blas(LR_Util::mat2ten_double(v_hxc_2d), *this->psi_ks, nocc, nvirt, AX); - return AX; + cal_AX_blas(LR_Util::mat2ten_double(v_hxc_2d), *this->psi_ks, nocc, nvirt, psi_out); #endif } } \ No newline at end of file diff --git a/source/module_beyonddft/potentials/kernel_base.h b/source/module_beyonddft/potentials/kernel_base.h index 002d8febce..ca928cc000 100644 --- a/source/module_beyonddft/potentials/kernel_base.h +++ b/source/module_beyonddft/potentials/kernel_base.h @@ -2,12 +2,13 @@ #include "module_basis/module_pw/pw_basis.h" #include "module_elecstate/module_charge/charge.h" #include "module_cell/unitcell.h" +// #include class KernelBase { public: virtual ~KernelBase() = default; virtual void cal_kernel(const Charge* chg_gs, const UnitCell* ucell, int& nspin) = 0; - virtual std::map& get_kernel() = 0; + virtual ModuleBase::matrix& get_kernel(const std::string& name) { return kernel_set_[name]; } protected: const ModulePW::PW_Basis* rho_basis_ = nullptr; std::map kernel_set_; // [kernel_type][nspin][nrxx] diff --git a/source/module_beyonddft/potentials/kernel_hartree.hpp b/source/module_beyonddft/potentials/kernel_hartree.hpp index d406ae7125..d69cfc2eca 100644 --- a/source/module_beyonddft/potentials/kernel_hartree.hpp +++ b/source/module_beyonddft/potentials/kernel_hartree.hpp @@ -12,15 +12,13 @@ namespace elecstate { this->rho_basis_ = rho_basis_in; } - std::map& get_kernel() override { return this->kernel_set_; } void cal_kernel(const Charge* chg_gs, const UnitCell* ucell, int& nspin) override { ModuleBase::TITLE("KernelHartree", "cal_v_eff"); ModuleBase::timer::tick("KernelHartree", "cal_v_eff"); - ModuleBase::matrix f_hartree(chg_gs->nrxx, nspin); - - //1. Coulomb kernel in reciprocal space + this->kernel_set_.emplace("Hartree", ModuleBase::matrix(nspin, chg_gs->nrxx)); + //1. Hartree kernel in reciprocal space std::vector> Porter(this->rho_basis_->nmaxgr, std::complex(0.0, 0.0)); #ifdef _OPENMP #pragma omp parallel for @@ -32,14 +30,14 @@ namespace elecstate //2. FFT to real space rho_basis_->recip2real(Porter.data(), Porter.data()); - //3. Add to f_hartree + //3. Add to kernel_set_ if (nspin == 4) { #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif for (int ir = 0;ir < this->rho_basis_->nrxx;ir++) - f_hartree(ir, 0) += Porter[ir].real(); + kernel_set_["Hartree"](0, ir) += Porter[ir].real(); } else { @@ -48,9 +46,9 @@ namespace elecstate #endif for (int is = 0;is < nspin;is++) for (int ir = 0;ir < this->rho_basis_->nrxx;ir++) - f_hartree(ir, is) += Porter[ir].real(); + kernel_set_["Hartree"](is, ir) += Porter[ir].real(); } - this->kernel_set_.insert({ "Hartree", f_hartree }); + ModuleBase::timer::tick("KernelHartree", "cal_v_eff"); return; } diff --git a/source/module_beyonddft/potentials/kernel_xc.hpp b/source/module_beyonddft/potentials/kernel_xc.hpp index 896e7b3d32..a5629530e2 100644 --- a/source/module_beyonddft/potentials/kernel_xc.hpp +++ b/source/module_beyonddft/potentials/kernel_xc.hpp @@ -10,7 +10,6 @@ namespace elecstate public: KernelXC() {}; ~KernelXC() {}; - std::map& get_kernel() override { return this->kernel_set_; } void cal_kernel(const Charge* chg_gs/* ground state*/, const UnitCell* ucell, int& nspin) override; @@ -23,8 +22,8 @@ namespace elecstate #include void KernelXC::f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const Charge* chg_gs) { - ModuleBase::TITLE("XC_Functional", "v_xc_libxc"); - ModuleBase::timer::tick("XC_Functional", "v_xc_libxc"); + ModuleBase::TITLE("XC_Functional", "f_xc_libxc"); + ModuleBase::timer::tick("XC_Functional", "f_xc_libxc"); // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ std::vector funcs = XC_Functional::init_func((1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); @@ -32,7 +31,7 @@ namespace elecstate // converting rho (extract it as a subfuntion in the future) // ----------------------------------------------------------------------------------- - ModuleBase::matrix rho(nrxx, nspin); //spin-major + std::vector rho(nspin * nrxx); //spin-major ModuleBase::matrix sigma; std::vector amag; if (1 == nspin || 2 == nspin) @@ -42,7 +41,7 @@ namespace elecstate #endif for (int is = 0; is < nspin; ++is) for (int ir = 0; ir < nrxx; ++ir) - rho(ir, is) = chg_gs->rho[is][ir] + 1.0 / nspin * chg_gs->rho_core[ir]; + rho[ir * nspin + is] = chg_gs->rho[is][ir] + 1.0 / nspin * chg_gs->rho_core[ir]; } else { @@ -55,8 +54,8 @@ namespace elecstate const double arhox = std::abs(chg_gs->rho[0][ir] + chg_gs->rho_core[ir]); amag[ir] = std::sqrt(std::pow(chg_gs->rho[1][ir], 2) + std::pow(chg_gs->rho[2][ir], 2) + std::pow(chg_gs->rho[3][ir], 2)); const double amag_clip = (amag[ir] < arhox) ? amag[ir] : arhox; - rho(ir, 0) = (arhox + amag_clip) / 2.0; - rho(ir, 1) = (arhox - amag_clip) / 2.0; + rho[ir * nspin + 0] = (arhox + amag_clip) / 2.0; + rho[ir * nspin + 1] = (arhox - amag_clip) / 2.0; } } @@ -67,10 +66,28 @@ namespace elecstate // ----------------------------------------------------------------------------------- //==================== XC Kernels (f_xc)============================= //LDA - ModuleBase::matrix v2rho2; //(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 + this->kernel_set_.emplace("v2rho2", ModuleBase::matrix(((1 == nspin) ? 1 : 3), nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 //GGA - ModuleBase::matrix v2rhosigma; //(nrxx* ((1 == nspin) ? 1 : 6)): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 - ModuleBase::matrix v2sigma2; //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 + const bool is_gga = [&funcs]() + { + for (xc_func_type& func : funcs) + { + switch (func.info->family) + { + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + return true; + } + } + return false; + }(); + + if (is_gga) + { + this->kernel_set_.emplace("v2rhosigma", ModuleBase::matrix(((1 == nspin) ? 1 : 6), nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 + this->kernel_set_.emplace("v2sigma2", ModuleBase::matrix(((1 == nspin) ? 1 : 6), nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 + } + //MetaGGA // ModuleBase::matrix v2rholapl; // ModuleBase::matrix v2rhotau; @@ -92,18 +109,14 @@ namespace elecstate { case XC_FAMILY_LDA: // call Libxc function: xc_lda_exc_vxc - xc_lda_fxc(&func, nrxx, rho.c, - v2rho2.c); - this->kernel_set_.insert({ "LDA-v2rho2", v2rho2 }); + xc_lda_fxc(&func, nrxx, rho.data(), + this->kernel_set_["v2rho2"].c); break; case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: // call Libxc function: xc_gga_exc_vxc - xc_gga_fxc(&func, nrxx, rho.c, sigma.c, - v2rho2.c, v2rhosigma.c, v2sigma2.c); - this->kernel_set_.insert({ "GGA-v2rho2", v2rho2 }); - this->kernel_set_.insert({ "GGA-v2rhosigma", v2rhosigma }); - this->kernel_set_.insert({ "GGA-v2sigma2", v2sigma2 }); + xc_gga_fxc(&func, nrxx, rho.data(), sigma.c, + this->kernel_set_["v2rho2"].c, this->kernel_set_["v2rhosigma"].c, this->kernel_set_["v2sigma2"].c); break; default: throw std::domain_error("func.info->family =" + std::to_string(func.info->family) diff --git a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp index 99bf8f3f4e..3773beace8 100644 --- a/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp +++ b/source/module_beyonddft/potentials/pot_hxc_lrtd.hpp @@ -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)) * rho[0][ir]; + v_eff(0, ir) += (this->kernel_hxc[0]->get_kernel("Hartree")(0, ir) + this->kernel_hxc[1]->get_kernel("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__)); diff --git a/source/module_hamilt_general/operator.h b/source/module_hamilt_general/operator.h index 187629907c..1a76ad46e8 100644 --- a/source/module_hamilt_general/operator.h +++ b/source/module_hamilt_general/operator.h @@ -62,6 +62,7 @@ class Operator /// an developer-friendly interface for act() function virtual psi::Psi act(const psi::Psi& psi_in) const { return psi_in; }; + virtual void act(const psi::Psi& psi_in, psi::Psi& psi_out) const {}; Operator* next_op = nullptr; diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 47c3bb0543..4040842222 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -593,6 +593,7 @@ void Input::Default(void) //========================================================== beyonddft_method = "none"; nstates = 0; + xc_kernel = "LDA"; return; } @@ -2136,7 +2137,11 @@ bool Input::Read(const std::string &fn) { read_value(ifs, nstates); } - //---------------------------------------------------------------------------------- + else if (strcmp("xc_kernel", word) == 0) + { + read_value(ifs, xc_kernel); + } + //---------------------------------------------------------------------------------- else { // xiaohui add 2015-09-15 diff --git a/source/module_io/input.h b/source/module_io/input.h index 95868731be..fd619328b2 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -562,6 +562,7 @@ class Input std::string beyonddft_method; // the method for solving excited state. e.g. lr-tddft int n_unocc; // the number of unoccupied bands to form the 2-particle states int nstates; // the number of 2-particle states to be solved + std::string xc_kernel; // xc kernel for LR-TDDFT private: //========================================================== diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 22bafab789..2920d3603f 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -362,6 +362,7 @@ TEST_F(InputTest, Default) EXPECT_TRUE(INPUT.mdp.dump_virial); EXPECT_EQ(INPUT.beyonddft_method, "none"); EXPECT_EQ(INPUT.nstates, 0); + EXPECT_EQ(INPUT.xc_kernel, "LDA"); } TEST_F(InputTest, Read) diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index c28fddc1c6..3fd400d68c 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -370,6 +370,7 @@ TEST_F(InputParaTest,Bcast) EXPECT_EQ(INPUT.beyonddft_method, "none"); EXPECT_EQ(INPUT.nstates, 0); + EXPECT_EQ(INPUT.xc_kernel, "LDA"); } } diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index b95ab8e378..6ff523186e 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -358,6 +358,7 @@ TEST_F(write_input,print) EXPECT_THAT(output, testing::HasSubstr("bessel_descriptor_sigma 0.1 #spherical bessel smearing_sigma")); EXPECT_THAT(output, testing::HasSubstr("beyonddft_method none #the method for solving excited state. e.g. lr-tddft")); EXPECT_THAT(output, testing::HasSubstr("nstates 0 #the number of 2-particle states to be solved")); + EXPECT_THAT(output, testing::HasSubstr("xc_kernel LDA #xc kernel for LR-TDDFT. default: LDA")); ifs.close(); remove("write_input_test.log"); } diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 68bde7653f..e94962a7c6 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -462,6 +462,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ofs << "\n#Parameters (22.beyond dft)" << std::endl; ModuleBase::GlobalFunc::OUTP(ofs, "beyonddft_method", beyonddft_method, "the method for solving excited state. e.g. lr-tddft"); ModuleBase::GlobalFunc::OUTP(ofs, "nstates", nstates, "the number of 2-particle states to be solved"); + ModuleBase::GlobalFunc::OUTP(ofs, "xc_kernel", xc_kernel, "xc kernel for LR-TDDFT. default: LDA"); ofs.close(); return; }