Skip to content

Feature: GGA XC kernel for LR-TDDFT at nspin=2 #5712

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Dec 10, 2024
Merged

Feature: GGA XC kernel for LR-TDDFT at nspin=2 #5712

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
037c85e
gga fxc->vxc for close-shell
maki49 Nov 5, 2024
788a887
spin2-gga (open shell)
maki49 Nov 30, 2024
c14b2c7
add gga cut-off
maki49 Dec 4, 2024
085a86f
add a test case
maki49 Dec 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 145 additions & 6 deletions source/module_lr/potentials/pot_hxc_lrtd.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ namespace LR
:xc_kernel_(xc_kernel), tpiba_(ucell.tpiba), spin_type_(st), rho_basis_(rho_basis), nrxx_(chg_gs.nrxx),
nspin_(PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z) ? 1 : 2),
pot_hartree_(LR_Util::make_unique<elecstate::PotHartree>(&rho_basis)),
xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel), //call XC_Functional::set_func_type and libxc
xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel, (st == SpinType::S2_updown)), //call XC_Functional::set_func_type and libxc
xc_type_(XCType(XC_Functional::get_func_type()))
{
if (std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); }
Expand Down Expand Up @@ -127,7 +127,7 @@ namespace LR
for (int ir = 0;ir < nrxx;++ir)
{
e_drho[ir] = -(fxc.v2rhosigma_2drho.at(ir) * rho[ir]
+ fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(ir) * drho.at(ir))
+ fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(0).at(ir) * drho.at(ir))
+ drho.at(ir) * fxc.vsigma.at(ir) * 2.);
}
XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
Expand All @@ -141,10 +141,149 @@ namespace LR
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
};
break;
// case SpinType::S2_singlet:
// break;
// case SpinType::S2_triplet:
// break;
case SpinType::S2_singlet:
funcs[s] = [this, &fxc](FXC_PARA_TYPE)-> void
{
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);

std::vector<double> vxc_tmp(nrxx, 0.0);

// 1. the terms in grad_dot int f_uu
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
for (int ir = 0;ir < nrxx;++ir)
{
// gdot terms in f_uu
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_singlet.at(ir)
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_singlet.at(ir)
+ drho.at(ir) * (fxc.vsigma.at(ir * 3) * 2. + fxc.vsigma.at(ir * 3 + 1)));
}
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);

// 2. terms not in grad_dot
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * (fxc.v2rho2.at(ir * 3) + fxc.v2rho2.at(ir * 3 + 1))
+ drho.at(ir) * fxc.v2rhosigma_drho_singlet.at(ir);
}
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
};
break;
case SpinType::S2_triplet:
funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
{
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);

std::vector<double> vxc_tmp(nrxx, 0.0);

// 1. the terms in grad_dot int f_uu
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
for (int ir = 0;ir < nrxx;++ir)
{
// gdot terms in f_uu
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_triplet.at(ir)
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_triplet.at(ir)
+ drho.at(ir) * (fxc.vsigma.at(ir * 3) * 2. - fxc.vsigma.at(ir * 3 + 1)));
}
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);

// 2. terms not in grad_dot
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * (fxc.v2rho2.at(ir * 3) - fxc.v2rho2.at(ir * 3 + 1))
+ drho.at(ir) * fxc.v2rhosigma_drho_triplet.at(ir);
}
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
};
break;
case SpinType::S2_updown:
funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
{
assert(ispin_op.size() >= 2);
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);

std::vector<double> vxc_tmp(nrxx, 0.0);

// 1. the terms in grad_dot int f_uu
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
switch (ispin_op[0] << 1 | ispin_op[1])
{
case 0: // (0,0)
for (int ir = 0;ir < nrxx;++ir)
{
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_uu.at(ir)
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_uu_u.at(ir)
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_uu_d.at(ir)
+ drho.at(ir) * fxc.vsigma.at(ir * 3) * 2.);
}
break;
case 1: // (0,1)
for (int ir = 0;ir < nrxx;++ir)
{
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_du.at(ir) // rho_d, drho_u
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_ud_u.at(ir)
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_ud_d.at(ir)
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 1));
}
break;
case 2: // (1,0)
for (int ir = 0;ir < nrxx;++ir)
{
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_ud.at(ir) // rho_u, drho_d
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_du_u.at(ir)
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_du_d.at(ir)
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 1));
}
break;
case 3: // (1,1)
for (int ir = 0;ir < nrxx;++ir)
{
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_dd.at(ir)
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_dd_u.at(ir)
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_dd_d.at(ir)
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 2) * 2.);
}
break;
default:
throw std::runtime_error("Invalid ispin_op");
}
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);

