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