Skip to content

Commit 2c109a0

Browse files
maki49maki49
authored
Feature: GGA XC kernel for LR-TDDFT at nspin=2 (#5712)
* gga fxc->vxc for close-shell modify gga fxc * spin2-gga (open shell) * add gga cut-off * add a test case
1 parent ab2a789 commit 2c109a0

File tree

8 files changed

+335
-14
lines changed

8 files changed

+335
-14
lines changed

‎source/module_lr/potentials/pot_hxc_lrtd.cpp‎

Lines changed: 145 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ namespace LR
1717
:xc_kernel_(xc_kernel), tpiba_(ucell.tpiba), spin_type_(st), rho_basis_(rho_basis), nrxx_(chg_gs.nrxx),
1818
nspin_(PARAM.inp.nspin == 1 || (PARAM.inp.nspin == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z) ? 1 : 2),
1919
pot_hartree_(LR_Util::make_unique<elecstate::PotHartree>(&rho_basis)),
20-
xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel), //call XC_Functional::set_func_type and libxc
20+
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
2121
xc_type_(XCType(XC_Functional::get_func_type()))
2222
{
2323
if (std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); }
@@ -127,7 +127,7 @@ namespace LR
127127
for (int ir = 0;ir < nrxx;++ir)
128128
{
129129
e_drho[ir] = -(fxc.v2rhosigma_2drho.at(ir) * rho[ir]
130-
+ fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(ir) * drho.at(ir))
130+
+ fxc.v2sigma2_4drho.at(ir) * (fxc.drho_gs.at(0).at(ir) * drho.at(ir))
131131
+ drho.at(ir) * fxc.vsigma.at(ir) * 2.);
132132
}
133133
XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
@@ -141,10 +141,149 @@ namespace LR
141141
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
142142
};
143143
break;
144-
// case SpinType::S2_singlet:
145-
// break;
146-
// case SpinType::S2_triplet:
147-
// break;
144+
case SpinType::S2_singlet:
145+
funcs[s] = [this, &fxc](FXC_PARA_TYPE)-> void
146+
{
147+
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
148+
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);
149+
150+
std::vector<double> vxc_tmp(nrxx, 0.0);
151+
152+
// 1. the terms in grad_dot int f_uu
153+
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
154+
for (int ir = 0;ir < nrxx;++ir)
155+
{
156+
// gdot terms in f_uu
157+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_singlet.at(ir)
158+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_singlet.at(ir)
159+
+ drho.at(ir) * (fxc.vsigma.at(ir * 3) * 2. + fxc.vsigma.at(ir * 3 + 1)));
160+
}
161+
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
162+
163+
// 2. terms not in grad_dot
164+
for (int ir = 0;ir < nrxx;++ir)
165+
{
166+
vxc_tmp[ir] += rho[ir] * (fxc.v2rho2.at(ir * 3) + fxc.v2rho2.at(ir * 3 + 1))
167+
+ drho.at(ir) * fxc.v2rhosigma_drho_singlet.at(ir);
168+
}
169+
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
170+
};
171+
break;
172+
case SpinType::S2_triplet:
173+
funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
174+
{
175+
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
176+
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);
177+
178+
std::vector<double> vxc_tmp(nrxx, 0.0);
179+
180+
// 1. the terms in grad_dot int f_uu
181+
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
182+
for (int ir = 0;ir < nrxx;++ir)
183+
{
184+
// gdot terms in f_uu
185+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_triplet.at(ir)
186+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_triplet.at(ir)
187+
+ drho.at(ir) * (fxc.vsigma.at(ir * 3) * 2. - fxc.vsigma.at(ir * 3 + 1)));
188+
}
189+
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
190+
191+
// 2. terms not in grad_dot
192+
for (int ir = 0;ir < nrxx;++ir)
193+
{
194+
vxc_tmp[ir] += rho[ir] * (fxc.v2rho2.at(ir * 3) - fxc.v2rho2.at(ir * 3 + 1))
195+
+ drho.at(ir) * fxc.v2rhosigma_drho_triplet.at(ir);
196+
}
197+
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
198+
};
199+
break;
200+
case SpinType::S2_updown:
201+
funcs[s] = [this, &fxc](FXC_PARA_TYPE)->void
202+
{
203+
assert(ispin_op.size() >= 2);
204+
std::vector<ModuleBase::Vector3<double>> drho(nrxx); // transition density gradient
205+
LR_Util::grad(rho, drho.data(), this->rho_basis_, this->tpiba_);
206+
207+
std::vector<double> vxc_tmp(nrxx, 0.0);
208+
209+
// 1. the terms in grad_dot int f_uu
210+
std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx);
211+
switch (ispin_op[0] << 1 | ispin_op[1])
212+
{
213+
case 0: // (0,0)
214+
for (int ir = 0;ir < nrxx;++ir)
215+
{
216+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_uu.at(ir)
217+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_uu_u.at(ir)
218+
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_uu_d.at(ir)
219+
+ drho.at(ir) * fxc.vsigma.at(ir * 3) * 2.);
220+
}
221+
break;
222+
case 1: // (0,1)
223+
for (int ir = 0;ir < nrxx;++ir)
224+
{
225+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_du.at(ir) // rho_d, drho_u
226+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_ud_u.at(ir)
227+
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_ud_d.at(ir)
228+
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 1));
229+
}
230+
break;
231+
case 2: // (1,0)
232+
for (int ir = 0;ir < nrxx;++ir)
233+
{
234+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_ud.at(ir) // rho_u, drho_d
235+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_du_u.at(ir)
236+
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_du_d.at(ir)
237+
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 1));
238+
}
239+
break;
240+
case 3: // (1,1)
241+
for (int ir = 0;ir < nrxx;++ir)
242+
{
243+
gdot_terms[ir] = -(rho[ir] * fxc.v2rhosigma_drho_dd.at(ir)
244+
+ (fxc.drho_gs.at(0).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_dd_u.at(ir)
245+
+ (fxc.drho_gs.at(1).at(ir) * drho.at(ir)) * fxc.v2sigma2_drho_dd_d.at(ir)
246+
+ drho.at(ir) * fxc.vsigma.at(ir * 3 + 2) * 2.);
247+
}
248+
break;
249+
default:
250+
throw std::runtime_error("Invalid ispin_op");
251+
}
252+
XC_Functional::grad_dot(gdot_terms.data(), vxc_tmp.data(), &this->rho_basis_, this->tpiba_);
253+
254+
// 2. terms not in grad_dot
255+
switch (ispin_op[0] << 1 | ispin_op[1])
256+
{
257+
case 0:
258+
for (int ir = 0;ir < nrxx;++ir)
259+
{
260+
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3) + drho.at(ir) * fxc.v2rhosigma_drho_uu.at(ir);
261+
}
262+
break;
263+
case 1:
264+
for (int ir = 0;ir < nrxx;++ir)
265+
{
266+
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 1) + drho.at(ir) * fxc.v2rhosigma_drho_ud.at(ir);
267+
}
268+
break;
269+
case 2:
270+
for (int ir = 0;ir < nrxx;++ir)
271+
{
272+
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 1) + drho.at(ir) * fxc.v2rhosigma_drho_du.at(ir);
273+
}
274+
break;
275+
case 3:
276+
for (int ir = 0;ir < nrxx;++ir)
277+
{
278+
vxc_tmp[ir] += rho[ir] * fxc.v2rho2.at(ir * 3 + 2) + drho.at(ir) * fxc.v2rhosigma_drho_dd.at(ir);
279+
}
280+
break;
281+
default:
282+
throw std::runtime_error("Invalid ispin_op");
283+
}
284+
BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1);
285+
};
286+
break;
148287
default:
149288
throw std::domain_error("SpinType =" + std::to_string(static_cast<int>(s)) + "for GGA or HYB_GGA is unfinished in "
150289
+ std::string(__FILE__) + " line " + std::to_string(__LINE__));