// 2. terms not in grad_dot
switch (ispin_op[0] << 1 | ispin_op[1])
{
case 0:
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3) + drho.at(ir) * fxc.v2rhosigma_drho_uu.at(ir);
}
break;
case 1:
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 1) + drho.at(ir) * fxc.v2rhosigma_drho_ud.at(ir);
}
break;
case 2:
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 1) + drho.at(ir) * fxc.v2rhosigma_drho_du.at(ir);
}
break;
case 3:
for (int ir = 0;ir < nrxx;++ir)
{
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 2) + drho.at(ir) * fxc.v2rhosigma_drho_dd.at(ir);
}
break;
default:
throw std::runtime_error("Invalid ispin_op");
}
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
};
break;
default:
throw std::domain_error("SpinType =" + std::to_string(static_cast<int>(s)) + "for GGA or HYB_GGA is unfinished in "
+ std::string(__FILE__) + " line " + std::to_string(__LINE__));
Expand Down
91 changes: 87 additions & 4 deletions source/module_lr/potentials/xc_kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ LR::KernelXC::KernelXC(const ModulePW::PW_Basis& rho_basis,
const Parallel_Grid& pgrid,
const int& nspin,
const std::string& kernel_name,
const std::vector<std::string>& lr_init_xc_kernel) :rho_basis_(rho_basis)
const std::vector<std::string>& lr_init_xc_kernel,
const bool openshell) :rho_basis_(rho_basis), openshell_(openshell)
{
if (!std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(kernel_name)) { return; }
XC_Functional::set_xc_type(kernel_name); // for hse, (1-alpha) and omega are set here
Expand Down Expand Up @@ -92,6 +93,23 @@ inline void add_assign_op(const std::vector<T>& src, std::vector<T>& dst)
{
add_op(src, dst, dst);
}
template<typename Telement, typename Tscalar>
inline void cutoff_grid_data_spin2(std::vector<Telement>& func, const std::vector<Tscalar>& mask)
{
const int& nrxx = mask.size() / 2;
assert(func.size() % nrxx == 0 && func.size() / nrxx > 1);
const int n_component = func.size() / nrxx;
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 4096)
#endif
for (int ir = 0;ir < nrxx;++ir)
{
const int& i2 = 2 * ir;
const int& istart = n_component * ir;
std::for_each(func.begin() + istart, func.begin() + istart + n_component - 1, [&](Telement& f) { f *= mask[i2]; }); //spin-up
std::for_each(func.begin() + istart + 1, func.begin() + istart + n_component, [&](Telement& f) { f *= mask[i2 + 1]; }); //spin-down
}
}

#ifdef USE_LIBXC
void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const double* const* const rho_gs, const double* const rho_core)
Expand Down Expand Up @@ -190,7 +208,8 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl

xc_func_set_dens_threshold(&func, rho_threshold);

//cut off grho if not LDA (future subfunc)
//cut off function
const std::vector<double> sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma);

