Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates a functional from libxc and its derivatives
10 : !> \note
11 : !> LibXC:
12 : !> (Marques, Oliveira, Burnus, CPC 183, 2272 (2012)).
13 : !>
14 : !> WARNING: In the subroutine libxc_spin_polarized_calc, it could be that the
15 : !> ordering for the 1st index of v2lapltau, v2rholapl, v2rhotau,
16 : !> v2sigmalapl and v2sigmatau is not correct. For the moment it does not
17 : !> matter since the calculation of the 2nd derivatives for meta-GGA
18 : !> functionals is not implemented in CP2K.
19 : !>
20 : !> \par History
21 : !> 01.2013 created [F. Tran]
22 : !> 07.2014 updates to versions 2.1 [JGH]
23 : !> 08.2015 refactoring [A. Gloess (agloess)]
24 : !> 01.2018 refactoring [A. Gloess (agloess)]
25 : !> 10.2018/04.2019 added hyb_mgga [S. Simko, included by F. Stein]
26 : !> \author F. Tran
27 : ! **************************************************************************************************
28 : #:include 'xc_gamma.fypp'
29 :
30 : #! Every staging buffer of a worker, so that the component count and the
31 : #! release below cannot fall out of step with the allocation.
32 : #:set worker_1d = ["exc", "nd", "nda", "ndb"]
33 : #:set worker_2d = ["rho", "sigma", "lapl", "tau", "vrho", "vsigma", "vlapl", "vtau", &
34 : & "v2rho2", "v2rhosigma", "v2sigma2", "v2rholapl", "v2rhotau", &
35 : & "v2sigmalapl", "v2sigmatau", "v2lapl2", "v2lapltau", "v2tau2", &
36 : & "v3rho3", "v3rho2sigma", "v3rhosigma2", "v3sigma3"] + v3_new
37 :
38 : MODULE xc_libxc
39 : USE bibliography, ONLY: Lehtola2018, &
40 : Marques2012, &
41 : cite_reference
42 : USE input_section_types, ONLY: section_add_keyword, &
43 : section_add_subsection, &
44 : section_create, &
45 : section_release, &
46 : section_type, &
47 : section_vals_type, &
48 : section_vals_val_get
49 : USE kinds, ONLY: default_string_length, &
50 : dp
51 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
52 : xc_dset_get_derivative
53 : USE xc_derivative_types, ONLY: xc_derivative_get, &
54 : xc_derivative_type
55 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
56 : USE xc_rho_set_types, ONLY: xc_rho_set_get, &
57 : xc_rho_set_type
58 : #if defined (__LIBXC)
59 : USE input_keyword_types, ONLY: keyword_create, &
60 : keyword_release, &
61 : keyword_type
62 : USE iso_c_binding, ONLY: C_SIZE_T, C_INT, C_DOUBLE
63 : USE xc_derivative_desc, ONLY: &
64 : deriv_rho, deriv_rhoa, deriv_rhob, &
65 : deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
66 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, &
67 : deriv_gamma, deriv_gamma_aa, deriv_gamma_ab, deriv_gamma_bb
68 : USE xc_libxc_wrap, ONLY: xc_f03_func_t, &
69 : xc_f03_func_init, &
70 : xc_f03_func_end, &
71 : xc_f03_func_info_t, &
72 : xc_f03_functional_get_name, &
73 : xc_f03_func_get_info, &
74 : xc_f03_func_info_get_family, &
75 : xc_f03_func_info_get_kind, &
76 : xc_f03_func_info_get_n_ext_params, &
77 : xc_f03_func_info_get_name, &
78 : xc_f03_available_functional_numbers, &
79 : xc_f03_available_functional_names, &
80 : xc_f03_maximum_name_length, &
81 : xc_f03_number_of_functionals, &
82 : xc_f03_func_info_get_ext_params_name, &
83 : xc_f03_func_info_get_ext_params_description, &
84 : xc_f03_func_info_get_ext_params_default_value, &
85 : xc_f03_gga_exc, &
86 : xc_f03_gga_exc_vxc, &
87 : xc_f03_gga_exc_vxc_fxc, &
88 : xc_f03_gga_exc_vxc_fxc_kxc, &
89 : xc_f03_gga_fxc, &
90 : xc_f03_gga_vxc, &
91 : xc_f03_gga_vxc_fxc, &
92 : xc_f03_gga_vxc_fxc_kxc, &
93 : xc_f03_mgga_exc_vxc_fxc_kxc, &
94 : xc_f03_mgga_vxc_fxc_kxc, &
95 : xc_f03_lda, &
96 : xc_f03_lda_exc, &
97 : xc_f03_lda_exc_vxc, &
98 : xc_f03_lda_exc_vxc_fxc, &
99 : xc_f03_lda_fxc, &
100 : xc_f03_lda_kxc, &
101 : xc_f03_lda_vxc, &
102 : xc_f03_mgga, &
103 : xc_f03_mgga_exc, &
104 : xc_f03_mgga_exc_vxc, &
105 : xc_f03_mgga_fxc, &
106 : xc_f03_mgga_vxc, &
107 : xc_f03_mgga_vxc_fxc, &
108 : XC_POLARIZED, &
109 : XC_UNPOLARIZED, &
110 : XC_FAMILY_LDA, &
111 : XC_FAMILY_GGA, &
112 : XC_FAMILY_MGGA, &
113 : XC_FAMILY_HYB_LDA, &
114 : XC_FAMILY_HYB_GGA, &
115 : XC_FAMILY_HYB_MGGA, &
116 : XC_CORRELATION, &
117 : XC_EXCHANGE, &
118 : XC_EXCHANGE_CORRELATION, &
119 : XC_KINETIC, &
120 : xc_libxc_wrap_info_refs, &
121 : xc_libxc_wrap_version, &
122 : xc_libxc_wrap_library_reference, &
123 : xc_libxc_wrap_functional_get_number, &
124 : xc_libxc_wrap_info_needs_laplace, &
125 : xc_libxc_wrap_info_no_exc, &
126 : xc_libxc_wrap_set_thresholds, &
127 : xc_f03_func_set_ext_params, &
128 : xc_libxc_wrap_is_under_development, &
129 : xc_libxc_get_reference_length, &
130 : xc_libxc_check_functional
131 : #endif
132 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
133 : !$ omp_get_num_threads, &
134 : !$ omp_get_thread_num
135 :
136 : #include "../base/base_uses.f90"
137 :
138 : IMPLICIT NONE
139 : PRIVATE
140 :
141 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_libxc'
142 :
143 : PUBLIC :: libxc_spin_unpolarized_info, libxc_spin_unpolarized_eval, &
144 : libxc_spin_polarized_info, libxc_spin_polarized_eval, &
145 : libxc_version_info, libxc_library_reference, &
146 : libxc_get_reference_length, libxc_add_sections, &
147 : libxc_check_existence_in_libxc, libxc_release_workers
148 :
149 : #if defined (__LIBXC)
150 : ! Number of grid points handed to LibXC in a single call. LibXC amortizes its
151 : ! per-call work (argument checks, output initialization and, for functionals
152 : ! built by mixing components, a full allocate/free cycle of the component
153 : ! buffers) over the points in the call, so evaluating point by point paid that
154 : ! cost once per point. Blocking rather than passing the whole grid keeps the
155 : ! staging buffers inside the cache. How wide a point is depends on the
156 : ! functional and on the derivative order: a spin-polarized meta-GGA stages 232
157 : ! doubles per point at third order against twelve for a polarized LDA, a span
158 : ! of twenty, so the block is chosen from a byte budget rather than fixed in
159 : ! points. The timing is flat over a wide range of block lengths, so the budget
160 : ! only has to keep the buffers cache-resident.
161 : INTEGER, PARAMETER, PRIVATE :: libxc_block_bytes = 256*1024
162 : INTEGER, PARAMETER, PRIVATE :: libxc_block_min = 64
163 : INTEGER, PARAMETER, PRIVATE :: libxc_block_max = 512
164 :
165 : ! **************************************************************************************************
166 : !> \brief One LibXC functional object plus the buffers used to hand it a block of
167 : !> grid points. One of these is kept per OpenMP thread so that neither the
168 : !> functional object nor the buffers are rebuilt between batches.
169 : ! **************************************************************************************************
170 : TYPE libxc_worker_type
171 : TYPE(xc_f03_func_t) :: func = xc_f03_func_t()
172 : TYPE(xc_f03_func_info_t) :: info = xc_f03_func_info_t()
173 : LOGICAL :: is_init = .FALSE.
174 : ! number of points the staging buffers below were allocated for
175 : INTEGER :: nb = 0
176 : ! staging buffers, dimensioned (component, point)
177 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rho, sigma, lapl, tau
178 : ! the floored gradient norms the spin-polarized chain rule is written in
179 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: nd, nda, ndb
180 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: exc
181 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: vrho, vsigma, vlapl, vtau
182 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: v2rho2, v2rhosigma, v2sigma2, &
183 : v2rholapl, v2rhotau, v2sigmalapl, &
184 : v2sigmatau, v2lapl2, v2lapltau, v2tau2
185 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: v3rho3, v3rho2sigma, &
186 : v3rhosigma2, v3sigma3
187 : #:for arr in v3_new
188 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: ${arr}$
189 : #:endfor
190 : END TYPE libxc_worker_type
191 :
192 : ! **************************************************************************************************
193 : !> \brief All threads' workers for one distinct functional setup, together with
194 : !> the key that identifies that setup.
195 : ! **************************************************************************************************
196 : TYPE libxc_worker_set_type
197 : INTEGER :: func_id = -1
198 : INTEGER :: nspin = -1
199 : INTEGER :: family = -1
200 : ! highest derivative order the workers' buffers were allocated for
201 : INTEGER :: grad_deriv = -1
202 : LOGICAL :: has_laplace = .FALSE.
203 : LOGICAL :: no_exc = .FALSE.
204 : LOGICAL :: has_params = .FALSE.
205 : REAL(KIND=dp) :: epsilon_rho = -1.0_dp
206 : REAL(KIND=dp) :: epsilon_tau = -1.0_dp
207 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: params
208 : CHARACTER(LEN=128), DIMENSION(:), ALLOCATABLE :: param_names
209 : TYPE(libxc_worker_type), DIMENSION(:), ALLOCATABLE :: worker
210 : END TYPE libxc_worker_set_type
211 :
212 : ! Workers are acquired from a serial region only (see libxc_get_workers), so
213 : ! the cache itself needs no locking.
214 : TYPE(libxc_worker_set_type), DIMENSION(:), ALLOCATABLE, TARGET, PRIVATE, SAVE :: libxc_cache
215 :
216 : ! Stands in for the derivative arrays a given evaluation does not produce.
217 : ! Those arguments still have to be associated with something, but they are
218 : ! never read or written; pointing them at the density instead would make an
219 : ! output argument alias a read-only input. Grown on demand and kept for the
220 : ! run rather than allocated per evaluation, and like the cache above it is
221 : ! only ever touched from a serial region.
222 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET, PRIVATE, SAVE :: libxc_unused
223 : #endif
224 :
225 : CONTAINS
226 :
227 : ! **************************************************************************************************
228 : !> \brief This function checks whether a functional name belongs to LibXC
229 : !> \param libxc_params (possible) LibXC input section
230 : !> \return exists whether the functional exists in LibXC
231 : ! **************************************************************************************************
232 4324 : FUNCTION libxc_check_existence_in_libxc(libxc_params) RESULT(exists)
233 :
234 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
235 : LOGICAL :: exists
236 :
237 : #if defined (__LIBXC)
238 :
239 4324 : exists = xc_libxc_check_functional(libxc_params%section%name)
240 : #else
241 : MARK_USED(libxc_params)
242 : exists = .FALSE.
243 : #endif
244 :
245 4324 : END FUNCTION libxc_check_existence_in_libxc
246 :
247 : ! **************************************************************************************************
248 : !> \brief This function returns the maximum length of the reference string for a given LibXC functional
249 : !> \param libxc_params LibXC input section
250 : !> \param lsd spin polarized calculation
251 : !> \return maximum length of the string
252 : ! **************************************************************************************************
253 124 : FUNCTION libxc_get_reference_length(libxc_params, lsd) RESULT(length)
254 :
255 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
256 : LOGICAL, INTENT(IN) :: lsd
257 : INTEGER :: length
258 :
259 : #if defined (__LIBXC)
260 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_get_reference_length'
261 :
262 : CHARACTER(LEN=default_string_length) :: func_name
263 : INTEGER :: func_id, handle
264 : TYPE(xc_f03_func_t) :: xc_func
265 : TYPE(xc_f03_func_info_t) :: xc_info
266 :
267 124 : CALL timeset(routineN, handle)
268 :
269 124 : func_name = libxc_params%section%name
270 :
271 124 : func_id = xc_libxc_wrap_functional_get_number(func_name)
272 124 : IF (lsd) THEN
273 60 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
274 : ELSE
275 64 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
276 : END IF
277 124 : xc_info = xc_f03_func_get_info(xc_func)
278 :
279 124 : length = xc_libxc_get_reference_length(xc_info)
280 :
281 124 : CALL xc_f03_func_end(xc_func)
282 :
283 124 : CALL timestop(handle)
284 : #else
285 : MARK_USED(libxc_params)
286 : MARK_USED(lsd)
287 : length = 0
288 : CPABORT("In order to use LibXC you have to download and install it!")
289 : #endif
290 :
291 124 : END FUNCTION libxc_get_reference_length
292 :
293 : ! **************************************************************************************************
294 : !> \brief ...
295 : !> \param section ...
296 : ! **************************************************************************************************
297 13675 : SUBROUTINE libxc_add_sections(section)
298 :
299 : TYPE(section_type), POINTER, INTENT(IN) :: section
300 :
301 : #if defined (__LIBXC)
302 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_add_sections'
303 :
304 : TYPE(section_type), POINTER :: subsection
305 : TYPE(keyword_type), POINTER :: keyword
306 : INTEGER :: handle, no_func, len_name, ii, func_id, n_param, iparam
307 : REAL(KIND=C_DOUBLE) :: default_val
308 : CHARACTER(LEN=128) :: func_name, param_name, param_descr, description
309 : CHARACTER(LEN=2*default_string_length) :: warning
310 13675 : INTEGER(KIND=C_INT), DIMENSION(:), ALLOCATABLE :: func_ids
311 : TYPE(xc_f03_func_t) :: xc_func
312 : TYPE(xc_f03_func_info_t) :: xc_info
313 :
314 13675 : CALL timeset(routineN, handle)
315 :
316 13675 : CPASSERT(ASSOCIATED(section))
317 13675 : NULLIFY (subsection, keyword)
318 :
319 13675 : no_func = xc_f03_number_of_functionals()
320 13675 : len_name = xc_f03_maximum_name_length()
321 :
322 41025 : ALLOCATE (func_ids(no_func))
323 :
324 13675 : CALL xc_f03_available_functional_numbers(func_ids)
325 :
326 9750275 : DO ii = 1, no_func
327 :
328 9736600 : func_id = func_ids(ii)
329 9736600 : IF (ii > 1) THEN
330 9722925 : IF (func_id == func_ids(ii - 1)) CYCLE
331 : END IF
332 9599850 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
333 9599850 : xc_info = xc_f03_func_get_info(xc_func)
334 :
335 9599850 : func_name = xc_f03_functional_get_name(func_id)
336 9599850 : description = xc_f03_func_info_get_name(xc_info)
337 9599850 : n_param = xc_f03_func_info_get_n_ext_params(xc_info)
338 :
339 9599850 : NULLIFY (subsection)
340 : CALL section_create(subsection, __LOCATION__, name=TRIM(func_name), description=TRIM(description), &
341 9599850 : n_keywords=2 + n_param, n_subsections=0, repeats=.FALSE.)
342 :
343 9599850 : IF (description(1:1) == "_") THEN
344 : warning = " This parameter is an internal parameter of the functional. Changing this "// &
345 0 : "parameter effectively changes the functional."
346 : ELSE
347 9599850 : warning = " "
348 : END IF
349 :
350 9599850 : NULLIFY (keyword)
351 : CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
352 : description="Activates the functional."//TRIM(warning), &
353 9599850 : lone_keyword_l_val=.TRUE., default_l_val=.FALSE.)
354 9599850 : CALL section_add_keyword(subsection, keyword)
355 9599850 : CALL keyword_release(keyword)
356 :
357 : CALL keyword_create(keyword, __LOCATION__, name="SCALE", description="Scales this functional", &
358 9599850 : default_r_val=1.0_dp)
359 9599850 : CALL section_add_keyword(subsection, keyword)
360 9599850 : CALL keyword_release(keyword)
361 :
362 58487975 : DO iparam = 1, n_param
363 48888125 : param_name = xc_f03_func_info_get_ext_params_name(xc_info, iparam - 1)
364 48888125 : param_descr = xc_f03_func_info_get_ext_params_description(xc_info, iparam - 1)
365 48888125 : default_val = xc_f03_func_info_get_ext_params_default_value(xc_info, iparam - 1)
366 48888125 : NULLIFY (keyword)
367 : CALL keyword_create(keyword, __LOCATION__, name=TRIM(param_name), &
368 48888125 : description=TRIM(param_descr), default_r_val=default_val)
369 48888125 : CALL section_add_keyword(subsection, keyword)
370 58487975 : CALL keyword_release(keyword)
371 : END DO
372 :
373 9599850 : CALL section_add_subsection(section, subsection)
374 9599850 : CALL section_release(subsection)
375 :
376 9750275 : CALL xc_f03_func_end(xc_func)
377 :
378 : END DO
379 :
380 13675 : DEALLOCATE (func_ids)
381 :
382 13675 : CALL timestop(handle)
383 : #else
384 : MARK_USED(section)
385 :
386 : #endif
387 :
388 13675 : END SUBROUTINE libxc_add_sections
389 :
390 : ! **************************************************************************************************
391 : !> \brief info about the functional from libxc
392 : !> \param libxc_params input parameter (functional name, scaling and parameters)
393 : !> \param reference string with the reference of the actual functional
394 : !> \param shortform string with the shortform of the functional name
395 : !> \param needs the components needed by this functional are set to
396 : !> true (does not set the unneeded components to false)
397 : !> \param max_deriv maximum implemented derivative of the xc functional
398 : !> \param print_warn whether to print warning about development status of a functional
399 : !> \param func_name_override optional LibXC functional name overriding the section name
400 : !> \author F. Tran
401 : ! **************************************************************************************************
402 14444 : SUBROUTINE libxc_spin_unpolarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
403 : func_name_override)
404 :
405 : TYPE(section_vals_type), POINTER :: libxc_params
406 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
407 : TYPE(xc_rho_cflags_type), &
408 : INTENT(inout), OPTIONAL :: needs
409 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
410 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
411 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
412 :
413 : #if defined (__LIBXC)
414 : CHARACTER(LEN=128) :: s1, s2
415 : CHARACTER(LEN=default_string_length) :: func_name
416 : INTEGER :: func_id
417 : REAL(KIND=dp) :: func_scale
418 : TYPE(xc_f03_func_t) :: xc_func
419 : TYPE(xc_f03_func_info_t) :: xc_info
420 :
421 28800 : IF (PRESENT(func_name_override)) THEN
422 88 : func_name = func_name_override
423 88 : func_scale = 1.0_dp
424 : ELSE
425 14356 : func_name = libxc_params%section%name
426 14356 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
427 : END IF
428 :
429 14444 : CALL cite_reference(Marques2012)
430 14444 : CALL cite_reference(Lehtola2018)
431 :
432 14444 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
433 :
434 14444 : func_id = xc_libxc_wrap_functional_get_number(func_name)
435 14444 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
436 14444 : xc_info = xc_f03_func_get_info(xc_func)
437 :
438 14444 : s1 = xc_f03_func_info_get_name(xc_info)
439 9738 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
440 9738 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
441 2928 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
442 1448 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
443 330 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
444 : CASE default
445 14444 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
446 : END SELECT
447 14444 : IF (PRESENT(shortform)) THEN
448 64 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
449 : END IF
450 14444 : IF (PRESENT(reference)) THEN
451 64 : CALL xc_libxc_wrap_info_refs(xc_info, XC_UNPOLARIZED, func_scale, reference)
452 : END IF
453 14444 : IF (PRESENT(needs)) THEN
454 6212 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
455 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
456 6212 : needs%rho = .TRUE.
457 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
458 4832 : needs%rho = .TRUE.
459 4832 : needs%norm_drho = .TRUE.
460 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
461 3328 : needs%rho = .TRUE.
462 3328 : needs%norm_drho = .TRUE.
463 3328 : needs%tau = .TRUE.
464 3328 : needs%laplace_rho = xc_libxc_wrap_info_needs_laplace(xc_info)
465 : CASE default
466 14372 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
467 : END SELECT
468 : END IF
469 14444 : IF (PRESENT(max_deriv)) THEN
470 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
471 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
472 0 : max_deriv = 3
473 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
474 88 : max_deriv = 2
475 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
476 0 : max_deriv = 2
477 : CASE default
478 88 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
479 : END SELECT
480 : END IF
481 14444 : IF (PRESENT(print_warn)) THEN
482 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
483 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
484 : END IF
485 : END IF
486 :
487 14444 : CALL xc_f03_func_end(xc_func)
488 : #else
489 : MARK_USED(libxc_params)
490 : MARK_USED(reference)
491 : MARK_USED(shortform)
492 : MARK_USED(needs)
493 : MARK_USED(max_deriv)
494 : MARK_USED(print_warn)
495 : MARK_USED(func_name_override)
496 :
497 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
498 : "for a functional of the LibXC library, "// &
499 : "you have to download and install the library!")
500 : #endif
501 :
502 14444 : END SUBROUTINE libxc_spin_unpolarized_info
503 :
504 : ! **************************************************************************************************
505 : !> \brief info about the functional from libxc
506 : !> \param libxc_params input parameter (functional name, scaling and parameters)
507 : !> \param reference string with the reference of the actual functional
508 : !> \param shortform string with the shortform of the functional name
509 : !> \param needs the components needed by this functional are set to
510 : !> true (does not set the unneeded components to false)
511 : !> \param max_deriv maximum implemented derivative of the xc functional
512 : !> \param print_warn whether to print warning about development status of a functional
513 : !> \param func_name_override optional LibXC functional name overriding the section name
514 : !> \author F. Tran
515 : ! **************************************************************************************************
516 3944 : SUBROUTINE libxc_spin_polarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
517 : func_name_override)
518 :
519 : TYPE(section_vals_type), POINTER :: libxc_params
520 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
521 : TYPE(xc_rho_cflags_type), &
522 : INTENT(inout), OPTIONAL :: needs
523 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
524 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
525 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
526 :
527 : #if defined (__LIBXC)
528 : CHARACTER(LEN=128) :: s1, s2
529 : CHARACTER(LEN=default_string_length) :: func_name
530 : INTEGER :: func_id
531 : REAL(KIND=dp) :: func_scale
532 : TYPE(xc_f03_func_t) :: xc_func
533 : TYPE(xc_f03_func_info_t) :: xc_info
534 :
535 7880 : IF (PRESENT(func_name_override)) THEN
536 8 : func_name = func_name_override
537 8 : func_scale = 1.0_dp
538 : ELSE
539 3936 : func_name = libxc_params%section%name
540 3936 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
541 : END IF
542 :
543 3944 : CALL cite_reference(Marques2012)
544 3944 : CALL cite_reference(Lehtola2018)
545 :
546 3944 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
547 :
548 3944 : func_id = xc_libxc_wrap_functional_get_number(func_name)
549 3944 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
550 3944 : xc_info = xc_f03_func_get_info(xc_func)
551 :
552 3944 : s1 = xc_f03_func_info_get_name(xc_info)
553 2014 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
554 2014 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
555 1684 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
556 246 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
557 0 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
558 : CASE default
559 3944 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
560 : END SELECT
561 3944 : IF (PRESENT(shortform)) THEN
562 60 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
563 : END IF
564 3944 : IF (PRESENT(reference)) THEN
565 60 : CALL xc_libxc_wrap_info_refs(xc_info, XC_POLARIZED, func_scale, reference)
566 : END IF
567 3944 : IF (PRESENT(needs)) THEN
568 1404 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
569 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
570 1404 : needs%rho_spin = .TRUE.
571 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
572 674 : needs%rho_spin = .TRUE.
573 674 : needs%norm_drho = .TRUE.
574 674 : needs%norm_drho_spin = .TRUE.
575 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
576 1806 : needs%rho_spin = .TRUE.
577 1806 : needs%norm_drho = .TRUE.
578 1806 : needs%norm_drho_spin = .TRUE.
579 1806 : needs%tau_spin = .TRUE.
580 1806 : needs%laplace_rho_spin = xc_libxc_wrap_info_needs_laplace(xc_info)
581 : CASE default
582 3884 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
583 : END SELECT
584 : END IF
585 3944 : IF (PRESENT(max_deriv)) THEN
586 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
587 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
588 0 : max_deriv = 3
589 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
590 8 : max_deriv = 2
591 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
592 0 : max_deriv = 2
593 : CASE default
594 8 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
595 : END SELECT
596 : END IF
597 3944 : IF (PRESENT(print_warn)) THEN
598 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
599 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
600 : END IF
601 : END IF
602 :
603 3944 : CALL xc_f03_func_end(xc_func)
604 : #else
605 : MARK_USED(libxc_params)
606 : MARK_USED(reference)
607 : MARK_USED(shortform)
608 : MARK_USED(needs)
609 : MARK_USED(max_deriv)
610 : MARK_USED(print_warn)
611 : MARK_USED(func_name_override)
612 :
613 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are "// &
614 : "asking for a functional of the LibXC library, "// &
615 : "you have to download and install the library!")
616 : #endif
617 :
618 3944 : END SUBROUTINE libxc_spin_polarized_info
619 :
620 : ! **************************************************************************************************
621 : !> \brief info about the LibXC version
622 : !> \param version version string reported by the library at run time
623 : !> \param compiled_version version of the LibXC headers CP2K was compiled against
624 : !> \author A. Gloess (agloess)
625 : ! **************************************************************************************************
626 78 : SUBROUTINE libxc_version_info(version, compiled_version)
627 : CHARACTER(LEN=*), INTENT(OUT) :: version ! the string that is output
628 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: compiled_version
629 :
630 : #if defined (__LIBXC)
631 78 : CALL xc_libxc_wrap_version(version, compiled_version)
632 : #else
633 : version = "none"
634 : IF (PRESENT(compiled_version)) compiled_version = "none"
635 : CPABORT("In order to use libxc you need to download and install it")
636 : #endif
637 :
638 78 : END SUBROUTINE libxc_version_info
639 :
640 : ! **************************************************************************************************
641 : !> \brief Returns the citation LibXC asks for the library itself.
642 : !> \param reference bibliographic reference
643 : !> \param doi digital object identifier of that reference
644 : !> \author S. Lehtola
645 : ! **************************************************************************************************
646 78 : SUBROUTINE libxc_library_reference(reference, doi)
647 : CHARACTER(LEN=*), INTENT(OUT) :: reference, doi
648 :
649 : #if defined (__LIBXC)
650 78 : CALL xc_libxc_wrap_library_reference(reference, doi)
651 : #else
652 : reference = "none"
653 : doi = "none"
654 : CPABORT("In order to use libxc you need to download and install it")
655 : #endif
656 :
657 78 : END SUBROUTINE libxc_library_reference
658 :
659 : ! **************************************************************************************************
660 : !> \brief Reads the external parameters of a LibXC functional from its input section.
661 : !> \param libxc_params LibXC input section
662 : !> \param param_names names of the external parameters, as reported by LibXC
663 : !> \param params values read from the input
664 : !> \author S. Lehtola
665 : ! **************************************************************************************************
666 : #if defined (__LIBXC)
667 13624 : SUBROUTINE libxc_read_params(libxc_params, param_names, params)
668 : TYPE(section_vals_type), INTENT(IN), POINTER :: libxc_params
669 : CHARACTER(LEN=128), DIMENSION(:), INTENT(IN) :: param_names
670 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: params
671 :
672 : INTEGER :: i
673 :
674 86442 : DO i = 1, SIZE(param_names)
675 86442 : CALL section_vals_val_get(libxc_params, TRIM(param_names(i)), r_val=params(i))
676 : END DO
677 :
678 13624 : END SUBROUTINE libxc_read_params
679 :
680 : ! **************************************************************************************************
681 : !> \brief Allocates the staging buffers of a single worker, sizing the block so
682 : !> that they stay within the cache budget.
683 : !> \param worker the worker
684 : !> \param family LibXC functional family
685 : !> \param nspin XC_UNPOLARIZED or XC_POLARIZED
686 : !> \param grad_deriv highest derivative order the buffers must hold
687 : !> \author S. Lehtola
688 : ! **************************************************************************************************
689 810 : SUBROUTINE libxc_worker_alloc(worker, family, nspin, grad_deriv)
690 : TYPE(libxc_worker_type), INTENT(INOUT) :: worker
691 : INTEGER, INTENT(IN) :: family, nspin, grad_deriv
692 :
693 : INTEGER :: nb, ncomp
694 :
695 : ! How many doubles a point occupies spans a factor of twenty over the
696 : ! families, spin cases and derivative orders, so a block length fixed in
697 : ! points would size the buffers anywhere between 49 KiB and 950 KiB. Derive
698 : ! it from the footprint instead: build a one-point worker, ask it how wide
699 : ! a point is, and cut the real block to the byte budget. Taking the count
700 : ! from the buffers themselves is what stops it drifting away from the
701 : ! allocation it is meant to describe.
702 34020 : BLOCK
703 810 : TYPE(libxc_worker_type) :: probe
704 810 : CALL libxc_worker_alloc_buffers(probe, family, nspin, grad_deriv, 1)
705 34032 : ncomp = libxc_worker_ncomp(probe)
706 : END BLOCK
707 : nb = MIN(libxc_block_max, &
708 810 : MAX(libxc_block_min, libxc_block_bytes/(ncomp*(STORAGE_SIZE(0.0_dp)/8))))
709 :
710 810 : CALL libxc_worker_alloc_buffers(worker, family, nspin, grad_deriv, nb)
711 :
712 810 : END SUBROUTINE libxc_worker_alloc
713 :
714 : ! **************************************************************************************************
715 : !> \brief Number of buffer components a single grid point occupies in a worker.
716 : !> \param worker the worker
717 : !> \return components per point, summed over the buffers it actually carries
718 : !> \author S. Lehtola
719 : ! **************************************************************************************************
720 810 : FUNCTION libxc_worker_ncomp(worker) RESULT(ncomp)
721 : TYPE(libxc_worker_type), INTENT(IN) :: worker
722 : INTEGER :: ncomp
723 :
724 810 : ncomp = 0
725 : #:for arr in worker_1d
726 810 : IF (ALLOCATED(worker%${arr}$)) ncomp = ncomp + 1
727 : #:endfor
728 : #:for arr in worker_2d
729 810 : IF (ALLOCATED(worker%${arr}$)) ncomp = ncomp + SIZE(worker%${arr}$, 1)
730 : #:endfor
731 :
732 810 : END FUNCTION libxc_worker_ncomp
733 :
734 : ! **************************************************************************************************
735 : !> \brief Releases the staging buffers of a worker, keeping its functional object.
736 : !> \param worker the worker
737 : !> \author S. Lehtola
738 : ! **************************************************************************************************
739 312 : SUBROUTINE libxc_worker_free_buffers(worker)
740 : TYPE(libxc_worker_type), INTENT(INOUT) :: worker
741 :
742 : #:for arr in worker_1d + worker_2d
743 312 : IF (ALLOCATED(worker%${arr}$)) DEALLOCATE (worker%${arr}$)
744 : #:endfor
745 312 : worker%nb = 0
746 :
747 312 : END SUBROUTINE libxc_worker_free_buffers
748 :
749 : ! **************************************************************************************************
750 : !> \brief Allocates the staging buffers of a single worker for a given block length.
751 : !> \param worker the worker
752 : !> \param family LibXC functional family
753 : !> \param nspin XC_UNPOLARIZED or XC_POLARIZED
754 : !> \param grad_deriv highest derivative order the buffers must hold
755 : !> \param nb number of grid points per block
756 : !> \author S. Lehtola
757 : !> \note The buffers follow LibXC's own layout, i.e. the leading dimension is the
758 : !> number of components of the quantity and the points run along the second
759 : !> dimension. They are allocated once per worker and reused for every batch.
760 : !> Only the orders this evaluation asks for are allocated: the third
761 : !> derivatives alone are 130 components per point for a spin-polarized
762 : !> meta-GGA, which a ground-state calculation would never fill.
763 : ! **************************************************************************************************
764 1620 : SUBROUTINE libxc_worker_alloc_buffers(worker, family, nspin, grad_deriv, nb)
765 : TYPE(libxc_worker_type), INTENT(INOUT) :: worker
766 : INTEGER, INTENT(IN) :: family, nspin, grad_deriv, nb
767 :
768 : INTEGER :: nl, np, ns
769 : LOGICAL :: is_gga, is_mgga, pol
770 :
771 1620 : pol = (nspin == XC_POLARIZED)
772 : ! number of components per point, cf. the LibXC manual
773 1620 : IF (pol) THEN
774 : np = 2 ! rho, vrho, lapl, tau, vlapl, vtau
775 : ns = 3 ! sigma, vsigma, v2rho2, v2lapl2, v2tau2
776 : nl = 4 ! v2rholapl, v2rhotau, v2lapltau, v3rho3
777 : ELSE
778 1164 : np = 1
779 1164 : ns = 1
780 1164 : nl = 1
781 : END IF
782 :
783 1620 : is_gga = .FALSE.
784 1620 : is_mgga = .FALSE.
785 2140 : SELECT CASE (family)
786 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
787 520 : is_gga = .TRUE.
788 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
789 484 : is_gga = .TRUE.
790 1620 : is_mgga = .TRUE.
791 : END SELECT
792 :
793 1620 : worker%nb = nb
794 :
795 : ! The energy density and the first derivatives are always staged: LibXC
796 : ! writes them whatever higher order was asked for.
797 6480 : ALLOCATE (worker%rho(np, nb))
798 4860 : ALLOCATE (worker%exc(nb))
799 4860 : ALLOCATE (worker%vrho(np, nb))
800 2940 : IF (grad_deriv >= 2) ALLOCATE (worker%v2rho2(ns, nb))
801 1788 : IF (grad_deriv >= 3) ALLOCATE (worker%v3rho3(nl, nb))
802 :
803 1620 : IF (is_gga) THEN
804 4016 : ALLOCATE (worker%sigma(ns, nb))
805 2008 : ALLOCATE (worker%vsigma(ns, nb))
806 1004 : IF (pol) THEN
807 1008 : ALLOCATE (worker%nd(nb), worker%nda(nb), worker%ndb(nb))
808 : END IF
809 1004 : IF (grad_deriv >= 2) THEN
810 1412 : ALLOCATE (worker%v2rhosigma(MERGE(6, 1, pol), nb))
811 924 : ALLOCATE (worker%v2sigma2(MERGE(6, 1, pol), nb))
812 : END IF
813 1004 : IF (grad_deriv >= 3) THEN
814 180 : ALLOCATE (worker%v3rho2sigma(MERGE(9, 1, pol), nb))
815 180 : ALLOCATE (worker%v3rhosigma2(MERGE(12, 1, pol), nb))
816 180 : ALLOCATE (worker%v3sigma3(MERGE(10, 1, pol), nb))
817 : END IF
818 : END IF
819 :
820 1620 : IF (is_mgga) THEN
821 1452 : ALLOCATE (worker%lapl(np, nb))
822 1452 : ALLOCATE (worker%tau(np, nb))
823 1452 : ALLOCATE (worker%vlapl(np, nb))
824 1452 : ALLOCATE (worker%vtau(np, nb))
825 484 : IF (grad_deriv >= 2) THEN
826 624 : ALLOCATE (worker%v2rholapl(nl, nb))
827 468 : ALLOCATE (worker%v2rhotau(nl, nb))
828 736 : ALLOCATE (worker%v2sigmalapl(MERGE(6, 1, pol), nb))
829 468 : ALLOCATE (worker%v2sigmatau(MERGE(6, 1, pol), nb))
830 624 : ALLOCATE (worker%v2lapl2(ns, nb))
831 468 : ALLOCATE (worker%v2lapltau(nl, nb))
832 312 : ALLOCATE (worker%v2tau2(ns, nb))
833 : END IF
834 484 : IF (grad_deriv >= 3) THEN
835 : #:for arr in v3_new
836 948 : ALLOCATE (worker%${arr}$ (MERGE(${mgga_ncomp[arr]}$, 1, pol), nb))
837 : #:endfor
838 : END IF
839 : ! LibXC leaves vlapl untouched for functionals that do not use the
840 : ! Laplacian, and CP2K then ignores it; zero it once so it never holds
841 : ! uninitialized memory.
842 276904 : worker%vlapl = 0.0_dp
843 : END IF
844 :
845 1620 : END SUBROUTINE libxc_worker_alloc_buffers
846 :
847 : ! **************************************************************************************************
848 : !> \brief Makes sure the placeholder for unproduced derivatives holds npoints values.
849 : !> \param npoints number of grid points in this evaluation
850 : !> \author S. Lehtola
851 : ! **************************************************************************************************
852 21842 : SUBROUTINE libxc_reserve_unused(npoints)
853 : INTEGER, INTENT(IN) :: npoints
854 :
855 21842 : IF (ALLOCATED(libxc_unused)) THEN
856 21546 : IF (SIZE(libxc_unused) >= npoints) RETURN
857 2 : DEALLOCATE (libxc_unused)
858 : END IF
859 894 : ALLOCATE (libxc_unused(npoints))
860 :
861 : END SUBROUTINE libxc_reserve_unused
862 :
863 : ! **************************************************************************************************
864 : !> \brief Returns the cache slot holding the per-thread workers for a functional.
865 : !> \param func_name LibXC functional name
866 : !> \param nspin XC_UNPOLARIZED or XC_POLARIZED
867 : !> \param libxc_params LibXC input section, or null when no parameters are to be set
868 : !> \param epsilon_rho density cutoff
869 : !> \param epsilon_tau kinetic energy density cutoff
870 : !> \param grad_deriv highest derivative order this evaluation needs
871 : !> \return index into libxc_cache
872 : !> \author S. Lehtola
873 : !> \note Must be called from a serial region: it is the only writer of the cache,
874 : !> and keeping it serial is what makes the cache lock-free. The workers it
875 : !> returns are then used concurrently, one per thread.
876 : ! **************************************************************************************************
877 21842 : FUNCTION libxc_get_workers(func_name, nspin, libxc_params, epsilon_rho, epsilon_tau, &
878 : grad_deriv) RESULT(idx)
879 : CHARACTER(LEN=*), INTENT(IN) :: func_name
880 : INTEGER, INTENT(IN) :: nspin
881 : TYPE(section_vals_type), INTENT(IN), POINTER :: libxc_params
882 : REAL(KIND=dp), INTENT(IN) :: epsilon_rho, epsilon_tau
883 : INTEGER, INTENT(IN) :: grad_deriv
884 : INTEGER :: idx
885 :
886 : CHARACTER(LEN=*), PARAMETER :: routineN = 'libxc_get_workers'
887 :
888 : INTEGER :: func_id, handle, i, ithread, n_params, &
889 : nthreads
890 : LOGICAL :: has_params
891 21842 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: params
892 : TYPE(libxc_worker_set_type), ALLOCATABLE, &
893 21842 : DIMENSION(:) :: tmp_cache
894 :
895 21842 : CALL timeset(routineN, handle)
896 :
897 21842 : func_id = xc_libxc_wrap_functional_get_number(func_name)
898 21842 : has_params = ASSOCIATED(libxc_params)
899 :
900 21842 : nthreads = 1
901 21842 : !$ nthreads = omp_get_max_threads()
902 :
903 21842 : IF (.NOT. ALLOCATED(libxc_cache)) ALLOCATE (libxc_cache(0))
904 :
905 21842 : idx = 0
906 30380 : DO i = 1, SIZE(libxc_cache)
907 : IF (libxc_cache(i)%func_id == func_id .AND. &
908 29882 : libxc_cache(i)%nspin == nspin .AND. &
909 498 : (libxc_cache(i)%has_params .EQV. has_params)) THEN
910 : IF (libxc_cache(i)%epsilon_rho == epsilon_rho .AND. &
911 21344 : libxc_cache(i)%epsilon_tau == epsilon_tau .AND. &
912 : SIZE(libxc_cache(i)%worker) >= nthreads) THEN
913 : idx = i
914 : EXIT
915 : END IF
916 : END IF
917 : END DO
918 :
919 : ! A set built for a lower derivative order carries buffers that are too
920 : ! narrow for this call. Widening them in place keeps one set per functional
921 : ! setup: the SCF fills the cache at first order and a later response
922 : ! calculation grows that same slot to third order.
923 21842 : IF (idx /= 0) THEN
924 21344 : IF (libxc_cache(idx)%grad_deriv < grad_deriv) THEN
925 624 : DO ithread = 1, SIZE(libxc_cache(idx)%worker)
926 312 : CALL libxc_worker_free_buffers(libxc_cache(idx)%worker(ithread))
927 : CALL libxc_worker_alloc(libxc_cache(idx)%worker(ithread), &
928 624 : libxc_cache(idx)%family, nspin, grad_deriv)
929 : END DO
930 312 : libxc_cache(idx)%grad_deriv = grad_deriv
931 : END IF
932 : END IF
933 :
934 : IF (idx == 0) THEN
935 : ! Grow the cache by one slot and build the workers for it.
936 498 : CALL MOVE_ALLOC(libxc_cache, tmp_cache)
937 2322 : ALLOCATE (libxc_cache(SIZE(tmp_cache) + 1))
938 828 : DO i = 1, SIZE(tmp_cache)
939 828 : CALL libxc_move_worker_set(tmp_cache(i), libxc_cache(i))
940 : END DO
941 1326 : DEALLOCATE (tmp_cache)
942 498 : idx = SIZE(libxc_cache)
943 :
944 : ASSOCIATE (ws => libxc_cache(idx))
945 498 : ws%func_id = func_id
946 498 : ws%nspin = nspin
947 498 : ws%has_params = has_params
948 498 : ws%epsilon_rho = epsilon_rho
949 498 : ws%epsilon_tau = epsilon_tau
950 498 : ws%grad_deriv = grad_deriv
951 1992 : ALLOCATE (ws%worker(nthreads))
952 :
953 996 : DO ithread = 1, nthreads
954 498 : CALL xc_f03_func_init(ws%worker(ithread)%func, func_id, nspin)
955 498 : ws%worker(ithread)%info = xc_f03_func_get_info(ws%worker(ithread)%func)
956 996 : ws%worker(ithread)%is_init = .TRUE.
957 : END DO
958 :
959 498 : ws%family = xc_f03_func_info_get_family(ws%worker(1)%info)
960 498 : ws%has_laplace = xc_libxc_wrap_info_needs_laplace(ws%worker(1)%info)
961 498 : ws%no_exc = xc_libxc_wrap_info_no_exc(ws%worker(1)%info)
962 :
963 : ! Remember the parameter names so that later evaluations can re-read
964 : ! the input without another functional object to ask for them.
965 498 : n_params = 0
966 498 : IF (has_params) n_params = xc_f03_func_info_get_n_ext_params(ws%worker(1)%info)
967 1392 : ALLOCATE (ws%params(n_params))
968 1392 : ALLOCATE (ws%param_names(n_params))
969 2562 : DO i = 1, n_params
970 2562 : ws%param_names(i) = xc_f03_func_info_get_ext_params_name(ws%worker(1)%info, i - 1)
971 : END DO
972 498 : IF (n_params > 0) THEN
973 396 : CALL libxc_read_params(libxc_params, ws%param_names, ws%params)
974 792 : DO ithread = 1, nthreads
975 792 : CALL xc_f03_func_set_ext_params(ws%worker(ithread)%func, ws%params)
976 : END DO
977 : END IF
978 :
979 1494 : DO ithread = 1, nthreads
980 : CALL xc_libxc_wrap_set_thresholds(ws%worker(ithread)%func, ws%worker(ithread)%info, &
981 498 : epsilon_rho, epsilon_tau)
982 996 : CALL libxc_worker_alloc(ws%worker(ithread), ws%family, nspin, grad_deriv)
983 : END DO
984 : END ASSOCIATE
985 : ELSE
986 : ! Cached slot: the input section may in principle carry different
987 : ! parameter values than the ones the workers were last set up with, so
988 : ! re-read them and push them through only when they actually changed.
989 : ASSOCIATE (ws => libxc_cache(idx))
990 21344 : n_params = SIZE(ws%params)
991 21344 : IF (n_params > 0) THEN
992 39684 : ALLOCATE (params(n_params))
993 13228 : CALL libxc_read_params(libxc_params, ws%param_names, params)
994 83982 : IF (ANY(params /= ws%params)) THEN
995 0 : ws%params(:) = params(:)
996 0 : DO ithread = 1, SIZE(ws%worker)
997 0 : CALL xc_f03_func_set_ext_params(ws%worker(ithread)%func, ws%params)
998 : ! re-assert the cutoffs: setting parameters runs the
999 : ! functional's own callback, which is free to rebuild the
1000 : ! component functionals a mixed functional is made of
1001 : CALL xc_libxc_wrap_set_thresholds(ws%worker(ithread)%func, &
1002 : ws%worker(ithread)%info, &
1003 0 : ws%epsilon_rho, ws%epsilon_tau)
1004 : END DO
1005 : END IF
1006 13228 : DEALLOCATE (params)
1007 : END IF
1008 : END ASSOCIATE
1009 : END IF
1010 :
1011 21842 : CALL timestop(handle)
1012 :
1013 21842 : END FUNCTION libxc_get_workers
1014 :
1015 : ! **************************************************************************************************
1016 : !> \brief Moves a worker set to a new cache slot without touching LibXC.
1017 : !> \param from source slot, left empty
1018 : !> \param to destination slot
1019 : !> \author S. Lehtola
1020 : !> \note Used when the cache array grows. The LibXC objects are opaque handles, so
1021 : !> they survive being moved; only the Fortran allocatables need transferring.
1022 : ! **************************************************************************************************
1023 330 : SUBROUTINE libxc_move_worker_set(from, to)
1024 : TYPE(libxc_worker_set_type), INTENT(INOUT) :: from, to
1025 :
1026 330 : to%func_id = from%func_id
1027 330 : to%nspin = from%nspin
1028 330 : to%family = from%family
1029 330 : to%has_laplace = from%has_laplace
1030 330 : to%no_exc = from%no_exc
1031 330 : to%has_params = from%has_params
1032 330 : to%epsilon_rho = from%epsilon_rho
1033 330 : to%epsilon_tau = from%epsilon_tau
1034 330 : IF (ALLOCATED(from%params)) CALL MOVE_ALLOC(from%params, to%params)
1035 330 : IF (ALLOCATED(from%param_names)) CALL MOVE_ALLOC(from%param_names, to%param_names)
1036 330 : IF (ALLOCATED(from%worker)) CALL MOVE_ALLOC(from%worker, to%worker)
1037 :
1038 330 : END SUBROUTINE libxc_move_worker_set
1039 : #endif
1040 :
1041 : ! **************************************************************************************************
1042 : !> \brief Destroys all cached LibXC functional objects.
1043 : !> \author S. Lehtola
1044 : !> \note Call once at the end of a run. The workers are deliberately kept alive
1045 : !> across SCF steps and across calculations, so nothing else releases them.
1046 : ! **************************************************************************************************
1047 11299 : SUBROUTINE libxc_release_workers()
1048 :
1049 : #if defined (__LIBXC)
1050 : INTEGER :: i, ithread
1051 :
1052 11299 : IF (.NOT. ALLOCATED(libxc_cache)) RETURN
1053 :
1054 794 : DO i = 1, SIZE(libxc_cache)
1055 498 : IF (.NOT. ALLOCATED(libxc_cache(i)%worker)) CYCLE
1056 1292 : DO ithread = 1, SIZE(libxc_cache(i)%worker)
1057 996 : IF (libxc_cache(i)%worker(ithread)%is_init) THEN
1058 498 : CALL xc_f03_func_end(libxc_cache(i)%worker(ithread)%func)
1059 498 : libxc_cache(i)%worker(ithread)%is_init = .FALSE.
1060 : END IF
1061 : END DO
1062 : END DO
1063 1292 : DEALLOCATE (libxc_cache)
1064 296 : IF (ALLOCATED(libxc_unused)) DEALLOCATE (libxc_unused)
1065 : #endif
1066 :
1067 : END SUBROUTINE libxc_release_workers
1068 :
1069 : ! **************************************************************************************************
1070 : !> \brief evaluates the functional from libxc
1071 : !> \param rho_set the density where you want to evaluate the functional
1072 : !> \param deriv_set place where to store the functional derivatives (they are
1073 : !> added to the derivatives)
1074 : !> \param grad_deriv degree of the derivative that should be evaluated;
1075 : !> all derivatives up to the given degree are evaluated, in a single
1076 : !> LibXC call per block of grid points
1077 : !> \param libxc_params input parameter (functional name, scaling and parameters)
1078 : !> \param func_name_override optional LibXC functional name overriding the section name
1079 : !> \author F. Tran
1080 : ! **************************************************************************************************
1081 18290 : SUBROUTINE libxc_spin_unpolarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
1082 :
1083 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1084 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
1085 : INTEGER, INTENT(in) :: grad_deriv
1086 : TYPE(section_vals_type), POINTER :: libxc_params
1087 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
1088 :
1089 : #if defined (__LIBXC)
1090 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_unpolarized_eval'
1091 :
1092 : CHARACTER(LEN=default_string_length) :: func_name
1093 : INTEGER :: handle, iw, npoints
1094 : INTEGER, DIMENSION(2, 3) :: bo
1095 : LOGICAL :: has_laplace
1096 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
1097 : TYPE(libxc_worker_set_type), POINTER :: workers
1098 : TYPE(section_vals_type), POINTER :: no_params
1099 18290 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rho, &
1100 18290 : e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_ndrho, &
1101 18290 : e_ndrho_laplace_rho, e_ndrho_ndrho, e_ndrho_rho, e_ndrho_tau, e_rho, &
1102 18290 : e_rho_laplace_rho, e_rho_rho, e_rho_rho_rho, e_rho_tau, e_tau, &
1103 18290 : e_tau_tau, laplace_rho, norm_drho, rho, tau
1104 : ! derivatives w.r.t. the reduced gradient gamma = |grad rho|^2
1105 : #:for descs, arr, idx in umgga_all
1106 18290 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_${'_'.join(descs)}$
1107 : #:endfor
1108 : TYPE(xc_derivative_type), POINTER :: deriv
1109 : TYPE(xc_f03_func_info_t) :: xc_info
1110 :
1111 18290 : CALL timeset(routineN, handle)
1112 :
1113 : ! Only "everything up to grad_deriv" is supported. Asking for a single
1114 : ! derivative order in isolation is rejected here rather than further down,
1115 : ! where it would quietly evaluate nothing at all.
1116 18290 : IF (grad_deriv < 0) THEN
1117 0 : CPABORT("Evaluating a single derivative order is not supported.")
1118 : END IF
1119 18290 : IF (grad_deriv > 3) THEN
1120 0 : CPABORT("derivatives larger than 3 not implemented")
1121 : END IF
1122 :
1123 18290 : has_laplace = .FALSE.
1124 18290 : NULLIFY (dummy)
1125 18290 : NULLIFY (workers)
1126 18290 : NULLIFY (rho, norm_drho, laplace_rho, tau)
1127 :
1128 18290 : IF (PRESENT(func_name_override)) THEN
1129 0 : func_name = func_name_override
1130 0 : func_scale = 1.0_dp
1131 : ELSE
1132 18290 : func_name = libxc_params%section%name
1133 18290 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
1134 : END IF
1135 :
1136 18290 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
1137 :
1138 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
1139 : rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, &
1140 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
1141 18290 : tau=tau, local_bounds=bo)
1142 :
1143 18290 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1144 :
1145 : ! One functional object per thread, kept across calls. The cutoffs are part
1146 : ! of the worker setup because LibXC does the screening itself.
1147 18290 : IF (PRESENT(func_name_override)) THEN
1148 : ! an overriding name carries no input section, so no parameters are set
1149 0 : NULLIFY (no_params)
1150 : iw = libxc_get_workers(func_name, XC_UNPOLARIZED, no_params, epsilon_rho, epsilon_tau, &
1151 0 : grad_deriv)
1152 : ELSE
1153 : iw = libxc_get_workers(func_name, XC_UNPOLARIZED, libxc_params, epsilon_rho, epsilon_tau, &
1154 18290 : grad_deriv)
1155 : END IF
1156 18290 : workers => libxc_cache(iw)
1157 18290 : xc_info = workers%worker(1)%info
1158 18290 : has_laplace = workers%has_laplace
1159 :
1160 : ! see libxc_unused: the arguments the requested order does not produce are
1161 : ! never touched, but they still have to point somewhere
1162 18290 : CALL libxc_reserve_unused(npoints)
1163 18290 : dummy(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) => libxc_unused(1:npoints)
1164 :
1165 : ! due to assumed shape array usage in next routine
1166 18290 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => rho
1167 18290 : IF (.NOT. ASSOCIATED(tau)) tau => rho
1168 :
1169 : ! only some MGGA functionals really need the Laplacian,
1170 : ! all others can work with rho (read-only) as dummy
1171 18290 : IF (.NOT. has_laplace) laplace_rho => rho
1172 :
1173 18290 : e_0 => dummy
1174 18290 : e_rho => dummy
1175 18290 : e_ndrho => dummy
1176 18290 : e_laplace_rho => dummy
1177 18290 : e_tau => dummy
1178 18290 : e_rho_rho => dummy
1179 18290 : e_ndrho_rho => dummy
1180 18290 : e_ndrho_ndrho => dummy
1181 18290 : e_rho_laplace_rho => dummy
1182 18290 : e_rho_tau => dummy
1183 18290 : e_ndrho_laplace_rho => dummy
1184 18290 : e_ndrho_tau => dummy
1185 18290 : e_laplace_rho_laplace_rho => dummy
1186 18290 : e_laplace_rho_tau => dummy
1187 18290 : e_tau_tau => dummy
1188 18290 : e_rho_rho_rho => dummy
1189 : #:for descs, arr, idx in umgga_all
1190 18290 : e_${'_'.join(descs)}$ => dummy
1191 : #:endfor
1192 :
1193 18290 : IF (grad_deriv >= 0) THEN
1194 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
1195 18290 : allocate_deriv=.TRUE.)
1196 18290 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1197 : END IF
1198 18290 : IF (grad_deriv >= 1) THEN
1199 10402 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1200 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1201 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1202 10402 : allocate_deriv=.TRUE.)
1203 10402 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1204 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1205 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1206 5314 : allocate_deriv=.TRUE.)
1207 5314 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1208 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1209 5314 : allocate_deriv=.TRUE.)
1210 5314 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1211 : #:for descs, arr, idx in ugga_fill_1
1212 : deriv => xc_dset_get_derivative(deriv_set, &
1213 : [${', '.join('deriv_' + d for d in descs)}$], &
1214 5314 : allocate_deriv=.TRUE.)
1215 5314 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1216 : #:endfor
1217 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1218 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1219 2320 : allocate_deriv=.TRUE.)
1220 2320 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1221 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1222 2320 : allocate_deriv=.TRUE.)
1223 2320 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1224 : #:for descs, arr, idx in ugga_fill_1
1225 : deriv => xc_dset_get_derivative(deriv_set, &
1226 : [${', '.join('deriv_' + d for d in descs)}$], &
1227 2320 : allocate_deriv=.TRUE.)
1228 2320 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1229 : #:endfor
1230 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
1231 2320 : allocate_deriv=.TRUE.)
1232 2320 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
1233 2320 : IF (has_laplace) THEN
1234 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
1235 616 : allocate_deriv=.TRUE.)
1236 616 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
1237 : END IF
1238 : CASE default
1239 18036 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1240 : END SELECT
1241 : END IF
1242 18290 : IF (grad_deriv >= 2) THEN
1243 1528 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1244 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1245 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1246 1528 : allocate_deriv=.TRUE.)
1247 1528 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1248 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1249 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1250 646 : allocate_deriv=.TRUE.)
1251 646 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1252 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
1253 646 : allocate_deriv=.TRUE.)
1254 646 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
1255 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1256 646 : allocate_deriv=.TRUE.)
1257 646 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1258 : #:for descs, arr, idx in ugga_fill_2
1259 : deriv => xc_dset_get_derivative(deriv_set, &
1260 : [${', '.join('deriv_' + d for d in descs)}$], &
1261 1292 : allocate_deriv=.TRUE.)
1262 1292 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1263 : #:endfor
1264 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1265 : ! not implemented ...
1266 :
1267 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1268 188 : allocate_deriv=.TRUE.)
1269 188 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1270 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
1271 188 : allocate_deriv=.TRUE.)
1272 188 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
1273 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1274 188 : allocate_deriv=.TRUE.)
1275 188 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1276 : #:for descs, arr, idx in ugga_fill_2
1277 : deriv => xc_dset_get_derivative(deriv_set, &
1278 : [${', '.join('deriv_' + d for d in descs)}$], &
1279 376 : allocate_deriv=.TRUE.)
1280 376 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1281 : #:endfor
1282 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_tau], &
1283 188 : allocate_deriv=.TRUE.)
1284 188 : CALL xc_derivative_get(deriv, deriv_data=e_rho_tau)
1285 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau], &
1286 188 : allocate_deriv=.TRUE.)
1287 188 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau)
1288 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau, deriv_tau], &
1289 188 : allocate_deriv=.TRUE.)
1290 188 : CALL xc_derivative_get(deriv, deriv_data=e_tau_tau)
1291 : #:for descs, arr, idx in umgga_extra_2_nl
1292 : deriv => xc_dset_get_derivative(deriv_set, &
1293 : [${', '.join('deriv_' + d for d in descs)}$], &
1294 188 : allocate_deriv=.TRUE.)
1295 188 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1296 : #:endfor
1297 188 : IF (has_laplace) THEN
1298 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_laplace_rho], &
1299 82 : allocate_deriv=.TRUE.)
1300 82 : CALL xc_derivative_get(deriv, deriv_data=e_rho_laplace_rho)
1301 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rho], &
1302 82 : allocate_deriv=.TRUE.)
1303 82 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rho)
1304 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_laplace_rho], &
1305 82 : allocate_deriv=.TRUE.)
1306 82 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_laplace_rho)
1307 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_tau], &
1308 82 : allocate_deriv=.TRUE.)
1309 82 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_tau)
1310 : END IF
1311 : CASE default
1312 2362 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1313 : END SELECT
1314 : END IF
1315 18290 : IF (grad_deriv >= 3) THEN
1316 8 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1317 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1318 8 : IF (has_laplace) THEN
1319 : #:for descs, arr, idx in umgga_fill_3_l
1320 : deriv => xc_dset_get_derivative(deriv_set, &
1321 : [${', '.join('deriv_' + d for d in descs)}$], &
1322 0 : allocate_deriv=.TRUE.)
1323 0 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1324 : #:endfor
1325 : END IF
1326 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
1327 8 : allocate_deriv=.TRUE.)
1328 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
1329 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1330 8 : IF (has_laplace) THEN
1331 : #:for descs, arr, idx in umgga_fill_3_l
1332 : deriv => xc_dset_get_derivative(deriv_set, &
1333 : [${', '.join('deriv_' + d for d in descs)}$], &
1334 0 : allocate_deriv=.TRUE.)
1335 0 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1336 : #:endfor
1337 : END IF
1338 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
1339 8 : allocate_deriv=.TRUE.)
1340 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
1341 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_gamma], &
1342 8 : allocate_deriv=.TRUE.)
1343 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_gamma)
1344 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_gamma, deriv_gamma], &
1345 8 : allocate_deriv=.TRUE.)
1346 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho_gamma_gamma)
1347 : deriv => xc_dset_get_derivative(deriv_set, [deriv_gamma, deriv_gamma, deriv_gamma], &
1348 8 : allocate_deriv=.TRUE.)
1349 8 : CALL xc_derivative_get(deriv, deriv_data=e_gamma_gamma_gamma)
1350 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1351 : #:for descs, arr, idx in umgga_fill_3_nl
1352 : deriv => xc_dset_get_derivative(deriv_set, &
1353 : [${', '.join('deriv_' + d for d in descs)}$], &
1354 72 : allocate_deriv=.TRUE.)
1355 72 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1356 : #:endfor
1357 8 : IF (has_laplace) THEN
1358 : #:for descs, arr, idx in umgga_fill_3_l
1359 : deriv => xc_dset_get_derivative(deriv_set, &
1360 : [${', '.join('deriv_' + d for d in descs)}$], &
1361 40 : allocate_deriv=.TRUE.)
1362 40 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1363 : #:endfor
1364 : END IF
1365 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
1366 8 : allocate_deriv=.TRUE.)
1367 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
1368 : CASE default
1369 24 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1370 : END SELECT
1371 : END IF
1372 :
1373 : !$OMP PARALLEL DEFAULT(NONE), &
1374 : !$OMP SHARED(rho,norm_drho,laplace_rho,tau,e_0,e_rho,e_ndrho,e_laplace_rho),&
1375 : !$OMP SHARED(e_tau,e_rho_rho,e_ndrho_rho,e_ndrho_ndrho,e_rho_laplace_rho),&
1376 : !$OMP SHARED(e_rho_tau,e_ndrho_laplace_rho,e_ndrho_tau,e_laplace_rho_laplace_rho),&
1377 : !$OMP SHARED(e_laplace_rho_tau,e_tau_tau,e_rho_rho_rho),&
1378 : #:for descs, arr, idx in umgga_all
1379 : !$OMP SHARED(e_${'_'.join(descs)}$),&
1380 : #:endfor
1381 : !$OMP SHARED(grad_deriv,npoints),&
1382 18290 : !$OMP SHARED(func_name,func_scale,workers)
1383 :
1384 : CALL libxc_spin_unpolarized_calc(rho=rho, norm_drho=norm_drho, &
1385 : laplace_rho=laplace_rho, tau=tau, &
1386 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_laplace_rho=e_laplace_rho, &
1387 : e_tau=e_tau, e_rho_rho=e_rho_rho, e_ndrho_rho=e_ndrho_rho, &
1388 : e_ndrho_ndrho=e_ndrho_ndrho, e_rho_laplace_rho=e_rho_laplace_rho, &
1389 : e_rho_tau=e_rho_tau, e_ndrho_laplace_rho=e_ndrho_laplace_rho, &
1390 : e_ndrho_tau=e_ndrho_tau, e_laplace_rho_laplace_rho=e_laplace_rho_laplace_rho, &
1391 : e_laplace_rho_tau=e_laplace_rho_tau, e_tau_tau=e_tau_tau, &
1392 : e_rho_rho_rho=e_rho_rho_rho, &
1393 : #:for descs, arr, idx in umgga_all
1394 : e_${'_'.join(descs)}$=e_${'_'.join(descs)}$, &
1395 : #:endfor
1396 : grad_deriv=grad_deriv, npoints=npoints, &
1397 : func_name=func_name, sc=func_scale, workers=workers)
1398 :
1399 : !$OMP END PARALLEL
1400 :
1401 18290 : NULLIFY (dummy)
1402 18290 : NULLIFY (workers)
1403 :
1404 18290 : CALL timestop(handle)
1405 : #else
1406 : MARK_USED(rho_set)
1407 : MARK_USED(deriv_set)
1408 : MARK_USED(grad_deriv)
1409 : MARK_USED(libxc_params)
1410 : MARK_USED(func_name_override)
1411 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
1412 : "for a functional of the LibXC library, "// &
1413 : "you have to download and install the library!")
1414 : #endif
1415 18290 : END SUBROUTINE libxc_spin_unpolarized_eval
1416 :
1417 : ! **************************************************************************************************
1418 : !> \brief evaluates the functional from libxc
1419 : !> \param rho_set the density where you want to evaluate the functional
1420 : !> \param deriv_set place where to store the functional derivatives (they are
1421 : !> added to the derivatives)
1422 : !> \param grad_deriv degree of the derivative that should be evaluated;
1423 : !> all derivatives up to the given degree are evaluated, in a single
1424 : !> LibXC call per block of grid points
1425 : !> \param libxc_params input parameter (functional name, scaling and parameters)
1426 : !> \param func_name_override optional LibXC functional name overriding the section name
1427 : !> \author F. Tran
1428 : ! **************************************************************************************************
1429 3552 : SUBROUTINE libxc_spin_polarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
1430 :
1431 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1432 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
1433 : INTEGER, INTENT(in) :: grad_deriv
1434 : TYPE(section_vals_type), POINTER :: libxc_params
1435 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
1436 :
1437 : #if defined (__LIBXC)
1438 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_polarized_eval'
1439 :
1440 : CHARACTER(LEN=default_string_length) :: func_name
1441 : INTEGER :: handle, iw, npoints
1442 : INTEGER, DIMENSION(2, 3) :: bo
1443 : LOGICAL :: has_laplace
1444 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
1445 : TYPE(libxc_worker_set_type), POINTER :: workers
1446 : TYPE(section_vals_type), POINTER :: no_params
1447 3552 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
1448 3552 : e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
1449 3552 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, e_laplace_rhob, &
1450 3552 : e_laplace_rhob_laplace_rhob, e_laplace_rhob_tau_a, &
1451 3552 : e_laplace_rhob_tau_b, e_ndrho, e_ndrho_laplace_rhoa, &
1452 3552 : e_ndrho_laplace_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
1453 3552 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa, &
1454 3552 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, e_ndrhoa_ndrhoa, &
1455 3552 : e_ndrhoa_ndrhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhoa_tau_a, &
1456 3552 : e_ndrhoa_tau_b, e_ndrhob
1457 3552 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_ndrhob_laplace_rhoa, &
1458 3552 : e_ndrhob_laplace_rhob, e_ndrhob_ndrhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
1459 3552 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_rhoa, e_rhoa_laplace_rhoa, &
1460 3552 : e_rhoa_laplace_rhob, e_rhoa_rhoa, e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
1461 3552 : e_rhoa_rhob, e_rhoa_rhob_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob, &
1462 3552 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhob_rhob, &
1463 3552 : e_rhob_rhob_rhob, e_rhob_tau_a, e_rhob_tau_b, e_tau_a, e_tau_a_tau_a, &
1464 3552 : e_tau_a_tau_b, e_tau_b, e_tau_b_tau_b, laplace_rhoa, laplace_rhob, &
1465 3552 : norm_drho, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
1466 : ! derivatives w.r.t. the reduced gradients gamma_ij
1467 : #:for descs, arr, idx in mgga_all
1468 3552 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_${'_'.join(descs)}$
1469 : #:endfor
1470 : TYPE(xc_derivative_type), POINTER :: deriv
1471 : TYPE(xc_f03_func_info_t) :: xc_info
1472 :
1473 3552 : CALL timeset(routineN, handle)
1474 :
1475 : ! Only "everything up to grad_deriv" is supported. Asking for a single
1476 : ! derivative order in isolation is rejected here rather than further down,
1477 : ! where it would quietly evaluate nothing at all.
1478 3552 : IF (grad_deriv < 0) THEN
1479 0 : CPABORT("Evaluating a single derivative order is not supported.")
1480 : END IF
1481 3552 : IF (grad_deriv > 3) THEN
1482 0 : CPABORT("derivatives larger than 3 not implemented")
1483 : END IF
1484 :
1485 3552 : NULLIFY (dummy)
1486 3552 : NULLIFY (workers)
1487 3552 : NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, &
1488 3552 : laplace_rhob, tau_a, tau_b)
1489 :
1490 3552 : IF (PRESENT(func_name_override)) THEN
1491 0 : func_name = func_name_override
1492 0 : func_scale = 1.0_dp
1493 : ELSE
1494 3552 : func_name = libxc_params%section%name
1495 3552 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
1496 : END IF
1497 :
1498 3552 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
1499 :
1500 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
1501 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1502 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
1503 : laplace_rhoa=laplace_rhoa, laplace_rhob=laplace_rhob, &
1504 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
1505 3552 : tau_a=tau_a, tau_b=tau_b, local_bounds=bo)
1506 :
1507 3552 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1508 :
1509 : ! One functional object per thread, kept across calls. The cutoffs are part
1510 : ! of the worker setup because LibXC does the screening itself.
1511 3552 : IF (PRESENT(func_name_override)) THEN
1512 : ! an overriding name carries no input section, so no parameters are set
1513 0 : NULLIFY (no_params)
1514 : iw = libxc_get_workers(func_name, XC_POLARIZED, no_params, epsilon_rho, epsilon_tau, &
1515 0 : grad_deriv)
1516 : ELSE
1517 : iw = libxc_get_workers(func_name, XC_POLARIZED, libxc_params, epsilon_rho, epsilon_tau, &
1518 3552 : grad_deriv)
1519 : END IF
1520 3552 : workers => libxc_cache(iw)
1521 3552 : xc_info = workers%worker(1)%info
1522 3552 : has_laplace = workers%has_laplace
1523 :
1524 : ! see libxc_unused: the arguments the requested order does not produce are
1525 : ! never touched, but they still have to point somewhere
1526 3552 : CALL libxc_reserve_unused(npoints)
1527 3552 : dummy(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) => libxc_unused(1:npoints)
1528 :
1529 : ! due to assumed shape array usage in next routine
1530 3552 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => rhoa
1531 3552 : IF (.NOT. ASSOCIATED(norm_drhoa)) norm_drhoa => rhoa
1532 3552 : IF (.NOT. ASSOCIATED(norm_drhob)) norm_drhob => rhoa
1533 3552 : IF (.NOT. ASSOCIATED(tau_a)) tau_a => rhoa
1534 3552 : IF (.NOT. ASSOCIATED(tau_b)) tau_b => rhoa
1535 :
1536 : ! only some MGGA functionals really need the Laplacian,
1537 : ! all others can work with rhoa (read-only) as dummy
1538 3552 : IF (.NOT. has_laplace) laplace_rhoa => rhoa
1539 3250 : IF (.NOT. has_laplace) laplace_rhob => rhoa
1540 :
1541 3552 : e_0 => dummy
1542 3552 : e_rhoa => dummy
1543 3552 : e_rhob => dummy
1544 3552 : e_ndrho => dummy
1545 3552 : e_ndrhoa => dummy
1546 3552 : e_ndrhob => dummy
1547 3552 : e_laplace_rhoa => dummy
1548 3552 : e_laplace_rhob => dummy
1549 3552 : e_tau_a => dummy
1550 3552 : e_tau_b => dummy
1551 3552 : e_rhoa_rhoa => dummy
1552 3552 : e_rhoa_rhob => dummy
1553 3552 : e_rhob_rhob => dummy
1554 3552 : e_ndrho_rhoa => dummy
1555 3552 : e_ndrho_rhob => dummy
1556 3552 : e_ndrhoa_rhoa => dummy
1557 3552 : e_ndrhoa_rhob => dummy
1558 3552 : e_ndrhob_rhoa => dummy
1559 3552 : e_ndrhob_rhob => dummy
1560 3552 : e_ndrho_ndrho => dummy
1561 3552 : e_ndrho_ndrhoa => dummy
1562 3552 : e_ndrho_ndrhob => dummy
1563 3552 : e_ndrhoa_ndrhoa => dummy
1564 3552 : e_ndrhoa_ndrhob => dummy
1565 3552 : e_ndrhob_ndrhob => dummy
1566 3552 : e_rhoa_laplace_rhoa => dummy
1567 3552 : e_rhoa_laplace_rhob => dummy
1568 3552 : e_rhob_laplace_rhoa => dummy
1569 3552 : e_rhob_laplace_rhob => dummy
1570 3552 : e_rhoa_tau_a => dummy
1571 3552 : e_rhoa_tau_b => dummy
1572 3552 : e_rhob_tau_a => dummy
1573 3552 : e_rhob_tau_b => dummy
1574 3552 : e_ndrho_laplace_rhoa => dummy
1575 3552 : e_ndrho_laplace_rhob => dummy
1576 3552 : e_ndrhoa_laplace_rhoa => dummy
1577 3552 : e_ndrhoa_laplace_rhob => dummy
1578 3552 : e_ndrhob_laplace_rhoa => dummy
1579 3552 : e_ndrhob_laplace_rhob => dummy
1580 3552 : e_ndrho_tau_a => dummy
1581 3552 : e_ndrho_tau_b => dummy
1582 3552 : e_ndrhoa_tau_a => dummy
1583 3552 : e_ndrhoa_tau_b => dummy
1584 3552 : e_ndrhob_tau_a => dummy
1585 3552 : e_ndrhob_tau_b => dummy
1586 3552 : e_laplace_rhoa_laplace_rhoa => dummy
1587 3552 : e_laplace_rhoa_laplace_rhob => dummy
1588 3552 : e_laplace_rhob_laplace_rhob => dummy
1589 3552 : e_laplace_rhoa_tau_a => dummy
1590 3552 : e_laplace_rhoa_tau_b => dummy
1591 3552 : e_laplace_rhob_tau_a => dummy
1592 3552 : e_laplace_rhob_tau_b => dummy
1593 3552 : e_tau_a_tau_a => dummy
1594 3552 : e_tau_a_tau_b => dummy
1595 3552 : e_tau_b_tau_b => dummy
1596 3552 : e_rhoa_rhoa_rhoa => dummy
1597 3552 : e_rhoa_rhoa_rhob => dummy
1598 3552 : e_rhoa_rhob_rhob => dummy
1599 3552 : e_rhob_rhob_rhob => dummy
1600 : #:for descs, arr, idx in mgga_all
1601 3552 : e_${'_'.join(descs)}$ => dummy
1602 : #:endfor
1603 :
1604 3552 : IF (grad_deriv >= 0) THEN
1605 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
1606 3552 : allocate_deriv=.TRUE.)
1607 3552 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1608 : END IF
1609 3552 : IF (grad_deriv >= 1) THEN
1610 1600 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1611 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1612 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1613 1600 : allocate_deriv=.TRUE.)
1614 1600 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1615 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1616 1600 : allocate_deriv=.TRUE.)
1617 1600 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1618 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1619 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1620 618 : allocate_deriv=.TRUE.)
1621 618 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1622 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1623 618 : allocate_deriv=.TRUE.)
1624 618 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1625 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1626 618 : allocate_deriv=.TRUE.)
1627 618 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1628 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1629 618 : allocate_deriv=.TRUE.)
1630 618 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
1631 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1632 618 : allocate_deriv=.TRUE.)
1633 618 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
1634 : #:for descs, arr, idx in gamma_only_1
1635 : deriv => xc_dset_get_derivative(deriv_set, &
1636 : [${', '.join('deriv_' + d for d in descs)}$], &
1637 1854 : allocate_deriv=.TRUE.)
1638 1854 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1639 : #:endfor
1640 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1641 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1642 1288 : allocate_deriv=.TRUE.)
1643 1288 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1644 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1645 1288 : allocate_deriv=.TRUE.)
1646 1288 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1647 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1648 1288 : allocate_deriv=.TRUE.)
1649 1288 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1650 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1651 1288 : allocate_deriv=.TRUE.)
1652 1288 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
1653 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1654 1288 : allocate_deriv=.TRUE.)
1655 1288 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
1656 : #:for descs, arr, idx in gamma_only_1
1657 : deriv => xc_dset_get_derivative(deriv_set, &
1658 : [${', '.join('deriv_' + d for d in descs)}$], &
1659 3864 : allocate_deriv=.TRUE.)
1660 3864 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1661 : #:endfor
1662 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
1663 1288 : allocate_deriv=.TRUE.)
1664 1288 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
1665 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
1666 1288 : allocate_deriv=.TRUE.)
1667 1288 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
1668 1288 : IF (has_laplace) THEN
1669 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
1670 294 : allocate_deriv=.TRUE.)
1671 294 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
1672 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
1673 294 : allocate_deriv=.TRUE.)
1674 294 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
1675 : END IF
1676 : CASE default
1677 3506 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1678 : END SELECT
1679 : END IF
1680 3552 : IF (grad_deriv >= 2) THEN
1681 78 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1682 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1683 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1684 78 : allocate_deriv=.TRUE.)
1685 78 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1686 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1687 78 : allocate_deriv=.TRUE.)
1688 78 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1689 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1690 78 : allocate_deriv=.TRUE.)
1691 78 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1692 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1693 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1694 74 : allocate_deriv=.TRUE.)
1695 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1696 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1697 74 : allocate_deriv=.TRUE.)
1698 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1699 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1700 74 : allocate_deriv=.TRUE.)
1701 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1702 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1703 74 : allocate_deriv=.TRUE.)
1704 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1705 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1706 74 : allocate_deriv=.TRUE.)
1707 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1708 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1709 74 : allocate_deriv=.TRUE.)
1710 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1711 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1712 74 : allocate_deriv=.TRUE.)
1713 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1714 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1715 74 : allocate_deriv=.TRUE.)
1716 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1717 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1718 74 : allocate_deriv=.TRUE.)
1719 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1720 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1721 74 : allocate_deriv=.TRUE.)
1722 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1723 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1724 74 : allocate_deriv=.TRUE.)
1725 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1726 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1727 74 : allocate_deriv=.TRUE.)
1728 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1729 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1730 74 : allocate_deriv=.TRUE.)
1731 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1732 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1733 74 : allocate_deriv=.TRUE.)
1734 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1735 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1736 74 : allocate_deriv=.TRUE.)
1737 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1738 : #:for descs, arr, idx in gamma_only_2
1739 : deriv => xc_dset_get_derivative(deriv_set, &
1740 : [${', '.join('deriv_' + d for d in descs)}$], &
1741 888 : allocate_deriv=.TRUE.)
1742 888 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1743 : #:endfor
1744 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1745 :
1746 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1747 74 : allocate_deriv=.TRUE.)
1748 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1749 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1750 74 : allocate_deriv=.TRUE.)
1751 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1752 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1753 74 : allocate_deriv=.TRUE.)
1754 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1755 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1756 74 : allocate_deriv=.TRUE.)
1757 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1758 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1759 74 : allocate_deriv=.TRUE.)
1760 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1761 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1762 74 : allocate_deriv=.TRUE.)
1763 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1764 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1765 74 : allocate_deriv=.TRUE.)
1766 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1767 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1768 74 : allocate_deriv=.TRUE.)
1769 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1770 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1771 74 : allocate_deriv=.TRUE.)
1772 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1773 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1774 74 : allocate_deriv=.TRUE.)
1775 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1776 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1777 74 : allocate_deriv=.TRUE.)
1778 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1779 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1780 74 : allocate_deriv=.TRUE.)
1781 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1782 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1783 74 : allocate_deriv=.TRUE.)
1784 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1785 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1786 74 : allocate_deriv=.TRUE.)
1787 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1788 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1789 74 : allocate_deriv=.TRUE.)
1790 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1791 : #:for descs, arr, idx in mgga_fill_2_nl
1792 : deriv => xc_dset_get_derivative(deriv_set, &
1793 : [${', '.join('deriv_' + d for d in descs)}$], &
1794 1332 : allocate_deriv=.TRUE.)
1795 1332 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1796 : #:endfor
1797 74 : IF (has_laplace) THEN
1798 : #:for descs, arr, idx in mgga_fill_2_l
1799 : deriv => xc_dset_get_derivative(deriv_set, &
1800 : [${', '.join('deriv_' + d for d in descs)}$], &
1801 156 : allocate_deriv=.TRUE.)
1802 156 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1803 : #:endfor
1804 : END IF
1805 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_a], &
1806 74 : allocate_deriv=.TRUE.)
1807 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_a)
1808 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_b], &
1809 74 : allocate_deriv=.TRUE.)
1810 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_b)
1811 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_a], &
1812 74 : allocate_deriv=.TRUE.)
1813 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_a)
1814 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_b], &
1815 74 : allocate_deriv=.TRUE.)
1816 74 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_b)
1817 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_a], &
1818 74 : allocate_deriv=.TRUE.)
1819 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_a)
1820 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_b], &
1821 74 : allocate_deriv=.TRUE.)
1822 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_b)
1823 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_a], &
1824 74 : allocate_deriv=.TRUE.)
1825 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_a)
1826 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_b], &
1827 74 : allocate_deriv=.TRUE.)
1828 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_b)
1829 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_a], &
1830 74 : allocate_deriv=.TRUE.)
1831 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_a)
1832 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_b], &
1833 74 : allocate_deriv=.TRUE.)
1834 74 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_b)
1835 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_a], &
1836 74 : allocate_deriv=.TRUE.)
1837 74 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_a)
1838 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_b], &
1839 74 : allocate_deriv=.TRUE.)
1840 74 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_b)
1841 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b, deriv_tau_b], &
1842 74 : allocate_deriv=.TRUE.)
1843 74 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b_tau_b)
1844 74 : IF (has_laplace) THEN
1845 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhoa], &
1846 26 : allocate_deriv=.TRUE.)
1847 26 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhoa)
1848 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhob], &
1849 26 : allocate_deriv=.TRUE.)
1850 26 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhob)
1851 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhoa], &
1852 26 : allocate_deriv=.TRUE.)
1853 26 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhoa)
1854 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhob], &
1855 26 : allocate_deriv=.TRUE.)
1856 26 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhob)
1857 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhoa], &
1858 26 : allocate_deriv=.TRUE.)
1859 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhoa)
1860 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhob], &
1861 26 : allocate_deriv=.TRUE.)
1862 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhob)
1863 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhoa], &
1864 26 : allocate_deriv=.TRUE.)
1865 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhoa)
1866 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhob], &
1867 26 : allocate_deriv=.TRUE.)
1868 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhob)
1869 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhoa], &
1870 26 : allocate_deriv=.TRUE.)
1871 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhoa)
1872 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhob], &
1873 26 : allocate_deriv=.TRUE.)
1874 26 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhob)
1875 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhoa], &
1876 26 : allocate_deriv=.TRUE.)
1877 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhoa)
1878 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhob], &
1879 26 : allocate_deriv=.TRUE.)
1880 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhob)
1881 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_laplace_rhob], &
1882 26 : allocate_deriv=.TRUE.)
1883 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_laplace_rhob)
1884 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_a], &
1885 26 : allocate_deriv=.TRUE.)
1886 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_a)
1887 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_b], &
1888 26 : allocate_deriv=.TRUE.)
1889 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_b)
1890 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_a], &
1891 26 : allocate_deriv=.TRUE.)
1892 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_a)
1893 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_b], &
1894 26 : allocate_deriv=.TRUE.)
1895 26 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_b)
1896 : END IF
1897 : CASE default
1898 226 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1899 : END SELECT
1900 : END IF
1901 3552 : IF (grad_deriv >= 3) THEN
1902 8 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1903 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1904 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
1905 8 : allocate_deriv=.TRUE.)
1906 8 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhoa)
1907 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
1908 8 : allocate_deriv=.TRUE.)
1909 8 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhob)
1910 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
1911 8 : allocate_deriv=.TRUE.)
1912 8 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob_rhob)
1913 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
1914 8 : allocate_deriv=.TRUE.)
1915 8 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob_rhob)
1916 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1917 : #:for descs, arr, idx in gamma_only_3 + rho3_entries
1918 : deriv => xc_dset_get_derivative(deriv_set, &
1919 : [${', '.join('deriv_' + d for d in descs)}$], &
1920 280 : allocate_deriv=.TRUE.)
1921 280 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1922 : #:endfor
1923 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1924 : #:for descs, arr, idx in mgga_fill_3_nl
1925 : deriv => xc_dset_get_derivative(deriv_set, &
1926 : [${', '.join('deriv_' + d for d in descs)}$], &
1927 960 : allocate_deriv=.TRUE.)
1928 960 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1929 : #:endfor
1930 12 : IF (has_laplace) THEN
1931 : #:for descs, arr, idx in mgga_fill_3_l
1932 : deriv => xc_dset_get_derivative(deriv_set, &
1933 : [${', '.join('deriv_' + d for d in descs)}$], &
1934 324 : allocate_deriv=.TRUE.)
1935 324 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1936 : #:endfor
1937 : END IF
1938 : #:for descs, arr, idx in rho3_entries
1939 : deriv => xc_dset_get_derivative(deriv_set, &
1940 : [${', '.join('deriv_' + d for d in descs)}$], &
1941 48 : allocate_deriv=.TRUE.)
1942 48 : CALL xc_derivative_get(deriv, deriv_data=e_${'_'.join(descs)}$)
1943 : #:endfor
1944 : CASE default
1945 28 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1946 : END SELECT
1947 : END IF
1948 :
1949 : !$OMP PARALLEL DEFAULT(NONE), &
1950 : #:for descs, arr, idx in mgga_all
1951 : !$OMP SHARED(e_${'_'.join(descs)}$),&
1952 : #:endfor
1953 : !$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob),&
1954 : !$OMP SHARED(laplace_rhoa,laplace_rhob,tau_a,tau_b),&
1955 : !$OMP SHARED(e_0,e_rhoa,e_rhob,e_ndrho,e_ndrhoa,e_ndrhob),&
1956 : !$OMP SHARED(e_laplace_rhoa,e_laplace_rhob,e_tau_a,e_tau_b),&
1957 : !$OMP SHARED(e_rhoa_rhoa,e_rhoa_rhob,e_rhob_rhob),&
1958 : !$OMP SHARED(e_ndrho_rhoa,e_ndrho_rhob),&
1959 : !$OMP SHARED(e_ndrhoa_rhoa,e_ndrhoa_rhob,e_ndrhob_rhoa,e_ndrhob_rhob),&
1960 : !$OMP SHARED(e_ndrho_ndrho,e_ndrho_ndrhoa,e_ndrho_ndrhob),&
1961 : !$OMP SHARED(e_ndrhoa_ndrhoa,e_ndrhoa_ndrhob,e_ndrhob_ndrhob),&
1962 : !$OMP SHARED(e_rhoa_laplace_rhoa,e_rhoa_laplace_rhob,e_rhob_laplace_rhoa,e_rhob_laplace_rhob),&
1963 : !$OMP SHARED(e_rhoa_tau_a,e_rhoa_tau_b,e_rhob_tau_a,e_rhob_tau_b),&
1964 : !$OMP SHARED(e_ndrho_laplace_rhoa,e_ndrho_laplace_rhob),&
1965 : !$OMP SHARED(e_ndrhoa_laplace_rhoa,e_ndrhoa_laplace_rhob,e_ndrhob_laplace_rhoa,e_ndrhob_laplace_rhob),&
1966 : !$OMP SHARED(e_ndrho_tau_a,e_ndrho_tau_b),&
1967 : !$OMP SHARED(e_ndrhoa_tau_a,e_ndrhoa_tau_b,e_ndrhob_tau_a,e_ndrhob_tau_b),&
1968 : !$OMP SHARED(e_laplace_rhoa_laplace_rhoa,e_laplace_rhoa_laplace_rhob,e_laplace_rhob_laplace_rhob),&
1969 : !$OMP SHARED(e_laplace_rhoa_tau_a,e_laplace_rhoa_tau_b,e_laplace_rhob_tau_a,e_laplace_rhob_tau_b),&
1970 : !$OMP SHARED(e_tau_a_tau_a,e_tau_a_tau_b,e_tau_b_tau_b),&
1971 : !$OMP SHARED(e_rhoa_rhoa_rhoa,e_rhoa_rhoa_rhob,e_rhoa_rhob_rhob,e_rhob_rhob_rhob),&
1972 : !$OMP SHARED(grad_deriv,npoints),&
1973 3552 : !$OMP SHARED(func_name,func_scale,workers)
1974 :
1975 : CALL libxc_spin_polarized_calc(rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1976 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, laplace_rhoa=laplace_rhoa, &
1977 : laplace_rhob=laplace_rhob, tau_a=tau_a, tau_b=tau_b, &
1978 : e_0=e_0, e_rhoa=e_rhoa, e_rhob=e_rhob, e_ndrho=e_ndrho, &
1979 : e_ndrhoa=e_ndrhoa, e_ndrhob=e_ndrhob, e_laplace_rhoa=e_laplace_rhoa, &
1980 : e_laplace_rhob=e_laplace_rhob, e_tau_a=e_tau_a, e_tau_b=e_tau_b, &
1981 : e_rhoa_rhoa=e_rhoa_rhoa, e_rhoa_rhob=e_rhoa_rhob, e_rhob_rhob=e_rhob_rhob, &
1982 : e_ndrho_rhoa=e_ndrho_rhoa, e_ndrho_rhob=e_ndrho_rhob, &
1983 : e_ndrhoa_rhoa=e_ndrhoa_rhoa, e_ndrhoa_rhob=e_ndrhoa_rhob, &
1984 : e_ndrhob_rhoa=e_ndrhob_rhoa, e_ndrhob_rhob=e_ndrhob_rhob, &
1985 : e_ndrho_ndrho=e_ndrho_ndrho, e_ndrho_ndrhoa=e_ndrho_ndrhoa, &
1986 : e_ndrho_ndrhob=e_ndrho_ndrhob, e_ndrhoa_ndrhoa=e_ndrhoa_ndrhoa, &
1987 : e_ndrhoa_ndrhob=e_ndrhoa_ndrhob, e_ndrhob_ndrhob=e_ndrhob_ndrhob, &
1988 : e_rhoa_laplace_rhoa=e_rhoa_laplace_rhoa, &
1989 : e_rhoa_laplace_rhob=e_rhoa_laplace_rhob, &
1990 : e_rhob_laplace_rhoa=e_rhob_laplace_rhoa, &
1991 : e_rhob_laplace_rhob=e_rhob_laplace_rhob, &
1992 : e_rhoa_tau_a=e_rhoa_tau_a, e_rhoa_tau_b=e_rhoa_tau_b, &
1993 : e_rhob_tau_a=e_rhob_tau_a, e_rhob_tau_b=e_rhob_tau_b, &
1994 : e_ndrho_laplace_rhoa=e_ndrho_laplace_rhoa, &
1995 : e_ndrho_laplace_rhob=e_ndrho_laplace_rhob, &
1996 : e_ndrhoa_laplace_rhoa=e_ndrhoa_laplace_rhoa, &
1997 : e_ndrhoa_laplace_rhob=e_ndrhoa_laplace_rhob, &
1998 : e_ndrhob_laplace_rhoa=e_ndrhob_laplace_rhoa, &
1999 : e_ndrhob_laplace_rhob=e_ndrhob_laplace_rhob, &
2000 : e_ndrho_tau_a=e_ndrho_tau_a, e_ndrho_tau_b=e_ndrho_tau_b, &
2001 : e_ndrhoa_tau_a=e_ndrhoa_tau_a, e_ndrhoa_tau_b=e_ndrhoa_tau_b, &
2002 : e_ndrhob_tau_a=e_ndrhob_tau_a, e_ndrhob_tau_b=e_ndrhob_tau_b, &
2003 : e_laplace_rhoa_laplace_rhoa=e_laplace_rhoa_laplace_rhoa, &
2004 : e_laplace_rhoa_laplace_rhob=e_laplace_rhoa_laplace_rhob, &
2005 : e_laplace_rhob_laplace_rhob=e_laplace_rhob_laplace_rhob, &
2006 : e_laplace_rhoa_tau_a=e_laplace_rhoa_tau_a, &
2007 : e_laplace_rhoa_tau_b=e_laplace_rhoa_tau_b, &
2008 : e_laplace_rhob_tau_a=e_laplace_rhob_tau_a, &
2009 : e_laplace_rhob_tau_b=e_laplace_rhob_tau_b, &
2010 : e_tau_a_tau_a=e_tau_a_tau_a, &
2011 : e_tau_a_tau_b=e_tau_a_tau_b, &
2012 : e_tau_b_tau_b=e_tau_b_tau_b, &
2013 : e_rhoa_rhoa_rhoa=e_rhoa_rhoa_rhoa, &
2014 : e_rhoa_rhoa_rhob=e_rhoa_rhoa_rhob, &
2015 : e_rhoa_rhob_rhob=e_rhoa_rhob_rhob, &
2016 : e_rhob_rhob_rhob=e_rhob_rhob_rhob, &
2017 : #:for descs, arr, idx in mgga_all
2018 : e_${'_'.join(descs)}$=e_${'_'.join(descs)}$, &
2019 : #:endfor
2020 : grad_deriv=grad_deriv, npoints=npoints, &
2021 : func_name=func_name, sc=func_scale, workers=workers)
2022 :
2023 : !$OMP END PARALLEL
2024 :
2025 3552 : NULLIFY (dummy)
2026 3552 : NULLIFY (workers)
2027 :
2028 3552 : CALL timestop(handle)
2029 : #else
2030 : MARK_USED(rho_set)
2031 : MARK_USED(deriv_set)
2032 : MARK_USED(grad_deriv)
2033 : MARK_USED(libxc_params)
2034 : MARK_USED(func_name_override)
2035 :
2036 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
2037 : "for a functional of the LibXC library, "// &
2038 : "you have to download and install the library!")
2039 : #endif
2040 3552 : END SUBROUTINE libxc_spin_polarized_eval
2041 :
2042 : ! **************************************************************************************************
2043 : !> \brief libxc exchange-correlation functionals
2044 : !> \param rho density
2045 : !> \param norm_drho norm of the gradient of the density
2046 : !> \param laplace_rho laplacian of the density
2047 : !> \param tau kinetic-energy density
2048 : !> \param e_0 energy density
2049 : !> \param e_rho derivative of the energy density with respect to rho
2050 : !> \param e_ndrho derivative of the energy density with respect to ndrho
2051 : !> \param e_laplace_rho derivative of the energy density with respect to laplace_rho
2052 : !> \param e_tau derivative of the energy density with respect to tau
2053 : !> \param e_rho_rho derivative of the energy density with respect to rho_rho
2054 : !> \param e_ndrho_rho derivative of the energy density with respect to ndrho_rho
2055 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
2056 : !> \param e_rho_laplace_rho derivative of the energy density with respect to rho_laplace_rho
2057 : !> \param e_rho_tau derivative of the energy density with respect to rho_tau
2058 : !> \param e_ndrho_laplace_rho derivative of the energy density with respect to ndrho_laplace_rho
2059 : !> \param e_ndrho_tau derivative of the energy density with respect to ndrho_tau
2060 : !> \param e_laplace_rho_laplace_rho derivative of the energy density with respect to laplace_rho_laplace_rho
2061 : !> \param e_laplace_rho_tau derivative of the energy density with respect to laplace_rho_tau
2062 : !> \param e_tau_tau derivative of the energy density with respect to tau_tau
2063 : !> \param e_rho_rho_rho derivative of the energy density with respect to rho_rho_rho
2064 : !> \param grad_deriv degree of the derivative that should be evaluated;
2065 : !> all derivatives up to the given degree are evaluated, in a single
2066 : !> LibXC call per block of grid points
2067 : !> \param npoints number of points on the grid
2068 : !> \param func_name name of the functional
2069 : !> \param sc scaling factor of the functional
2070 : !> \param workers cached LibXC functional objects and staging buffers, one per thread
2071 : !> \author F. Tran
2072 : ! **************************************************************************************************
2073 : #if defined (__LIBXC)
2074 18290 : SUBROUTINE libxc_spin_unpolarized_calc(rho, norm_drho, laplace_rho, tau, &
2075 : e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, e_rho_rho, e_ndrho_rho, &
2076 : e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
2077 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, &
2078 : e_tau_tau, e_rho_rho_rho, &
2079 : #:for descs, arr, idx in umgga_all
2080 : e_${'_'.join(descs)}$, &
2081 : #:endfor
2082 : grad_deriv, npoints, func_name, sc, workers)
2083 :
2084 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho, laplace_rho, tau
2085 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, &
2086 : e_rho_rho, e_ndrho_rho, e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
2087 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_tau_tau, e_rho_rho_rho
2088 : ! Derivatives with respect to the reduced gradient gamma = |grad rho|^2.
2089 : ! These are LibXC's own sigma derivatives, copied through unchanged: unlike
2090 : ! the norm_drho ones they need no chain rule and carry no 1/|grad rho|.
2091 : #:for descs, arr, idx in umgga_all
2092 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_${'_'.join(descs)}$
2093 : #:endfor
2094 : INTEGER, INTENT(in) :: grad_deriv, npoints
2095 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
2096 : REAL(KIND=dp), INTENT(in) :: sc
2097 : TYPE(libxc_worker_set_type), INTENT(INOUT), TARGET :: workers
2098 :
2099 : INTEGER :: bsize, family, i, i0, ib, ii, ithread, &
2100 : nb, nblocks, nthreads
2101 : INTEGER(C_SIZE_T) :: np
2102 : LOGICAL :: is_gga, is_mgga
2103 :
2104 18290 : ithread = 0
2105 18290 : nthreads = 1
2106 18290 : !$ ithread = omp_get_thread_num()
2107 18290 : !$ nthreads = omp_get_num_threads()
2108 18290 : CPASSERT(ithread < SIZE(workers%worker))
2109 :
2110 : ! Blocks are capped by what the staging buffers were actually allocated for,
2111 : ! but shrunk when there are too few points to give every thread a block of its
2112 : ! own: small atomic grids would otherwise all be evaluated by thread 0.
2113 18290 : bsize = MAX(1, MIN(workers%worker(1)%nb, (npoints + nthreads - 1)/nthreads))
2114 18290 : nblocks = (npoints + bsize - 1)/bsize
2115 :
2116 18290 : family = workers%family
2117 : SELECT CASE (family)
2118 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
2119 : is_gga = .FALSE.
2120 5456 : is_mgga = .FALSE.
2121 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
2122 5456 : is_gga = .TRUE.
2123 5456 : is_mgga = .FALSE.
2124 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
2125 2348 : is_gga = .TRUE.
2126 2348 : is_mgga = .TRUE.
2127 : CASE default
2128 0 : is_gga = .FALSE.
2129 0 : is_mgga = .FALSE.
2130 18290 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
2131 : END SELECT
2132 :
2133 : ! Points below the density cutoff need no test here. The cutoffs were handed
2134 : ! to the functional object in xc_libxc_wrap_set_thresholds, and LibXC leaves
2135 : ! the outputs of the points it screens at zero, so accumulating them is a
2136 : ! no-op.
2137 : BLOCK
2138 : TYPE(libxc_worker_type), POINTER :: w
2139 : LOGICAL :: has_laplace, no_exc
2140 : REAL(KIND=dp) :: eps_rho, eps_tau
2141 :
2142 : ! Avoid ASSOCIATE aliases across the orphaned OpenMP worksharing loop:
2143 : ! ifx can lose their bindings even at -O0.
2144 18290 : w => workers%worker(ithread + 1)
2145 18290 : has_laplace = workers%has_laplace
2146 18290 : no_exc = workers%no_exc
2147 18290 : eps_rho = workers%epsilon_rho
2148 18290 : eps_tau = workers%epsilon_tau
2149 : !$OMP DO
2150 : DO ib = 1, nblocks
2151 1083874 : i0 = (ib - 1)*bsize
2152 1083874 : nb = MIN(bsize, npoints - i0)
2153 1083874 : np = INT(nb, KIND=C_SIZE_T)
2154 :
2155 : ! stage this block's inputs in LibXC's layout
2156 554048992 : DO i = 1, nb
2157 554048992 : w%rho(1, i) = rho(i0 + i)
2158 : END DO
2159 1083874 : IF (is_gga) THEN
2160 298335134 : DO i = 1, nb
2161 298335134 : w%sigma(1, i) = norm_drho(i0 + i)**2
2162 : END DO
2163 : END IF
2164 1083874 : IF (is_mgga) THEN
2165 118937092 : DO i = 1, nb
2166 118704312 : ii = i0 + i
2167 118704312 : w%lapl(1, i) = laplace_rho(ii)
2168 118937092 : w%tau(1, i) = tau(ii)
2169 : END DO
2170 : ! Meta-GGAs are screened on the kinetic energy density as well as on the
2171 : ! density, and LibXC only screens on the latter. A point failing the tau
2172 : ! test is handed over with zero density, which makes LibXC screen it out
2173 : ! and leave its outputs at zero: what skipping it used to achieve. The
2174 : ! Fermi hole curvature bound tau >= sigma/(8*rho) follows in the same pass,
2175 : ! with the denominator floored at the cutoff so the division stays defined
2176 : ! for the points that are about to be screened.
2177 118937092 : DO i = 1, nb
2178 118704312 : IF (w%tau(1, i) <= eps_tau) w%rho(1, i) = 0.0_dp
2179 : w%tau(1, i) = MAX(w%tau(1, i), &
2180 118937092 : w%sigma(1, i)/(8.0_dp*MAX(w%rho(1, i), eps_rho)))
2181 : END DO
2182 : END IF
2183 :
2184 : ! one LibXC call for the whole block
2185 : IF (is_mgga) THEN
2186 233824 : SELECT CASE (grad_deriv)
2187 : CASE (0)
2188 1044 : CALL xc_f03_mgga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), w%exc(1))
2189 : CASE (1)
2190 217026 : IF (no_exc) THEN
2191 : CALL xc_f03_mgga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2192 0 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2193 0 : w%exc(1:nb) = 0.0_dp
2194 : ELSE
2195 : CALL xc_f03_mgga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2196 217026 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2197 : END IF
2198 : CASE (2)
2199 13508 : IF (no_exc) THEN
2200 : CALL xc_f03_mgga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2201 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2202 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2203 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2204 0 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2205 0 : w%exc(1:nb) = 0.0_dp
2206 : ELSE
2207 : CALL xc_f03_mgga(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2208 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2209 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2210 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2211 13508 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2212 : END IF
2213 : CASE (3)
2214 232780 : IF (no_exc) THEN
2215 : CALL xc_f03_mgga_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2216 : w%lapl(1, 1), w%tau(1, 1), &
2217 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2218 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2219 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2220 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1), &
2221 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), w%v3rho2lapl(1, 1), w%v3rho2tau(1, 1), &
2222 : w%v3rhosigma2(1, 1), w%v3rhosigmalapl(1, 1), w%v3rhosigmatau(1, 1), &
2223 : w%v3rholapl2(1, 1), w%v3rholapltau(1, 1), w%v3rhotau2(1, 1), &
2224 : w%v3sigma3(1, 1), w%v3sigma2lapl(1, 1), w%v3sigma2tau(1, 1), &
2225 : w%v3sigmalapl2(1, 1), w%v3sigmalapltau(1, 1), w%v3sigmatau2(1, 1), &
2226 0 : w%v3lapl3(1, 1), w%v3lapl2tau(1, 1), w%v3lapltau2(1, 1), w%v3tau3(1, 1))
2227 0 : w%exc(1:nb) = 0.0_dp
2228 : ELSE
2229 : CALL xc_f03_mgga_exc_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2230 : w%lapl(1, 1), w%tau(1, 1), w%exc(1), &
2231 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2232 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2233 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2234 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1), &
2235 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), w%v3rho2lapl(1, 1), w%v3rho2tau(1, 1), &
2236 : w%v3rhosigma2(1, 1), w%v3rhosigmalapl(1, 1), w%v3rhosigmatau(1, 1), &
2237 : w%v3rholapl2(1, 1), w%v3rholapltau(1, 1), w%v3rhotau2(1, 1), &
2238 : w%v3sigma3(1, 1), w%v3sigma2lapl(1, 1), w%v3sigma2tau(1, 1), &
2239 : w%v3sigmalapl2(1, 1), w%v3sigmalapltau(1, 1), w%v3sigmatau2(1, 1), &
2240 1202 : w%v3lapl3(1, 1), w%v3lapl2tau(1, 1), w%v3lapltau2(1, 1), w%v3tau3(1, 1))
2241 : END IF
2242 : END SELECT
2243 851094 : ELSE IF (is_gga) THEN
2244 362958 : SELECT CASE (grad_deriv)
2245 : CASE (0)
2246 11962 : CALL xc_f03_gga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%exc(1))
2247 : CASE (1)
2248 329470 : IF (no_exc) THEN
2249 0 : CALL xc_f03_gga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%vrho(1, 1), w%vsigma(1, 1))
2250 0 : w%exc(1:nb) = 0.0_dp
2251 : ELSE
2252 : CALL xc_f03_gga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2253 329470 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1))
2254 : END IF
2255 : CASE (2)
2256 8848 : IF (no_exc) THEN
2257 : CALL xc_f03_gga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2258 : w%vrho(1, 1), w%vsigma(1, 1), &
2259 0 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2260 0 : w%exc(1:nb) = 0.0_dp
2261 : ELSE
2262 : CALL xc_f03_gga_exc_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2263 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
2264 8848 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2265 : END IF
2266 : CASE (3)
2267 350996 : IF (no_exc) THEN
2268 : CALL xc_f03_gga_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2269 : w%vrho(1, 1), w%vsigma(1, 1), &
2270 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1), &
2271 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), &
2272 0 : w%v3rhosigma2(1, 1), w%v3sigma3(1, 1))
2273 0 : w%exc(1:nb) = 0.0_dp
2274 : ELSE
2275 : CALL xc_f03_gga_exc_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2276 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
2277 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1), &
2278 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), &
2279 716 : w%v3rhosigma2(1, 1), w%v3sigma3(1, 1))
2280 : END IF
2281 : END SELECT
2282 : ELSE
2283 512640 : SELECT CASE (grad_deriv)
2284 : CASE (0)
2285 12542 : CALL xc_f03_lda_exc(w%func, np, w%rho(1, 1), w%exc(1))
2286 : CASE (1)
2287 464244 : CALL xc_f03_lda_exc_vxc(w%func, np, w%rho(1, 1), w%exc(1), w%vrho(1, 1))
2288 : CASE (2)
2289 : CALL xc_f03_lda_exc_vxc_fxc(w%func, np, w%rho(1, 1), &
2290 22596 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1))
2291 : CASE (3)
2292 : CALL xc_f03_lda(w%func, np, w%rho(1, 1), &
2293 500098 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1), w%v3rho3(1, 1))
2294 : END SELECT
2295 : END IF
2296 :
2297 : ! accumulate; each derivative is added by exactly one loop, selected by
2298 : ! the order requested and the family, not by a per-branch copy of the whole
2299 : ! block scaffolding
2300 1083874 : IF (grad_deriv >= 0) THEN
2301 554048992 : DO i = 1, nb
2302 552965118 : ii = i0 + i
2303 554048992 : e_0(ii) = e_0(ii) + sc*w%exc(i)*rho(ii)
2304 : END DO
2305 : END IF
2306 1083874 : IF (grad_deriv >= 1) THEN
2307 540981496 : DO i = 1, nb
2308 539923170 : ii = i0 + i
2309 540981496 : e_rho(ii) = e_rho(ii) + sc*w%vrho(1, i)
2310 : END DO
2311 1058326 : IF (is_gga) THEN
2312 291691134 : DO i = 1, nb
2313 291120364 : ii = i0 + i
2314 291120364 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*w%vsigma(1, i)*norm_drho(ii)
2315 291691134 : e_gamma(ii) = e_gamma(ii) + sc*w%vsigma(1, i)
2316 : END DO
2317 : END IF
2318 1058326 : IF (is_mgga) THEN
2319 118410048 : DO i = 1, nb
2320 118178312 : ii = i0 + i
2321 118410048 : e_tau(ii) = e_tau(ii) + sc*w%vtau(1, i)
2322 : END DO
2323 : END IF
2324 231736 : IF (is_mgga .AND. has_laplace) THEN
2325 14209482 : DO i = 1, nb
2326 14181507 : ii = i0 + i
2327 14209482 : e_laplace_rho(ii) = e_laplace_rho(ii) + sc*w%vlapl(1, i)
2328 : END DO
2329 : END IF
2330 : END IF
2331 1083874 : IF (grad_deriv >= 2) THEN
2332 24207320 : DO i = 1, nb
2333 24159734 : ii = i0 + i
2334 24207320 : e_rho_rho(ii) = e_rho_rho(ii) + sc*w%v2rho2(1, i)
2335 : END DO
2336 47586 : IF (is_gga) THEN
2337 12365008 : DO i = 1, nb
2338 12340734 : ii = i0 + i
2339 12340734 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*w%v2rhosigma(1, i)*norm_drho(ii)
2340 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2341 12340734 : sc*2.0_dp*(2.0_dp*w%sigma(1, i)*w%v2sigma2(1, i) + w%vsigma(1, i))
2342 : #:for descs, arr, idx in ugga_fill_2
2343 12365008 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2344 : #:endfor
2345 : END DO
2346 : END IF
2347 47586 : IF (is_mgga) THEN
2348 7511256 : DO i = 1, nb
2349 7496546 : ii = i0 + i
2350 : #:for descs, arr, idx in umgga_extra_2
2351 7496546 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2352 : #:endfor
2353 7496546 : e_rho_tau(ii) = e_rho_tau(ii) + sc*w%v2rhotau(1, i)
2354 7496546 : e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*w%v2sigmatau(1, i)*norm_drho(ii)
2355 7511256 : e_tau_tau(ii) = e_tau_tau(ii) + sc*w%v2tau2(1, i)
2356 : END DO
2357 : END IF
2358 14710 : IF (is_mgga .AND. has_laplace) THEN
2359 2234462 : DO i = 1, nb
2360 2230074 : ii = i0 + i
2361 2230074 : e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*w%v2rholapl(1, i)
2362 : e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
2363 2230074 : sc*2.0_dp*w%v2sigmalapl(1, i)*norm_drho(ii)
2364 2230074 : e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*w%v2lapl2(1, i)
2365 2234462 : e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*w%v2lapltau(1, i)
2366 : END DO
2367 : END IF
2368 : END IF
2369 1083874 : IF (grad_deriv >= 3) THEN
2370 1345884 : DO i = 1, nb
2371 1343250 : ii = i0 + i
2372 1345884 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*w%v3rho3(1, i)
2373 : END DO
2374 2634 : IF (is_mgga) THEN
2375 615452 : DO i = 1, nb
2376 614250 : ii = i0 + i
2377 : #:for descs, arr, idx in umgga_fill_3
2378 615452 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2379 : #:endfor
2380 : END DO
2381 1432 : ELSE IF (is_gga) THEN
2382 365216 : DO i = 1, nb
2383 364500 : ii = i0 + i
2384 : #:for descs, arr, idx in ugga_fill_3
2385 365216 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2386 : #:endfor
2387 : END DO
2388 : END IF
2389 : END IF
2390 : END DO
2391 : !$OMP END DO
2392 : END BLOCK
2393 :
2394 18290 : END SUBROUTINE libxc_spin_unpolarized_calc
2395 : #endif
2396 :
2397 : ! **************************************************************************************************
2398 : !> \brief libxc exchange-correlation functionals
2399 : !> \param rhoa alpha density
2400 : !> \param rhob beta density
2401 : !> \param norm_drho ...
2402 : !> \param norm_drhoa norm of the gradient of the alpha density
2403 : !> \param norm_drhob norm of the gradient of the beta density
2404 : !> \param laplace_rhoa laplacian of the alpha density
2405 : !> \param laplace_rhob laplacian of the beta density
2406 : !> \param tau_a alpha kinetic-energy density
2407 : !> \param tau_b beta kinetic-energy density
2408 : !> \param e_0 energy density
2409 : !> \param e_rhoa derivative of the energy density with respect to rhoa
2410 : !> \param e_rhob derivative of the energy density with respect to rhob
2411 : !> \param e_ndrho derivative of the energy density with respect to ndrho
2412 : !> \param e_ndrhoa derivative of the energy density with respect to ndrhoa
2413 : !> \param e_ndrhob derivative of the energy density with respect to ndrhob
2414 : !> \param e_laplace_rhoa derivative of the energy density with respect to laplace_rhoa
2415 : !> \param e_laplace_rhob derivative of the energy density with respect to laplace_rhob
2416 : !> \param e_tau_a derivative of the energy density with respect to tau_a
2417 : !> \param e_tau_b derivative of the energy density with respect to tau_b
2418 : !> \param e_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa
2419 : !> \param e_rhoa_rhob derivative of the energy density with respect to rhoa_rhob
2420 : !> \param e_rhob_rhob derivative of the energy density with respect to rhob_rhob
2421 : !> \param e_ndrho_rhoa derivative of the energy density with respect to ndrho_rhoa
2422 : !> \param e_ndrho_rhob derivative of the energy density with respect to ndrho_rhob
2423 : !> \param e_ndrhoa_rhoa derivative of the energy density with respect to ndrhoa_rhoa
2424 : !> \param e_ndrhoa_rhob derivative of the energy density with respect to ndrhoa_rhob
2425 : !> \param e_ndrhob_rhoa derivative of the energy density with respect to ndrhob_rhoa
2426 : !> \param e_ndrhob_rhob derivative of the energy density with respect to ndrhob_rhob
2427 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
2428 : !> \param e_ndrho_ndrhoa derivative of the energy density with respect to ndrho_ndrhoa
2429 : !> \param e_ndrho_ndrhob derivative of the energy density with respect to ndrho_ndrhob
2430 : !> \param e_ndrhoa_ndrhoa derivative of the energy density with respect to ndrhoa_ndrhoa
2431 : !> \param e_ndrhoa_ndrhob derivative of the energy density with respect to ndrhoa_ndrhob
2432 : !> \param e_ndrhob_ndrhob derivative of the energy density with respect to ndrhob_ndrhob
2433 : !> \param e_rhoa_laplace_rhoa derivative of the energy density with respect to rhoa_laplace_rhoa
2434 : !> \param e_rhoa_laplace_rhob derivative of the energy density with respect to rhoa_laplace_rhob
2435 : !> \param e_rhob_laplace_rhoa derivative of the energy density with respect to rhob_laplace_rhoa
2436 : !> \param e_rhob_laplace_rhob derivative of the energy density with respect to rhob_laplace_rhob
2437 : !> \param e_rhoa_tau_a derivative of the energy density with respect to rhoa_tau_a
2438 : !> \param e_rhoa_tau_b derivative of the energy density with respect to rhoa_tau_b
2439 : !> \param e_rhob_tau_a derivative of the energy density with respect to rhob_tau_a
2440 : !> \param e_rhob_tau_b derivative of the energy density with respect to rhob_tau_b
2441 : !> \param e_ndrho_laplace_rhoa derivative of the energy density with respect to ndrho_laplace_rhoa
2442 : !> \param e_ndrho_laplace_rhob derivative of the energy density with respect to ndrho_laplace_rhob
2443 : !> \param e_ndrhoa_laplace_rhoa derivative of the energy density with respect to ndrhoa_laplace_rhoa
2444 : !> \param e_ndrhoa_laplace_rhob derivative of the energy density with respect to ndrhoa_laplace_rhob
2445 : !> \param e_ndrhob_laplace_rhoa derivative of the energy density with respect to ndrhob_laplace_rhoa
2446 : !> \param e_ndrhob_laplace_rhob derivative of the energy density with respect to ndrhob_laplace_rhob
2447 : !> \param e_ndrho_tau_a derivative of the energy density with respect to ndrho_tau_a
2448 : !> \param e_ndrho_tau_b derivative of the energy density with respect to ndrho_tau_b
2449 : !> \param e_ndrhoa_tau_a derivative of the energy density with respect to ndrhoa_tau_a
2450 : !> \param e_ndrhoa_tau_b derivative of the energy density with respect to ndrhoa_tau_b
2451 : !> \param e_ndrhob_tau_a derivative of the energy density with respect to ndrhob_tau_a
2452 : !> \param e_ndrhob_tau_b derivative of the energy density with respect to ndrhob_tau_b
2453 : !> \param e_laplace_rhoa_laplace_rhoa derivative of the energy density with respect to laplace_rhoa_laplace_rhoa
2454 : !> \param e_laplace_rhoa_laplace_rhob derivative of the energy density with respect to laplace_rhoa_laplace_rhob
2455 : !> \param e_laplace_rhob_laplace_rhob derivative of the energy density with respect to laplace_rhob_laplace_rhob
2456 : !> \param e_laplace_rhoa_tau_a derivative of the energy density with respect to laplace_rhoa_tau_a
2457 : !> \param e_laplace_rhoa_tau_b derivative of the energy density with respect to laplace_rhoa_tau_b
2458 : !> \param e_laplace_rhob_tau_a derivative of the energy density with respect to laplace_rhob_tau_a
2459 : !> \param e_laplace_rhob_tau_b derivative of the energy density with respect to laplace_rhob_tau_b
2460 : !> \param e_tau_a_tau_a derivative of the energy density with respect to tau_a_tau_a
2461 : !> \param e_tau_a_tau_b derivative of the energy density with respect to tau_a_tau_b
2462 : !> \param e_tau_b_tau_b derivative of the energy density with respect to tau_b_tau_b
2463 : !> \param e_rhoa_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa_rhoa
2464 : !> \param e_rhoa_rhoa_rhob derivative of the energy density with respect to rhoa_rhoa_rhob
2465 : !> \param e_rhoa_rhob_rhob derivative of the energy density with respect to rhoa_rhob_rhob
2466 : !> \param e_rhob_rhob_rhob derivative of the energy density with respect to rhob_rhob_rhob
2467 : !> \param grad_deriv degree of the derivative that should be evaluated;
2468 : !> all derivatives up to the given degree are evaluated, in a single
2469 : !> LibXC call per block of grid points
2470 : !> \param npoints number of points on the grid
2471 : !> \param func_name name of the functional
2472 : !> \param sc scaling factor of the functional
2473 : !> \param workers cached LibXC functional objects and staging buffers, one per thread
2474 : !> \author F. Tran
2475 : ! **************************************************************************************************
2476 : #if defined (__LIBXC)
2477 3552 : SUBROUTINE libxc_spin_polarized_calc(rhoa, rhob, norm_drho, norm_drhoa, &
2478 : norm_drhob, laplace_rhoa, laplace_rhob, tau_a, tau_b, &
2479 : e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, e_ndrhob, &
2480 : e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, &
2481 : e_rhoa_rhoa, e_rhoa_rhob, e_rhob_rhob, &
2482 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, &
2483 : e_ndrhoa_rhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
2484 : e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
2485 : e_ndrhoa_ndrhoa, e_ndrhoa_ndrhob, e_ndrhob_ndrhob, &
2486 : e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
2487 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, &
2488 : e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, e_rhob_tau_b, &
2489 : e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, &
2490 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, &
2491 : e_ndrhob_laplace_rhoa, e_ndrhob_laplace_rhob, &
2492 : e_ndrho_tau_a, e_ndrho_tau_b, &
2493 : e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
2494 : e_ndrhob_tau_a, e_ndrhob_tau_b, &
2495 : e_laplace_rhoa_laplace_rhoa, &
2496 : e_laplace_rhoa_laplace_rhob, &
2497 : e_laplace_rhob_laplace_rhob, &
2498 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
2499 : e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, &
2500 : e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
2501 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
2502 : e_rhoa_rhob_rhob, e_rhob_rhob_rhob, &
2503 : #:for descs, arr, idx in mgga_all
2504 : e_${'_'.join(descs)}$, &
2505 : #:endfor
2506 : grad_deriv, npoints, func_name, sc, workers)
2507 :
2508 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, norm_drho, norm_drhoa, &
2509 : norm_drhob, laplace_rhoa, &
2510 : laplace_rhob, tau_a, tau_b
2511 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, &
2512 : e_ndrhob, e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, e_rhoa_rhoa, e_rhoa_rhob, &
2513 : e_rhob_rhob, e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhob_rhoa, &
2514 : e_ndrhob_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, e_ndrhoa_ndrhoa, &
2515 : e_ndrhoa_ndrhob, e_ndrhob_ndrhob, e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
2516 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, &
2517 : e_rhob_tau_b, e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, e_ndrhoa_laplace_rhoa
2518 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_ndrhoa_laplace_rhob, e_ndrhob_laplace_rhoa, &
2519 : e_ndrhob_laplace_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
2520 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
2521 : e_laplace_rhob_laplace_rhob, e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
2522 : e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
2523 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, e_rhoa_rhob_rhob, e_rhob_rhob_rhob
2524 : ! Derivatives with respect to the reduced gradients gamma_ij; these are
2525 : ! LibXC's own sigma derivatives, copied through without a chain rule.
2526 : #:for descs, arr, idx in mgga_all
2527 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_${'_'.join(descs)}$
2528 : #:endfor
2529 : INTEGER, INTENT(in) :: grad_deriv, npoints
2530 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
2531 : REAL(KIND=dp), INTENT(in) :: sc
2532 : TYPE(libxc_worker_set_type), INTENT(INOUT), TARGET :: workers
2533 :
2534 : INTEGER :: bsize, family, i, i0, ib, ii, ithread, &
2535 : nb, nblocks, nthreads
2536 : INTEGER(C_SIZE_T) :: np
2537 : LOGICAL :: is_gga, is_mgga
2538 :
2539 3552 : ithread = 0
2540 3552 : nthreads = 1
2541 3552 : !$ ithread = omp_get_thread_num()
2542 3552 : !$ nthreads = omp_get_num_threads()
2543 3552 : CPASSERT(ithread < SIZE(workers%worker))
2544 :
2545 3552 : bsize = MAX(1, MIN(workers%worker(1)%nb, (npoints + nthreads - 1)/nthreads))
2546 3552 : nblocks = (npoints + bsize - 1)/bsize
2547 :
2548 3552 : family = workers%family
2549 : SELECT CASE (family)
2550 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
2551 : is_gga = .FALSE.
2552 634 : is_mgga = .FALSE.
2553 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
2554 634 : is_gga = .TRUE.
2555 634 : is_mgga = .FALSE.
2556 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
2557 1318 : is_gga = .TRUE.
2558 1318 : is_mgga = .TRUE.
2559 : CASE default
2560 0 : is_gga = .FALSE.
2561 0 : is_mgga = .FALSE.
2562 3552 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
2563 : END SELECT
2564 :
2565 : ! As in the spin-unpolarized case the density and kinetic energy density
2566 : ! cutoffs are applied by LibXC itself.
2567 : BLOCK
2568 : TYPE(libxc_worker_type), POINTER :: w
2569 : LOGICAL :: has_laplace, no_exc
2570 : REAL(KIND=dp) :: eps_rho, eps_tau
2571 :
2572 : ! Keep the per-thread bindings explicit, as in the spin-unpolarized routine.
2573 3552 : w => workers%worker(ithread + 1)
2574 3552 : has_laplace = workers%has_laplace
2575 3552 : no_exc = workers%no_exc
2576 3552 : eps_rho = workers%epsilon_rho
2577 3552 : eps_tau = workers%epsilon_tau
2578 : !$OMP DO
2579 : DO ib = 1, nblocks
2580 242430 : i0 = (ib - 1)*bsize
2581 242430 : nb = MIN(bsize, npoints - i0)
2582 242430 : np = INT(nb, KIND=C_SIZE_T)
2583 :
2584 : ! stage this block's inputs in LibXC's layout
2585 110186815 : DO i = 1, nb
2586 109944385 : ii = i0 + i
2587 109944385 : w%rho(1, i) = MAX(rhoa(ii), 0.0_dp)
2588 110186815 : w%rho(2, i) = MAX(rhob(ii), 0.0_dp)
2589 : END DO
2590 242430 : IF (is_gga) THEN
2591 : ! CP2K works with the norms of the gradients, LibXC with their contractions
2592 39322809 : DO i = 1, nb
2593 39218797 : ii = i0 + i
2594 39218797 : w%nda(i) = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2595 39218797 : w%ndb(i) = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2596 39218797 : w%nd(i) = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2597 39218797 : w%sigma(1, i) = w%nda(i)**2
2598 39218797 : w%sigma(3, i) = w%ndb(i)**2
2599 39322809 : w%sigma(2, i) = 0.5_dp*(w%nd(i)**2 - w%sigma(1, i) - w%sigma(3, i))
2600 : END DO
2601 : END IF
2602 242430 : IF (is_mgga) THEN
2603 20149719 : DO i = 1, nb
2604 20083269 : ii = i0 + i
2605 20083269 : w%lapl(1, i) = laplace_rhoa(ii)
2606 20083269 : w%lapl(2, i) = laplace_rhob(ii)
2607 20083269 : w%tau(1, i) = MAX(tau_a(ii), 0.0_dp)
2608 20149719 : w%tau(2, i) = MAX(tau_b(ii), 0.0_dp)
2609 : END DO
2610 : ! screening on tau and the Fermi hole curvature bound, per spin channel;
2611 : ! see the spin-unpolarized routine. The tau test is on the sum of the two
2612 : ! channels, as the density test is.
2613 20149719 : DO i = 1, nb
2614 20083269 : IF (w%tau(1, i) + w%tau(2, i) <= eps_tau) THEN
2615 10732 : w%rho(1, i) = 0.0_dp
2616 10732 : w%rho(2, i) = 0.0_dp
2617 : END IF
2618 : w%tau(1, i) = MAX(w%tau(1, i), &
2619 20083269 : w%sigma(1, i)/(8.0_dp*MAX(w%rho(1, i), eps_rho)))
2620 : w%tau(2, i) = MAX(w%tau(2, i), &
2621 20149719 : w%sigma(3, i)/(8.0_dp*MAX(w%rho(2, i), eps_rho)))
2622 : END DO
2623 : END IF
2624 :
2625 : ! one LibXC call for the whole block
2626 : IF (is_mgga) THEN
2627 67092 : SELECT CASE (grad_deriv)
2628 : CASE (0)
2629 642 : CALL xc_f03_mgga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), w%exc(1))
2630 : CASE (1)
2631 54350 : IF (no_exc) THEN
2632 : CALL xc_f03_mgga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2633 0 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2634 0 : w%exc(1:nb) = 0.0_dp
2635 : ELSE
2636 : CALL xc_f03_mgga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2637 54350 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2638 : END IF
2639 : CASE (2)
2640 8836 : IF (no_exc) THEN
2641 : CALL xc_f03_mgga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2642 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2643 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2644 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2645 0 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2646 0 : w%exc(1:nb) = 0.0_dp
2647 : ELSE
2648 : CALL xc_f03_mgga(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2649 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2650 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2651 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2652 8836 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2653 : END IF
2654 : CASE (3)
2655 66450 : IF (no_exc) THEN
2656 : CALL xc_f03_mgga_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2657 : w%lapl(1, 1), w%tau(1, 1), &
2658 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2659 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2660 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2661 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1), &
2662 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), w%v3rho2lapl(1, 1), w%v3rho2tau(1, 1), &
2663 : w%v3rhosigma2(1, 1), w%v3rhosigmalapl(1, 1), w%v3rhosigmatau(1, 1), &
2664 : w%v3rholapl2(1, 1), w%v3rholapltau(1, 1), w%v3rhotau2(1, 1), &
2665 : w%v3sigma3(1, 1), w%v3sigma2lapl(1, 1), w%v3sigma2tau(1, 1), &
2666 : w%v3sigmalapl2(1, 1), w%v3sigmalapltau(1, 1), w%v3sigmatau2(1, 1), &
2667 0 : w%v3lapl3(1, 1), w%v3lapl2tau(1, 1), w%v3lapltau2(1, 1), w%v3tau3(1, 1))
2668 0 : w%exc(1:nb) = 0.0_dp
2669 : ELSE
2670 : CALL xc_f03_mgga_exc_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2671 : w%lapl(1, 1), w%tau(1, 1), w%exc(1), &
2672 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2673 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2674 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2675 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1), &
2676 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), w%v3rho2lapl(1, 1), w%v3rho2tau(1, 1), &
2677 : w%v3rhosigma2(1, 1), w%v3rhosigmalapl(1, 1), w%v3rhosigmatau(1, 1), &
2678 : w%v3rholapl2(1, 1), w%v3rholapltau(1, 1), w%v3rhotau2(1, 1), &
2679 : w%v3sigma3(1, 1), w%v3sigma2lapl(1, 1), w%v3sigma2tau(1, 1), &
2680 : w%v3sigmalapl2(1, 1), w%v3sigmalapltau(1, 1), w%v3sigmatau2(1, 1), &
2681 2622 : w%v3lapl3(1, 1), w%v3lapl2tau(1, 1), w%v3lapltau2(1, 1), w%v3tau3(1, 1))
2682 : END IF
2683 : END SELECT
2684 175980 : ELSE IF (is_gga) THEN
2685 38074 : SELECT CASE (grad_deriv)
2686 : CASE (0)
2687 512 : CALL xc_f03_gga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%exc(1))
2688 : CASE (1)
2689 32152 : IF (no_exc) THEN
2690 0 : CALL xc_f03_gga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%vrho(1, 1), w%vsigma(1, 1))
2691 0 : w%exc(1:nb) = 0.0_dp
2692 : ELSE
2693 : CALL xc_f03_gga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2694 32152 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1))
2695 : END IF
2696 : CASE (2)
2697 4182 : IF (no_exc) THEN
2698 : CALL xc_f03_gga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2699 : w%vrho(1, 1), w%vsigma(1, 1), &
2700 0 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2701 0 : w%exc(1:nb) = 0.0_dp
2702 : ELSE
2703 : CALL xc_f03_gga_exc_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2704 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
2705 4182 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2706 : END IF
2707 : CASE (3)
2708 37562 : IF (no_exc) THEN
2709 : CALL xc_f03_gga_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2710 : w%vrho(1, 1), w%vsigma(1, 1), &
2711 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1), &
2712 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), &
2713 0 : w%v3rhosigma2(1, 1), w%v3sigma3(1, 1))
2714 0 : w%exc(1:nb) = 0.0_dp
2715 : ELSE
2716 : CALL xc_f03_gga_exc_vxc_fxc_kxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2717 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
2718 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1), &
2719 : w%v3rho3(1, 1), w%v3rho2sigma(1, 1), &
2720 716 : w%v3rhosigma2(1, 1), w%v3sigma3(1, 1))
2721 : END IF
2722 : END SELECT
2723 : ELSE
2724 138418 : SELECT CASE (grad_deriv)
2725 : CASE (0)
2726 0 : CALL xc_f03_lda_exc(w%func, np, w%rho(1, 1), w%exc(1))
2727 : CASE (1)
2728 133102 : CALL xc_f03_lda_exc_vxc(w%func, np, w%rho(1, 1), w%exc(1), w%vrho(1, 1))
2729 : CASE (2)
2730 : CALL xc_f03_lda_exc_vxc_fxc(w%func, np, w%rho(1, 1), &
2731 4600 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1))
2732 : CASE (3)
2733 : CALL xc_f03_lda(w%func, np, w%rho(1, 1), &
2734 138418 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1), w%v3rho3(1, 1))
2735 : END SELECT
2736 : END IF
2737 :
2738 : ! accumulate; one loop per group of derivatives, selected by the order
2739 : ! requested and the family
2740 242430 : IF (grad_deriv >= 0) THEN
2741 110186815 : DO i = 1, nb
2742 109944385 : ii = i0 + i
2743 110186815 : e_0(ii) = e_0(ii) + sc*w%exc(i)*(w%rho(1, i) + w%rho(2, i))
2744 : END DO
2745 : END IF
2746 242430 : IF (grad_deriv >= 1) THEN
2747 109608601 : DO i = 1, nb
2748 109367325 : ii = i0 + i
2749 109367325 : e_rhoa(ii) = e_rhoa(ii) + sc*w%vrho(1, i)
2750 109608601 : e_rhob(ii) = e_rhob(ii) + sc*w%vrho(2, i)
2751 : END DO
2752 241276 : IF (is_gga) THEN
2753 38744595 : DO i = 1, nb
2754 38641737 : ii = i0 + i
2755 38641737 : e_ndrho(ii) = e_ndrho(ii) + sc*w%vsigma(2, i)*w%nd(i)
2756 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2757 38641737 : sc*(2.0_dp*w%vsigma(1, i) - w%vsigma(2, i))*w%nda(i)
2758 : e_ndrhob(ii) = e_ndrhob(ii) + &
2759 38744595 : sc*(2.0_dp*w%vsigma(3, i) - w%vsigma(2, i))*w%ndb(i)
2760 : END DO
2761 : END IF
2762 241276 : IF (is_mgga) THEN
2763 19834161 : DO i = 1, nb
2764 19768353 : ii = i0 + i
2765 19768353 : e_tau_a(ii) = e_tau_a(ii) + sc*w%vtau(1, i)
2766 19834161 : e_tau_b(ii) = e_tau_b(ii) + sc*w%vtau(2, i)
2767 : END DO
2768 : END IF
2769 65808 : IF (is_mgga .AND. has_laplace) THEN
2770 6463239 : DO i = 1, nb
2771 6438285 : ii = i0 + i
2772 6438285 : e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*w%vlapl(1, i)
2773 6463239 : e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*w%vlapl(2, i)
2774 : END DO
2775 : END IF
2776 : END IF
2777 242430 : IF (grad_deriv >= 2) THEN
2778 7144194 : DO i = 1, nb
2779 7122522 : ii = i0 + i
2780 7122522 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*w%v2rho2(1, i)
2781 7122522 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*w%v2rho2(2, i)
2782 7144194 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*w%v2rho2(3, i)
2783 : END DO
2784 21672 : IF (is_gga) THEN
2785 4443030 : DO i = 1, nb
2786 4426674 : ii = i0 + i
2787 4426674 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*w%v2rhosigma(2, i)*w%nd(i)
2788 4426674 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*w%v2rhosigma(5, i)*w%nd(i)
2789 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2790 4426674 : sc*(2.0_dp*w%v2rhosigma(1, i) - w%v2rhosigma(2, i))*w%nda(i)
2791 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2792 4426674 : sc*(2.0_dp*w%v2rhosigma(4, i) - w%v2rhosigma(5, i))*w%nda(i)
2793 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2794 4426674 : sc*(2.0_dp*w%v2rhosigma(3, i) - w%v2rhosigma(2, i))*w%ndb(i)
2795 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2796 4426674 : sc*(2.0_dp*w%v2rhosigma(6, i) - w%v2rhosigma(5, i))*w%ndb(i)
2797 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2798 4426674 : sc*(w%vsigma(2, i) + w%nd(i)**2*w%v2sigma2(4, i))
2799 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2800 4426674 : sc*(2.0_dp*w%v2sigma2(2, i) - w%v2sigma2(4, i))*w%nd(i)*w%nda(i)
2801 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2802 4426674 : sc*(2.0_dp*w%v2sigma2(5, i) - w%v2sigma2(4, i))*w%nd(i)*w%ndb(i)
2803 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2804 : sc*(2.0_dp*w%vsigma(1, i) - w%vsigma(2, i) + w%nda(i)**2*( &
2805 4426674 : 4.0_dp*w%v2sigma2(1, i) - 4.0_dp*w%v2sigma2(2, i) + w%v2sigma2(4, i)))
2806 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2807 : sc*(4.0_dp*w%v2sigma2(3, i) - 2.0_dp*w%v2sigma2(2, i) - &
2808 4426674 : 2.0_dp*w%v2sigma2(5, i) + w%v2sigma2(4, i))*w%nda(i)*w%ndb(i)
2809 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2810 : sc*(2.0_dp*w%vsigma(3, i) - w%vsigma(2, i) + w%ndb(i)**2*( &
2811 4443030 : 4.0_dp*w%v2sigma2(6, i) - 4.0_dp*w%v2sigma2(5, i) + w%v2sigma2(4, i)))
2812 : END DO
2813 : END IF
2814 21672 : IF (is_mgga) THEN
2815 1952596 : DO i = 1, nb
2816 1941138 : ii = i0 + i
2817 1941138 : e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*w%v2rhotau(1, i)
2818 1941138 : e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*w%v2rhotau(2, i)
2819 1941138 : e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*w%v2rhotau(3, i)
2820 1941138 : e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*w%v2rhotau(4, i)
2821 1941138 : e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*w%v2sigmatau(3, i)*w%nd(i)
2822 1941138 : e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*w%v2sigmatau(4, i)*w%nd(i)
2823 : e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2824 1941138 : sc*(2.0_dp*w%v2sigmatau(1, i) - w%v2sigmatau(3, i))*w%nda(i)
2825 : e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2826 1941138 : sc*(2.0_dp*w%v2sigmatau(2, i) - w%v2sigmatau(4, i))*w%nda(i)
2827 : e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2828 1941138 : sc*(2.0_dp*w%v2sigmatau(5, i) - w%v2sigmatau(3, i))*w%ndb(i)
2829 : e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2830 1941138 : sc*(2.0_dp*w%v2sigmatau(6, i) - w%v2sigmatau(4, i))*w%ndb(i)
2831 1941138 : e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*w%v2tau2(1, i)
2832 1941138 : e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*w%v2tau2(2, i)
2833 1952596 : e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*w%v2tau2(3, i)
2834 : END DO
2835 : END IF
2836 11458 : IF (is_mgga .AND. has_laplace) THEN
2837 958364 : DO i = 1, nb
2838 952722 : ii = i0 + i
2839 952722 : e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*w%v2rholapl(1, i)
2840 952722 : e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*w%v2rholapl(2, i)
2841 952722 : e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*w%v2rholapl(3, i)
2842 952722 : e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*w%v2rholapl(4, i)
2843 952722 : e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*w%v2sigmalapl(3, i)*w%nd(i)
2844 952722 : e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*w%v2sigmalapl(4, i)*w%nd(i)
2845 : e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2846 952722 : sc*(2.0_dp*w%v2sigmalapl(1, i) - w%v2sigmalapl(3, i))*w%nda(i)
2847 : e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2848 952722 : sc*(2.0_dp*w%v2sigmalapl(2, i) - w%v2sigmalapl(4, i))*w%nda(i)
2849 : e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2850 952722 : sc*(2.0_dp*w%v2sigmalapl(5, i) - w%v2sigmalapl(3, i))*w%ndb(i)
2851 : e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2852 952722 : sc*(2.0_dp*w%v2sigmalapl(6, i) - w%v2sigmalapl(4, i))*w%ndb(i)
2853 952722 : e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*w%v2lapl2(1, i)
2854 952722 : e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*w%v2lapl2(2, i)
2855 952722 : e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*w%v2lapl2(3, i)
2856 952722 : e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*w%v2lapltau(1, i)
2857 952722 : e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*w%v2lapltau(2, i)
2858 952722 : e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*w%v2lapltau(3, i)
2859 958364 : e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*w%v2lapltau(4, i)
2860 : END DO
2861 : END IF
2862 : END IF
2863 242430 : IF (grad_deriv >= 1 .AND. is_gga) THEN
2864 38744595 : DO i = 1, nb
2865 38641737 : ii = i0 + i
2866 : #:for descs, arr, idx in gamma_only_1
2867 38744595 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2868 : #:endfor
2869 : END DO
2870 : END IF
2871 242430 : IF (grad_deriv >= 2 .AND. is_gga) THEN
2872 4443030 : DO i = 1, nb
2873 4426674 : ii = i0 + i
2874 : #:for descs, arr, idx in gamma_only_2
2875 4443030 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2876 : #:endfor
2877 : END DO
2878 : END IF
2879 242430 : IF (grad_deriv >= 2 .AND. is_mgga) THEN
2880 1952596 : DO i = 1, nb
2881 1941138 : ii = i0 + i
2882 : #:for descs, arr, idx in [e for e in mgga_fill_2 if e not in gga_fill_2]
2883 1952596 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2884 : #:endfor
2885 : END DO
2886 : END IF
2887 242430 : IF (grad_deriv >= 3 .AND. is_mgga) THEN
2888 371496 : DO i = 1, nb
2889 368874 : ii = i0 + i
2890 : #:for descs, arr, idx in mgga_fill_3
2891 371496 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2892 : #:endfor
2893 : END DO
2894 : END IF
2895 242430 : IF (grad_deriv >= 3 .AND. is_gga .AND. .NOT. is_mgga) THEN
2896 365216 : DO i = 1, nb
2897 364500 : ii = i0 + i
2898 : #! the pure-density third derivatives are filled by the block below, which
2899 : #! runs for every family; adding them here as well would double-count them
2900 : #:for descs, arr, idx in gamma_only_3
2901 365216 : e_${'_'.join(descs)}$ (ii) = e_${'_'.join(descs)}$ (ii) + sc*w%${arr}$ (${idx}$, i)
2902 : #:endfor
2903 : END DO
2904 : END IF
2905 242430 : IF (grad_deriv >= 3) THEN
2906 1101928 : DO i = 1, nb
2907 1097874 : ii = i0 + i
2908 1097874 : e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*w%v3rho3(1, i)
2909 1097874 : e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*w%v3rho3(2, i)
2910 1097874 : e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*w%v3rho3(3, i)
2911 1101928 : e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*w%v3rho3(4, i)
2912 : END DO
2913 : END IF
2914 : END DO
2915 : !$OMP END DO
2916 : END BLOCK
2917 :
2918 3552 : END SUBROUTINE libxc_spin_polarized_calc
2919 : #endif
2920 :
2921 0 : END MODULE xc_libxc
|