‎source/module_lr/potentials/xc_kernel.cpp‎

Lines changed: 87 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,8 @@ LR::KernelXC::KernelXC(const ModulePW::PW_Basis& rho_basis,
2121
const Parallel_Grid& pgrid,
2222
const int& nspin,
2323
const std::string& kernel_name,
24-
const std::vector<std::string>& lr_init_xc_kernel) :rho_basis_(rho_basis)
24+
const std::vector<std::string>& lr_init_xc_kernel,
25+
const bool openshell) :rho_basis_(rho_basis), openshell_(openshell)
2526
{
2627
if (!std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(kernel_name)) { return; }
2728
XC_Functional::set_xc_type(kernel_name); // for hse, (1-alpha) and omega are set here
@@ -92,6 +93,23 @@ inline void add_assign_op(const std::vector<T>& src, std::vector<T>& dst)
9293
{
9394
add_op(src, dst, dst);
9495
}
96+
template<typename Telement, typename Tscalar>
97+
inline void cutoff_grid_data_spin2(std::vector<Telement>& func, const std::vector<Tscalar>& mask)
98+
{
99+
const int& nrxx = mask.size() / 2;
100+
assert(func.size() % nrxx == 0 && func.size() / nrxx > 1);
101+
const int n_component = func.size() / nrxx;
102+
#ifdef _OPENMP
103+
#pragma omp parallel for schedule(static, 4096)
104+
#endif
105+
for (int ir = 0;ir < nrxx;++ir)
106+
{
107+
const int& i2 = 2 * ir;
108+
const int& istart = n_component * ir;
109+
std::for_each(func.begin() + istart, func.begin() + istart + n_component - 1, [&](Telement& f) { f *= mask[i2]; }); //spin-up
110+
std::for_each(func.begin() + istart + 1, func.begin() + istart + n_component, [&](Telement& f) { f *= mask[i2 + 1]; }); //spin-down
111+
}
112+
}
95113

96114
#ifdef USE_LIBXC
97115
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)
@@ -190,7 +208,8 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
190208

