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 : MODULE xc_libxc
29 : USE bibliography, ONLY: Lehtola2018, &
30 : Marques2012, &
31 : cite_reference
32 : USE input_section_types, ONLY: section_add_keyword, &
33 : section_add_subsection, &
34 : section_create, &
35 : section_release, &
36 : section_type, &
37 : section_vals_type, &
38 : section_vals_val_get
39 : USE kinds, ONLY: default_string_length, &
40 : dp
41 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
42 : xc_dset_get_derivative
43 : USE xc_derivative_types, ONLY: xc_derivative_get, &
44 : xc_derivative_type
45 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
46 : USE xc_rho_set_types, ONLY: xc_rho_set_get, &
47 : xc_rho_set_type
48 : #if defined (__LIBXC)
49 : USE input_keyword_types, ONLY: keyword_create, &
50 : keyword_release, &
51 : keyword_type
52 : USE iso_c_binding, ONLY: C_SIZE_T, C_INT, C_DOUBLE
53 : USE xc_derivative_desc, ONLY: &
54 : deriv_rho, deriv_rhoa, deriv_rhob, &
55 : deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
56 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob
57 : USE xc_libxc_wrap, ONLY: xc_f03_func_t, &
58 : xc_f03_func_init, &
59 : xc_f03_func_end, &
60 : xc_f03_func_info_t, &
61 : xc_f03_functional_get_name, &
62 : xc_f03_func_get_info, &
63 : xc_f03_func_info_get_family, &
64 : xc_f03_func_info_get_kind, &
65 : xc_f03_func_info_get_n_ext_params, &
66 : xc_f03_func_info_get_name, &
67 : xc_f03_available_functional_numbers, &
68 : xc_f03_available_functional_names, &
69 : xc_f03_maximum_name_length, &
70 : xc_f03_number_of_functionals, &
71 : xc_f03_func_info_get_ext_params_name, &
72 : xc_f03_func_info_get_ext_params_description, &
73 : xc_f03_func_info_get_ext_params_default_value, &
74 : xc_f03_gga_exc, &
75 : xc_f03_gga_exc_vxc, &
76 : xc_f03_gga_exc_vxc_fxc, &
77 : xc_f03_gga_fxc, &
78 : xc_f03_gga_vxc, &
79 : xc_f03_gga_vxc_fxc, &
80 : xc_f03_lda, &
81 : xc_f03_lda_exc, &
82 : xc_f03_lda_exc_vxc, &
83 : xc_f03_lda_exc_vxc_fxc, &
84 : xc_f03_lda_fxc, &
85 : xc_f03_lda_kxc, &
86 : xc_f03_lda_vxc, &
87 : xc_f03_mgga, &
88 : xc_f03_mgga_exc, &
89 : xc_f03_mgga_exc_vxc, &
90 : xc_f03_mgga_fxc, &
91 : xc_f03_mgga_vxc, &
92 : xc_f03_mgga_vxc_fxc, &
93 : XC_POLARIZED, &
94 : XC_UNPOLARIZED, &
95 : XC_FAMILY_LDA, &
96 : XC_FAMILY_GGA, &
97 : XC_FAMILY_MGGA, &
98 : XC_FAMILY_HYB_LDA, &
99 : XC_FAMILY_HYB_GGA, &
100 : XC_FAMILY_HYB_MGGA, &
101 : XC_CORRELATION, &
102 : XC_EXCHANGE, &
103 : XC_EXCHANGE_CORRELATION, &
104 : XC_KINETIC, &
105 : xc_libxc_wrap_info_refs, &
106 : xc_libxc_wrap_version, &
107 : xc_libxc_wrap_library_reference, &
108 : xc_libxc_wrap_functional_get_number, &
109 : xc_libxc_wrap_info_needs_laplace, &
110 : xc_libxc_wrap_info_no_exc, &
111 : xc_libxc_wrap_set_thresholds, &
112 : xc_f03_func_set_ext_params, &
113 : xc_libxc_wrap_is_under_development, &
114 : xc_libxc_get_reference_length, &
115 : xc_libxc_check_functional
116 : #endif
117 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
118 : !$ omp_get_num_threads, &
119 : !$ omp_get_thread_num
120 :
121 : #include "../base/base_uses.f90"
122 :
123 : IMPLICIT NONE
124 : PRIVATE
125 :
126 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_libxc'
127 :
128 : PUBLIC :: libxc_spin_unpolarized_info, libxc_spin_unpolarized_eval, &
129 : libxc_spin_polarized_info, libxc_spin_polarized_eval, &
130 : libxc_version_info, libxc_library_reference, &
131 : libxc_get_reference_length, libxc_add_sections, &
132 : libxc_check_existence_in_libxc, libxc_release_workers
133 :
134 : #if defined (__LIBXC)
135 : ! Number of grid points handed to LibXC in a single call. LibXC amortizes its
136 : ! per-call work (argument checks, output initialization and, for functionals
137 : ! built by mixing components, a full allocate/free cycle of the component
138 : ! buffers) over the points in the call, so evaluating point by point paid that
139 : ! cost once per point. Blocking rather than passing the whole grid keeps the
140 : ! staging buffers inside the cache; the timing is flat between roughly 128 and
141 : ! 2048 points, so the exact value is not critical.
142 : INTEGER, PARAMETER, PRIVATE :: libxc_block_size = 512
143 :
144 : ! **************************************************************************************************
145 : !> \brief One LibXC functional object plus the buffers used to hand it a block of
146 : !> grid points. One of these is kept per OpenMP thread so that neither the
147 : !> functional object nor the buffers are rebuilt between batches.
148 : ! **************************************************************************************************
149 : TYPE libxc_worker_type
150 : TYPE(xc_f03_func_t) :: func = xc_f03_func_t()
151 : TYPE(xc_f03_func_info_t) :: info = xc_f03_func_info_t()
152 : LOGICAL :: is_init = .FALSE.
153 : ! staging buffers, dimensioned (component, point)
154 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rho, sigma, lapl, tau
155 : ! the floored gradient norms the spin-polarized chain rule is written in
156 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: nd, nda, ndb
157 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: exc
158 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: vrho, vsigma, vlapl, vtau
159 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: v2rho2, v2rhosigma, v2sigma2, &
160 : v2rholapl, v2rhotau, v2sigmalapl, &
161 : v2sigmatau, v2lapl2, v2lapltau, v2tau2
162 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: v3rho3
163 : END TYPE libxc_worker_type
164 :
165 : ! **************************************************************************************************
166 : !> \brief All threads' workers for one distinct functional setup, together with
167 : !> the key that identifies that setup.
168 : ! **************************************************************************************************
169 : TYPE libxc_worker_set_type
170 : INTEGER :: func_id = -1
171 : INTEGER :: nspin = -1
172 : INTEGER :: family = -1
173 : LOGICAL :: has_laplace = .FALSE.
174 : LOGICAL :: no_exc = .FALSE.
175 : LOGICAL :: has_params = .FALSE.
176 : REAL(KIND=dp) :: epsilon_rho = -1.0_dp
177 : REAL(KIND=dp) :: epsilon_tau = -1.0_dp
178 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: params
179 : CHARACTER(LEN=128), DIMENSION(:), ALLOCATABLE :: param_names
180 : TYPE(libxc_worker_type), DIMENSION(:), ALLOCATABLE :: worker
181 : END TYPE libxc_worker_set_type
182 :
183 : ! Workers are acquired from a serial region only (see libxc_get_workers), so
184 : ! the cache itself needs no locking.
185 : TYPE(libxc_worker_set_type), DIMENSION(:), ALLOCATABLE, TARGET, PRIVATE, SAVE :: libxc_cache
186 :
187 : ! Stands in for the derivative arrays a given evaluation does not produce.
188 : ! Those arguments still have to be associated with something, but they are
189 : ! never read or written; pointing them at the density instead would make an
190 : ! output argument alias a read-only input. Grown on demand and kept for the
191 : ! run rather than allocated per evaluation, and like the cache above it is
192 : ! only ever touched from a serial region.
193 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET, PRIVATE, SAVE :: libxc_unused
194 : #endif
195 :
196 : CONTAINS
197 :
198 : ! **************************************************************************************************
199 : !> \brief This function checks whether a functional name belongs to LibXC
200 : !> \param libxc_params (possible) LibXC input section
201 : !> \return exists whether the functional exists in LibXC
202 : ! **************************************************************************************************
203 4266 : FUNCTION libxc_check_existence_in_libxc(libxc_params) RESULT(exists)
204 :
205 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
206 : LOGICAL :: exists
207 :
208 : #if defined (__LIBXC)
209 :
210 4266 : exists = xc_libxc_check_functional(libxc_params%section%name)
211 : #else
212 : MARK_USED(libxc_params)
213 : exists = .FALSE.
214 : #endif
215 :
216 4266 : END FUNCTION libxc_check_existence_in_libxc
217 :
218 : ! **************************************************************************************************
219 : !> \brief This function returns the maximum length of the reference string for a given LibXC functional
220 : !> \param libxc_params LibXC input section
221 : !> \param lsd spin polarized calculation
222 : !> \return maximum length of the string
223 : ! **************************************************************************************************
224 98 : FUNCTION libxc_get_reference_length(libxc_params, lsd) RESULT(length)
225 :
226 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
227 : LOGICAL, INTENT(IN) :: lsd
228 : INTEGER :: length
229 :
230 : #if defined (__LIBXC)
231 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_get_reference_length'
232 :
233 : CHARACTER(LEN=default_string_length) :: func_name
234 : INTEGER :: func_id, handle
235 : TYPE(xc_f03_func_t) :: xc_func
236 : TYPE(xc_f03_func_info_t) :: xc_info
237 :
238 98 : CALL timeset(routineN, handle)
239 :
240 98 : func_name = libxc_params%section%name
241 :
242 98 : func_id = xc_libxc_wrap_functional_get_number(func_name)
243 98 : IF (lsd) THEN
244 46 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
245 : ELSE
246 52 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
247 : END IF
248 98 : xc_info = xc_f03_func_get_info(xc_func)
249 :
250 98 : length = xc_libxc_get_reference_length(xc_info)
251 :
252 98 : CALL xc_f03_func_end(xc_func)
253 :
254 98 : CALL timestop(handle)
255 : #else
256 : MARK_USED(libxc_params)
257 : MARK_USED(lsd)
258 : length = 0
259 : CPABORT("In order to use LibXC you have to download and install it!")
260 : #endif
261 :
262 98 : END FUNCTION libxc_get_reference_length
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param section ...
267 : ! **************************************************************************************************
268 108595 : SUBROUTINE libxc_add_sections(section)
269 :
270 : TYPE(section_type), POINTER, INTENT(IN) :: section
271 :
272 : #if defined (__LIBXC)
273 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_add_sections'
274 :
275 : TYPE(section_type), POINTER :: subsection
276 : TYPE(keyword_type), POINTER :: keyword
277 : INTEGER :: handle, no_func, len_name, ii, func_id, n_param, iparam
278 : REAL(KIND=C_DOUBLE) :: default_val
279 : CHARACTER(LEN=128) :: func_name, param_name, param_descr, description
280 : CHARACTER(LEN=2*default_string_length) :: warning
281 108595 : INTEGER(KIND=C_INT), DIMENSION(:), ALLOCATABLE :: func_ids
282 : TYPE(xc_f03_func_t) :: xc_func
283 : TYPE(xc_f03_func_info_t) :: xc_info
284 :
285 108595 : CALL timeset(routineN, handle)
286 :
287 108595 : CPASSERT(ASSOCIATED(section))
288 108595 : NULLIFY (subsection, keyword)
289 :
290 108595 : no_func = xc_f03_number_of_functionals()
291 108595 : len_name = xc_f03_maximum_name_length()
292 :
293 325785 : ALLOCATE (func_ids(no_func))
294 :
295 108595 : CALL xc_f03_available_functional_numbers(func_ids)
296 :
297 77428235 : DO ii = 1, no_func
298 :
299 77319640 : func_id = func_ids(ii)
300 77319640 : IF (ii > 1) THEN
301 77211045 : IF (func_id == func_ids(ii - 1)) CYCLE
302 : END IF
303 76233690 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
304 76233690 : xc_info = xc_f03_func_get_info(xc_func)
305 :
306 76233690 : func_name = xc_f03_functional_get_name(func_id)
307 76233690 : description = xc_f03_func_info_get_name(xc_info)
308 76233690 : n_param = xc_f03_func_info_get_n_ext_params(xc_info)
309 :
310 76233690 : NULLIFY (subsection)
311 : CALL section_create(subsection, __LOCATION__, name=TRIM(func_name), description=TRIM(description), &
312 76233690 : n_keywords=2 + n_param, n_subsections=0, repeats=.FALSE.)
313 :
314 76233690 : IF (description(1:1) == "_") THEN
315 : warning = " This parameter is an internal parameter of the functional. Changing this "// &
316 0 : "parameter effectively changes the functional."
317 : ELSE
318 76233690 : warning = " "
319 : END IF
320 :
321 76233690 : NULLIFY (keyword)
322 : CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
323 : description="Activates the functional."//TRIM(warning), &
324 76233690 : lone_keyword_l_val=.TRUE., default_l_val=.FALSE.)
325 76233690 : CALL section_add_keyword(subsection, keyword)
326 76233690 : CALL keyword_release(keyword)
327 :
328 : CALL keyword_create(keyword, __LOCATION__, name="SCALE", description="Scales this functional", &
329 76233690 : default_r_val=1.0_dp)
330 76233690 : CALL section_add_keyword(subsection, keyword)
331 76233690 : CALL keyword_release(keyword)
332 :
333 464460815 : DO iparam = 1, n_param
334 388227125 : param_name = xc_f03_func_info_get_ext_params_name(xc_info, iparam - 1)
335 388227125 : param_descr = xc_f03_func_info_get_ext_params_description(xc_info, iparam - 1)
336 388227125 : default_val = xc_f03_func_info_get_ext_params_default_value(xc_info, iparam - 1)
337 388227125 : NULLIFY (keyword)
338 : CALL keyword_create(keyword, __LOCATION__, name=TRIM(param_name), &
339 388227125 : description=TRIM(param_descr), default_r_val=default_val)
340 388227125 : CALL section_add_keyword(subsection, keyword)
341 464460815 : CALL keyword_release(keyword)
342 : END DO
343 :
344 76233690 : CALL section_add_subsection(section, subsection)
345 76233690 : CALL section_release(subsection)
346 :
347 77428235 : CALL xc_f03_func_end(xc_func)
348 :
349 : END DO
350 :
351 108595 : DEALLOCATE (func_ids)
352 :
353 108595 : CALL timestop(handle)
354 : #else
355 : MARK_USED(section)
356 :
357 : #endif
358 :
359 108595 : END SUBROUTINE libxc_add_sections
360 :
361 : ! **************************************************************************************************
362 : !> \brief info about the functional from libxc
363 : !> \param libxc_params input parameter (functional name, scaling and parameters)
364 : !> \param reference string with the reference of the actual functional
365 : !> \param shortform string with the shortform of the functional name
366 : !> \param needs the components needed by this functional are set to
367 : !> true (does not set the unneeded components to false)
368 : !> \param max_deriv maximum implemented derivative of the xc functional
369 : !> \param print_warn whether to print warning about development status of a functional
370 : !> \param func_name_override optional LibXC functional name overriding the section name
371 : !> \author F. Tran
372 : ! **************************************************************************************************
373 13390 : SUBROUTINE libxc_spin_unpolarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
374 : func_name_override)
375 :
376 : TYPE(section_vals_type), POINTER :: libxc_params
377 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
378 : TYPE(xc_rho_cflags_type), &
379 : INTENT(inout), OPTIONAL :: needs
380 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
381 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
382 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
383 :
384 : #if defined (__LIBXC)
385 : CHARACTER(LEN=128) :: s1, s2
386 : CHARACTER(LEN=default_string_length) :: func_name
387 : INTEGER :: func_id
388 : REAL(KIND=dp) :: func_scale
389 : TYPE(xc_f03_func_t) :: xc_func
390 : TYPE(xc_f03_func_info_t) :: xc_info
391 :
392 26692 : IF (PRESENT(func_name_override)) THEN
393 88 : func_name = func_name_override
394 88 : func_scale = 1.0_dp
395 : ELSE
396 13302 : func_name = libxc_params%section%name
397 13302 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
398 : END IF
399 :
400 13390 : CALL cite_reference(Marques2012)
401 13390 : CALL cite_reference(Lehtola2018)
402 :
403 13390 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
404 :
405 13390 : func_id = xc_libxc_wrap_functional_get_number(func_name)
406 13390 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
407 13390 : xc_info = xc_f03_func_get_info(xc_func)
408 :
409 13390 : s1 = xc_f03_func_info_get_name(xc_info)
410 9121 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
411 9121 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
412 2491 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
413 1448 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
414 330 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
415 : CASE default
416 13390 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
417 : END SELECT
418 13390 : IF (PRESENT(shortform)) THEN
419 52 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
420 : END IF
421 13390 : IF (PRESENT(reference)) THEN
422 52 : CALL xc_libxc_wrap_info_refs(xc_info, XC_UNPOLARIZED, func_scale, reference)
423 : END IF
424 13390 : IF (PRESENT(needs)) THEN
425 5880 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
426 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
427 5880 : needs%rho = .TRUE.
428 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
429 4504 : needs%rho = .TRUE.
430 4504 : needs%norm_drho = .TRUE.
431 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
432 2946 : needs%rho = .TRUE.
433 2946 : needs%norm_drho = .TRUE.
434 2946 : needs%tau = .TRUE.
435 2946 : needs%laplace_rho = xc_libxc_wrap_info_needs_laplace(xc_info)
436 : CASE default
437 13330 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
438 : END SELECT
439 : END IF
440 13390 : IF (PRESENT(max_deriv)) THEN
441 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
442 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
443 0 : max_deriv = 3
444 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
445 88 : max_deriv = 2
446 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
447 0 : max_deriv = 2
448 : CASE default
449 88 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
450 : END SELECT
451 : END IF
452 13390 : IF (PRESENT(print_warn)) THEN
453 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
454 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
455 : END IF
456 : END IF
457 :
458 13390 : CALL xc_f03_func_end(xc_func)
459 : #else
460 : MARK_USED(libxc_params)
461 : MARK_USED(reference)
462 : MARK_USED(shortform)
463 : MARK_USED(needs)
464 : MARK_USED(max_deriv)
465 : MARK_USED(print_warn)
466 : MARK_USED(func_name_override)
467 :
468 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
469 : "for a functional of the LibXC library, "// &
470 : "you have to download and install the library!")
471 : #endif
472 :
473 13390 : END SUBROUTINE libxc_spin_unpolarized_info
474 :
475 : ! **************************************************************************************************
476 : !> \brief info about the functional from libxc
477 : !> \param libxc_params input parameter (functional name, scaling and parameters)
478 : !> \param reference string with the reference of the actual functional
479 : !> \param shortform string with the shortform of the functional name
480 : !> \param needs the components needed by this functional are set to
481 : !> true (does not set the unneeded components to false)
482 : !> \param max_deriv maximum implemented derivative of the xc functional
483 : !> \param print_warn whether to print warning about development status of a functional
484 : !> \param func_name_override optional LibXC functional name overriding the section name
485 : !> \author F. Tran
486 : ! **************************************************************************************************
487 2476 : SUBROUTINE libxc_spin_polarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
488 : func_name_override)
489 :
490 : TYPE(section_vals_type), POINTER :: libxc_params
491 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
492 : TYPE(xc_rho_cflags_type), &
493 : INTENT(inout), OPTIONAL :: needs
494 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
495 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
496 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
497 :
498 : #if defined (__LIBXC)
499 : CHARACTER(LEN=128) :: s1, s2
500 : CHARACTER(LEN=default_string_length) :: func_name
501 : INTEGER :: func_id
502 : REAL(KIND=dp) :: func_scale
503 : TYPE(xc_f03_func_t) :: xc_func
504 : TYPE(xc_f03_func_info_t) :: xc_info
505 :
506 4944 : IF (PRESENT(func_name_override)) THEN
507 8 : func_name = func_name_override
508 8 : func_scale = 1.0_dp
509 : ELSE
510 2468 : func_name = libxc_params%section%name
511 2468 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
512 : END IF
513 :
514 2476 : CALL cite_reference(Marques2012)
515 2476 : CALL cite_reference(Lehtola2018)
516 :
517 2476 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
518 :
519 2476 : func_id = xc_libxc_wrap_functional_get_number(func_name)
520 2476 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
521 2476 : xc_info = xc_f03_func_get_info(xc_func)
522 :
523 2476 : s1 = xc_f03_func_info_get_name(xc_info)
524 1158 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
525 1158 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
526 1072 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
527 246 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
528 0 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
529 : CASE default
530 2476 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
531 : END SELECT
532 2476 : IF (PRESENT(shortform)) THEN
533 46 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
534 : END IF
535 2476 : IF (PRESENT(reference)) THEN
536 46 : CALL xc_libxc_wrap_info_refs(xc_info, XC_POLARIZED, func_scale, reference)
537 : END IF
538 2476 : IF (PRESENT(needs)) THEN
539 1016 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
540 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
541 1016 : needs%rho_spin = .TRUE.
542 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
543 262 : needs%rho_spin = .TRUE.
544 262 : needs%norm_drho = .TRUE.
545 262 : needs%norm_drho_spin = .TRUE.
546 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
547 1152 : needs%rho_spin = .TRUE.
548 1152 : needs%norm_drho = .TRUE.
549 1152 : needs%norm_drho_spin = .TRUE.
550 1152 : needs%tau_spin = .TRUE.
551 1152 : needs%laplace_rho_spin = xc_libxc_wrap_info_needs_laplace(xc_info)
552 : CASE default
553 2430 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
554 : END SELECT
555 : END IF
556 2476 : IF (PRESENT(max_deriv)) THEN
557 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
558 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
559 0 : max_deriv = 3
560 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
561 8 : max_deriv = 2
562 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
563 0 : max_deriv = 2
564 : CASE default
565 8 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
566 : END SELECT
567 : END IF
568 2476 : IF (PRESENT(print_warn)) THEN
569 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
570 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
571 : END IF
572 : END IF
573 :
574 2476 : CALL xc_f03_func_end(xc_func)
575 : #else
576 : MARK_USED(libxc_params)
577 : MARK_USED(reference)
578 : MARK_USED(shortform)
579 : MARK_USED(needs)
580 : MARK_USED(max_deriv)
581 : MARK_USED(print_warn)
582 : MARK_USED(func_name_override)
583 :
584 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are "// &
585 : "asking for a functional of the LibXC library, "// &
586 : "you have to download and install the library!")
587 : #endif
588 :
589 2476 : END SUBROUTINE libxc_spin_polarized_info
590 :
591 : ! **************************************************************************************************
592 : !> \brief info about the LibXC version
593 : !> \param version version string reported by the library at run time
594 : !> \param compiled_version version of the LibXC headers CP2K was compiled against
595 : !> \author A. Gloess (agloess)
596 : ! **************************************************************************************************
597 63 : SUBROUTINE libxc_version_info(version, compiled_version)
598 : CHARACTER(LEN=*), INTENT(OUT) :: version ! the string that is output
599 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: compiled_version
600 :
601 : #if defined (__LIBXC)
602 63 : CALL xc_libxc_wrap_version(version, compiled_version)
603 : #else
604 : version = "none"
605 : IF (PRESENT(compiled_version)) compiled_version = "none"
606 : CPABORT("In order to use libxc you need to download and install it")
607 : #endif
608 :
609 63 : END SUBROUTINE libxc_version_info
610 :
611 : ! **************************************************************************************************
612 : !> \brief Returns the citation LibXC asks for the library itself.
613 : !> \param reference bibliographic reference
614 : !> \param doi digital object identifier of that reference
615 : !> \author S. Lehtola
616 : ! **************************************************************************************************
617 63 : SUBROUTINE libxc_library_reference(reference, doi)
618 : CHARACTER(LEN=*), INTENT(OUT) :: reference, doi
619 :
620 : #if defined (__LIBXC)
621 63 : CALL xc_libxc_wrap_library_reference(reference, doi)
622 : #else
623 : reference = "none"
624 : doi = "none"
625 : CPABORT("In order to use libxc you need to download and install it")
626 : #endif
627 :
628 63 : END SUBROUTINE libxc_library_reference
629 :
630 : ! **************************************************************************************************
631 : !> \brief Reads the external parameters of a LibXC functional from its input section.
632 : !> \param libxc_params LibXC input section
633 : !> \param param_names names of the external parameters, as reported by LibXC
634 : !> \param params values read from the input
635 : !> \author S. Lehtola
636 : ! **************************************************************************************************
637 : #if defined (__LIBXC)
638 12440 : SUBROUTINE libxc_read_params(libxc_params, param_names, params)
639 : TYPE(section_vals_type), INTENT(IN), POINTER :: libxc_params
640 : CHARACTER(LEN=128), DIMENSION(:), INTENT(IN) :: param_names
641 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: params
642 :
643 : INTEGER :: i
644 :
645 76952 : DO i = 1, SIZE(param_names)
646 76952 : CALL section_vals_val_get(libxc_params, TRIM(param_names(i)), r_val=params(i))
647 : END DO
648 :
649 12440 : END SUBROUTINE libxc_read_params
650 :
651 : ! **************************************************************************************************
652 : !> \brief Allocates the staging buffers of a single worker.
653 : !> \param worker the worker
654 : !> \param family LibXC functional family
655 : !> \param nspin XC_UNPOLARIZED or XC_POLARIZED
656 : !> \author S. Lehtola
657 : !> \note The buffers follow LibXC's own layout, i.e. the leading dimension is the
658 : !> number of components of the quantity and the points run along the second
659 : !> dimension. They are allocated once per worker and reused for every batch.
660 : ! **************************************************************************************************
661 470 : SUBROUTINE libxc_worker_alloc(worker, family, nspin)
662 : TYPE(libxc_worker_type), INTENT(INOUT) :: worker
663 : INTEGER, INTENT(IN) :: family, nspin
664 :
665 : INTEGER :: nb, nl, np, ns
666 :
667 470 : nb = libxc_block_size
668 : ! number of components per point, cf. the LibXC manual
669 470 : IF (nspin == XC_POLARIZED) THEN
670 : np = 2 ! rho, vrho, lapl, tau, vlapl, vtau
671 : ns = 3 ! sigma, vsigma, v2rho2, v2lapl2, v2tau2
672 : nl = 4 ! v2rholapl, v2rhotau, v2lapltau, v3rho3
673 : ELSE
674 328 : np = 1
675 328 : ns = 1
676 328 : nl = 1
677 : END IF
678 :
679 940 : ALLOCATE (worker%rho(np, nb))
680 470 : ALLOCATE (worker%exc(nb))
681 940 : ALLOCATE (worker%vrho(np, nb))
682 940 : ALLOCATE (worker%v2rho2(ns, nb))
683 940 : ALLOCATE (worker%v3rho3(nl, nb))
684 :
685 770 : SELECT CASE (family)
686 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
687 300 : ALLOCATE (worker%sigma(ns, nb))
688 300 : ALLOCATE (worker%vsigma(ns, nb))
689 300 : IF (nspin == XC_POLARIZED) THEN
690 74 : ALLOCATE (worker%nd(nb), worker%nda(nb), worker%ndb(nb))
691 : END IF
692 826 : ALLOCATE (worker%v2rhosigma(MERGE(6, 1, nspin == XC_POLARIZED), nb))
693 1070 : ALLOCATE (worker%v2sigma2(MERGE(6, 1, nspin == XC_POLARIZED), nb))
694 : END SELECT
695 :
696 126 : SELECT CASE (family)
697 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
698 252 : ALLOCATE (worker%lapl(np, nb))
699 252 : ALLOCATE (worker%tau(np, nb))
700 252 : ALLOCATE (worker%vlapl(np, nb))
701 252 : ALLOCATE (worker%vtau(np, nb))
702 252 : ALLOCATE (worker%v2rholapl(nl, nb))
703 252 : ALLOCATE (worker%v2rhotau(nl, nb))
704 346 : ALLOCATE (worker%v2sigmalapl(MERGE(6, 1, nspin == XC_POLARIZED), nb))
705 252 : ALLOCATE (worker%v2sigmatau(MERGE(6, 1, nspin == XC_POLARIZED), nb))
706 126 : ALLOCATE (worker%v2lapl2(ns, nb))
707 252 : ALLOCATE (worker%v2lapltau(nl, nb))
708 126 : ALLOCATE (worker%v2tau2(ns, nb))
709 : ! LibXC leaves vlapl untouched for functionals that do not use the
710 : ! Laplacian, and CP2K then ignores it; zero it once so it never holds
711 : ! uninitialized memory.
712 146004 : worker%vlapl = 0.0_dp
713 : END SELECT
714 :
715 470 : END SUBROUTINE libxc_worker_alloc
716 :
717 : ! **************************************************************************************************
718 : !> \brief Makes sure the placeholder for unproduced derivatives holds npoints values.
719 : !> \param npoints number of grid points in this evaluation
720 : !> \author S. Lehtola
721 : ! **************************************************************************************************
722 20450 : SUBROUTINE libxc_reserve_unused(npoints)
723 : INTEGER, INTENT(IN) :: npoints
724 :
725 20450 : IF (ALLOCATED(libxc_unused)) THEN
726 20170 : IF (SIZE(libxc_unused) >= npoints) RETURN
727 2 : DEALLOCATE (libxc_unused)
728 : END IF
729 846 : ALLOCATE (libxc_unused(npoints))
730 :
731 : END SUBROUTINE libxc_reserve_unused
732 :
733 : ! **************************************************************************************************
734 : !> \brief Returns the cache slot holding the per-thread workers for a functional.
735 : !> \param func_name LibXC functional name
736 : !> \param nspin XC_UNPOLARIZED or XC_POLARIZED
737 : !> \param libxc_params LibXC input section, or null when no parameters are to be set
738 : !> \param epsilon_rho density cutoff
739 : !> \param epsilon_tau kinetic energy density cutoff
740 : !> \return index into libxc_cache
741 : !> \author S. Lehtola
742 : !> \note Must be called from a serial region: it is the only writer of the cache,
743 : !> and keeping it serial is what makes the cache lock-free. The workers it
744 : !> returns are then used concurrently, one per thread.
745 : ! **************************************************************************************************
746 20450 : FUNCTION libxc_get_workers(func_name, nspin, libxc_params, epsilon_rho, epsilon_tau) RESULT(idx)
747 : CHARACTER(LEN=*), INTENT(IN) :: func_name
748 : INTEGER, INTENT(IN) :: nspin
749 : TYPE(section_vals_type), INTENT(IN), POINTER :: libxc_params
750 : REAL(KIND=dp), INTENT(IN) :: epsilon_rho, epsilon_tau
751 : INTEGER :: idx
752 :
753 : CHARACTER(LEN=*), PARAMETER :: routineN = 'libxc_get_workers'
754 :
755 : INTEGER :: func_id, handle, i, ithread, n_params, &
756 : nthreads
757 : LOGICAL :: has_params
758 20450 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: params
759 : TYPE(libxc_worker_set_type), ALLOCATABLE, &
760 20450 : DIMENSION(:) :: tmp_cache
761 :
762 20450 : CALL timeset(routineN, handle)
763 :
764 20450 : func_id = xc_libxc_wrap_functional_get_number(func_name)
765 20450 : has_params = ASSOCIATED(libxc_params)
766 :
767 20450 : nthreads = 1
768 20450 : !$ nthreads = omp_get_max_threads()
769 :
770 20450 : IF (.NOT. ALLOCATED(libxc_cache)) ALLOCATE (libxc_cache(0))
771 :
772 20450 : idx = 0
773 28396 : DO i = 1, SIZE(libxc_cache)
774 : IF (libxc_cache(i)%func_id == func_id .AND. &
775 27926 : libxc_cache(i)%nspin == nspin .AND. &
776 470 : (libxc_cache(i)%has_params .EQV. has_params)) THEN
777 : IF (libxc_cache(i)%epsilon_rho == epsilon_rho .AND. &
778 19980 : libxc_cache(i)%epsilon_tau == epsilon_tau .AND. &
779 : SIZE(libxc_cache(i)%worker) >= nthreads) THEN
780 : idx = i
781 : EXIT
782 : END IF
783 : END IF
784 : END DO
785 :
786 20450 : IF (idx == 0) THEN
787 : ! Grow the cache by one slot and build the workers for it.
788 470 : CALL MOVE_ALLOC(libxc_cache, tmp_cache)
789 2198 : ALLOCATE (libxc_cache(SIZE(tmp_cache) + 1))
790 788 : DO i = 1, SIZE(tmp_cache)
791 788 : CALL libxc_move_worker_set(tmp_cache(i), libxc_cache(i))
792 : END DO
793 1258 : DEALLOCATE (tmp_cache)
794 470 : idx = SIZE(libxc_cache)
795 :
796 : ASSOCIATE (ws => libxc_cache(idx))
797 470 : ws%func_id = func_id
798 470 : ws%nspin = nspin
799 470 : ws%has_params = has_params
800 470 : ws%epsilon_rho = epsilon_rho
801 470 : ws%epsilon_tau = epsilon_tau
802 1880 : ALLOCATE (ws%worker(nthreads))
803 :
804 940 : DO ithread = 1, nthreads
805 470 : CALL xc_f03_func_init(ws%worker(ithread)%func, func_id, nspin)
806 470 : ws%worker(ithread)%info = xc_f03_func_get_info(ws%worker(ithread)%func)
807 940 : ws%worker(ithread)%is_init = .TRUE.
808 : END DO
809 :
810 470 : ws%family = xc_f03_func_info_get_family(ws%worker(1)%info)
811 470 : ws%has_laplace = xc_libxc_wrap_info_needs_laplace(ws%worker(1)%info)
812 470 : ws%no_exc = xc_libxc_wrap_info_no_exc(ws%worker(1)%info)
813 :
814 : ! Remember the parameter names so that later evaluations can re-read
815 : ! the input without another functional object to ask for them.
816 470 : n_params = 0
817 470 : IF (has_params) n_params = xc_f03_func_info_get_n_ext_params(ws%worker(1)%info)
818 1312 : ALLOCATE (ws%params(n_params))
819 1312 : ALLOCATE (ws%param_names(n_params))
820 2366 : DO i = 1, n_params
821 2366 : ws%param_names(i) = xc_f03_func_info_get_ext_params_name(ws%worker(1)%info, i - 1)
822 : END DO
823 470 : IF (n_params > 0) THEN
824 372 : CALL libxc_read_params(libxc_params, ws%param_names, ws%params)
825 744 : DO ithread = 1, nthreads
826 744 : CALL xc_f03_func_set_ext_params(ws%worker(ithread)%func, ws%params)
827 : END DO
828 : END IF
829 :
830 1410 : DO ithread = 1, nthreads
831 : CALL xc_libxc_wrap_set_thresholds(ws%worker(ithread)%func, ws%worker(ithread)%info, &
832 470 : epsilon_rho, epsilon_tau)
833 940 : CALL libxc_worker_alloc(ws%worker(ithread), ws%family, nspin)
834 : END DO
835 : END ASSOCIATE
836 : ELSE
837 : ! Cached slot: the input section may in principle carry different
838 : ! parameter values than the ones the workers were last set up with, so
839 : ! re-read them and push them through only when they actually changed.
840 : ASSOCIATE (ws => libxc_cache(idx))
841 19980 : n_params = SIZE(ws%params)
842 19980 : IF (n_params > 0) THEN
843 36204 : ALLOCATE (params(n_params))
844 12068 : CALL libxc_read_params(libxc_params, ws%param_names, params)
845 74684 : IF (ANY(params /= ws%params)) THEN
846 0 : ws%params(:) = params(:)
847 0 : DO ithread = 1, SIZE(ws%worker)
848 0 : CALL xc_f03_func_set_ext_params(ws%worker(ithread)%func, ws%params)
849 : ! re-assert the cutoffs: setting parameters runs the
850 : ! functional's own callback, which is free to rebuild the
851 : ! component functionals a mixed functional is made of
852 : CALL xc_libxc_wrap_set_thresholds(ws%worker(ithread)%func, &
853 : ws%worker(ithread)%info, &
854 0 : ws%epsilon_rho, ws%epsilon_tau)
855 : END DO
856 : END IF
857 12068 : DEALLOCATE (params)
858 : END IF
859 : END ASSOCIATE
860 : END IF
861 :
862 20450 : CALL timestop(handle)
863 :
864 20450 : END FUNCTION libxc_get_workers
865 :
866 : ! **************************************************************************************************
867 : !> \brief Moves a worker set to a new cache slot without touching LibXC.
868 : !> \param from source slot, left empty
869 : !> \param to destination slot
870 : !> \author S. Lehtola
871 : !> \note Used when the cache array grows. The LibXC objects are opaque handles, so
872 : !> they survive being moved; only the Fortran allocatables need transferring.
873 : ! **************************************************************************************************
874 318 : SUBROUTINE libxc_move_worker_set(from, to)
875 : TYPE(libxc_worker_set_type), INTENT(INOUT) :: from, to
876 :
877 318 : to%func_id = from%func_id
878 318 : to%nspin = from%nspin
879 318 : to%family = from%family
880 318 : to%has_laplace = from%has_laplace
881 318 : to%no_exc = from%no_exc
882 318 : to%has_params = from%has_params
883 318 : to%epsilon_rho = from%epsilon_rho
884 318 : to%epsilon_tau = from%epsilon_tau
885 318 : IF (ALLOCATED(from%params)) CALL MOVE_ALLOC(from%params, to%params)
886 318 : IF (ALLOCATED(from%param_names)) CALL MOVE_ALLOC(from%param_names, to%param_names)
887 318 : IF (ALLOCATED(from%worker)) CALL MOVE_ALLOC(from%worker, to%worker)
888 :
889 318 : END SUBROUTINE libxc_move_worker_set
890 : #endif
891 :
892 : ! **************************************************************************************************
893 : !> \brief Destroys all cached LibXC functional objects.
894 : !> \author S. Lehtola
895 : !> \note Call once at the end of a run. The workers are deliberately kept alive
896 : !> across SCF steps and across calculations, so nothing else releases them.
897 : ! **************************************************************************************************
898 11241 : SUBROUTINE libxc_release_workers()
899 :
900 : #if defined (__LIBXC)
901 : INTEGER :: i, ithread
902 :
903 11241 : IF (.NOT. ALLOCATED(libxc_cache)) RETURN
904 :
905 750 : DO i = 1, SIZE(libxc_cache)
906 470 : IF (.NOT. ALLOCATED(libxc_cache(i)%worker)) CYCLE
907 1220 : DO ithread = 1, SIZE(libxc_cache(i)%worker)
908 940 : IF (libxc_cache(i)%worker(ithread)%is_init) THEN
909 470 : CALL xc_f03_func_end(libxc_cache(i)%worker(ithread)%func)
910 470 : libxc_cache(i)%worker(ithread)%is_init = .FALSE.
911 : END IF
912 : END DO
913 : END DO
914 1220 : DEALLOCATE (libxc_cache)
915 280 : IF (ALLOCATED(libxc_unused)) DEALLOCATE (libxc_unused)
916 : #endif
917 :
918 : END SUBROUTINE libxc_release_workers
919 :
920 : ! **************************************************************************************************
921 : !> \brief evaluates the functional from libxc
922 : !> \param rho_set the density where you want to evaluate the functional
923 : !> \param deriv_set place where to store the functional derivatives (they are
924 : !> added to the derivatives)
925 : !> \param grad_deriv degree of the derivative that should be evaluated;
926 : !> all derivatives up to the given degree are evaluated, in a single
927 : !> LibXC call per block of grid points
928 : !> \param libxc_params input parameter (functional name, scaling and parameters)
929 : !> \param func_name_override optional LibXC functional name overriding the section name
930 : !> \author F. Tran
931 : ! **************************************************************************************************
932 17700 : SUBROUTINE libxc_spin_unpolarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
933 :
934 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
935 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
936 : INTEGER, INTENT(in) :: grad_deriv
937 : TYPE(section_vals_type), POINTER :: libxc_params
938 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
939 :
940 : #if defined (__LIBXC)
941 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_unpolarized_eval'
942 :
943 : CHARACTER(LEN=default_string_length) :: func_name
944 : INTEGER :: handle, iw, npoints
945 : INTEGER, DIMENSION(2, 3) :: bo
946 : LOGICAL :: has_laplace
947 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
948 : TYPE(libxc_worker_set_type), POINTER :: workers
949 : TYPE(section_vals_type), POINTER :: no_params
950 17700 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rho, &
951 17700 : e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_ndrho, &
952 17700 : e_ndrho_laplace_rho, e_ndrho_ndrho, e_ndrho_rho, e_ndrho_tau, e_rho, &
953 17700 : e_rho_laplace_rho, e_rho_rho, e_rho_rho_rho, e_rho_tau, e_tau, &
954 17700 : e_tau_tau, laplace_rho, norm_drho, rho, tau
955 : TYPE(xc_derivative_type), POINTER :: deriv
956 : TYPE(xc_f03_func_info_t) :: xc_info
957 :
958 17700 : CALL timeset(routineN, handle)
959 :
960 : ! Only "everything up to grad_deriv" is supported. Asking for a single
961 : ! derivative order in isolation is rejected here rather than further down,
962 : ! where it would quietly evaluate nothing at all.
963 17700 : IF (grad_deriv < 0) THEN
964 0 : CPABORT("Evaluating a single derivative order is not supported.")
965 : END IF
966 17700 : IF (grad_deriv > 3) THEN
967 0 : CPABORT("derivatives larger than 3 not implemented")
968 : END IF
969 :
970 17700 : has_laplace = .FALSE.
971 17700 : NULLIFY (dummy)
972 17700 : NULLIFY (workers)
973 17700 : NULLIFY (rho, norm_drho, laplace_rho, tau)
974 :
975 17700 : IF (PRESENT(func_name_override)) THEN
976 0 : func_name = func_name_override
977 0 : func_scale = 1.0_dp
978 : ELSE
979 17700 : func_name = libxc_params%section%name
980 17700 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
981 : END IF
982 :
983 17700 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
984 :
985 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
986 : rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, &
987 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
988 17700 : tau=tau, local_bounds=bo)
989 :
990 17700 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
991 :
992 : ! One functional object per thread, kept across calls. The cutoffs are part
993 : ! of the worker setup because LibXC does the screening itself.
994 17700 : IF (PRESENT(func_name_override)) THEN
995 : ! an overriding name carries no input section, so no parameters are set
996 0 : NULLIFY (no_params)
997 0 : iw = libxc_get_workers(func_name, XC_UNPOLARIZED, no_params, epsilon_rho, epsilon_tau)
998 : ELSE
999 17700 : iw = libxc_get_workers(func_name, XC_UNPOLARIZED, libxc_params, epsilon_rho, epsilon_tau)
1000 : END IF
1001 17700 : workers => libxc_cache(iw)
1002 17700 : xc_info = workers%worker(1)%info
1003 17700 : has_laplace = workers%has_laplace
1004 :
1005 : ! see libxc_unused: the arguments the requested order does not produce are
1006 : ! never touched, but they still have to point somewhere
1007 17700 : CALL libxc_reserve_unused(npoints)
1008 17700 : dummy(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) => libxc_unused(1:npoints)
1009 :
1010 : ! due to assumed shape array usage in next routine
1011 17700 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => rho
1012 17700 : IF (.NOT. ASSOCIATED(tau)) tau => rho
1013 :
1014 : ! only some MGGA functionals really need the Laplacian,
1015 : ! all others can work with rho (read-only) as dummy
1016 17700 : IF (.NOT. has_laplace) laplace_rho => rho
1017 :
1018 17700 : e_0 => dummy
1019 17700 : e_rho => dummy
1020 17700 : e_ndrho => dummy
1021 17700 : e_laplace_rho => dummy
1022 17700 : e_tau => dummy
1023 17700 : e_rho_rho => dummy
1024 17700 : e_ndrho_rho => dummy
1025 17700 : e_ndrho_ndrho => dummy
1026 17700 : e_rho_laplace_rho => dummy
1027 17700 : e_rho_tau => dummy
1028 17700 : e_ndrho_laplace_rho => dummy
1029 17700 : e_ndrho_tau => dummy
1030 17700 : e_laplace_rho_laplace_rho => dummy
1031 17700 : e_laplace_rho_tau => dummy
1032 17700 : e_tau_tau => dummy
1033 17700 : e_rho_rho_rho => dummy
1034 :
1035 17700 : IF (grad_deriv >= 0) THEN
1036 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
1037 17700 : allocate_deriv=.TRUE.)
1038 17700 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1039 : END IF
1040 17700 : IF (grad_deriv >= 1) THEN
1041 10210 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1042 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1043 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1044 10210 : allocate_deriv=.TRUE.)
1045 10210 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1046 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1047 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1048 5126 : allocate_deriv=.TRUE.)
1049 5126 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1050 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1051 5126 : allocate_deriv=.TRUE.)
1052 5126 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1053 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1054 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
1055 2110 : allocate_deriv=.TRUE.)
1056 2110 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
1057 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1058 2110 : allocate_deriv=.TRUE.)
1059 2110 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1060 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
1061 2110 : allocate_deriv=.TRUE.)
1062 2110 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
1063 2110 : IF (has_laplace) THEN
1064 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
1065 522 : allocate_deriv=.TRUE.)
1066 522 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
1067 : END IF
1068 : CASE default
1069 17446 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1070 : END SELECT
1071 : END IF
1072 17700 : IF (grad_deriv >= 2) THEN
1073 1488 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1074 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1075 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1076 1488 : allocate_deriv=.TRUE.)
1077 1488 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1078 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1079 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1080 606 : allocate_deriv=.TRUE.)
1081 606 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1082 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
1083 606 : allocate_deriv=.TRUE.)
1084 606 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
1085 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1086 606 : allocate_deriv=.TRUE.)
1087 606 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1088 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1089 : ! not implemented ...
1090 :
1091 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
1092 148 : allocate_deriv=.TRUE.)
1093 148 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
1094 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
1095 148 : allocate_deriv=.TRUE.)
1096 148 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
1097 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1098 148 : allocate_deriv=.TRUE.)
1099 148 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1100 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_tau], &
1101 148 : allocate_deriv=.TRUE.)
1102 148 : CALL xc_derivative_get(deriv, deriv_data=e_rho_tau)
1103 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau], &
1104 148 : allocate_deriv=.TRUE.)
1105 148 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau)
1106 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau, deriv_tau], &
1107 148 : allocate_deriv=.TRUE.)
1108 148 : CALL xc_derivative_get(deriv, deriv_data=e_tau_tau)
1109 148 : IF (has_laplace) THEN
1110 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_laplace_rho], &
1111 62 : allocate_deriv=.TRUE.)
1112 62 : CALL xc_derivative_get(deriv, deriv_data=e_rho_laplace_rho)
1113 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rho], &
1114 62 : allocate_deriv=.TRUE.)
1115 62 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rho)
1116 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_laplace_rho], &
1117 62 : allocate_deriv=.TRUE.)
1118 62 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_laplace_rho)
1119 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_tau], &
1120 62 : allocate_deriv=.TRUE.)
1121 62 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_tau)
1122 : END IF
1123 : CASE default
1124 2242 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1125 : END SELECT
1126 : END IF
1127 17700 : IF (grad_deriv >= 3) THEN
1128 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1129 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1130 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
1131 0 : allocate_deriv=.TRUE.)
1132 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
1133 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1134 0 : CPABORT("derivatives larger than 2 not implemented")
1135 : CASE default
1136 0 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1137 : END SELECT
1138 : END IF
1139 :
1140 : !$OMP PARALLEL DEFAULT(NONE), &
1141 : !$OMP SHARED(rho,norm_drho,laplace_rho,tau,e_0,e_rho,e_ndrho,e_laplace_rho),&
1142 : !$OMP SHARED(e_tau,e_rho_rho,e_ndrho_rho,e_ndrho_ndrho,e_rho_laplace_rho),&
1143 : !$OMP SHARED(e_rho_tau,e_ndrho_laplace_rho,e_ndrho_tau,e_laplace_rho_laplace_rho),&
1144 : !$OMP SHARED(e_laplace_rho_tau,e_tau_tau,e_rho_rho_rho),&
1145 : !$OMP SHARED(grad_deriv,npoints),&
1146 17700 : !$OMP SHARED(func_name,func_scale,workers)
1147 :
1148 : CALL libxc_spin_unpolarized_calc(rho=rho, norm_drho=norm_drho, &
1149 : laplace_rho=laplace_rho, tau=tau, &
1150 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_laplace_rho=e_laplace_rho, &
1151 : e_tau=e_tau, e_rho_rho=e_rho_rho, e_ndrho_rho=e_ndrho_rho, &
1152 : e_ndrho_ndrho=e_ndrho_ndrho, e_rho_laplace_rho=e_rho_laplace_rho, &
1153 : e_rho_tau=e_rho_tau, e_ndrho_laplace_rho=e_ndrho_laplace_rho, &
1154 : e_ndrho_tau=e_ndrho_tau, e_laplace_rho_laplace_rho=e_laplace_rho_laplace_rho, &
1155 : e_laplace_rho_tau=e_laplace_rho_tau, e_tau_tau=e_tau_tau, &
1156 : e_rho_rho_rho=e_rho_rho_rho, &
1157 : grad_deriv=grad_deriv, npoints=npoints, &
1158 : func_name=func_name, sc=func_scale, workers=workers)
1159 :
1160 : !$OMP END PARALLEL
1161 :
1162 17700 : NULLIFY (dummy)
1163 17700 : NULLIFY (workers)
1164 :
1165 17700 : CALL timestop(handle)
1166 : #else
1167 : MARK_USED(rho_set)
1168 : MARK_USED(deriv_set)
1169 : MARK_USED(grad_deriv)
1170 : MARK_USED(libxc_params)
1171 : MARK_USED(func_name_override)
1172 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
1173 : "for a functional of the LibXC library, "// &
1174 : "you have to download and install the library!")
1175 : #endif
1176 17700 : END SUBROUTINE libxc_spin_unpolarized_eval
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief evaluates the functional from libxc
1180 : !> \param rho_set the density where you want to evaluate the functional
1181 : !> \param deriv_set place where to store the functional derivatives (they are
1182 : !> added to the derivatives)
1183 : !> \param grad_deriv degree of the derivative that should be evaluated;
1184 : !> all derivatives up to the given degree are evaluated, in a single
1185 : !> LibXC call per block of grid points
1186 : !> \param libxc_params input parameter (functional name, scaling and parameters)
1187 : !> \param func_name_override optional LibXC functional name overriding the section name
1188 : !> \author F. Tran
1189 : ! **************************************************************************************************
1190 2750 : SUBROUTINE libxc_spin_polarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
1191 :
1192 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1193 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
1194 : INTEGER, INTENT(in) :: grad_deriv
1195 : TYPE(section_vals_type), POINTER :: libxc_params
1196 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
1197 :
1198 : #if defined (__LIBXC)
1199 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_polarized_eval'
1200 :
1201 : CHARACTER(LEN=default_string_length) :: func_name
1202 : INTEGER :: handle, iw, npoints
1203 : INTEGER, DIMENSION(2, 3) :: bo
1204 : LOGICAL :: has_laplace
1205 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
1206 : TYPE(libxc_worker_set_type), POINTER :: workers
1207 : TYPE(section_vals_type), POINTER :: no_params
1208 2750 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
1209 2750 : e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
1210 2750 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, e_laplace_rhob, &
1211 2750 : e_laplace_rhob_laplace_rhob, e_laplace_rhob_tau_a, &
1212 2750 : e_laplace_rhob_tau_b, e_ndrho, e_ndrho_laplace_rhoa, &
1213 2750 : e_ndrho_laplace_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
1214 2750 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa, &
1215 2750 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, e_ndrhoa_ndrhoa, &
1216 2750 : e_ndrhoa_ndrhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhoa_tau_a, &
1217 2750 : e_ndrhoa_tau_b, e_ndrhob
1218 2750 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_ndrhob_laplace_rhoa, &
1219 2750 : e_ndrhob_laplace_rhob, e_ndrhob_ndrhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
1220 2750 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_rhoa, e_rhoa_laplace_rhoa, &
1221 2750 : e_rhoa_laplace_rhob, e_rhoa_rhoa, e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
1222 2750 : e_rhoa_rhob, e_rhoa_rhob_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob, &
1223 2750 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhob_rhob, &
1224 2750 : e_rhob_rhob_rhob, e_rhob_tau_a, e_rhob_tau_b, e_tau_a, e_tau_a_tau_a, &
1225 2750 : e_tau_a_tau_b, e_tau_b, e_tau_b_tau_b, laplace_rhoa, laplace_rhob, &
1226 2750 : norm_drho, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
1227 : TYPE(xc_derivative_type), POINTER :: deriv
1228 : TYPE(xc_f03_func_info_t) :: xc_info
1229 :
1230 2750 : CALL timeset(routineN, handle)
1231 :
1232 : ! Only "everything up to grad_deriv" is supported. Asking for a single
1233 : ! derivative order in isolation is rejected here rather than further down,
1234 : ! where it would quietly evaluate nothing at all.
1235 2750 : IF (grad_deriv < 0) THEN
1236 0 : CPABORT("Evaluating a single derivative order is not supported.")
1237 : END IF
1238 2750 : IF (grad_deriv > 3) THEN
1239 0 : CPABORT("derivatives larger than 3 not implemented")
1240 : END IF
1241 :
1242 2750 : NULLIFY (dummy)
1243 2750 : NULLIFY (workers)
1244 2750 : NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, &
1245 2750 : laplace_rhob, tau_a, tau_b)
1246 :
1247 2750 : IF (PRESENT(func_name_override)) THEN
1248 0 : func_name = func_name_override
1249 0 : func_scale = 1.0_dp
1250 : ELSE
1251 2750 : func_name = libxc_params%section%name
1252 2750 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
1253 : END IF
1254 :
1255 2750 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
1256 :
1257 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
1258 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1259 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
1260 : laplace_rhoa=laplace_rhoa, laplace_rhob=laplace_rhob, &
1261 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
1262 2750 : tau_a=tau_a, tau_b=tau_b, local_bounds=bo)
1263 :
1264 2750 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1265 :
1266 : ! One functional object per thread, kept across calls. The cutoffs are part
1267 : ! of the worker setup because LibXC does the screening itself.
1268 2750 : IF (PRESENT(func_name_override)) THEN
1269 : ! an overriding name carries no input section, so no parameters are set
1270 0 : NULLIFY (no_params)
1271 0 : iw = libxc_get_workers(func_name, XC_POLARIZED, no_params, epsilon_rho, epsilon_tau)
1272 : ELSE
1273 2750 : iw = libxc_get_workers(func_name, XC_POLARIZED, libxc_params, epsilon_rho, epsilon_tau)
1274 : END IF
1275 2750 : workers => libxc_cache(iw)
1276 2750 : xc_info = workers%worker(1)%info
1277 2750 : has_laplace = workers%has_laplace
1278 :
1279 : ! see libxc_unused: the arguments the requested order does not produce are
1280 : ! never touched, but they still have to point somewhere
1281 2750 : CALL libxc_reserve_unused(npoints)
1282 2750 : dummy(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) => libxc_unused(1:npoints)
1283 :
1284 : ! due to assumed shape array usage in next routine
1285 2750 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => rhoa
1286 2750 : IF (.NOT. ASSOCIATED(norm_drhoa)) norm_drhoa => rhoa
1287 2750 : IF (.NOT. ASSOCIATED(norm_drhob)) norm_drhob => rhoa
1288 2750 : IF (.NOT. ASSOCIATED(tau_a)) tau_a => rhoa
1289 2750 : IF (.NOT. ASSOCIATED(tau_b)) tau_b => rhoa
1290 :
1291 : ! only some MGGA functionals really need the Laplacian,
1292 : ! all others can work with rhoa (read-only) as dummy
1293 2750 : IF (.NOT. has_laplace) laplace_rhoa => rhoa
1294 2562 : IF (.NOT. has_laplace) laplace_rhob => rhoa
1295 :
1296 2750 : e_0 => dummy
1297 2750 : e_rhoa => dummy
1298 2750 : e_rhob => dummy
1299 2750 : e_ndrho => dummy
1300 2750 : e_ndrhoa => dummy
1301 2750 : e_ndrhob => dummy
1302 2750 : e_laplace_rhoa => dummy
1303 2750 : e_laplace_rhob => dummy
1304 2750 : e_tau_a => dummy
1305 2750 : e_tau_b => dummy
1306 2750 : e_rhoa_rhoa => dummy
1307 2750 : e_rhoa_rhob => dummy
1308 2750 : e_rhob_rhob => dummy
1309 2750 : e_ndrho_rhoa => dummy
1310 2750 : e_ndrho_rhob => dummy
1311 2750 : e_ndrhoa_rhoa => dummy
1312 2750 : e_ndrhoa_rhob => dummy
1313 2750 : e_ndrhob_rhoa => dummy
1314 2750 : e_ndrhob_rhob => dummy
1315 2750 : e_ndrho_ndrho => dummy
1316 2750 : e_ndrho_ndrhoa => dummy
1317 2750 : e_ndrho_ndrhob => dummy
1318 2750 : e_ndrhoa_ndrhoa => dummy
1319 2750 : e_ndrhoa_ndrhob => dummy
1320 2750 : e_ndrhob_ndrhob => dummy
1321 2750 : e_rhoa_laplace_rhoa => dummy
1322 2750 : e_rhoa_laplace_rhob => dummy
1323 2750 : e_rhob_laplace_rhoa => dummy
1324 2750 : e_rhob_laplace_rhob => dummy
1325 2750 : e_rhoa_tau_a => dummy
1326 2750 : e_rhoa_tau_b => dummy
1327 2750 : e_rhob_tau_a => dummy
1328 2750 : e_rhob_tau_b => dummy
1329 2750 : e_ndrho_laplace_rhoa => dummy
1330 2750 : e_ndrho_laplace_rhob => dummy
1331 2750 : e_ndrhoa_laplace_rhoa => dummy
1332 2750 : e_ndrhoa_laplace_rhob => dummy
1333 2750 : e_ndrhob_laplace_rhoa => dummy
1334 2750 : e_ndrhob_laplace_rhob => dummy
1335 2750 : e_ndrho_tau_a => dummy
1336 2750 : e_ndrho_tau_b => dummy
1337 2750 : e_ndrhoa_tau_a => dummy
1338 2750 : e_ndrhoa_tau_b => dummy
1339 2750 : e_ndrhob_tau_a => dummy
1340 2750 : e_ndrhob_tau_b => dummy
1341 2750 : e_laplace_rhoa_laplace_rhoa => dummy
1342 2750 : e_laplace_rhoa_laplace_rhob => dummy
1343 2750 : e_laplace_rhob_laplace_rhob => dummy
1344 2750 : e_laplace_rhoa_tau_a => dummy
1345 2750 : e_laplace_rhoa_tau_b => dummy
1346 2750 : e_laplace_rhob_tau_a => dummy
1347 2750 : e_laplace_rhob_tau_b => dummy
1348 2750 : e_tau_a_tau_a => dummy
1349 2750 : e_tau_a_tau_b => dummy
1350 2750 : e_tau_b_tau_b => dummy
1351 2750 : e_rhoa_rhoa_rhoa => dummy
1352 2750 : e_rhoa_rhoa_rhob => dummy
1353 2750 : e_rhoa_rhob_rhob => dummy
1354 2750 : e_rhob_rhob_rhob => dummy
1355 :
1356 2750 : IF (grad_deriv >= 0) THEN
1357 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
1358 2750 : allocate_deriv=.TRUE.)
1359 2750 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1360 : END IF
1361 2750 : IF (grad_deriv >= 1) THEN
1362 1376 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1363 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1364 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1365 1376 : allocate_deriv=.TRUE.)
1366 1376 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1367 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1368 1376 : allocate_deriv=.TRUE.)
1369 1376 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1370 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1371 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1372 386 : allocate_deriv=.TRUE.)
1373 386 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1374 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1375 386 : allocate_deriv=.TRUE.)
1376 386 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1377 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1378 386 : allocate_deriv=.TRUE.)
1379 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1380 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1381 386 : allocate_deriv=.TRUE.)
1382 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
1383 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1384 386 : allocate_deriv=.TRUE.)
1385 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
1386 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1387 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
1388 942 : allocate_deriv=.TRUE.)
1389 942 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
1390 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
1391 942 : allocate_deriv=.TRUE.)
1392 942 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
1393 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
1394 942 : allocate_deriv=.TRUE.)
1395 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1396 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1397 942 : allocate_deriv=.TRUE.)
1398 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
1399 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1400 942 : allocate_deriv=.TRUE.)
1401 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
1402 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
1403 942 : allocate_deriv=.TRUE.)
1404 942 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
1405 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
1406 942 : allocate_deriv=.TRUE.)
1407 942 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
1408 942 : IF (has_laplace) THEN
1409 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
1410 180 : allocate_deriv=.TRUE.)
1411 180 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
1412 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
1413 180 : allocate_deriv=.TRUE.)
1414 180 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
1415 : END IF
1416 : CASE default
1417 2704 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1418 : END SELECT
1419 : END IF
1420 2750 : IF (grad_deriv >= 2) THEN
1421 38 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1422 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1423 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1424 38 : allocate_deriv=.TRUE.)
1425 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1426 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1427 38 : allocate_deriv=.TRUE.)
1428 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1429 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1430 38 : allocate_deriv=.TRUE.)
1431 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1432 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1433 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1434 34 : allocate_deriv=.TRUE.)
1435 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1436 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1437 34 : allocate_deriv=.TRUE.)
1438 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1439 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1440 34 : allocate_deriv=.TRUE.)
1441 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1442 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1443 34 : allocate_deriv=.TRUE.)
1444 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1445 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1446 34 : allocate_deriv=.TRUE.)
1447 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1448 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1449 34 : allocate_deriv=.TRUE.)
1450 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1451 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1452 34 : allocate_deriv=.TRUE.)
1453 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1454 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1455 34 : allocate_deriv=.TRUE.)
1456 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1457 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1458 34 : allocate_deriv=.TRUE.)
1459 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1460 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1461 34 : allocate_deriv=.TRUE.)
1462 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1463 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1464 34 : allocate_deriv=.TRUE.)
1465 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1466 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1467 34 : allocate_deriv=.TRUE.)
1468 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1469 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1470 34 : allocate_deriv=.TRUE.)
1471 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1472 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1473 34 : allocate_deriv=.TRUE.)
1474 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1475 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1476 34 : allocate_deriv=.TRUE.)
1477 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1478 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1479 :
1480 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1481 14 : allocate_deriv=.TRUE.)
1482 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1483 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1484 14 : allocate_deriv=.TRUE.)
1485 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1486 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1487 14 : allocate_deriv=.TRUE.)
1488 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1489 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1490 14 : allocate_deriv=.TRUE.)
1491 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1492 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1493 14 : allocate_deriv=.TRUE.)
1494 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1495 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1496 14 : allocate_deriv=.TRUE.)
1497 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1498 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1499 14 : allocate_deriv=.TRUE.)
1500 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1501 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1502 14 : allocate_deriv=.TRUE.)
1503 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1504 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1505 14 : allocate_deriv=.TRUE.)
1506 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1507 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1508 14 : allocate_deriv=.TRUE.)
1509 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1510 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1511 14 : allocate_deriv=.TRUE.)
1512 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1513 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1514 14 : allocate_deriv=.TRUE.)
1515 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1516 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1517 14 : allocate_deriv=.TRUE.)
1518 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1519 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1520 14 : allocate_deriv=.TRUE.)
1521 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1522 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1523 14 : allocate_deriv=.TRUE.)
1524 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1525 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_a], &
1526 14 : allocate_deriv=.TRUE.)
1527 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_a)
1528 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_b], &
1529 14 : allocate_deriv=.TRUE.)
1530 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_b)
1531 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_a], &
1532 14 : allocate_deriv=.TRUE.)
1533 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_a)
1534 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_b], &
1535 14 : allocate_deriv=.TRUE.)
1536 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_b)
1537 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_a], &
1538 14 : allocate_deriv=.TRUE.)
1539 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_a)
1540 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_b], &
1541 14 : allocate_deriv=.TRUE.)
1542 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_b)
1543 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_a], &
1544 14 : allocate_deriv=.TRUE.)
1545 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_a)
1546 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_b], &
1547 14 : allocate_deriv=.TRUE.)
1548 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_b)
1549 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_a], &
1550 14 : allocate_deriv=.TRUE.)
1551 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_a)
1552 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_b], &
1553 14 : allocate_deriv=.TRUE.)
1554 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_b)
1555 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_a], &
1556 14 : allocate_deriv=.TRUE.)
1557 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_a)
1558 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_b], &
1559 14 : allocate_deriv=.TRUE.)
1560 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_b)
1561 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b, deriv_tau_b], &
1562 14 : allocate_deriv=.TRUE.)
1563 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b_tau_b)
1564 14 : IF (has_laplace) THEN
1565 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhoa], &
1566 6 : allocate_deriv=.TRUE.)
1567 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhoa)
1568 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhob], &
1569 6 : allocate_deriv=.TRUE.)
1570 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhob)
1571 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhoa], &
1572 6 : allocate_deriv=.TRUE.)
1573 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhoa)
1574 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhob], &
1575 6 : allocate_deriv=.TRUE.)
1576 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhob)
1577 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhoa], &
1578 6 : allocate_deriv=.TRUE.)
1579 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhoa)
1580 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhob], &
1581 6 : allocate_deriv=.TRUE.)
1582 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhob)
1583 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhoa], &
1584 6 : allocate_deriv=.TRUE.)
1585 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhoa)
1586 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhob], &
1587 6 : allocate_deriv=.TRUE.)
1588 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhob)
1589 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhoa], &
1590 6 : allocate_deriv=.TRUE.)
1591 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhoa)
1592 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhob], &
1593 6 : allocate_deriv=.TRUE.)
1594 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhob)
1595 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhoa], &
1596 6 : allocate_deriv=.TRUE.)
1597 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhoa)
1598 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhob], &
1599 6 : allocate_deriv=.TRUE.)
1600 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhob)
1601 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_laplace_rhob], &
1602 6 : allocate_deriv=.TRUE.)
1603 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_laplace_rhob)
1604 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_a], &
1605 6 : allocate_deriv=.TRUE.)
1606 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_a)
1607 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_b], &
1608 6 : allocate_deriv=.TRUE.)
1609 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_b)
1610 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_a], &
1611 6 : allocate_deriv=.TRUE.)
1612 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_a)
1613 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_b], &
1614 6 : allocate_deriv=.TRUE.)
1615 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_b)
1616 : END IF
1617 : CASE default
1618 86 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1619 : END SELECT
1620 : END IF
1621 2750 : IF (grad_deriv >= 3) THEN
1622 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1623 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1624 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
1625 0 : allocate_deriv=.TRUE.)
1626 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhoa)
1627 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
1628 0 : allocate_deriv=.TRUE.)
1629 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhob)
1630 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
1631 0 : allocate_deriv=.TRUE.)
1632 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob_rhob)
1633 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
1634 0 : allocate_deriv=.TRUE.)
1635 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob_rhob)
1636 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1637 0 : CPABORT("derivatives larger than 2 not implemented")
1638 : CASE default
1639 0 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1640 : END SELECT
1641 : END IF
1642 :
1643 : !$OMP PARALLEL DEFAULT(NONE), &
1644 : !$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob),&
1645 : !$OMP SHARED(laplace_rhoa,laplace_rhob,tau_a,tau_b),&
1646 : !$OMP SHARED(e_0,e_rhoa,e_rhob,e_ndrho,e_ndrhoa,e_ndrhob),&
1647 : !$OMP SHARED(e_laplace_rhoa,e_laplace_rhob,e_tau_a,e_tau_b),&
1648 : !$OMP SHARED(e_rhoa_rhoa,e_rhoa_rhob,e_rhob_rhob),&
1649 : !$OMP SHARED(e_ndrho_rhoa,e_ndrho_rhob),&
1650 : !$OMP SHARED(e_ndrhoa_rhoa,e_ndrhoa_rhob,e_ndrhob_rhoa,e_ndrhob_rhob),&
1651 : !$OMP SHARED(e_ndrho_ndrho,e_ndrho_ndrhoa,e_ndrho_ndrhob),&
1652 : !$OMP SHARED(e_ndrhoa_ndrhoa,e_ndrhoa_ndrhob,e_ndrhob_ndrhob),&
1653 : !$OMP SHARED(e_rhoa_laplace_rhoa,e_rhoa_laplace_rhob,e_rhob_laplace_rhoa,e_rhob_laplace_rhob),&
1654 : !$OMP SHARED(e_rhoa_tau_a,e_rhoa_tau_b,e_rhob_tau_a,e_rhob_tau_b),&
1655 : !$OMP SHARED(e_ndrho_laplace_rhoa,e_ndrho_laplace_rhob),&
1656 : !$OMP SHARED(e_ndrhoa_laplace_rhoa,e_ndrhoa_laplace_rhob,e_ndrhob_laplace_rhoa,e_ndrhob_laplace_rhob),&
1657 : !$OMP SHARED(e_ndrho_tau_a,e_ndrho_tau_b),&
1658 : !$OMP SHARED(e_ndrhoa_tau_a,e_ndrhoa_tau_b,e_ndrhob_tau_a,e_ndrhob_tau_b),&
1659 : !$OMP SHARED(e_laplace_rhoa_laplace_rhoa,e_laplace_rhoa_laplace_rhob,e_laplace_rhob_laplace_rhob),&
1660 : !$OMP SHARED(e_laplace_rhoa_tau_a,e_laplace_rhoa_tau_b,e_laplace_rhob_tau_a,e_laplace_rhob_tau_b),&
1661 : !$OMP SHARED(e_tau_a_tau_a,e_tau_a_tau_b,e_tau_b_tau_b),&
1662 : !$OMP SHARED(e_rhoa_rhoa_rhoa,e_rhoa_rhoa_rhob,e_rhoa_rhob_rhob,e_rhob_rhob_rhob),&
1663 : !$OMP SHARED(grad_deriv,npoints),&
1664 2750 : !$OMP SHARED(func_name,func_scale,workers)
1665 :
1666 : CALL libxc_spin_polarized_calc(rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1667 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, laplace_rhoa=laplace_rhoa, &
1668 : laplace_rhob=laplace_rhob, tau_a=tau_a, tau_b=tau_b, &
1669 : e_0=e_0, e_rhoa=e_rhoa, e_rhob=e_rhob, e_ndrho=e_ndrho, &
1670 : e_ndrhoa=e_ndrhoa, e_ndrhob=e_ndrhob, e_laplace_rhoa=e_laplace_rhoa, &
1671 : e_laplace_rhob=e_laplace_rhob, e_tau_a=e_tau_a, e_tau_b=e_tau_b, &
1672 : e_rhoa_rhoa=e_rhoa_rhoa, e_rhoa_rhob=e_rhoa_rhob, e_rhob_rhob=e_rhob_rhob, &
1673 : e_ndrho_rhoa=e_ndrho_rhoa, e_ndrho_rhob=e_ndrho_rhob, &
1674 : e_ndrhoa_rhoa=e_ndrhoa_rhoa, e_ndrhoa_rhob=e_ndrhoa_rhob, &
1675 : e_ndrhob_rhoa=e_ndrhob_rhoa, e_ndrhob_rhob=e_ndrhob_rhob, &
1676 : e_ndrho_ndrho=e_ndrho_ndrho, e_ndrho_ndrhoa=e_ndrho_ndrhoa, &
1677 : e_ndrho_ndrhob=e_ndrho_ndrhob, e_ndrhoa_ndrhoa=e_ndrhoa_ndrhoa, &
1678 : e_ndrhoa_ndrhob=e_ndrhoa_ndrhob, e_ndrhob_ndrhob=e_ndrhob_ndrhob, &
1679 : e_rhoa_laplace_rhoa=e_rhoa_laplace_rhoa, &
1680 : e_rhoa_laplace_rhob=e_rhoa_laplace_rhob, &
1681 : e_rhob_laplace_rhoa=e_rhob_laplace_rhoa, &
1682 : e_rhob_laplace_rhob=e_rhob_laplace_rhob, &
1683 : e_rhoa_tau_a=e_rhoa_tau_a, e_rhoa_tau_b=e_rhoa_tau_b, &
1684 : e_rhob_tau_a=e_rhob_tau_a, e_rhob_tau_b=e_rhob_tau_b, &
1685 : e_ndrho_laplace_rhoa=e_ndrho_laplace_rhoa, &
1686 : e_ndrho_laplace_rhob=e_ndrho_laplace_rhob, &
1687 : e_ndrhoa_laplace_rhoa=e_ndrhoa_laplace_rhoa, &
1688 : e_ndrhoa_laplace_rhob=e_ndrhoa_laplace_rhob, &
1689 : e_ndrhob_laplace_rhoa=e_ndrhob_laplace_rhoa, &
1690 : e_ndrhob_laplace_rhob=e_ndrhob_laplace_rhob, &
1691 : e_ndrho_tau_a=e_ndrho_tau_a, e_ndrho_tau_b=e_ndrho_tau_b, &
1692 : e_ndrhoa_tau_a=e_ndrhoa_tau_a, e_ndrhoa_tau_b=e_ndrhoa_tau_b, &
1693 : e_ndrhob_tau_a=e_ndrhob_tau_a, e_ndrhob_tau_b=e_ndrhob_tau_b, &
1694 : e_laplace_rhoa_laplace_rhoa=e_laplace_rhoa_laplace_rhoa, &
1695 : e_laplace_rhoa_laplace_rhob=e_laplace_rhoa_laplace_rhob, &
1696 : e_laplace_rhob_laplace_rhob=e_laplace_rhob_laplace_rhob, &
1697 : e_laplace_rhoa_tau_a=e_laplace_rhoa_tau_a, &
1698 : e_laplace_rhoa_tau_b=e_laplace_rhoa_tau_b, &
1699 : e_laplace_rhob_tau_a=e_laplace_rhob_tau_a, &
1700 : e_laplace_rhob_tau_b=e_laplace_rhob_tau_b, &
1701 : e_tau_a_tau_a=e_tau_a_tau_a, &
1702 : e_tau_a_tau_b=e_tau_a_tau_b, &
1703 : e_tau_b_tau_b=e_tau_b_tau_b, &
1704 : e_rhoa_rhoa_rhoa=e_rhoa_rhoa_rhoa, &
1705 : e_rhoa_rhoa_rhob=e_rhoa_rhoa_rhob, &
1706 : e_rhoa_rhob_rhob=e_rhoa_rhob_rhob, &
1707 : e_rhob_rhob_rhob=e_rhob_rhob_rhob, &
1708 : grad_deriv=grad_deriv, npoints=npoints, &
1709 : func_name=func_name, sc=func_scale, workers=workers)
1710 :
1711 : !$OMP END PARALLEL
1712 :
1713 2750 : NULLIFY (dummy)
1714 2750 : NULLIFY (workers)
1715 :
1716 2750 : CALL timestop(handle)
1717 : #else
1718 : MARK_USED(rho_set)
1719 : MARK_USED(deriv_set)
1720 : MARK_USED(grad_deriv)
1721 : MARK_USED(libxc_params)
1722 : MARK_USED(func_name_override)
1723 :
1724 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
1725 : "for a functional of the LibXC library, "// &
1726 : "you have to download and install the library!")
1727 : #endif
1728 2750 : END SUBROUTINE libxc_spin_polarized_eval
1729 :
1730 : ! **************************************************************************************************
1731 : !> \brief libxc exchange-correlation functionals
1732 : !> \param rho density
1733 : !> \param norm_drho norm of the gradient of the density
1734 : !> \param laplace_rho laplacian of the density
1735 : !> \param tau kinetic-energy density
1736 : !> \param e_0 energy density
1737 : !> \param e_rho derivative of the energy density with respect to rho
1738 : !> \param e_ndrho derivative of the energy density with respect to ndrho
1739 : !> \param e_laplace_rho derivative of the energy density with respect to laplace_rho
1740 : !> \param e_tau derivative of the energy density with respect to tau
1741 : !> \param e_rho_rho derivative of the energy density with respect to rho_rho
1742 : !> \param e_ndrho_rho derivative of the energy density with respect to ndrho_rho
1743 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
1744 : !> \param e_rho_laplace_rho derivative of the energy density with respect to rho_laplace_rho
1745 : !> \param e_rho_tau derivative of the energy density with respect to rho_tau
1746 : !> \param e_ndrho_laplace_rho derivative of the energy density with respect to ndrho_laplace_rho
1747 : !> \param e_ndrho_tau derivative of the energy density with respect to ndrho_tau
1748 : !> \param e_laplace_rho_laplace_rho derivative of the energy density with respect to laplace_rho_laplace_rho
1749 : !> \param e_laplace_rho_tau derivative of the energy density with respect to laplace_rho_tau
1750 : !> \param e_tau_tau derivative of the energy density with respect to tau_tau
1751 : !> \param e_rho_rho_rho derivative of the energy density with respect to rho_rho_rho
1752 : !> \param grad_deriv degree of the derivative that should be evaluated;
1753 : !> all derivatives up to the given degree are evaluated, in a single
1754 : !> LibXC call per block of grid points
1755 : !> \param npoints number of points on the grid
1756 : !> \param func_name name of the functional
1757 : !> \param sc scaling factor of the functional
1758 : !> \param workers cached LibXC functional objects and staging buffers, one per thread
1759 : !> \author F. Tran
1760 : ! **************************************************************************************************
1761 : #if defined (__LIBXC)
1762 17700 : SUBROUTINE libxc_spin_unpolarized_calc(rho, norm_drho, laplace_rho, tau, &
1763 : e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, e_rho_rho, e_ndrho_rho, &
1764 : e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1765 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, &
1766 : e_tau_tau, e_rho_rho_rho, &
1767 : grad_deriv, npoints, func_name, sc, workers)
1768 :
1769 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho, laplace_rho, tau
1770 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, &
1771 : e_rho_rho, e_ndrho_rho, e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1772 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_tau_tau, e_rho_rho_rho
1773 : INTEGER, INTENT(in) :: grad_deriv, npoints
1774 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
1775 : REAL(KIND=dp), INTENT(in) :: sc
1776 : TYPE(libxc_worker_set_type), INTENT(INOUT) :: workers
1777 :
1778 : INTEGER :: bsize, family, i, i0, ib, ii, ithread, &
1779 : nb, nblocks, nthreads
1780 : INTEGER(C_SIZE_T) :: np
1781 : LOGICAL :: is_gga, is_mgga
1782 :
1783 17700 : ithread = 0
1784 17700 : nthreads = 1
1785 17700 : !$ ithread = omp_get_thread_num()
1786 17700 : !$ nthreads = omp_get_num_threads()
1787 17700 : CPASSERT(ithread < SIZE(workers%worker))
1788 :
1789 : ! Blocks are capped by libxc_block_size, which is what the staging buffers
1790 : ! were allocated for, but shrunk when there are too few points to give every
1791 : ! thread a block of its own: small atomic grids would otherwise all be
1792 : ! evaluated by thread 0.
1793 17700 : bsize = MAX(1, MIN(libxc_block_size, (npoints + nthreads - 1)/nthreads))
1794 17700 : nblocks = (npoints + bsize - 1)/bsize
1795 :
1796 17700 : family = workers%family
1797 : SELECT CASE (family)
1798 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1799 : is_gga = .FALSE.
1800 5268 : is_mgga = .FALSE.
1801 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1802 5268 : is_gga = .TRUE.
1803 5268 : is_mgga = .FALSE.
1804 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1805 2138 : is_gga = .TRUE.
1806 2138 : is_mgga = .TRUE.
1807 : CASE default
1808 0 : is_gga = .FALSE.
1809 0 : is_mgga = .FALSE.
1810 17700 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1811 : END SELECT
1812 :
1813 : ! Points below the density cutoff need no test here. The cutoffs were handed
1814 : ! to the functional object in xc_libxc_wrap_set_thresholds, and LibXC leaves
1815 : ! the outputs of the points it screens at zero, so accumulating them is a
1816 : ! no-op.
1817 : ASSOCIATE (w => workers%worker(ithread + 1), has_laplace => workers%has_laplace, &
1818 : no_exc => workers%no_exc, eps_rho => workers%epsilon_rho, &
1819 : eps_tau => workers%epsilon_tau)
1820 17700 : !$OMP DO
1821 : DO ib = 1, nblocks
1822 1018131 : i0 = (ib - 1)*bsize
1823 1018131 : nb = MIN(bsize, npoints - i0)
1824 1018131 : np = INT(nb, KIND=C_SIZE_T)
1825 :
1826 : ! stage this block's inputs in LibXC's layout
1827 520450624 : DO i = 1, nb
1828 520450624 : w%rho(1, i) = rho(i0 + i)
1829 : END DO
1830 1018131 : IF (is_gga) THEN
1831 273501950 : DO i = 1, nb
1832 273501950 : w%sigma(1, i) = norm_drho(i0 + i)**2
1833 : END DO
1834 : END IF
1835 1018131 : IF (is_mgga) THEN
1836 102686484 : DO i = 1, nb
1837 102485437 : ii = i0 + i
1838 102485437 : w%lapl(1, i) = laplace_rho(ii)
1839 102686484 : w%tau(1, i) = tau(ii)
1840 : END DO
1841 : ! Meta-GGAs are screened on the kinetic energy density as well as on the
1842 : ! density, and LibXC only screens on the latter. A point failing the tau
1843 : ! test is handed over with zero density, which makes LibXC screen it out
1844 : ! and leave its outputs at zero: what skipping it used to achieve. The
1845 : ! Fermi hole curvature bound tau >= sigma/(8*rho) follows in the same pass,
1846 : ! with the denominator floored at the cutoff so the division stays defined
1847 : ! for the points that are about to be screened.
1848 102686484 : DO i = 1, nb
1849 102485437 : IF (w%tau(1, i) <= eps_tau) w%rho(1, i) = 0.0_dp
1850 : w%tau(1, i) = MAX(w%tau(1, i), &
1851 102686484 : w%sigma(1, i)/(8.0_dp*MAX(w%rho(1, i), eps_rho)))
1852 : END DO
1853 : END IF
1854 :
1855 : ! one LibXC call for the whole block
1856 : IF (is_mgga) THEN
1857 202091 : SELECT CASE (grad_deriv)
1858 : CASE (0)
1859 1044 : CALL xc_f03_mgga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), w%exc(1))
1860 : CASE (1)
1861 191303 : IF (no_exc) THEN
1862 : CALL xc_f03_mgga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
1863 0 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
1864 0 : w%exc(1:nb) = 0.0_dp
1865 : ELSE
1866 : CALL xc_f03_mgga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
1867 191303 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
1868 : END IF
1869 : CASE (2)
1870 201047 : IF (no_exc) THEN
1871 : CALL xc_f03_mgga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
1872 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
1873 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
1874 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
1875 0 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
1876 0 : w%exc(1:nb) = 0.0_dp
1877 : ELSE
1878 : CALL xc_f03_mgga(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
1879 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
1880 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
1881 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
1882 8700 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
1883 : END IF
1884 : END SELECT
1885 817084 : ELSE IF (is_gga) THEN
1886 346132 : SELECT CASE (grad_deriv)
1887 : CASE (0)
1888 11962 : CALL xc_f03_gga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%exc(1))
1889 : CASE (1)
1890 316224 : IF (no_exc) THEN
1891 0 : CALL xc_f03_gga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%vrho(1, 1), w%vsigma(1, 1))
1892 0 : w%exc(1:nb) = 0.0_dp
1893 : ELSE
1894 : CALL xc_f03_gga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
1895 316224 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1))
1896 : END IF
1897 : CASE (2)
1898 334170 : IF (no_exc) THEN
1899 : CALL xc_f03_gga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
1900 : w%vrho(1, 1), w%vsigma(1, 1), &
1901 0 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
1902 0 : w%exc(1:nb) = 0.0_dp
1903 : ELSE
1904 : CALL xc_f03_gga_exc_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
1905 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
1906 5984 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
1907 : END IF
1908 : END SELECT
1909 : ELSE
1910 495456 : SELECT CASE (grad_deriv)
1911 : CASE (0)
1912 12542 : CALL xc_f03_lda_exc(w%func, np, w%rho(1, 1), w%exc(1))
1913 : CASE (1)
1914 450640 : CALL xc_f03_lda_exc_vxc(w%func, np, w%rho(1, 1), w%exc(1), w%vrho(1, 1))
1915 : CASE (2)
1916 : CALL xc_f03_lda_exc_vxc_fxc(w%func, np, w%rho(1, 1), &
1917 19732 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1))
1918 : CASE (3)
1919 : CALL xc_f03_lda(w%func, np, w%rho(1, 1), &
1920 482914 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1), w%v3rho3(1, 1))
1921 : END SELECT
1922 : END IF
1923 :
1924 : ! accumulate; each derivative is added by exactly one loop, selected by
1925 : ! the order requested and the family, not by a per-branch copy of the whole
1926 : ! block scaffolding
1927 1018131 : IF (grad_deriv >= 0) THEN
1928 520450624 : DO i = 1, nb
1929 519432493 : ii = i0 + i
1930 520450624 : e_0(ii) = e_0(ii) + sc*w%exc(i)*rho(ii)
1931 : END DO
1932 : END IF
1933 1018131 : IF (grad_deriv >= 1) THEN
1934 507383128 : DO i = 1, nb
1935 506390545 : ii = i0 + i
1936 507383128 : e_rho(ii) = e_rho(ii) + sc*w%vrho(1, i)
1937 : END DO
1938 992583 : IF (is_gga) THEN
1939 266857950 : DO i = 1, nb
1940 266335739 : ii = i0 + i
1941 266857950 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*w%vsigma(1, i)*norm_drho(ii)
1942 : END DO
1943 : END IF
1944 992583 : IF (is_mgga) THEN
1945 102159440 : DO i = 1, nb
1946 101959437 : ii = i0 + i
1947 102159440 : e_tau(ii) = e_tau(ii) + sc*w%vtau(1, i)
1948 : END DO
1949 : END IF
1950 200003 : IF (is_mgga .AND. has_laplace) THEN
1951 9918194 : DO i = 1, nb
1952 9898632 : ii = i0 + i
1953 9918194 : e_laplace_rho(ii) = e_laplace_rho(ii) + sc*w%vlapl(1, i)
1954 : END DO
1955 : END IF
1956 : END IF
1957 1018131 : IF (grad_deriv >= 2) THEN
1958 17477900 : DO i = 1, nb
1959 17443484 : ii = i0 + i
1960 17477900 : e_rho_rho(ii) = e_rho_rho(ii) + sc*w%v2rho2(1, i)
1961 : END DO
1962 34416 : IF (is_gga) THEN
1963 7461668 : DO i = 1, nb
1964 7446984 : ii = i0 + i
1965 7446984 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*w%v2rhosigma(1, i)*norm_drho(ii)
1966 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1967 7461668 : sc*2.0_dp*(2.0_dp*w%sigma(1, i)*w%v2sigma2(1, i) + w%vsigma(1, i))
1968 : END DO
1969 : END IF
1970 34416 : IF (is_mgga) THEN
1971 4433996 : DO i = 1, nb
1972 4425296 : ii = i0 + i
1973 4425296 : e_rho_tau(ii) = e_rho_tau(ii) + sc*w%v2rhotau(1, i)
1974 4425296 : e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*w%v2sigmatau(1, i)*norm_drho(ii)
1975 4433996 : e_tau_tau(ii) = e_tau_tau(ii) + sc*w%v2tau2(1, i)
1976 : END DO
1977 : END IF
1978 8700 : IF (is_mgga .AND. has_laplace) THEN
1979 1321422 : DO i = 1, nb
1980 1318824 : ii = i0 + i
1981 1318824 : e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*w%v2rholapl(1, i)
1982 : e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
1983 1318824 : sc*2.0_dp*w%v2sigmalapl(1, i)*norm_drho(ii)
1984 1318824 : e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*w%v2lapl2(1, i)
1985 1321422 : e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*w%v2lapltau(1, i)
1986 : END DO
1987 : END IF
1988 : END IF
1989 1018131 : IF (grad_deriv >= 3) THEN
1990 0 : DO i = 1, nb
1991 0 : ii = i0 + i
1992 0 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*w%v3rho3(1, i)
1993 : END DO
1994 : END IF
1995 : END DO
1996 : !$OMP END DO
1997 : END ASSOCIATE
1998 :
1999 17700 : END SUBROUTINE libxc_spin_unpolarized_calc
2000 : #endif
2001 :
2002 : ! **************************************************************************************************
2003 : !> \brief libxc exchange-correlation functionals
2004 : !> \param rhoa alpha density
2005 : !> \param rhob beta density
2006 : !> \param norm_drho ...
2007 : !> \param norm_drhoa norm of the gradient of the alpha density
2008 : !> \param norm_drhob norm of the gradient of the beta density
2009 : !> \param laplace_rhoa laplacian of the alpha density
2010 : !> \param laplace_rhob laplacian of the beta density
2011 : !> \param tau_a alpha kinetic-energy density
2012 : !> \param tau_b beta kinetic-energy density
2013 : !> \param e_0 energy density
2014 : !> \param e_rhoa derivative of the energy density with respect to rhoa
2015 : !> \param e_rhob derivative of the energy density with respect to rhob
2016 : !> \param e_ndrho derivative of the energy density with respect to ndrho
2017 : !> \param e_ndrhoa derivative of the energy density with respect to ndrhoa
2018 : !> \param e_ndrhob derivative of the energy density with respect to ndrhob
2019 : !> \param e_laplace_rhoa derivative of the energy density with respect to laplace_rhoa
2020 : !> \param e_laplace_rhob derivative of the energy density with respect to laplace_rhob
2021 : !> \param e_tau_a derivative of the energy density with respect to tau_a
2022 : !> \param e_tau_b derivative of the energy density with respect to tau_b
2023 : !> \param e_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa
2024 : !> \param e_rhoa_rhob derivative of the energy density with respect to rhoa_rhob
2025 : !> \param e_rhob_rhob derivative of the energy density with respect to rhob_rhob
2026 : !> \param e_ndrho_rhoa derivative of the energy density with respect to ndrho_rhoa
2027 : !> \param e_ndrho_rhob derivative of the energy density with respect to ndrho_rhob
2028 : !> \param e_ndrhoa_rhoa derivative of the energy density with respect to ndrhoa_rhoa
2029 : !> \param e_ndrhoa_rhob derivative of the energy density with respect to ndrhoa_rhob
2030 : !> \param e_ndrhob_rhoa derivative of the energy density with respect to ndrhob_rhoa
2031 : !> \param e_ndrhob_rhob derivative of the energy density with respect to ndrhob_rhob
2032 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
2033 : !> \param e_ndrho_ndrhoa derivative of the energy density with respect to ndrho_ndrhoa
2034 : !> \param e_ndrho_ndrhob derivative of the energy density with respect to ndrho_ndrhob
2035 : !> \param e_ndrhoa_ndrhoa derivative of the energy density with respect to ndrhoa_ndrhoa
2036 : !> \param e_ndrhoa_ndrhob derivative of the energy density with respect to ndrhoa_ndrhob
2037 : !> \param e_ndrhob_ndrhob derivative of the energy density with respect to ndrhob_ndrhob
2038 : !> \param e_rhoa_laplace_rhoa derivative of the energy density with respect to rhoa_laplace_rhoa
2039 : !> \param e_rhoa_laplace_rhob derivative of the energy density with respect to rhoa_laplace_rhob
2040 : !> \param e_rhob_laplace_rhoa derivative of the energy density with respect to rhob_laplace_rhoa
2041 : !> \param e_rhob_laplace_rhob derivative of the energy density with respect to rhob_laplace_rhob
2042 : !> \param e_rhoa_tau_a derivative of the energy density with respect to rhoa_tau_a
2043 : !> \param e_rhoa_tau_b derivative of the energy density with respect to rhoa_tau_b
2044 : !> \param e_rhob_tau_a derivative of the energy density with respect to rhob_tau_a
2045 : !> \param e_rhob_tau_b derivative of the energy density with respect to rhob_tau_b
2046 : !> \param e_ndrho_laplace_rhoa derivative of the energy density with respect to ndrho_laplace_rhoa
2047 : !> \param e_ndrho_laplace_rhob derivative of the energy density with respect to ndrho_laplace_rhob
2048 : !> \param e_ndrhoa_laplace_rhoa derivative of the energy density with respect to ndrhoa_laplace_rhoa
2049 : !> \param e_ndrhoa_laplace_rhob derivative of the energy density with respect to ndrhoa_laplace_rhob
2050 : !> \param e_ndrhob_laplace_rhoa derivative of the energy density with respect to ndrhob_laplace_rhoa
2051 : !> \param e_ndrhob_laplace_rhob derivative of the energy density with respect to ndrhob_laplace_rhob
2052 : !> \param e_ndrho_tau_a derivative of the energy density with respect to ndrho_tau_a
2053 : !> \param e_ndrho_tau_b derivative of the energy density with respect to ndrho_tau_b
2054 : !> \param e_ndrhoa_tau_a derivative of the energy density with respect to ndrhoa_tau_a
2055 : !> \param e_ndrhoa_tau_b derivative of the energy density with respect to ndrhoa_tau_b
2056 : !> \param e_ndrhob_tau_a derivative of the energy density with respect to ndrhob_tau_a
2057 : !> \param e_ndrhob_tau_b derivative of the energy density with respect to ndrhob_tau_b
2058 : !> \param e_laplace_rhoa_laplace_rhoa derivative of the energy density with respect to laplace_rhoa_laplace_rhoa
2059 : !> \param e_laplace_rhoa_laplace_rhob derivative of the energy density with respect to laplace_rhoa_laplace_rhob
2060 : !> \param e_laplace_rhob_laplace_rhob derivative of the energy density with respect to laplace_rhob_laplace_rhob
2061 : !> \param e_laplace_rhoa_tau_a derivative of the energy density with respect to laplace_rhoa_tau_a
2062 : !> \param e_laplace_rhoa_tau_b derivative of the energy density with respect to laplace_rhoa_tau_b
2063 : !> \param e_laplace_rhob_tau_a derivative of the energy density with respect to laplace_rhob_tau_a
2064 : !> \param e_laplace_rhob_tau_b derivative of the energy density with respect to laplace_rhob_tau_b
2065 : !> \param e_tau_a_tau_a derivative of the energy density with respect to tau_a_tau_a
2066 : !> \param e_tau_a_tau_b derivative of the energy density with respect to tau_a_tau_b
2067 : !> \param e_tau_b_tau_b derivative of the energy density with respect to tau_b_tau_b
2068 : !> \param e_rhoa_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa_rhoa
2069 : !> \param e_rhoa_rhoa_rhob derivative of the energy density with respect to rhoa_rhoa_rhob
2070 : !> \param e_rhoa_rhob_rhob derivative of the energy density with respect to rhoa_rhob_rhob
2071 : !> \param e_rhob_rhob_rhob derivative of the energy density with respect to rhob_rhob_rhob
2072 : !> \param grad_deriv degree of the derivative that should be evaluated;
2073 : !> all derivatives up to the given degree are evaluated, in a single
2074 : !> LibXC call per block of grid points
2075 : !> \param npoints number of points on the grid
2076 : !> \param func_name name of the functional
2077 : !> \param sc scaling factor of the functional
2078 : !> \param workers cached LibXC functional objects and staging buffers, one per thread
2079 : !> \author F. Tran
2080 : ! **************************************************************************************************
2081 : #if defined (__LIBXC)
2082 2750 : SUBROUTINE libxc_spin_polarized_calc(rhoa, rhob, norm_drho, norm_drhoa, &
2083 : norm_drhob, laplace_rhoa, laplace_rhob, tau_a, tau_b, &
2084 : e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, e_ndrhob, &
2085 : e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, &
2086 : e_rhoa_rhoa, e_rhoa_rhob, e_rhob_rhob, &
2087 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, &
2088 : e_ndrhoa_rhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
2089 : e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
2090 : e_ndrhoa_ndrhoa, e_ndrhoa_ndrhob, e_ndrhob_ndrhob, &
2091 : e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
2092 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, &
2093 : e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, e_rhob_tau_b, &
2094 : e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, &
2095 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, &
2096 : e_ndrhob_laplace_rhoa, e_ndrhob_laplace_rhob, &
2097 : e_ndrho_tau_a, e_ndrho_tau_b, &
2098 : e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
2099 : e_ndrhob_tau_a, e_ndrhob_tau_b, &
2100 : e_laplace_rhoa_laplace_rhoa, &
2101 : e_laplace_rhoa_laplace_rhob, &
2102 : e_laplace_rhob_laplace_rhob, &
2103 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
2104 : e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, &
2105 : e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
2106 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
2107 : e_rhoa_rhob_rhob, e_rhob_rhob_rhob, &
2108 : grad_deriv, npoints, func_name, sc, workers)
2109 :
2110 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, norm_drho, norm_drhoa, &
2111 : norm_drhob, laplace_rhoa, &
2112 : laplace_rhob, tau_a, tau_b
2113 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, &
2114 : e_ndrhob, e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, e_rhoa_rhoa, e_rhoa_rhob, &
2115 : e_rhob_rhob, e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhob_rhoa, &
2116 : e_ndrhob_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, e_ndrhoa_ndrhoa, &
2117 : e_ndrhoa_ndrhob, e_ndrhob_ndrhob, e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
2118 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, &
2119 : e_rhob_tau_b, e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, e_ndrhoa_laplace_rhoa
2120 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_ndrhoa_laplace_rhob, e_ndrhob_laplace_rhoa, &
2121 : e_ndrhob_laplace_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
2122 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
2123 : e_laplace_rhob_laplace_rhob, e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
2124 : 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, &
2125 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, e_rhoa_rhob_rhob, e_rhob_rhob_rhob
2126 : INTEGER, INTENT(in) :: grad_deriv, npoints
2127 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
2128 : REAL(KIND=dp), INTENT(in) :: sc
2129 : TYPE(libxc_worker_set_type), INTENT(INOUT) :: workers
2130 :
2131 : INTEGER :: bsize, family, i, i0, ib, ii, ithread, &
2132 : nb, nblocks, nthreads
2133 : INTEGER(C_SIZE_T) :: np
2134 : LOGICAL :: is_gga, is_mgga
2135 :
2136 2750 : ithread = 0
2137 2750 : nthreads = 1
2138 2750 : !$ ithread = omp_get_thread_num()
2139 2750 : !$ nthreads = omp_get_num_threads()
2140 2750 : CPASSERT(ithread < SIZE(workers%worker))
2141 :
2142 2750 : bsize = MAX(1, MIN(libxc_block_size, (npoints + nthreads - 1)/nthreads))
2143 2750 : nblocks = (npoints + bsize - 1)/bsize
2144 :
2145 2750 : family = workers%family
2146 : SELECT CASE (family)
2147 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
2148 : is_gga = .FALSE.
2149 402 : is_mgga = .FALSE.
2150 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
2151 402 : is_gga = .TRUE.
2152 402 : is_mgga = .FALSE.
2153 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
2154 972 : is_gga = .TRUE.
2155 972 : is_mgga = .TRUE.
2156 : CASE default
2157 0 : is_gga = .FALSE.
2158 0 : is_mgga = .FALSE.
2159 2750 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
2160 : END SELECT
2161 :
2162 : ! As in the spin-unpolarized case the density and kinetic energy density
2163 : ! cutoffs are applied by LibXC itself.
2164 : ASSOCIATE (w => workers%worker(ithread + 1), has_laplace => workers%has_laplace, &
2165 : no_exc => workers%no_exc, eps_rho => workers%epsilon_rho, &
2166 : eps_tau => workers%epsilon_tau)
2167 2750 : !$OMP DO
2168 : DO ib = 1, nblocks
2169 154032 : i0 = (ib - 1)*bsize
2170 154032 : nb = MIN(bsize, npoints - i0)
2171 154032 : np = INT(nb, KIND=C_SIZE_T)
2172 :
2173 : ! stage this block's inputs in LibXC's layout
2174 78715696 : DO i = 1, nb
2175 78561664 : ii = i0 + i
2176 78561664 : w%rho(1, i) = MAX(rhoa(ii), 0.0_dp)
2177 78715696 : w%rho(2, i) = MAX(rhob(ii), 0.0_dp)
2178 : END DO
2179 154032 : IF (is_gga) THEN
2180 : ! CP2K works with the norms of the gradients, LibXC with their contractions
2181 18077738 : DO i = 1, nb
2182 18042076 : ii = i0 + i
2183 18042076 : w%nda(i) = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2184 18042076 : w%ndb(i) = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2185 18042076 : w%nd(i) = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2186 18042076 : w%sigma(1, i) = w%nda(i)**2
2187 18042076 : w%sigma(3, i) = w%ndb(i)**2
2188 18077738 : w%sigma(2, i) = 0.5_dp*(w%nd(i)**2 - w%sigma(1, i) - w%sigma(3, i))
2189 : END DO
2190 : END IF
2191 154032 : IF (is_mgga) THEN
2192 9495912 : DO i = 1, nb
2193 9477048 : ii = i0 + i
2194 9477048 : w%lapl(1, i) = laplace_rhoa(ii)
2195 9477048 : w%lapl(2, i) = laplace_rhob(ii)
2196 9477048 : w%tau(1, i) = MAX(tau_a(ii), 0.0_dp)
2197 9495912 : w%tau(2, i) = MAX(tau_b(ii), 0.0_dp)
2198 : END DO
2199 : ! screening on tau and the Fermi hole curvature bound, per spin channel;
2200 : ! see the spin-unpolarized routine. The tau test is on the sum of the two
2201 : ! channels, as the density test is.
2202 9495912 : DO i = 1, nb
2203 9477048 : IF (w%tau(1, i) + w%tau(2, i) <= eps_tau) THEN
2204 10732 : w%rho(1, i) = 0.0_dp
2205 10732 : w%rho(2, i) = 0.0_dp
2206 : END IF
2207 : w%tau(1, i) = MAX(w%tau(1, i), &
2208 9477048 : w%sigma(1, i)/(8.0_dp*MAX(w%rho(1, i), eps_rho)))
2209 : w%tau(2, i) = MAX(w%tau(2, i), &
2210 9495912 : w%sigma(3, i)/(8.0_dp*MAX(w%rho(2, i), eps_rho)))
2211 : END DO
2212 : END IF
2213 :
2214 : ! one LibXC call for the whole block
2215 : IF (is_mgga) THEN
2216 19492 : SELECT CASE (grad_deriv)
2217 : CASE (0)
2218 628 : CALL xc_f03_mgga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), w%exc(1))
2219 : CASE (1)
2220 18040 : IF (no_exc) THEN
2221 : CALL xc_f03_mgga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2222 0 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2223 0 : w%exc(1:nb) = 0.0_dp
2224 : ELSE
2225 : CALL xc_f03_mgga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2226 18040 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1))
2227 : END IF
2228 : CASE (2)
2229 18864 : IF (no_exc) THEN
2230 : CALL xc_f03_mgga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2231 : w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2232 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2233 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2234 0 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2235 0 : w%exc(1:nb) = 0.0_dp
2236 : ELSE
2237 : CALL xc_f03_mgga(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%lapl(1, 1), w%tau(1, 1), &
2238 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), w%vlapl(1, 1), w%vtau(1, 1), &
2239 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2rholapl(1, 1), &
2240 : w%v2rhotau(1, 1), w%v2sigma2(1, 1), w%v2sigmalapl(1, 1), &
2241 196 : w%v2sigmatau(1, 1), w%v2lapl2(1, 1), w%v2lapltau(1, 1), w%v2tau2(1, 1))
2242 : END IF
2243 : END SELECT
2244 135168 : ELSE IF (is_gga) THEN
2245 17310 : SELECT CASE (grad_deriv)
2246 : CASE (0)
2247 512 : CALL xc_f03_gga_exc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%exc(1))
2248 : CASE (1)
2249 14968 : IF (no_exc) THEN
2250 0 : CALL xc_f03_gga_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), w%vrho(1, 1), w%vsigma(1, 1))
2251 0 : w%exc(1:nb) = 0.0_dp
2252 : ELSE
2253 : CALL xc_f03_gga_exc_vxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2254 14968 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1))
2255 : END IF
2256 : CASE (2)
2257 16798 : IF (no_exc) THEN
2258 : CALL xc_f03_gga_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2259 : w%vrho(1, 1), w%vsigma(1, 1), &
2260 0 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2261 0 : w%exc(1:nb) = 0.0_dp
2262 : ELSE
2263 : CALL xc_f03_gga_exc_vxc_fxc(w%func, np, w%rho(1, 1), w%sigma(1, 1), &
2264 : w%exc(1), w%vrho(1, 1), w%vsigma(1, 1), &
2265 1318 : w%v2rho2(1, 1), w%v2rhosigma(1, 1), w%v2sigma2(1, 1))
2266 : END IF
2267 : END SELECT
2268 : ELSE
2269 118370 : SELECT CASE (grad_deriv)
2270 : CASE (0)
2271 0 : CALL xc_f03_lda_exc(w%func, np, w%rho(1, 1), w%exc(1))
2272 : CASE (1)
2273 116634 : CALL xc_f03_lda_exc_vxc(w%func, np, w%rho(1, 1), w%exc(1), w%vrho(1, 1))
2274 : CASE (2)
2275 : CALL xc_f03_lda_exc_vxc_fxc(w%func, np, w%rho(1, 1), &
2276 1736 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1))
2277 : CASE (3)
2278 : CALL xc_f03_lda(w%func, np, w%rho(1, 1), &
2279 118370 : w%exc(1), w%vrho(1, 1), w%v2rho2(1, 1), w%v3rho3(1, 1))
2280 : END SELECT
2281 : END IF
2282 :
2283 : ! accumulate; one loop per group of derivatives, selected by the order
2284 : ! requested and the family
2285 154032 : IF (grad_deriv >= 0) THEN
2286 78715696 : DO i = 1, nb
2287 78561664 : ii = i0 + i
2288 78715696 : e_0(ii) = e_0(ii) + sc*w%exc(i)*(w%rho(1, i) + w%rho(2, i))
2289 : END DO
2290 : END IF
2291 154032 : IF (grad_deriv >= 1) THEN
2292 78137496 : DO i = 1, nb
2293 77984604 : ii = i0 + i
2294 77984604 : e_rhoa(ii) = e_rhoa(ii) + sc*w%vrho(1, i)
2295 78137496 : e_rhob(ii) = e_rhob(ii) + sc*w%vrho(2, i)
2296 : END DO
2297 152892 : IF (is_gga) THEN
2298 17499538 : DO i = 1, nb
2299 17465016 : ii = i0 + i
2300 17465016 : e_ndrho(ii) = e_ndrho(ii) + sc*w%vsigma(2, i)*w%nd(i)
2301 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2302 17465016 : sc*(2.0_dp*w%vsigma(1, i) - w%vsigma(2, i))*w%nda(i)
2303 : e_ndrhob(ii) = e_ndrhob(ii) + &
2304 17499538 : sc*(2.0_dp*w%vsigma(3, i) - w%vsigma(2, i))*w%ndb(i)
2305 : END DO
2306 : END IF
2307 152892 : IF (is_mgga) THEN
2308 9180368 : DO i = 1, nb
2309 9162132 : ii = i0 + i
2310 9162132 : e_tau_a(ii) = e_tau_a(ii) + sc*w%vtau(1, i)
2311 9180368 : e_tau_b(ii) = e_tau_b(ii) + sc*w%vtau(2, i)
2312 : END DO
2313 : END IF
2314 18236 : IF (is_mgga .AND. has_laplace) THEN
2315 1246680 : DO i = 1, nb
2316 1244160 : ii = i0 + i
2317 1244160 : e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*w%vlapl(1, i)
2318 1246680 : e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*w%vlapl(2, i)
2319 : END DO
2320 : END IF
2321 : END IF
2322 154032 : IF (grad_deriv >= 2) THEN
2323 1636402 : DO i = 1, nb
2324 1633152 : ii = i0 + i
2325 1633152 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*w%v2rho2(1, i)
2326 1633152 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*w%v2rho2(2, i)
2327 1636402 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*w%v2rho2(3, i)
2328 : END DO
2329 3250 : IF (is_gga) THEN
2330 761318 : DO i = 1, nb
2331 759804 : ii = i0 + i
2332 759804 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*w%v2rhosigma(2, i)*w%nd(i)
2333 759804 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*w%v2rhosigma(5, i)*w%nd(i)
2334 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2335 759804 : sc*(2.0_dp*w%v2rhosigma(1, i) - w%v2rhosigma(2, i))*w%nda(i)
2336 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2337 759804 : sc*(2.0_dp*w%v2rhosigma(4, i) - w%v2rhosigma(5, i))*w%nda(i)
2338 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2339 759804 : sc*(2.0_dp*w%v2rhosigma(3, i) - w%v2rhosigma(2, i))*w%ndb(i)
2340 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2341 759804 : sc*(2.0_dp*w%v2rhosigma(6, i) - w%v2rhosigma(5, i))*w%ndb(i)
2342 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2343 759804 : sc*(w%vsigma(2, i) + w%nd(i)**2*w%v2sigma2(4, i))
2344 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2345 759804 : sc*(2.0_dp*w%v2sigma2(2, i) - w%v2sigma2(4, i))*w%nd(i)*w%nda(i)
2346 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2347 759804 : sc*(2.0_dp*w%v2sigma2(5, i) - w%v2sigma2(4, i))*w%nd(i)*w%ndb(i)
2348 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2349 : sc*(2.0_dp*w%vsigma(1, i) - w%vsigma(2, i) + w%nda(i)**2*( &
2350 759804 : 4.0_dp*w%v2sigma2(1, i) - 4.0_dp*w%v2sigma2(2, i) + w%v2sigma2(4, i)))
2351 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2352 : sc*(4.0_dp*w%v2sigma2(3, i) - 2.0_dp*w%v2sigma2(2, i) - &
2353 759804 : 2.0_dp*w%v2sigma2(5, i) + w%v2sigma2(4, i))*w%nda(i)*w%ndb(i)
2354 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2355 : sc*(2.0_dp*w%vsigma(3, i) - w%vsigma(2, i) + w%ndb(i)**2*( &
2356 761318 : 4.0_dp*w%v2sigma2(6, i) - 4.0_dp*w%v2sigma2(5, i) + w%v2sigma2(4, i)))
2357 : END DO
2358 : END IF
2359 3250 : IF (is_mgga) THEN
2360 96964 : DO i = 1, nb
2361 96768 : ii = i0 + i
2362 96768 : e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*w%v2rhotau(1, i)
2363 96768 : e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*w%v2rhotau(2, i)
2364 96768 : e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*w%v2rhotau(3, i)
2365 96768 : e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*w%v2rhotau(4, i)
2366 96768 : e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*w%v2sigmatau(3, i)*w%nd(i)
2367 96768 : e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*w%v2sigmatau(4, i)*w%nd(i)
2368 : e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2369 96768 : sc*(2.0_dp*w%v2sigmatau(1, i) - w%v2sigmatau(3, i))*w%nda(i)
2370 : e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2371 96768 : sc*(2.0_dp*w%v2sigmatau(2, i) - w%v2sigmatau(4, i))*w%nda(i)
2372 : e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2373 96768 : sc*(2.0_dp*w%v2sigmatau(5, i) - w%v2sigmatau(3, i))*w%ndb(i)
2374 : e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2375 96768 : sc*(2.0_dp*w%v2sigmatau(6, i) - w%v2sigmatau(4, i))*w%ndb(i)
2376 96768 : e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*w%v2tau2(1, i)
2377 96768 : e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*w%v2tau2(2, i)
2378 96964 : e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*w%v2tau2(3, i)
2379 : END DO
2380 : END IF
2381 196 : IF (is_mgga .AND. has_laplace) THEN
2382 41556 : DO i = 1, nb
2383 41472 : ii = i0 + i
2384 41472 : e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*w%v2rholapl(1, i)
2385 41472 : e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*w%v2rholapl(2, i)
2386 41472 : e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*w%v2rholapl(3, i)
2387 41472 : e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*w%v2rholapl(4, i)
2388 41472 : e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*w%v2sigmalapl(3, i)*w%nd(i)
2389 41472 : e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*w%v2sigmalapl(4, i)*w%nd(i)
2390 : e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2391 41472 : sc*(2.0_dp*w%v2sigmalapl(1, i) - w%v2sigmalapl(3, i))*w%nda(i)
2392 : e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2393 41472 : sc*(2.0_dp*w%v2sigmalapl(2, i) - w%v2sigmalapl(4, i))*w%nda(i)
2394 : e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2395 41472 : sc*(2.0_dp*w%v2sigmalapl(5, i) - w%v2sigmalapl(3, i))*w%ndb(i)
2396 : e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2397 41472 : sc*(2.0_dp*w%v2sigmalapl(6, i) - w%v2sigmalapl(4, i))*w%ndb(i)
2398 41472 : e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*w%v2lapl2(1, i)
2399 41472 : e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*w%v2lapl2(2, i)
2400 41472 : e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*w%v2lapl2(3, i)
2401 41472 : e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*w%v2lapltau(1, i)
2402 41472 : e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*w%v2lapltau(2, i)
2403 41472 : e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*w%v2lapltau(3, i)
2404 41556 : e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*w%v2lapltau(4, i)
2405 : END DO
2406 : END IF
2407 : END IF
2408 154032 : IF (grad_deriv >= 3) THEN
2409 0 : DO i = 1, nb
2410 0 : ii = i0 + i
2411 0 : e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*w%v3rho3(1, i)
2412 0 : e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*w%v3rho3(2, i)
2413 0 : e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*w%v3rho3(3, i)
2414 0 : e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*w%v3rho3(4, i)
2415 : END DO
2416 : END IF
2417 : END DO
2418 : !$OMP END DO
2419 : END ASSOCIATE
2420 :
2421 2750 : END SUBROUTINE libxc_spin_polarized_calc
2422 : #endif
2423 :
2424 0 : END MODULE xc_libxc
|
