Skip to content

Commit d57aa19

Browse files
author
dyzheng
committed
Feat(xc/psi): libxc gga_grad=3 SF, PAGED_GPU to_range/fix_k fixes
P0 - libxc gga_grad=3 SF transform (b5b2898): - Implement compute_mag_part_nspin4(), cal_gdr_sf(), cal_dh_sf(), convert_v_nspin4_sf() in xc_functional_libxc_tools.cpp - Route gga_grad>=2 nspin=4 through SF path in vxc and convert_vtxc_v P1 - PAGED_GPU pointer offset fixes: - fix to_range() for PAGED_GPU: use base offset 0 instead of k*nbands*nbasis - fix fix_k() for PAGED_GPU: set psi_bias=0, psi_current=psi - fix copy constructors: separate psi_bias/psi_current for PAGED_GPU vs ALL_GPU branches, set current_k_gpu_=-1 after copy P2 - GPU transfer (already handled by PagingTransfer template specialization)
1 parent 4dfc525 commit d57aa19

3 files changed

Lines changed: 272 additions & 22 deletions

File tree

‎source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp‎

Lines changed: 214 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,28 @@ XC_Functional_Libxc::convert_rho_amag_nspin4(
5151
return std::make_tuple(std::move(rho), std::move(amag));
5252
}
5353

54+
std::vector<double> XC_Functional_Libxc::compute_mag_part_nspin4(
55+
const std::size_t nrxx,
56+
const Charge* const chr)
57+
{
58+
std::vector<double> mag_part(3 * nrxx, 0.0);
59+
#ifdef _OPENMP
60+
#pragma omp parallel for schedule(static, 1024)
61+
#endif
62+
for (std::size_t ir = 0; ir < nrxx; ++ir)
63+
{
64+
double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir];
65+
double amag = std::sqrt(mx * mx + my * my + mz * mz);
66+
if (amag > 1e-12)
67+
{
68+
mag_part[ir] = mx / amag;
69+
mag_part[ir + nrxx] = my / amag;
70+
mag_part[ir + 2 * nrxx] = mz / amag;
71+
}
72+
}
73+
return mag_part;
74+
}
75+
5476
// calculating grho
5577
std::vector<std::vector<ModuleBase::Vector3<double>>>
5678
XC_Functional_Libxc::cal_gdr(
@@ -213,22 +235,48 @@ std::pair<double,ModuleBase::matrix> XC_Functional_Libxc::convert_vtxc_v(
213235

214236
if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
215237
{
216-
const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr);
217-
218-
double rvtxc = 0.0;
219-
#ifdef _OPENMP
220-
#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256)
221-
#endif
222-
for( int is=0; is<nspin; ++is )
238+
if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
223239
{
224-
for( std::size_t ir=0; ir<nrxx; ++ir )
240+
std::vector<double> mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);
241+
const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh_sf(nspin, nrxx, sgn, gdr, vsigma, mag_part_tmp, tpiba, chr);
242+
243+
constexpr int nspin4 = 4;
244+
double rvtxc = 0.0;
245+
#ifdef _OPENMP
246+
#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256)
247+
#endif
248+
for (int is = 0; is < nspin4; ++is)
225249
{
226-
rvtxc += dh[is][ir] * rho[ir*nspin+is];
227-
v(is,ir) -= dh[is][ir];
250+
for (std::size_t ir = 0; ir < nrxx; ++ir)
251+
{
252+
double rho_ir = 0.0;
253+
if (is == 0) { rho_ir = rho[ir * nspin + 0]; }
254+
else { rho_ir = rho[ir * nspin + 0] * mag_part_tmp[ir + (is - 1) * nrxx]; }
255+
rvtxc += dh[is][ir] * rho_ir;
256+
v(is, ir) -= dh[is][ir];
257+
}
228258
}
259+
vtxc -= rvtxc;
229260
}
261+
else
262+
{
263+
const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr);
264+
265+
double rvtxc = 0.0;
266+
#ifdef _OPENMP
267+
#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256)
268+
#endif
269+
for( int is=0; is<nspin; ++is )
270+
{
271+
for( std::size_t ir=0; ir<nrxx; ++ir )
272+
{
273+
rvtxc += dh[is][ir] * rho[ir*nspin+is];
274+
v(is,ir) -= dh[is][ir];
275+
}
276+
}
230277

231-
vtxc -= rvtxc;
278+
vtxc -= rvtxc;
279+
}
232280
} // end if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA))
233281

