1+ #ifdef USE_LIBXC
2+
3+ #include " xc_functional_libxc.h"
4+ #include " xc_functional.h"
5+ #include " module_elecstate/module_charge/charge.h"
6+
7+ // converting rho (abacus=>libxc)
8+ std::vector<double > XC_Functional_Libxc::convert_rho (
9+ const int nspin,
10+ const std::size_t nrxx,
11+ const Charge* const chr)
12+ {
13+ std::vector<double > rho (nrxx*nspin);
14+ #ifdef _OPENMP
15+ #pragma omp parallel for collapse(2) schedule(static, 1024)
16+ #endif
17+ for ( int is=0 ; is<nspin; ++is )
18+ for ( int ir=0 ; ir<nrxx; ++ir )
19+ rho[ir*nspin+is] = chr->rho [is][ir] + 1.0 /nspin*chr->rho_core [ir];
20+ return rho;
21+ }
22+
23+ // converting rho (abacus=>libxc)
24+ std::tuple<std::vector<double >, std::vector<double >>
25+ XC_Functional_Libxc::convert_rho_amag_nspin4 (
26+ const int nspin,
27+ const std::size_t nrxx,
28+ const Charge* const chr)
29+ {
30+ assert (GlobalV::NSPIN==4 );
31+ std::vector<double > rho (nrxx*nspin);
32+ std::vector<double > amag (nrxx);
33+ #ifdef _OPENMP
34+ #pragma omp parallel for
35+ #endif
36+ for ( int ir=0 ; ir<nrxx; ++ir )
37+ {
38+ const double arhox = std::abs ( chr->rho [0 ][ir] + chr->rho_core [ir] );
39+ amag[ir] = std::sqrt ( std::pow (chr->rho [1 ][ir],2 )
40+ + std::pow (chr->rho [2 ][ir],2 )
41+ + std::pow (chr->rho [3 ][ir],2 ) );
42+ const double amag_clip = (amag[ir]<arhox) ? amag[ir] : arhox;
43+ rho[ir*nspin+0 ] = (arhox + amag_clip) / 2.0 ;
44+ rho[ir*nspin+1 ] = (arhox - amag_clip) / 2.0 ;
45+ }
46+ return std::make_tuple (std::move (rho), std::move (amag));
47+ }
48+
49+ // calculating grho
50+ std::vector<std::vector<ModuleBase::Vector3<double >>>
51+ XC_Functional_Libxc::cal_gdr (
52+ const int nspin,
53+ const std::vector<double > &rho,
54+ const double tpiba,
55+ const Charge* const chr)
56+ {
57+ std::vector<std::vector<ModuleBase::Vector3<double >>> gdr (nspin);
58+ const std::size_t nrxx = rho.size ();
59+ for ( int is=0 ; is!=nspin; ++is )
60+ {
61+ std::vector<double > rhor (nrxx);
62+ #ifdef _OPENMP
63+ #pragma omp parallel for schedule(static, 1024)
64+ #endif
65+ for (std::size_t ir=0 ; ir<nrxx; ++ir)
66+ rhor[ir] = rho[ir*nspin+is];
67+ // ------------------------------------------
68+ // initialize the charge density array in reciprocal space
69+ // bring electron charge density from real space to reciprocal space
70+ // ------------------------------------------
71+ std::vector<std::complex <double >> rhog (chr->rhopw ->npw );
72+ chr->rhopw ->real2recip (rhor.data (), rhog.data ());
73+
74+ // -------------------------------------------
75+ // compute the gradient of charge density and
76+ // store the gradient in gdr[is]
77+ // -------------------------------------------
78+ gdr[is].resize (nrxx);
79+ XC_Functional::grad_rho (rhog.data (), gdr[is].data (), chr->rhopw , tpiba);
80+ } // end for(is)
81+ return gdr;
82+ }
83+
84+ // converting grho (abacus=>libxc)
85+ std::vector<double > XC_Functional_Libxc::convert_sigma (
86+ const std::vector<std::vector<ModuleBase::Vector3<double >>> &gdr)
87+ {
88+ const std::size_t nspin = gdr.size ();
89+ assert (nspin>0 );
90+ const std::size_t nrxx = gdr[0 ].size ();
91+ for (std::size_t is=1 ; is<nspin; ++is)
92+ assert (nrxx==gdr[is].size ());
93+
94+ std::vector<double > sigma ( nrxx * ((1 ==nspin)?1 :3 ) );
95+ if ( 1 ==nspin )
96+ {
97+ #ifdef _OPENMP
98+ #pragma omp parallel for schedule(static, 1024)
99+ #endif
100+ for ( std::size_t ir=0 ; ir<nrxx; ++ir )
101+ sigma[ir] = gdr[0 ][ir]*gdr[0 ][ir];
102+ }
103+ else
104+ {
105+ #ifdef _OPENMP
106+ #pragma omp parallel for schedule(static, 256)
107+ #endif
108+ for ( std::size_t ir=0 ; ir<nrxx; ++ir )
109+ {
110+ sigma[ir*3 ] = gdr[0 ][ir]*gdr[0 ][ir];
111+ sigma[ir*3 +1 ] = gdr[0 ][ir]*gdr[1 ][ir];
112+ sigma[ir*3 +2 ] = gdr[1 ][ir]*gdr[1 ][ir];
113+ }
114+ }
115+ return sigma;
116+ }
117+
118+ // sgn for threshold mask
119+ std::vector<double > XC_Functional_Libxc::cal_sgn (
120+ const double rho_threshold,
121+ const double grho_threshold,
122+ const xc_func_type &func,
123+ const int nspin,
124+ const std::size_t nrxx,
125+ const std::vector<double > &rho,
126+ const std::vector<double > &sigma)
127+ {
128+ std::vector<double > sgn (nrxx*nspin, 1.0 );
129+ // in the case of GGA correlation for polarized case,
130+ // a cutoff for grho is required to ensure that libxc gives reasonable results
131+ if (nspin==2 && func.info ->family != XC_FAMILY_LDA && func.info ->kind ==XC_CORRELATION)
132+ {
133+ #ifdef _OPENMP
134+ #pragma omp parallel for schedule(static, 512)
135+ #endif
136+ for ( int ir=0 ; ir<nrxx; ++ir )
137+ {
138+ if ( rho[ir*2 ]<rho_threshold || std::sqrt (std::abs (sigma[ir*3 ]))<grho_threshold )
139+ sgn[ir*2 ] = 0.0 ;
140+ if ( rho[ir*2 +1 ]<rho_threshold || std::sqrt (std::abs (sigma[ir*3 +2 ]))<grho_threshold )
141+ sgn[ir*2 +1 ] = 0.0 ;
142+ }
143+ }
144+ return sgn;
145+ }
146+
147+ // converting etxc from exc (libxc=>abacus)
148+ double XC_Functional_Libxc::convert_etxc (
149+ const int nspin,
150+ const std::size_t nrxx,
151+ const std::vector<double > &sgn,
152+ const std::vector<double > &rho,
153+ std::vector<double > exc)
154+ {
155+ double etxc = 0.0 ;
156+ #ifdef _OPENMP
157+ #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256)
158+ #endif
159+ for ( int is=0 ; is<nspin; ++is )
160+ for ( int ir=0 ; ir<nrxx; ++ir )
161+ etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is];
162+ return etxc;
163+ }
164+
165+ // converting etxc from exc (libxc=>abacus)
166+ double XC_Functional_Libxc::convert_vtxc_v (
167+ const xc_func_type &func,
168+ const int nspin,
169+ const std::size_t nrxx,
170+ const std::vector<double > &sgn,
171+ const std::vector<double > &rho,
172+ const std::vector<std::vector<ModuleBase::Vector3<double >>> &gdr,
173+ const std::vector<double > &vrho,
174+ const std::vector<double > &vsigma,
175+ const double tpiba,
176+ const Charge* const chr,
177+ ModuleBase::matrix &v)
178+ {
179+ assert (v.nr ==nspin);
180+ assert (v.nc ==nrxx);
181+
182+ double vtxc = 0.0 ;
183+
184+ #ifdef _OPENMP
185+ #pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256)
186+ #endif
187+ for ( int is=0 ; is<nspin; ++is )
188+ {
189+ for ( std::size_t ir=0 ; ir<nrxx; ++ir )
190+ {
191+ const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is];
192+ v (is,ir) += v_tmp;
193+ vtxc += v_tmp * rho[ir*nspin+is];
194+ }
195+ }
196+
197+ if (func.info ->family == XC_FAMILY_GGA || func.info ->family == XC_FAMILY_HYB_GGA)
198+ {
199+ const std::vector<std::vector<double >> dh = XC_Functional_Libxc::cal_dh (nspin, nrxx, sgn, gdr, vsigma, tpiba, chr);
200+
201+ double rvtxc = 0.0 ;
202+ #ifdef _OPENMP
203+ #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256)
204+ #endif
205+ for ( int is=0 ; is<nspin; ++is )
206+ {
207+ for ( std::size_t ir=0 ; ir<nrxx; ++ir )
208+ {
209+ rvtxc += dh[is][ir] * rho[ir*nspin+is];
210+ v (is,ir) -= dh[is][ir];
211+ }
212+ }
213+
214+ vtxc -= rvtxc;
215+ } // end if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA))
216+
217+ return vtxc;
218+ }
219+
220+
221+ // dh for gga v
222+ std::vector<std::vector<double >> XC_Functional_Libxc::cal_dh (
223+ const int nspin,
224+ const std::size_t nrxx,
225+ const std::vector<double > &sgn,
226+ const std::vector<std::vector<ModuleBase::Vector3<double >>> &gdr,
227+ const std::vector<double > &vsigma,
228+ const double tpiba,
229+ const Charge* const chr)
230+ {
231+ std::vector<std::vector<ModuleBase::Vector3<double >>> h (
232+ nspin,
233+ std::vector<ModuleBase::Vector3<double >>(nrxx) );
234+ if ( 1 ==nspin )
235+ {
236+ #ifdef _OPENMP
237+ #pragma omp parallel for schedule(static, 1024)
238+ #endif
239+ for ( std::size_t ir=0 ; ir<nrxx; ++ir )
240+ h[0 ][ir] = 2.0 * gdr[0 ][ir] * vsigma[ir] * 2.0 * sgn[ir];
241+ }
242+ else
243+ {
244+ #ifdef _OPENMP
245+ #pragma omp parallel for schedule(static, 1024)
246+ #endif
247+ for ( std::size_t ir=0 ; ir< nrxx; ++ir )
248+ {
249+ h[0 ][ir] = 2.0 * (gdr[0 ][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0
250+ + gdr[1 ][ir] * vsigma[ir*3 +1 ] * sgn[ir*2 ] * sgn[ir*2 +1 ]);
251+ h[1 ][ir] = 2.0 * (gdr[1 ][ir] * vsigma[ir*3 +2 ] * sgn[ir*2 +1 ] * 2.0
252+ + gdr[0 ][ir] * vsigma[ir*3 +1 ] * sgn[ir*2 ] * sgn[ir*2 +1 ]);
253+ }
254+ }
255+
256+ // define two dimensional array dh [ nspin, nrxx ]
257+ std::vector<std::vector<double >> dh (nspin, std::vector<double >(nrxx));
258+ for ( int is=0 ; is!=nspin; ++is )
259+ XC_Functional::grad_dot ( h[is].data (), dh[is].data (), chr->rhopw , tpiba);
260+
261+ return dh;
262+ }
263+
264+
265+ // convert v for NSPIN=4
266+ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4 (
267+ const std::size_t nrxx,
268+ const Charge* const chr,
269+ const std::vector<double > &amag,
270+ const ModuleBase::matrix &v)
271+ {
272+ assert (GlobalV::NSPIN==4 );
273+ constexpr double vanishing_charge = 1.0e-10 ;
274+ ModuleBase::matrix v_nspin4 (GlobalV::NSPIN, nrxx);
275+ for ( int ir=0 ; ir<nrxx; ++ir )
276+ v_nspin4 (0 ,ir) = 0.5 * (v (0 ,ir)+v (1 ,ir));
277+ if (GlobalV::DOMAG || GlobalV::DOMAG_Z)
278+ {
279+ for ( int ir=0 ; ir<nrxx; ++ir )
280+ {
281+ if ( amag[ir] > vanishing_charge )
282+ {
283+ const double vs = 0.5 * (v (0 ,ir)-v (1 ,ir));
284+ for (int ipol=1 ; ipol<GlobalV::NSPIN; ++ipol)
285+ v_nspin4 (ipol,ir) = vs * chr->rho [ipol][ir] / amag[ir];
286+ }
287+ }
288+ }
289+ return v_nspin4;
290+ }
291+
292+ #endif
0 commit comments