// Libxc interfaces overwrite (instead of add onto) the output arrays, so we need temporary copies
std::vector<double> vrho_tmp(this->vrho_.size());
Expand All @@ -206,9 +225,19 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
break;
case XC_FAMILY_GGA:
case XC_FAMILY_HYB_GGA:
{
xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho_tmp.data(), vsigma_tmp.data());
xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_tmp.data(), v2rhosigma_tmp.data(), v2sigma2_tmp.data());
// std::cout << "max element of v2sigma2_tmp: " << *std::max_element(v2sigma2_tmp.begin(), v2sigma2_tmp.end()) << std::endl;
// std::cout << "rho corresponding to max element of v2sigma2_tmp: " << rho[(std::max_element(v2sigma2_tmp.begin(), v2sigma2_tmp.end()) - v2sigma2_tmp.begin()) / 6] << std::endl;
// cut off by sgn
cutoff_grid_data_spin2(vrho_tmp, sgn);
cutoff_grid_data_spin2(vsigma_tmp, sgn);
cutoff_grid_data_spin2(v2rho2_tmp, sgn);
cutoff_grid_data_spin2(v2rhosigma_tmp, sgn);
cutoff_grid_data_spin2(v2sigma2_tmp, sgn);
break;
}
default:
throw std::domain_error("func.info->family =" + std::to_string(func.info->family)
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
Expand Down Expand Up @@ -239,8 +268,6 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl

if (nspin == 1)
{
// 0. drho
this->drho_gs_ = gradrho[0];
// 1. $2f^{\rho\sigma}*\nabla\rho$
this->v2rhosigma_2drho_.resize(nrxx);
#ifdef _OPENMP
Expand All @@ -261,11 +288,67 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
this->v2sigma2_4drho_[i] = gradrho[0][i] * v2s2[i] * 4.;
}
}
else if (2 == nspin) //close-shell
{
if (!openshell_)
{
this->v2rhosigma_drho_singlet_.resize(nrxx);
this->v2sigma2_drho_singlet_.resize(nrxx);
this->v2rhosigma_drho_triplet_.resize(nrxx);
this->v2sigma2_drho_triplet_.resize(nrxx);
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 4096)
#endif
for (int i = 0;i < nrxx;++i)
{
const int istart = i * 6;
this->v2rhosigma_drho_singlet_[i] = gradrho[0][i] * (v2rs[istart] + v2rs[istart + 1] + v2rs[istart + 2]) * 2.;
this->v2sigma2_drho_singlet_[i] = gradrho[0][i] * (v2s2[istart] * 2. + v2s2[istart + 1] * 3. + v2s2[istart + 2] * 2. + v2s2[istart + 3] + v2s2[istart + 4]) * 2.;
this->v2rhosigma_drho_triplet_[i] = gradrho[0][i] * (v2rs[istart] - v2rs[istart + 2]) * 2.;
this->v2sigma2_drho_triplet_[i] = gradrho[0][i] * (v2s2[istart] * 2. + v2s2[istart + 1] - v2s2[istart + 2] * 2. - v2s2[istart + 4]) * 2.;
}
}
else
{
this->v2rhosigma_drho_uu_.resize(nrxx);
this->v2rhosigma_drho_ud_.resize(nrxx);
this->v2rhosigma_drho_du_.resize(nrxx);
this->v2rhosigma_drho_dd_.resize(nrxx);
this->v2sigma2_drho_uu_u_.resize(nrxx);
this->v2sigma2_drho_uu_d_.resize(nrxx);
this->v2sigma2_drho_ud_u_.resize(nrxx);
this->v2sigma2_drho_ud_d_.resize(nrxx);
this->v2sigma2_drho_du_u_.resize(nrxx);
this->v2sigma2_drho_du_d_.resize(nrxx);
this->v2sigma2_drho_dd_u_.resize(nrxx);
this->v2sigma2_drho_dd_d_.resize(nrxx);
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 4096)
#endif
for (int i = 0;i < nrxx;++i)
{
const int istart = i * 6;
this->v2rhosigma_drho_uu_[i] = gradrho[0][i] * v2rs[istart] * 2. + gradrho[1][i] * v2rs[istart + 1];
this->v2rhosigma_drho_ud_[i] = gradrho[0][i] * v2rs[istart + 1] + gradrho[1][i] * v2rs[istart + 2] * 2.;
this->v2rhosigma_drho_du_[i] = gradrho[0][i] * v2rs[istart + 3] * 2. + gradrho[1][i] * v2rs[istart + 4];
this->v2rhosigma_drho_dd_[i] = gradrho[0][i] * v2rs[istart + 4] + gradrho[1][i] * v2rs[istart + 5] * 2.;
this->v2sigma2_drho_uu_u_[i] = gradrho[0][i] * v2s2[istart] * 4. + gradrho[1][i] * v2s2[istart + 1] * 2.;
this->v2sigma2_drho_uu_d_[i] = gradrho[0][i] * v2s2[istart + 1] * 2. + gradrho[1][i] * v2s2[istart + 3];
this->v2sigma2_drho_ud_u_[i] = gradrho[0][i] * v2s2[istart + 1] * 2. + gradrho[1][i] * v2s2[istart + 3];
this->v2sigma2_drho_ud_d_[i] = gradrho[0][i] * v2s2[istart + 2] * 4. + gradrho[1][i] * v2s2[istart + 4] * 2.;
this->v2sigma2_drho_du_u_[i] = gradrho[1][i] * v2s2[istart + 2] * 4. + gradrho[0][i] * v2s2[istart + 1] * 2.;
this->v2sigma2_drho_du_d_[i] = gradrho[1][i] * v2s2[istart + 4] * 2. + gradrho[0][i] * v2s2[istart + 3];
this->v2sigma2_drho_dd_u_[i] = gradrho[1][i] * v2s2[istart + 4] * 2. + gradrho[0][i] * v2s2[istart + 3];
this->v2sigma2_drho_dd_d_[i] = gradrho[1][i] * v2s2[istart + 5] * 4. + gradrho[0][i] * v2s2[istart + 4] * 2.;
}
}
}
else
{
throw std::domain_error("nspin =" + std::to_string(nspin)
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
}
this->drho_gs_ = std::move(gradrho);
}
if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) {
return;
Expand Down
Loading
Loading