234282
return std::make_pair(vtxc, std::move(v));
@@ -317,4 +365,159 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
317365
return v_nspin4;
318366
}
319367

368+
std::vector<std::vector<ModuleBase::Vector3<double>>> XC_Functional_Libxc::cal_gdr_sf(
369+
const int nspin,
370+
const std::size_t nrxx,
371+
const std::vector<double> &rho,
372+
const std::vector<double> &mag_part,
373+
const double tpiba,
374+
const Charge* const chr)
375+
{
376+
std::vector<std::vector<ModuleBase::Vector3<double>>> gdr(nspin);
377+
std::vector<double> rhor(nrxx);
378+
std::vector<std::complex<double>> rhog(chr->rhopw->npw);
379+
std::vector<ModuleBase::Vector3<double>> gdr_tmp(nrxx);
380+
381+
#ifdef _OPENMP
382+
#pragma omp parallel for schedule(static, 1024)
383+
#endif
384+
for (std::size_t ir = 0; ir < nrxx; ++ir)
385+
{
386+
rhor[ir] = rho[ir * nspin + 0];
387+
}
388+
chr->rhopw->real2recip(rhor.data(), rhog.data());
389+
gdr[0].resize(nrxx);
390+
XC_Functional::grad_rho(rhog.data(), gdr[0].data(), chr->rhopw, tpiba);
391+
392+
gdr[1].resize(nrxx);
393+
#ifdef _OPENMP
394+
#pragma omp parallel for schedule(static, 1024)
395+
#endif
396+
for (std::size_t ir = 0; ir < nrxx; ++ir)
397+
{
398+
gdr_tmp[ir] = gdr[0][ir];
399+
gdr[0][ir] = 0.5 * gdr_tmp[ir];
400+
gdr[1][ir] = 0.5 * gdr_tmp[ir];
401+
}
402+
403+
for (int is = 1; is <= 3; ++is)
404+
{
405+
chr->rhopw->real2recip(chr->rho[is], rhog.data());
406+
XC_Functional::grad_rho(rhog.data(), gdr_tmp.data(), chr->rhopw, tpiba);
407+
const double* mp = mag_part.data() + (is - 1) * nrxx;
408+
#ifdef _OPENMP
409+
#pragma omp parallel for schedule(static, 1024)
410+
#endif
411+
for (std::size_t ir = 0; ir < nrxx; ++ir)
412+
{
413+
const ModuleBase::Vector3<double> g = 0.5 * gdr_tmp[ir] * mp[ir];
414+
gdr[0][ir] += g;
415+
gdr[1][ir] -= g;
416+
}
417+
}
418+
419+
return gdr;
420+
}
421+
422+
std::vector<std::vector<double>> XC_Functional_Libxc::cal_dh_sf(
423+
const int nspin,
424+
const std::size_t nrxx,
425+
const std::vector<double> &sgn,
426+
const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
427+
const std::vector<double> &vsigma,
428+
const std::vector<double> &mag_part,
429+
const double tpiba,
430+
const Charge* const chr)
431+
{
432+
std::vector<ModuleBase::Vector3<double>> h1(nrxx), h2(nrxx);
433+
#ifdef _OPENMP
434+
#pragma omp parallel for schedule(static, 1024)
435+
#endif
436+
for (std::size_t ir = 0; ir < nrxx; ++ir)
437+
{
438+
h1[ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0
439+
+ gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]);
440+
h2[ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0
441+
+ gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]);
442+
}
443+
444+
std::vector<ModuleBase::Vector3<double>> tmp_h(nrxx);
445+
std::vector<double> dh0(nrxx), dh_mu(nrxx);
446+
447+
#ifdef _OPENMP
448+
#pragma omp parallel for schedule(static, 1024)
449+
#endif
450+
for (std::size_t ir = 0; ir < nrxx; ++ir)
451+
{
452+
tmp_h[ir] = 0.5 * (h1[ir] + h2[ir]);
453+
}
454+
XC_Functional::grad_dot(tmp_h.data(), dh0.data(), chr->rhopw, tpiba);
455+
456+
std::vector<std::vector<double>> dh_total(4, std::vector<double>(nrxx, 0.0));
457+
#ifdef _OPENMP
458+
#pragma omp parallel for schedule(static, 1024)
459+
#endif
460+
for (std::size_t ir = 0; ir < nrxx; ++ir)
461+
{
462+
dh_total[0][ir] = dh0[ir];
463+
}
464+
465+
for (int mu = 1; mu < 4; ++mu)
466+
{
467+
const double* mp_mu = mag_part.data() + (mu - 1) * nrxx;
468+
#ifdef _OPENMP
469+
#pragma omp parallel for schedule(static, 1024)
470+
#endif
471+
for (std::size_t ir = 0; ir < nrxx; ++ir)
472+
{
473+
tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]) * mp_mu[ir];
474+
}
475+
XC_Functional::grad_dot(tmp_h.data(), dh_mu.data(), chr->rhopw, tpiba);
476+
#ifdef _OPENMP
477+
#pragma omp parallel for schedule(static, 1024)
478+
#endif
479+
for (std::size_t ir = 0; ir < nrxx; ++ir)
480+
{
481+
dh_total[mu][ir] = dh_mu[ir];
482+
}
483+
}
484+
485+
return dh_total;
486+
}
487+
488+
ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4_sf(
489+
const std::size_t nrxx,
490+
const Charge* const chr,
491+
const std::vector<double> &mag_part,
492+
const ModuleBase::matrix &v)
493+
{
494+
assert(PARAM.inp.nspin==4);
495+
constexpr double vanishing_charge = 1.0e-10;
496+
ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx);
497+
498+
for (std::size_t ir = 0; ir < nrxx; ++ir)
499+
{
500+
v_nspin4(0, ir) = 0.5 * (v(0, ir) + v(1, ir));
501+
}
502+
503+
if (PARAM.globalv.domag || PARAM.globalv.domag_z)
504+
{
505+
for (std::size_t ir = 0; ir < nrxx; ++ir)
506+
{
507+
double amag = std::sqrt(std::pow(chr->rho[1][ir], 2)
508+
+ std::pow(chr->rho[2][ir], 2)
509+
+ std::pow(chr->rho[3][ir], 2));
510+
if (amag > vanishing_charge)
511+
{
512+
const double vs = 0.5 * (v(0, ir) - v(1, ir));
513+
for (int ipol = 1; ipol < PARAM.inp.nspin; ++ipol)
514+
{
515+
v_nspin4(ipol, ir) = vs * mag_part[ir + (ipol - 1) * nrxx];
516+
}
517+
}
518+
}
519+
}
520+
return v_nspin4;
521+
}
522+
320523
#endif

