diff --git a/source/module_beyonddft/potentials/kernel.h b/source/module_beyonddft/potentials/kernel.h index 1e0f93486c..c3e188bb75 100644 --- a/source/module_beyonddft/potentials/kernel.h +++ b/source/module_beyonddft/potentials/kernel.h @@ -6,27 +6,32 @@ namespace elecstate { - class KernelBase - { - public: - virtual ~KernelBase() = default; - virtual void cal_kernel(const Charge* chg_gs, const UnitCell* ucell, int& nspin) = 0; - virtual std::vector& get_kernel(const std::string& name) { return kernel_set_[name]; } - protected: - const ModulePW::PW_Basis* rho_basis_ = nullptr; - std::map> kernel_set_; // [kernel_type][nrxx][nspin] - }; - - class KernelXC : public KernelBase + class KernelXC { public: KernelXC() {}; ~KernelXC() {}; - void cal_kernel(const Charge* chg_gs/* ground state*/, const UnitCell* ucell, int& nspin) override; + void cal_kernel(const Charge* chg_gs/* ground state*/, const UnitCell* ucell, int& nspin); + const std::vector& get_kernel(const std::string& name) { return kernel_set_[name]; } + const std::vector& get_factor_rho() { return to_mul_rho_; } + const double& get_factor_rho(const int& index) { return to_mul_rho_.at(index); } + const std::vector>& get_factor_drho() { return to_mul_drho_; } + const ModuleBase::Vector3& get_factor_drho(const int& index) { return to_mul_drho_.at(index); } + const std::vector& get_factor_d2rho() { return to_mul_d2rho_; } + const double& get_factor_d2rho(const int& index) { return to_mul_d2rho_.at(index); } + + + protected: // xc kernel for LR-TDDFT void f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const Charge* chg_gs); + + const ModulePW::PW_Basis* rho_basis_ = nullptr; + std::map> kernel_set_; // [kernel_type][nrxx][nspin] + std::vector to_mul_rho_; + std::vector> to_mul_drho_; + std::vector to_mul_d2rho_; }; } diff --git a/source/module_beyonddft/potentials/kernel_xc.cpp b/source/module_beyonddft/potentials/kernel_xc.cpp index b5f7f819e9..b7d7d7e91e 100644 --- a/source/module_beyonddft/potentials/kernel_xc.cpp +++ b/source/module_beyonddft/potentials/kernel_xc.cpp @@ -60,8 +60,13 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons }(); std::vector>> gdr; // \nabla \rho std::vector sigma; // |\nabla\rho|^2 + std::vector sgn; // sgn for threshold mask if (is_gga) { + // 0. set up sgn for threshold mask + // in the case of GGA correlation for polarized case, + // a cutoff for grho is required to ensure that libxc gives reasonable results + // 1. \nabla \rho gdr.resize(nspin); for (int is = 0; is < nspin; ++is) @@ -99,11 +104,11 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons } // ----------------------------------------------------------------------------------- //==================== XC Kernels (f_xc)============================= - //LDA + this->kernel_set_.emplace("vrho", std::vector(nspin * nrxx)); this->kernel_set_.emplace("v2rho2", std::vector(((1 == nspin) ? 1 : 3) * nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11 - //GGA if (is_gga) { + this->kernel_set_.emplace("vsigma", std::vector(((1 == nspin) ? 1 : 3) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 this->kernel_set_.emplace("v2rhosigma", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12 this->kernel_set_.emplace("v2sigma2", std::vector(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22 } @@ -121,13 +126,13 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons switch (func.info->family) { case XC_FAMILY_LDA: - // call Libxc function: xc_lda_exc_vxc - xc_lda_fxc(&func, nrxx, rho.data(), - this->kernel_set_["v2rho2"].data()); + xc_lda_vxc(&func, nrxx, rho.data(), this->kernel_set_["vrho"].data()); + xc_lda_fxc(&func, nrxx, rho.data(), this->kernel_set_["v2rho2"].data()); break; case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: - // call Libxc function: xc_gga_exc_vxc + xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), + this->kernel_set_["vrho"].data(), this->kernel_set_["vsigma"].data()); xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), this->kernel_set_["v2rho2"].data(), this->kernel_set_["v2rhosigma"].data(), this->kernel_set_["v2sigma2"].data()); break; @@ -139,29 +144,70 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons // some formulas for GGA if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - std::vector& v2r2 = this->kernel_set_["v2rho2"]; + this->to_mul_rho_.resize(nrxx); + this->to_mul_drho_.resize(nrxx); + this->to_mul_d2rho_.resize(nrxx); + const std::vector& v2r2 = this->kernel_set_["v2rho2"]; const std::vector& v2rs = this->kernel_set_["v2rhosigma"]; const std::vector& v2s2 = this->kernel_set_["v2sigma2"]; + const std::vector& vs = this->kernel_set_["vsigma"]; const double tpiba2 = tpiba * tpiba; + // method 1 + // if (1 == nspin) + // { + // // ============= to be multiplied by transition density =========== + // // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$ + // std::vector div_v2rhosigma_gdrho_r(nrxx); + // std::vector> v2rhosigma_gdrho_r(nrxx); + // for (int ir = 0; ir < nrxx; ++ir)v2rhosigma_gdrho_r[ir] = gdr[0][ir] * v2rs.at(ir); + // XC_Functional::grad_dot(v2rhosigma_gdrho_r.data(), div_v2rhosigma_gdrho_r.data(), chg_gs->rhopw, tpiba); + // // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$ + // std::vector v2sigma2_sigma_r(nrxx); + // for (int ir = 0; ir < nrxx; ++ir) v2sigma2_sigma_r[ir] = v2s2.at(ir) * sigma[ir]; + // gdr[0].resize(nrxx); + // LR_Util::laplace(v2sigma2_sigma_r.data(), v2sigma2_sigma_r.data(), *(chg_gs->rhopw), tpiba2); + // // 3. $\nabla^2(v^\sigma)$ + // std::vector lap_vsigma(nrxx); + // LR_Util::laplace(vs.data(), lap_vsigma.data(), *(chg_gs->rhopw), tpiba2); + // // add up + // BlasConnector::axpy(nrxx, 1.0, v2r2.data(), 1, to_mul_rho_.data(), 1); + // BlasConnector::axpy(nrxx, -4.0, div_v2rhosigma_gdrho_r.data(), 1, to_mul_rho_.data(), 1); + // BlasConnector::axpy(nrxx, 4.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1); + // BlasConnector::axpy(nrxx, 2.0, lap_vsigma.data(), 1, to_mul_rho_.data(), 1); + + // // ============= to be multiplied by transition density gradient=========== + // // 1. $-2 f^{\rho\sigma}*\nabla\rho$ + // std::vector> v2rhosigma_gdrho(nrxx); + // for (int ir = 0; ir < nrxx; ++ir)v2rhosigma_gdrho[ir] = -2.0 * gdr[0][ir] * v2rs.at(ir); + // // 2. $\nabla(4f^{\sigma\sigma}*sigma+2v^{\sigma})$ + // std::vector v2sigma2_sigma_vsigma(nrxx); + // for (int ir = 0; ir < nrxx; ++ir)v2sigma2_sigma_vsigma[ir] = 4.0 * sigma.at(ir) * v2s2.at(ir) + 2.0 * vs.at(ir); + // std::vector> grad_v2sigma2_sigma_vsigma(nrxx); + // LR_Util::grad(v2sigma2_sigma_vsigma.data(), grad_v2sigma2_sigma_vsigma.data(), *(chg_gs->rhopw), tpiba); + // // add up + // for (int ir = 0; ir < nrxx; ++ir) + // this->to_mul_drho_[ir] = v2rhosigma_gdrho[ir] + grad_v2sigma2_sigma_vsigma[ir]; + // } if (1 == nspin) { - // \div(v2rhosigma*gdrho) + // ============= to be multiplied by transition rho =========== + // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$ std::vector div_v2rhosigma_gdrho_r(nrxx); std::vector> v2rhosigma_gdrho_r(nrxx); for (int ir = 0; ir < nrxx; ++ir)v2rhosigma_gdrho_r[ir] = gdr[0][ir] * v2rs.at(ir); XC_Functional::grad_dot(v2rhosigma_gdrho_r.data(), div_v2rhosigma_gdrho_r.data(), chg_gs->rhopw, tpiba); - // \lap(v2sigma2*sigma) - std::vector v2sigma2_sigma_r(nrxx); - for (int ir = 0; ir < nrxx; ++ir) v2sigma2_sigma_r[ir] = v2s2.at(ir) * sigma[ir]; - gdr[0].resize(nrxx); - LR_Util::laplace(v2sigma2_sigma_r.data(), v2sigma2_sigma_r.data(), *(chg_gs->rhopw), tpiba2); - // add to v2rho2 - BlasConnector::axpy(nrxx, -4.0, div_v2rhosigma_gdrho_r.data(), 1, v2r2.data(), 1); - BlasConnector::axpy(nrxx, 4.0, v2sigma2_sigma_r.data(), 1, v2r2.data(), 1); + BlasConnector::axpy(nrxx, 1.0, v2r2.data(), 1, to_mul_rho_.data(), 1); + // BlasConnector::axpy(nrxx, -2.0, div_v2rhosigma_gdrho_r.data(), 1, to_mul_rho_.data(), 1); + + // ============= to be multiplied by transition drho and d2rho=========== + // 1. $\nabla(4f^{\sigma\sigma}*\sigma+2v^{\sigma})$ and its gradient + for (int ir = 0; ir < nrxx; ++ir) to_mul_d2rho_[ir] = -(4.0 * sigma.at(ir) * v2s2.at(ir) + 2.0 * vs.at(ir)); + LR_Util::grad(to_mul_d2rho_.data(), to_mul_drho_.data(), *(chg_gs->rhopw), tpiba); + } - else if (2 == nspin) + else if (2 == nspin) // wrong, to be fixed { - //\div(v2rhosigma*gdrho) + // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$ std::vector div_v2rhosigma_gdrho_r(3 * nrxx); std::vector> v2rhosigma_gdrho_r(3 * nrxx); for (int ir = 0; ir < nrxx; ++ir) @@ -175,7 +221,7 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons } for (int isig = 0;isig < 3;++isig) XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs->rhopw, tpiba); - // \lap(v2sigma2*sigma) + // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$ std::vector v2sigma2_sigma_r(3 * nrxx); for (int ir = 0; ir < nrxx; ++ir) { @@ -191,10 +237,21 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons } for (int isig = 0;isig < 3;++isig) LR_Util::laplace(v2sigma2_sigma_r.data() + isig * nrxx, v2sigma2_sigma_r.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2); + // 3. $\nabla^2(v^\sigma)$ + std::vector lap_vsigma(3 * nrxx); + for (int ir = 0;ir < nrxx;++ir) + { + lap_vsigma[ir] = vs.at(ir * 3) * 2.0; + lap_vsigma[nrxx + ir] = vs.at(ir * 3 + 1); + lap_vsigma[2 * nrxx + ir] = vs.at(ir * 3 + 2) * 2.0; + } + for (int isig = 0;isig < 3;++isig) + LR_Util::laplace(lap_vsigma.data() + isig * nrxx, lap_vsigma.data() + isig * nrxx, *(chg_gs->rhopw), tpiba2); // add to v2rho2 - BlasConnector::axpy(3 * nrxx, -1.0, div_v2rhosigma_gdrho_r.data(), 1, v2r2.data(), 1); - BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, v2r2.data(), 1); - + BlasConnector::axpy(3 * nrxx, 1.0, v2r2.data(), 1, to_mul_rho_.data(), 1); + BlasConnector::axpy(3 * nrxx, -1.0, div_v2rhosigma_gdrho_r.data(), 1, to_mul_rho_.data(), 1); + BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1); + BlasConnector::axpy(3 * nrxx, 1.0, lap_vsigma.data(), 1, to_mul_rho_.data(), 1); } } } // end for( xc_func_type &func : funcs ) diff --git a/source/module_beyonddft/potentials/pot_hxc_lrtd.cpp b/source/module_beyonddft/potentials/pot_hxc_lrtd.cpp index b3157f5d73..486866269c 100644 --- a/source/module_beyonddft/potentials/pot_hxc_lrtd.cpp +++ b/source/module_beyonddft/potentials/pot_hxc_lrtd.cpp @@ -4,6 +4,7 @@ #include "module_base/timer.h" #include "module_hamilt_general/module_xc/xc_functional.h" #include +#include "module_beyonddft/utils/lr_util.h" namespace elecstate { // constructor for exchange-correlation kernel @@ -26,8 +27,7 @@ namespace elecstate if (local_xc.find(this->xc_kernel) != local_xc.end()) { XC_Functional::set_xc_type(this->xc_kernel); - this->xc_kernel_components_ = new KernelXC(); - this->xc_kernel_components_->cal_kernel(chg_gs, ucell_in, this->nspin); + this->xc_kernel_components_.cal_kernel(chg_gs, ucell_in, this->nspin); } } @@ -35,28 +35,59 @@ namespace elecstate { ModuleBase::TITLE("PotHxcLR", "cal_v_eff"); ModuleBase::timer::tick("PotHxcLR", "cal_v_eff"); +#ifdef USE_LIBXC 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" || xc_kernel == "hf") return; - else if (XC_Functional::get_func_type() == 1 || XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4)//LDA or GGA or HYBGGA + else if (XC_Functional::get_func_type() == 1)//LDA if (1 == nspin)// for LDA-spin0, just fxc*rho where fxc=v2rho2; for GGA, v2rho2 has been replaced by the true fxc for (int ir = 0;ir < nrxx;++ir) - v_eff(0, ir) += this->xc_kernel_components_->get_kernel("v2rho2").at(ir) * rho[0][ir]; + v_eff(0, ir) += ModuleBase::e2 * this->xc_kernel_components_.get_kernel("v2rho2").at(ir) * rho[0][ir]; else if (2 == nspin) for (int ir = 0;ir < nrxx;++ir) { const int irs0 = 2 * ir; const int irs1 = irs0 + 1; const int irs2 = irs0 + 2; - v_eff(0, ir) += this->xc_kernel_components_->get_kernel("v2rho2").at(irs0) * rho[0][ir] - + this->xc_kernel_components_->get_kernel("v2rho2").at(irs1) * rho[1][ir]; - v_eff(1, ir) += this->xc_kernel_components_->get_kernel("v2rho2").at(irs1) * rho[0][ir] - + this->xc_kernel_components_->get_kernel("v2rho2").at(irs2) * rho[1][ir]; + v_eff(0, ir) += ModuleBase::e2 * this->xc_kernel_components_.get_kernel("v2rho2").at(irs0) * rho[0][ir] + + this->xc_kernel_components_.get_kernel("v2rho2").at(irs1) * rho[1][ir]; + v_eff(1, ir) += ModuleBase::e2 * this->xc_kernel_components_.get_kernel("v2rho2").at(irs1) * rho[0][ir] + + this->xc_kernel_components_.get_kernel("v2rho2").at(irs2) * rho[1][ir]; } else //remain for spin 4 throw std::domain_error("nspin =" + std::to_string(nspin) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + else if (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4) // GGA or HYB_GGA + { + if (1 == nspin) + { + std::vector> drho(nrxx); + LR_Util::grad(rho[0], drho.data(), *(this->rho_basis_), ucell->tpiba); + std::vector d2rho(nrxx); + LR_Util::laplace(rho[0], d2rho.data(), *(this->rho_basis_), ucell->tpiba); + for (int ir = 0;ir < nrxx;++ir) + v_eff(0, ir) += ModuleBase::e2 * + (this->xc_kernel_components_.get_factor_rho(ir) * rho[0][ir]); + // + this->xc_kernel_components_.get_factor_drho(ir) * drho.at(ir)); + // + this->xc_kernel_components_.get_factor_d2rho(ir) * d2rho.at(ir)); + } + else if (2 == nspin) // wrong, to be fixed + for (int ir = 0;ir < nrxx;++ir) + { + const int irs0 = 2 * ir; + const int irs1 = irs0 + 1; + const int irs2 = irs0 + 2; + v_eff(0, ir) += ModuleBase::e2 * this->xc_kernel_components_.get_kernel("v2rho2").at(irs0) * rho[0][ir] + + this->xc_kernel_components_.get_kernel("v2rho2").at(irs1) * rho[1][ir]; + v_eff(1, ir) += ModuleBase::e2 * this->xc_kernel_components_.get_kernel("v2rho2").at(irs1) * rho[0][ir] + + this->xc_kernel_components_.get_kernel("v2rho2").at(irs2) * rho[1][ir]; + } + else //remain for spin 4 + throw std::domain_error("nspin =" + std::to_string(nspin) + + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + } else +#endif throw std::domain_error("GlobalV::XC_Functional::get_func_type() =" + std::to_string(XC_Functional::get_func_type()) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); diff --git a/source/module_beyonddft/potentials/pot_hxc_lrtd.h b/source/module_beyonddft/potentials/pot_hxc_lrtd.h index d3132727a4..15148e6ac4 100644 --- a/source/module_beyonddft/potentials/pot_hxc_lrtd.h +++ b/source/module_beyonddft/potentials/pot_hxc_lrtd.h @@ -10,11 +10,7 @@ namespace elecstate public: // constructor for exchange-correlation kernel PotHxcLR(const std::string& xc_kernel_in, const ModulePW::PW_Basis* rho_basis_in, const UnitCell* ucell_in, const Charge* chg_gs/*ground state*/); - ~PotHxcLR() - { - if (this->xc_kernel == "lda" /*|| pbe... */) - delete this->xc_kernel_components_; - } + ~PotHxcLR() {} 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); int nrxx; @@ -24,7 +20,7 @@ namespace elecstate /// LDA: v2rho2 /// GGA: v2rho2, v2rhosigma, v2sigma2 /// meta-GGA: v2rho2, v2rhosigma, v2sigma2, v2rholap, v2rhotau, v2sigmalap, v2sigmatau, v2laptau, v2lap2, v2tau2 - KernelBase* xc_kernel_components_; + KernelXC xc_kernel_components_; const std::string xc_kernel; };