Skip to content

Commit

Permalink
refactor XC-kernel and fix the missing e2 (GGA still bug)
Browse files Browse the repository at this point in the history
  • Loading branch information
@maki49
maki49 committed May 26, 2024
1 parent 354af82 commit 5fea1e3
Show file tree
Hide file tree
Showing 4 changed files with 138 additions and 49 deletions.
31 changes: 18 additions & 13 deletions source/module_beyonddft/potentials/kernel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& get_kernel(const std::string& name) { return kernel_set_[name]; }
protected:
const ModulePW::PW_Basis* rho_basis_ = nullptr;
std::map<std::string, std::vector<double>> 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<double>& get_kernel(const std::string& name) { return kernel_set_[name]; }
const std::vector<double>& get_factor_rho() { return to_mul_rho_; }
const double& get_factor_rho(const int& index) { return to_mul_rho_.at(index); }
const std::vector<ModuleBase::Vector3<double>>& get_factor_drho() { return to_mul_drho_; }
const ModuleBase::Vector3<double>& get_factor_drho(const int& index) { return to_mul_drho_.at(index); }
const std::vector<double>& 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<std::string, std::vector<double>> kernel_set_; // [kernel_type][nrxx][nspin]
std::vector<double> to_mul_rho_;
std::vector<ModuleBase::Vector3<double>> to_mul_drho_;
std::vector<double> to_mul_d2rho_;
};
}

101 changes: 79 additions & 22 deletions source/module_beyonddft/potentials/kernel_xc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,13 @@ void elecstate::KernelXC::f_xc_libxc(const int& nspin, const double& omega, cons
}();
std::vector<std::vector<ModuleBase::Vector3<double>>> gdr; // \nabla \rho
std::vector<double> sigma; // |\nabla\rho|^2
std::vector<double> 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)
Expand Down Expand Up @@ -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<double>(nspin * nrxx));
this->kernel_set_.emplace("v2rho2", std::vector<double>(((1 == nspin) ? 1 : 3) * nrxx));//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11
//GGA
if (is_gga)
{
this->kernel_set_.emplace("vsigma", std::vector<double>(((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<double>(((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<double>(((1 == nspin) ? 1 : 6) * nrxx)); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22
}
Expand All @@ -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;
Expand All @@ -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<double>& 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<double>& v2r2 = this->kernel_set_["v2rho2"];
const std::vector<double>& v2rs = this->kernel_set_["v2rhosigma"];
const std::vector<double>& v2s2 = this->kernel_set_["v2sigma2"];
const std::vector<double>& 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<double> div_v2rhosigma_gdrho_r(nrxx);
// std::vector<ModuleBase::Vector3<double>> 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<double> 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<double> 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<ModuleBase::Vector3<double>> 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<double> 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<ModuleBase::Vector3<double>> 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<double> div_v2rhosigma_gdrho_r(nrxx);
std::vector<ModuleBase::Vector3<double>> 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<double> 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<double> div_v2rhosigma_gdrho_r(3 * nrxx);
std::vector<ModuleBase::Vector3<double>> v2rhosigma_gdrho_r(3 * nrxx);
for (int ir = 0; ir < nrxx; ++ir)
Expand All @@ -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<double> v2sigma2_sigma_r(3 * nrxx);
for (int ir = 0; ir < nrxx; ++ir)
{
Expand All @@ -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<double> 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 )
Expand Down
47 changes: 39 additions & 8 deletions source/module_beyonddft/potentials/pot_hxc_lrtd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "module_base/timer.h"
#include "module_hamilt_general/module_xc/xc_functional.h"
#include <set>
#include "module_beyonddft/utils/lr_util.h"
namespace elecstate
{
// constructor for exchange-correlation kernel
Expand All @@ -26,37 +27,67 @@ 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);
}
}

void PotHxcLR::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");
#ifdef USE_LIBXC
const int nspin = v_eff.nr;
v_eff += H_Hartree_pw::v_hartree(*ucell, const_cast<ModulePW::PW_Basis*>(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<ModuleBase::Vector3<double>> drho(nrxx);
LR_Util::grad(rho[0], drho.data(), *(this->rho_basis_), ucell->tpiba);
std::vector<double> 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__));

Expand Down
8 changes: 2 additions & 6 deletions source/module_beyonddft/potentials/pot_hxc_lrtd.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
};

Expand Down

0 comments on commit 5fea1e3

Please sign in to comment.