‎source/source_hamilt/module_xc/xc_functional_libxc_vxc.cpp‎

Lines changed: 21 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( /
5656
// converting rho
5757
std::vector<double> rho;
5858
std::vector<double> amag;
59+
std::vector<double> mag_part;
5960
if(1==nspin || 2==PARAM.inp.nspin)
6061
{
6162
rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
@@ -65,13 +66,24 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( /
6566
std::tuple<std::vector<double>,std::vector<double>> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr);
6667
rho = std::get<0>(std::move(rho_amag));
6768
amag = std::get<1>(std::move(rho_amag));
69+
if(PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
70+
{
71+
mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);
72+
}
6873
}
6974

7075
std::vector<std::vector<ModuleBase::Vector3<double>>> gdr;
7176
std::vector<double> sigma;
7277
if(is_gga)
7378
{
74-
gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
79+
if(PARAM.inp.nspin==4 && PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
80+
{
81+
gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr);
82+
}
83+
else
84+
{
85+
gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
86+
}
7587
sigma = XC_Functional_Libxc::convert_sigma(gdr);
7688
}
7789

@@ -169,7 +181,14 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( /
169181

170182
if(4==PARAM.inp.nspin)
171183
{
172-
v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v);
184+
if(PARAM.inp.gga_grad >= 2 && (PARAM.globalv.domag || PARAM.globalv.domag_z))
185+
{
186+
v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v);
187+
}
188+
else
189+
{
190+
v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v);
191+
}
173192
}
174193

175194
//-------------------------------------------------

0 commit comments

Comments
 (0)