191209
xc_func_set_dens_threshold(&func, rho_threshold);
192210

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

195214
// Libxc interfaces overwrite (instead of add onto) the output arrays, so we need temporary copies
196215
std::vector<double> vrho_tmp(this->vrho_.size());
@@ -206,9 +225,19 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
206225
break;
207226
case XC_FAMILY_GGA:
208227
case XC_FAMILY_HYB_GGA:
228+
{
209229
xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho_tmp.data(), vsigma_tmp.data());
210230
xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_tmp.data(), v2rhosigma_tmp.data(), v2sigma2_tmp.data());
231+
// std::cout << "max element of v2sigma2_tmp: " << *std::max_element(v2sigma2_tmp.begin(), v2sigma2_tmp.end()) << std::endl;
232+
// 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;
233+
// cut off by sgn
234+
cutoff_grid_data_spin2(vrho_tmp, sgn);
235+
cutoff_grid_data_spin2(vsigma_tmp, sgn);
236+
cutoff_grid_data_spin2(v2rho2_tmp, sgn);
237+
cutoff_grid_data_spin2(v2rhosigma_tmp, sgn);
238+
cutoff_grid_data_spin2(v2sigma2_tmp, sgn);
211239
break;
240+
}
212241
default:
213242
throw std::domain_error("func.info->family =" + std::to_string(func.info->family)
214243
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
@@ -239,8 +268,6 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
239268

240269
if (nspin == 1)
241270
{
242-
// 0. drho
243-
this->drho_gs_ = gradrho[0];
244271
// 1. $2f^{\rho\sigma}*\nabla\rho$
245272
this->v2rhosigma_2drho_.resize(nrxx);
246273
#ifdef _OPENMP
@@ -261,11 +288,67 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
261288
this->v2sigma2_4drho_[i] = gradrho[0][i] * v2s2[i] * 4.;
262289
}
263290
}
291+
else if (2 == nspin) //close-shell
292+
{
293+
if (!openshell_)
294+
{
295+
this->v2rhosigma_drho_singlet_.resize(nrxx);
296+
this->v2sigma2_drho_singlet_.resize(nrxx);
297+
this->v2rhosigma_drho_triplet_.resize(nrxx);
298+
this->v2sigma2_drho_triplet_.resize(nrxx);
299+
#ifdef _OPENMP
300+
#pragma omp parallel for schedule(static, 4096)
301+
#endif
302+
for (int i = 0;i < nrxx;++i)
303+
{
304+
const int istart = i * 6;
305+
this->v2rhosigma_drho_singlet_[i] = gradrho[0][i] * (v2rs[istart] + v2rs[istart + 1] + v2rs[istart + 2]) * 2.;
306+
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.;
307+
this->v2rhosigma_drho_triplet_[i] = gradrho[0][i] * (v2rs[istart] - v2rs[istart + 2]) * 2.;
308+
this->v2sigma2_drho_triplet_[i] = gradrho[0][i] * (v2s2[istart] * 2. + v2s2[istart + 1] - v2s2[istart + 2] * 2. - v2s2[istart + 4]) * 2.;
309+
}
310+
}
311+
else
312+
{
313+
this->v2rhosigma_drho_uu_.resize(nrxx);
314+
this->v2rhosigma_drho_ud_.resize(nrxx);
315+
this->v2rhosigma_drho_du_.resize(nrxx);
316+
this->v2rhosigma_drho_dd_.resize(nrxx);
317+
this->v2sigma2_drho_uu_u_.resize(nrxx);
318+
this->v2sigma2_drho_uu_d_.resize(nrxx);
319+
this->v2sigma2_drho_ud_u_.resize(nrxx);
320+
this->v2sigma2_drho_ud_d_.resize(nrxx);
321+
this->v2sigma2_drho_du_u_.resize(nrxx);
322+
this->v2sigma2_drho_du_d_.resize(nrxx);
323+
this->v2sigma2_drho_dd_u_.resize(nrxx);
324+
this->v2sigma2_drho_dd_d_.resize(nrxx);
325+
#ifdef _OPENMP
326+
#pragma omp parallel for schedule(static, 4096)
327+
#endif
328+
for (int i = 0;i < nrxx;++i)
329+
{
330+
const int istart = i * 6;
331+
this->v2rhosigma_drho_uu_[i] = gradrho[0][i] * v2rs[istart] * 2. + gradrho[1][i] * v2rs[istart + 1];
332+
this->v2rhosigma_drho_ud_[i] = gradrho[0][i] * v2rs[istart + 1] + gradrho[1][i] * v2rs[istart + 2] * 2.;
333+
this->v2rhosigma_drho_du_[i] = gradrho[0][i] * v2rs[istart + 3] * 2. + gradrho[1][i] * v2rs[istart + 4];
334+
this->v2rhosigma_drho_dd_[i] = gradrho[0][i] * v2rs[istart + 4] + gradrho[1][i] * v2rs[istart + 5] * 2.;
335+
this->v2sigma2_drho_uu_u_[i] = gradrho[0][i] * v2s2[istart] * 4. + gradrho[1][i] * v2s2[istart + 1] * 2.;
336+
this->v2sigma2_drho_uu_d_[i] = gradrho[0][i] * v2s2[istart + 1] * 2. + gradrho[1][i] * v2s2[istart + 3];
337+
this->v2sigma2_drho_ud_u_[i] = gradrho[0][i] * v2s2[istart + 1] * 2. + gradrho[1][i] * v2s2[istart + 3];
338+
this->v2sigma2_drho_ud_d_[i] = gradrho[0][i] * v2s2[istart + 2] * 4. + gradrho[1][i] * v2s2[istart + 4] * 2.;
339+
this->v2sigma2_drho_du_u_[i] = gradrho[1][i] * v2s2[istart + 2] * 4. + gradrho[0][i] * v2s2[istart + 1] * 2.;
340+
this->v2sigma2_drho_du_d_[i] = gradrho[1][i] * v2s2[istart + 4] * 2. + gradrho[0][i] * v2s2[istart + 3];
341+
this->v2sigma2_drho_dd_u_[i] = gradrho[1][i] * v2s2[istart + 4] * 2. + gradrho[0][i] * v2s2[istart + 3];
342+
this->v2sigma2_drho_dd_d_[i] = gradrho[1][i] * v2s2[istart + 5] * 4. + gradrho[0][i] * v2s2[istart + 4] * 2.;
343+
}
344+
}
345+
}
264346
else
265347
{
266348
throw std::domain_error("nspin =" + std::to_string(nspin)
267349
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
268350
}
351+
this->drho_gs_ = std::move(gradrho);
269352
}
270353
if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) {
271354
return;

0 commit comments

Comments
 (0)