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_functional_get_number, &
108 : xc_libxc_wrap_needs_laplace, &
109 : xc_libxc_wrap_functional_set_params, &
110 : xc_libxc_wrap_is_under_development, &
111 : xc_libxc_get_reference_length, &
112 : xc_libxc_check_functional
113 : #endif
114 :
115 : #include "../base/base_uses.f90"
116 :
117 : IMPLICIT NONE
118 : PRIVATE
119 :
120 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_libxc'
121 :
122 : PUBLIC :: libxc_spin_unpolarized_info, libxc_spin_unpolarized_eval, &
123 : libxc_spin_polarized_info, libxc_spin_polarized_eval, &
124 : libxc_version_info, libxc_get_reference_length, libxc_add_sections, &
125 : libxc_check_existence_in_libxc
126 :
127 : #if defined (__LIBXC)
128 : INTEGER(C_SIZE_T), PARAMETER, PRIVATE :: one = 1
129 : #endif
130 :
131 : CONTAINS
132 :
133 : ! **************************************************************************************************
134 : !> \brief This function checks whether a functional name belongs to LibXC
135 : !> \param libxc_params (possible) LibXC input section
136 : !> \return exists whether the functional exists in LibXC
137 : ! **************************************************************************************************
138 2255 : FUNCTION libxc_check_existence_in_libxc(libxc_params) RESULT(exists)
139 :
140 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
141 : LOGICAL :: exists
142 :
143 : #if defined (__LIBXC)
144 :
145 2255 : exists = xc_libxc_check_functional(libxc_params%section%name)
146 : #else
147 : MARK_USED(libxc_params)
148 : exists = .FALSE.
149 : #endif
150 :
151 2255 : END FUNCTION libxc_check_existence_in_libxc
152 :
153 : ! **************************************************************************************************
154 : !> \brief This function returns the maximum length of the reference string for a given LibXC functional
155 : !> \param libxc_params LibXC input section
156 : !> \param lsd spin polarized calculation
157 : !> \return maximum length of the string
158 : ! **************************************************************************************************
159 98 : FUNCTION libxc_get_reference_length(libxc_params, lsd) RESULT(length)
160 :
161 : TYPE(section_vals_type), POINTER, INTENT(IN) :: libxc_params
162 : LOGICAL, INTENT(IN) :: lsd
163 : INTEGER :: length
164 :
165 : #if defined (__LIBXC)
166 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_get_reference_length'
167 :
168 : CHARACTER(LEN=default_string_length) :: func_name
169 : INTEGER :: func_id, handle
170 : TYPE(xc_f03_func_t) :: xc_func
171 : TYPE(xc_f03_func_info_t) :: xc_info
172 :
173 98 : CALL timeset(routineN, handle)
174 :
175 98 : func_name = libxc_params%section%name
176 :
177 98 : func_id = xc_libxc_wrap_functional_get_number(func_name)
178 196 : !$OMP CRITICAL(libxc_init)
179 98 : IF (lsd) THEN
180 46 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
181 : ELSE
182 52 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
183 : END IF
184 98 : xc_info = xc_f03_func_get_info(xc_func)
185 : !$OMP END CRITICAL(libxc_init)
186 98 : !$OMP BARRIER
187 :
188 98 : length = xc_libxc_get_reference_length(xc_info)
189 :
190 98 : CALL xc_f03_func_end(xc_func)
191 :
192 98 : CALL timestop(handle)
193 : #else
194 : MARK_USED(libxc_params)
195 : MARK_USED(lsd)
196 : length = 0
197 : CPABORT("In order to use LibXC you have to download and install it!")
198 : #endif
199 :
200 98 : END FUNCTION libxc_get_reference_length
201 :
202 : ! **************************************************************************************************
203 : !> \brief ...
204 : !> \param section ...
205 : ! **************************************************************************************************
206 106335 : SUBROUTINE libxc_add_sections(section)
207 :
208 : TYPE(section_type), POINTER, INTENT(IN) :: section
209 :
210 : #if defined (__LIBXC)
211 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_add_sections'
212 :
213 : TYPE(section_type), POINTER :: subsection
214 : TYPE(keyword_type), POINTER :: keyword
215 : INTEGER :: handle, no_func, len_name, ii, func_id, n_param, iparam
216 : REAL(KIND=C_DOUBLE) :: default_val
217 : CHARACTER(LEN=128) :: func_name, param_name, param_descr, description
218 : CHARACTER(LEN=2*default_string_length) :: warning
219 106335 : INTEGER(KIND=C_INT), DIMENSION(:), ALLOCATABLE :: func_ids
220 : TYPE(xc_f03_func_t) :: xc_func
221 : TYPE(xc_f03_func_info_t) :: xc_info
222 :
223 106335 : CALL timeset(routineN, handle)
224 :
225 106335 : CPASSERT(ASSOCIATED(section))
226 106335 : NULLIFY (subsection, keyword)
227 :
228 106335 : no_func = xc_f03_number_of_functionals()
229 106335 : len_name = xc_f03_maximum_name_length()
230 :
231 319005 : ALLOCATE (func_ids(no_func))
232 :
233 106335 : CALL xc_f03_available_functional_numbers(func_ids)
234 :
235 75816855 : DO ii = 1, no_func
236 :
237 75710520 : func_id = func_ids(ii)
238 75710520 : IF (ii > 1) THEN
239 75604185 : IF (func_id == func_ids(ii - 1)) CYCLE
240 : END IF
241 149294340 : !$OMP CRITICAL(libxc_init)
242 74647170 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
243 74647170 : xc_info = xc_f03_func_get_info(xc_func)
244 : !$OMP END CRITICAL(libxc_init)
245 74647170 : !$OMP BARRIER
246 :
247 74647170 : func_name = xc_f03_functional_get_name(func_id)
248 74647170 : description = xc_f03_func_info_get_name(xc_info)
249 74647170 : n_param = xc_f03_func_info_get_n_ext_params(xc_info)
250 :
251 74647170 : NULLIFY (subsection)
252 : CALL section_create(subsection, __LOCATION__, name=TRIM(func_name), description=TRIM(description), &
253 74647170 : n_keywords=2 + n_param, n_subsections=0, repeats=.FALSE.)
254 :
255 74647170 : IF (description(1:1) == "_") THEN
256 : warning = " This parameter is an internal parameter of the functional. Changing this "// &
257 0 : "parameter effectively changes the functional."
258 : ELSE
259 74647170 : warning = " "
260 : END IF
261 :
262 74647170 : NULLIFY (keyword)
263 : CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
264 : description="Activates the functional."//TRIM(warning), &
265 74647170 : lone_keyword_l_val=.TRUE., default_l_val=.FALSE.)
266 74647170 : CALL section_add_keyword(subsection, keyword)
267 74647170 : CALL keyword_release(keyword)
268 :
269 : CALL keyword_create(keyword, __LOCATION__, name="SCALE", description="Scales this functional", &
270 74647170 : default_r_val=1.0_dp)
271 74647170 : CALL section_add_keyword(subsection, keyword)
272 74647170 : CALL keyword_release(keyword)
273 :
274 454794795 : DO iparam = 1, n_param
275 380147625 : param_name = xc_f03_func_info_get_ext_params_name(xc_info, iparam - 1)
276 380147625 : param_descr = xc_f03_func_info_get_ext_params_description(xc_info, iparam - 1)
277 380147625 : default_val = xc_f03_func_info_get_ext_params_default_value(xc_info, iparam - 1)
278 380147625 : NULLIFY (keyword)
279 : CALL keyword_create(keyword, __LOCATION__, name=TRIM(param_name), &
280 380147625 : description=TRIM(param_descr), default_r_val=default_val)
281 380147625 : CALL section_add_keyword(subsection, keyword)
282 454794795 : CALL keyword_release(keyword)
283 : END DO
284 :
285 74647170 : CALL section_add_subsection(section, subsection)
286 74647170 : CALL section_release(subsection)
287 :
288 75816855 : CALL xc_f03_func_end(xc_func)
289 :
290 : END DO
291 :
292 106335 : DEALLOCATE (func_ids)
293 :
294 106335 : CALL timestop(handle)
295 : #else
296 : MARK_USED(section)
297 :
298 : #endif
299 :
300 106335 : END SUBROUTINE libxc_add_sections
301 :
302 : ! **************************************************************************************************
303 : !> \brief info about the functional from libxc
304 : !> \param libxc_params input parameter (functional name, scaling and parameters)
305 : !> \param reference string with the reference of the actual functional
306 : !> \param shortform string with the shortform of the functional name
307 : !> \param needs the components needed by this functional are set to
308 : !> true (does not set the unneeded components to false)
309 : !> \param max_deriv maximum implemented derivative of the xc functional
310 : !> \param print_warn whether to print warning about development status of a functional
311 : !> \param func_name_override optional LibXC functional name overriding the section name
312 : !> \author F. Tran
313 : ! **************************************************************************************************
314 13666 : SUBROUTINE libxc_spin_unpolarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
315 : func_name_override)
316 :
317 : TYPE(section_vals_type), POINTER :: libxc_params
318 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
319 : TYPE(xc_rho_cflags_type), &
320 : INTENT(inout), OPTIONAL :: needs
321 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
322 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
323 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
324 :
325 : #if defined (__LIBXC)
326 : CHARACTER(LEN=128) :: s1, s2
327 : CHARACTER(LEN=default_string_length) :: func_name
328 : INTEGER :: func_id
329 : REAL(KIND=dp) :: func_scale
330 : TYPE(xc_f03_func_t) :: xc_func
331 : TYPE(xc_f03_func_info_t) :: xc_info
332 :
333 27244 : IF (PRESENT(func_name_override)) THEN
334 88 : func_name = func_name_override
335 88 : func_scale = 1.0_dp
336 : ELSE
337 13578 : func_name = libxc_params%section%name
338 13578 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
339 : END IF
340 :
341 13666 : CALL cite_reference(Marques2012)
342 13666 : CALL cite_reference(Lehtola2018)
343 :
344 13666 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
345 :
346 13666 : func_id = xc_libxc_wrap_functional_get_number(func_name)
347 27332 : !$OMP CRITICAL(libxc_init)
348 13666 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
349 13666 : xc_info = xc_f03_func_get_info(xc_func)
350 : !$OMP END CRITICAL(libxc_init)
351 13666 : !$OMP BARRIER
352 :
353 13666 : s1 = xc_f03_func_info_get_name(xc_info)
354 9585 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
355 9585 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
356 2331 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
357 1296 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
358 454 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
359 : CASE default
360 13666 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
361 : END SELECT
362 13666 : IF (PRESENT(shortform)) THEN
363 52 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
364 : END IF
365 13666 : IF (PRESENT(reference)) THEN
366 52 : CALL xc_libxc_wrap_info_refs(xc_info, XC_UNPOLARIZED, func_scale, reference)
367 : END IF
368 13666 : IF (PRESENT(needs)) THEN
369 6100 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
370 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
371 6100 : needs%rho = .TRUE.
372 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
373 4544 : needs%rho = .TRUE.
374 4544 : needs%norm_drho = .TRUE.
375 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
376 2962 : needs%rho = .TRUE.
377 2962 : needs%norm_drho = .TRUE.
378 2962 : needs%tau = .TRUE.
379 2962 : needs%laplace_rho = xc_libxc_wrap_needs_laplace(func_id)
380 : CASE default
381 13606 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
382 : END SELECT
383 : END IF
384 13666 : IF (PRESENT(max_deriv)) THEN
385 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
386 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
387 0 : max_deriv = 3
388 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
389 88 : max_deriv = 2
390 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
391 0 : max_deriv = 2
392 : CASE default
393 88 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
394 : END SELECT
395 : END IF
396 13666 : IF (PRESENT(print_warn)) THEN
397 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
398 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
399 : END IF
400 : END IF
401 :
402 13666 : CALL xc_f03_func_end(xc_func)
403 : #else
404 : MARK_USED(libxc_params)
405 : MARK_USED(reference)
406 : MARK_USED(shortform)
407 : MARK_USED(needs)
408 : MARK_USED(max_deriv)
409 : MARK_USED(print_warn)
410 : MARK_USED(func_name_override)
411 :
412 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
413 : "for a functional of the LibXC library, "// &
414 : "you have to download and install the library!")
415 : #endif
416 :
417 13666 : END SUBROUTINE libxc_spin_unpolarized_info
418 :
419 : ! **************************************************************************************************
420 : !> \brief info about the functional from libxc
421 : !> \param libxc_params input parameter (functional name, scaling and parameters)
422 : !> \param reference string with the reference of the actual functional
423 : !> \param shortform string with the shortform of the functional name
424 : !> \param needs the components needed by this functional are set to
425 : !> true (does not set the unneeded components to false)
426 : !> \param max_deriv maximum implemented derivative of the xc functional
427 : !> \param print_warn whether to print warning about development status of a functional
428 : !> \param func_name_override optional LibXC functional name overriding the section name
429 : !> \author F. Tran
430 : ! **************************************************************************************************
431 2476 : SUBROUTINE libxc_spin_polarized_info(libxc_params, reference, shortform, needs, max_deriv, print_warn, &
432 : func_name_override)
433 :
434 : TYPE(section_vals_type), POINTER :: libxc_params
435 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
436 : TYPE(xc_rho_cflags_type), &
437 : INTENT(inout), OPTIONAL :: needs
438 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
439 : LOGICAL, INTENT(IN), OPTIONAL :: print_warn
440 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
441 :
442 : #if defined (__LIBXC)
443 : CHARACTER(LEN=128) :: s1, s2
444 : CHARACTER(LEN=default_string_length) :: func_name
445 : INTEGER :: func_id
446 : REAL(KIND=dp) :: func_scale
447 : TYPE(xc_f03_func_t) :: xc_func
448 : TYPE(xc_f03_func_info_t) :: xc_info
449 :
450 4944 : IF (PRESENT(func_name_override)) THEN
451 8 : func_name = func_name_override
452 8 : func_scale = 1.0_dp
453 : ELSE
454 2468 : func_name = libxc_params%section%name
455 2468 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
456 : END IF
457 :
458 2476 : CALL cite_reference(Marques2012)
459 2476 : CALL cite_reference(Lehtola2018)
460 :
461 2476 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
462 :
463 2476 : func_id = xc_libxc_wrap_functional_get_number(func_name)
464 4952 : !$OMP CRITICAL(libxc_init)
465 2476 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
466 2476 : xc_info = xc_f03_func_get_info(xc_func)
467 : !$OMP END CRITICAL(libxc_init)
468 2476 : !$OMP BARRIER
469 :
470 2476 : s1 = xc_f03_func_info_get_name(xc_info)
471 1158 : SELECT CASE (xc_f03_func_info_get_kind(xc_info))
472 1158 : CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange"
473 1072 : CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation"
474 246 : CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation"
475 0 : CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic"
476 : CASE default
477 2476 : CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.")
478 : END SELECT
479 2476 : IF (PRESENT(shortform)) THEN
480 46 : shortform = TRIM(s1)//' ('//TRIM(s2)//')'
481 : END IF
482 2476 : IF (PRESENT(reference)) THEN
483 46 : CALL xc_libxc_wrap_info_refs(xc_info, XC_POLARIZED, func_scale, reference)
484 : END IF
485 2476 : IF (PRESENT(needs)) THEN
486 1016 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
487 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
488 1016 : needs%rho_spin = .TRUE.
489 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
490 262 : needs%rho_spin = .TRUE.
491 262 : needs%norm_drho = .TRUE.
492 262 : needs%norm_drho_spin = .TRUE.
493 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
494 1152 : needs%rho_spin = .TRUE.
495 1152 : needs%norm_drho = .TRUE.
496 1152 : needs%norm_drho_spin = .TRUE.
497 1152 : needs%tau_spin = .TRUE.
498 1152 : needs%laplace_rho_spin = xc_libxc_wrap_needs_laplace(func_id)
499 : CASE default
500 2430 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
501 : END SELECT
502 : END IF
503 2476 : IF (PRESENT(max_deriv)) THEN
504 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
505 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
506 0 : max_deriv = 3
507 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
508 8 : max_deriv = 2
509 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
510 0 : max_deriv = 2
511 : CASE default
512 8 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
513 : END SELECT
514 : END IF
515 2476 : IF (PRESENT(print_warn)) THEN
516 0 : IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN
517 0 : CPWARN(TRIM(func_name)//" is under development. Use with caution.")
518 : END IF
519 : END IF
520 :
521 2476 : CALL xc_f03_func_end(xc_func)
522 : #else
523 : MARK_USED(libxc_params)
524 : MARK_USED(reference)
525 : MARK_USED(shortform)
526 : MARK_USED(needs)
527 : MARK_USED(max_deriv)
528 : MARK_USED(print_warn)
529 : MARK_USED(func_name_override)
530 :
531 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are "// &
532 : "asking for a functional of the LibXC library, "// &
533 : "you have to download and install the library!")
534 : #endif
535 :
536 2476 : END SUBROUTINE libxc_spin_polarized_info
537 :
538 : ! **************************************************************************************************
539 : !> \brief info about the LibXC version
540 : !> \param version ...
541 : !> \author A. Gloess (agloess)
542 : ! **************************************************************************************************
543 0 : SUBROUTINE libxc_version_info(version)
544 : CHARACTER(LEN=*), INTENT(OUT) :: version ! the string that is output
545 :
546 : #if defined (__LIBXC)
547 0 : CALL xc_libxc_wrap_version(version)
548 : #else
549 : version = "none"
550 : CPABORT("In order to use libxc you need to download and install it")
551 : #endif
552 :
553 0 : END SUBROUTINE libxc_version_info
554 :
555 : ! **************************************************************************************************
556 : !> \brief evaluates the functional from libxc
557 : !> \param rho_set the density where you want to evaluate the functional
558 : !> \param deriv_set place where to store the functional derivatives (they are
559 : !> added to the derivatives)
560 : !> \param grad_deriv degree of the derivative that should be evaluated,
561 : !> if positive all the derivatives up to the given degree are evaluated,
562 : !> if negative only the given degree is calculated
563 : !> \param libxc_params input parameter (functional name, scaling and parameters)
564 : !> \param func_name_override optional LibXC functional name overriding the section name
565 : !> \author F. Tran
566 : ! **************************************************************************************************
567 18110 : SUBROUTINE libxc_spin_unpolarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
568 :
569 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
570 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
571 : INTEGER, INTENT(in) :: grad_deriv
572 : TYPE(section_vals_type), POINTER :: libxc_params
573 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
574 :
575 : #if defined (__LIBXC)
576 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_unpolarized_eval'
577 :
578 : CHARACTER(LEN=default_string_length) :: func_name
579 : INTEGER :: func_id, handle, npoints
580 : INTEGER, DIMENSION(2, 3) :: bo
581 : LOGICAL :: has_laplace, no_exc
582 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
583 18110 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rho, &
584 18110 : e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_ndrho, &
585 18110 : e_ndrho_laplace_rho, e_ndrho_ndrho, e_ndrho_rho, e_ndrho_tau, e_rho, &
586 18110 : e_rho_laplace_rho, e_rho_rho, e_rho_rho_rho, e_rho_tau, e_tau, &
587 18110 : e_tau_tau, laplace_rho, norm_drho, rho, tau
588 : TYPE(xc_derivative_type), POINTER :: deriv
589 : TYPE(xc_f03_func_t) :: xc_func
590 : TYPE(xc_f03_func_info_t) :: xc_info
591 :
592 18110 : CALL timeset(routineN, handle)
593 :
594 18110 : has_laplace = .FALSE.
595 18110 : NULLIFY (dummy)
596 18110 : NULLIFY (rho, norm_drho, laplace_rho, tau)
597 :
598 18110 : IF (PRESENT(func_name_override)) THEN
599 0 : func_name = func_name_override
600 0 : func_scale = 1.0_dp
601 : ELSE
602 18110 : func_name = libxc_params%section%name
603 18110 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
604 : END IF
605 :
606 18110 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
607 :
608 18110 : func_id = xc_libxc_wrap_functional_get_number(func_name)
609 18110 : CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
610 18110 : xc_info = xc_f03_func_get_info(xc_func)
611 18110 : no_exc = .FALSE.
612 18110 : IF (.NOT. PRESENT(func_name_override)) THEN
613 18110 : CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
614 : END IF
615 :
616 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
617 : rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, &
618 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
619 18110 : tau=tau, local_bounds=bo)
620 :
621 18110 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
622 :
623 18110 : dummy => rho
624 :
625 : ! due to assumed shape array usage in next routine
626 18110 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => dummy
627 18110 : IF (.NOT. ASSOCIATED(tau)) tau => dummy
628 :
629 : ! only some MGGA functionals really need the Laplacian,
630 : ! all others can work with rho (read-only) as dummy
631 18110 : has_laplace = xc_libxc_wrap_needs_laplace(func_id)
632 18110 : IF (.NOT. has_laplace) laplace_rho => dummy
633 :
634 18110 : e_0 => dummy
635 18110 : e_rho => dummy
636 18110 : e_ndrho => dummy
637 18110 : e_laplace_rho => dummy
638 18110 : e_tau => dummy
639 18110 : e_rho_rho => dummy
640 18110 : e_ndrho_rho => dummy
641 18110 : e_ndrho_ndrho => dummy
642 18110 : e_rho_laplace_rho => dummy
643 18110 : e_rho_tau => dummy
644 18110 : e_ndrho_laplace_rho => dummy
645 18110 : e_ndrho_tau => dummy
646 18110 : e_laplace_rho_laplace_rho => dummy
647 18110 : e_laplace_rho_tau => dummy
648 18110 : e_tau_tau => dummy
649 18110 : e_rho_rho_rho => dummy
650 :
651 18110 : IF (grad_deriv >= 0) THEN
652 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
653 18110 : allocate_deriv=.TRUE.)
654 18110 : CALL xc_derivative_get(deriv, deriv_data=e_0)
655 : END IF
656 18110 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
657 10306 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
658 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
659 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
660 10306 : allocate_deriv=.TRUE.)
661 10306 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
662 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
663 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
664 5340 : allocate_deriv=.TRUE.)
665 5340 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
666 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
667 5340 : allocate_deriv=.TRUE.)
668 5340 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
669 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
670 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
671 2210 : allocate_deriv=.TRUE.)
672 2210 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
673 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
674 2210 : allocate_deriv=.TRUE.)
675 2210 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
676 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
677 2210 : allocate_deriv=.TRUE.)
678 2210 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
679 2210 : IF (has_laplace) THEN
680 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
681 568 : allocate_deriv=.TRUE.)
682 568 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
683 : END IF
684 : CASE default
685 17856 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
686 : END SELECT
687 : END IF
688 18110 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
689 1584 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
690 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
691 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
692 1584 : allocate_deriv=.TRUE.)
693 1584 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
694 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
695 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
696 700 : allocate_deriv=.TRUE.)
697 700 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
698 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
699 700 : allocate_deriv=.TRUE.)
700 700 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
701 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
702 700 : allocate_deriv=.TRUE.)
703 700 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
704 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
705 : ! not implemented ...
706 :
707 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
708 308 : allocate_deriv=.TRUE.)
709 308 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
710 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
711 308 : allocate_deriv=.TRUE.)
712 308 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
713 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
714 308 : allocate_deriv=.TRUE.)
715 308 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
716 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_tau], &
717 308 : allocate_deriv=.TRUE.)
718 308 : CALL xc_derivative_get(deriv, deriv_data=e_rho_tau)
719 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau], &
720 308 : allocate_deriv=.TRUE.)
721 308 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau)
722 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau, deriv_tau], &
723 308 : allocate_deriv=.TRUE.)
724 308 : CALL xc_derivative_get(deriv, deriv_data=e_tau_tau)
725 308 : IF (has_laplace) THEN
726 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_laplace_rho], &
727 108 : allocate_deriv=.TRUE.)
728 108 : CALL xc_derivative_get(deriv, deriv_data=e_rho_laplace_rho)
729 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rho], &
730 108 : allocate_deriv=.TRUE.)
731 108 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rho)
732 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_laplace_rho], &
733 108 : allocate_deriv=.TRUE.)
734 108 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_laplace_rho)
735 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho, deriv_tau], &
736 108 : allocate_deriv=.TRUE.)
737 108 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho_tau)
738 : END IF
739 : CASE default
740 2592 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
741 : END SELECT
742 : END IF
743 18110 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
744 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
745 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
746 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
747 0 : allocate_deriv=.TRUE.)
748 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
749 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
750 0 : CPABORT("derivatives larger than 2 not implemented")
751 : CASE default
752 0 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
753 : END SELECT
754 : END IF
755 18110 : IF (grad_deriv >= 4 .OR. grad_deriv <= -4) THEN
756 0 : CPABORT("derivatives larger than 3 not implemented")
757 : END IF
758 :
759 : !$OMP PARALLEL DEFAULT(NONE), &
760 : !$OMP SHARED(rho,norm_drho,laplace_rho,tau,e_0,e_rho,e_ndrho,e_laplace_rho),&
761 : !$OMP SHARED(e_tau,e_rho_rho,e_ndrho_rho,e_ndrho_ndrho,e_rho_laplace_rho),&
762 : !$OMP SHARED(e_rho_tau,e_ndrho_laplace_rho,e_ndrho_tau,e_laplace_rho_laplace_rho),&
763 : !$OMP SHARED(e_laplace_rho_tau,e_tau_tau,e_rho_rho_rho),&
764 : !$OMP SHARED(grad_deriv,npoints),&
765 : !$OMP SHARED(epsilon_rho,epsilon_tau),&
766 18110 : !$OMP SHARED(func_name,func_scale,xc_func,xc_info,no_exc,has_laplace)
767 :
768 : CALL libxc_spin_unpolarized_calc(rho=rho, norm_drho=norm_drho, &
769 : laplace_rho=laplace_rho, tau=tau, &
770 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_laplace_rho=e_laplace_rho, &
771 : e_tau=e_tau, e_rho_rho=e_rho_rho, e_ndrho_rho=e_ndrho_rho, &
772 : e_ndrho_ndrho=e_ndrho_ndrho, e_rho_laplace_rho=e_rho_laplace_rho, &
773 : e_rho_tau=e_rho_tau, e_ndrho_laplace_rho=e_ndrho_laplace_rho, &
774 : e_ndrho_tau=e_ndrho_tau, e_laplace_rho_laplace_rho=e_laplace_rho_laplace_rho, &
775 : e_laplace_rho_tau=e_laplace_rho_tau, e_tau_tau=e_tau_tau, &
776 : e_rho_rho_rho=e_rho_rho_rho, &
777 : grad_deriv=grad_deriv, npoints=npoints, &
778 : epsilon_rho=epsilon_rho, &
779 : epsilon_tau=epsilon_tau, func_name=func_name, &
780 : sc=func_scale, xc_func=xc_func, xc_info=xc_info, no_exc=no_exc, has_laplace=has_laplace)
781 :
782 : !$OMP END PARALLEL
783 :
784 18110 : NULLIFY (dummy)
785 :
786 18110 : CALL xc_f03_func_end(xc_func)
787 :
788 18110 : CALL timestop(handle)
789 : #else
790 : MARK_USED(rho_set)
791 : MARK_USED(deriv_set)
792 : MARK_USED(grad_deriv)
793 : MARK_USED(libxc_params)
794 : MARK_USED(func_name_override)
795 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
796 : "for a functional of the LibXC library, "// &
797 : "you have to download and install the library!")
798 : #endif
799 18110 : END SUBROUTINE libxc_spin_unpolarized_eval
800 :
801 : ! **************************************************************************************************
802 : !> \brief evaluates the functional from libxc
803 : !> \param rho_set the density where you want to evaluate the functional
804 : !> \param deriv_set place where to store the functional derivatives (they are
805 : !> added to the derivatives)
806 : !> \param grad_deriv degree of the derivative that should be evaluated,
807 : !> if positive all the derivatives up to the given degree are evaluated,
808 : !> if negative only the given degree is calculated
809 : !> \param libxc_params input parameter (functional name, scaling and parameters)
810 : !> \param func_name_override optional LibXC functional name overriding the section name
811 : !> \author F. Tran
812 : ! **************************************************************************************************
813 2750 : SUBROUTINE libxc_spin_polarized_eval(rho_set, deriv_set, grad_deriv, libxc_params, func_name_override)
814 :
815 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
816 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
817 : INTEGER, INTENT(in) :: grad_deriv
818 : TYPE(section_vals_type), POINTER :: libxc_params
819 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: func_name_override
820 :
821 : #if defined (__LIBXC)
822 : CHARACTER(len=*), PARAMETER :: routineN = 'libxc_spin_polarized_eval'
823 :
824 : CHARACTER(LEN=default_string_length) :: func_name
825 : INTEGER :: func_id, handle, npoints
826 : INTEGER, DIMENSION(2, 3) :: bo
827 : LOGICAL :: has_laplace, no_exc
828 : REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale
829 2750 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
830 2750 : e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
831 2750 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, e_laplace_rhob, &
832 2750 : e_laplace_rhob_laplace_rhob, e_laplace_rhob_tau_a, &
833 2750 : e_laplace_rhob_tau_b, e_ndrho, e_ndrho_laplace_rhoa, &
834 2750 : e_ndrho_laplace_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
835 2750 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa, &
836 2750 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, e_ndrhoa_ndrhoa, &
837 2750 : e_ndrhoa_ndrhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhoa_tau_a, &
838 2750 : e_ndrhoa_tau_b, e_ndrhob
839 2750 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_ndrhob_laplace_rhoa, &
840 2750 : e_ndrhob_laplace_rhob, e_ndrhob_ndrhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
841 2750 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_rhoa, e_rhoa_laplace_rhoa, &
842 2750 : e_rhoa_laplace_rhob, e_rhoa_rhoa, e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
843 2750 : e_rhoa_rhob, e_rhoa_rhob_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob, &
844 2750 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhob_rhob, &
845 2750 : e_rhob_rhob_rhob, e_rhob_tau_a, e_rhob_tau_b, e_tau_a, e_tau_a_tau_a, &
846 2750 : e_tau_a_tau_b, e_tau_b, e_tau_b_tau_b, laplace_rhoa, laplace_rhob, &
847 2750 : norm_drho, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
848 : TYPE(xc_derivative_type), POINTER :: deriv
849 : TYPE(xc_f03_func_t) :: xc_func
850 : TYPE(xc_f03_func_info_t) :: xc_info
851 :
852 2750 : CALL timeset(routineN, handle)
853 :
854 2750 : NULLIFY (dummy)
855 2750 : NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, &
856 2750 : laplace_rhob, tau_a, tau_b)
857 :
858 2750 : IF (PRESENT(func_name_override)) THEN
859 0 : func_name = func_name_override
860 0 : func_scale = 1.0_dp
861 : ELSE
862 2750 : func_name = libxc_params%section%name
863 2750 : CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale)
864 : END IF
865 :
866 2750 : IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
867 :
868 2750 : func_id = xc_libxc_wrap_functional_get_number(func_name)
869 2750 : CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED)
870 2750 : xc_info = xc_f03_func_get_info(xc_func)
871 2750 : no_exc = .FALSE.
872 2750 : IF (.NOT. PRESENT(func_name_override)) THEN
873 2750 : CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
874 : END IF
875 :
876 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., &
877 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
878 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
879 : laplace_rhoa=laplace_rhoa, laplace_rhob=laplace_rhob, &
880 : rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
881 2750 : tau_a=tau_a, tau_b=tau_b, local_bounds=bo)
882 :
883 2750 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
884 :
885 2750 : dummy => rhoa
886 :
887 : ! due to assumed shape array usage in next routine
888 2750 : IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => dummy
889 2750 : IF (.NOT. ASSOCIATED(norm_drhoa)) norm_drhoa => dummy
890 2750 : IF (.NOT. ASSOCIATED(norm_drhob)) norm_drhob => dummy
891 2750 : IF (.NOT. ASSOCIATED(tau_a)) tau_a => dummy
892 2750 : IF (.NOT. ASSOCIATED(tau_b)) tau_b => dummy
893 :
894 : ! only some MGGA functionals really need the Laplacian,
895 : ! all others can work with rhoa (read-only) as dummy
896 2750 : has_laplace = xc_libxc_wrap_needs_laplace(func_id)
897 2750 : IF (.NOT. has_laplace) laplace_rhoa => dummy
898 2750 : IF (.NOT. has_laplace) laplace_rhob => dummy
899 :
900 2750 : e_0 => dummy
901 2750 : e_rhoa => dummy
902 2750 : e_rhob => dummy
903 2750 : e_ndrho => dummy
904 2750 : e_ndrhoa => dummy
905 2750 : e_ndrhob => dummy
906 2750 : e_laplace_rhoa => dummy
907 2750 : e_laplace_rhob => dummy
908 2750 : e_tau_a => dummy
909 2750 : e_tau_b => dummy
910 2750 : e_rhoa_rhoa => dummy
911 2750 : e_rhoa_rhob => dummy
912 2750 : e_rhob_rhob => dummy
913 2750 : e_ndrho_rhoa => dummy
914 2750 : e_ndrho_rhob => dummy
915 2750 : e_ndrhoa_rhoa => dummy
916 2750 : e_ndrhoa_rhob => dummy
917 2750 : e_ndrhob_rhoa => dummy
918 2750 : e_ndrhob_rhob => dummy
919 2750 : e_ndrho_ndrho => dummy
920 2750 : e_ndrho_ndrhoa => dummy
921 2750 : e_ndrho_ndrhob => dummy
922 2750 : e_ndrhoa_ndrhoa => dummy
923 2750 : e_ndrhoa_ndrhob => dummy
924 2750 : e_ndrhob_ndrhob => dummy
925 2750 : e_rhoa_laplace_rhoa => dummy
926 2750 : e_rhoa_laplace_rhob => dummy
927 2750 : e_rhob_laplace_rhoa => dummy
928 2750 : e_rhob_laplace_rhob => dummy
929 2750 : e_rhoa_tau_a => dummy
930 2750 : e_rhoa_tau_b => dummy
931 2750 : e_rhob_tau_a => dummy
932 2750 : e_rhob_tau_b => dummy
933 2750 : e_ndrho_laplace_rhoa => dummy
934 2750 : e_ndrho_laplace_rhob => dummy
935 2750 : e_ndrhoa_laplace_rhoa => dummy
936 2750 : e_ndrhoa_laplace_rhob => dummy
937 2750 : e_ndrhob_laplace_rhoa => dummy
938 2750 : e_ndrhob_laplace_rhob => dummy
939 2750 : e_ndrho_tau_a => dummy
940 2750 : e_ndrho_tau_b => dummy
941 2750 : e_ndrhoa_tau_a => dummy
942 2750 : e_ndrhoa_tau_b => dummy
943 2750 : e_ndrhob_tau_a => dummy
944 2750 : e_ndrhob_tau_b => dummy
945 2750 : e_laplace_rhoa_laplace_rhoa => dummy
946 2750 : e_laplace_rhoa_laplace_rhob => dummy
947 2750 : e_laplace_rhob_laplace_rhob => dummy
948 2750 : e_laplace_rhoa_tau_a => dummy
949 2750 : e_laplace_rhoa_tau_b => dummy
950 2750 : e_laplace_rhob_tau_a => dummy
951 2750 : e_laplace_rhob_tau_b => dummy
952 2750 : e_tau_a_tau_a => dummy
953 2750 : e_tau_a_tau_b => dummy
954 2750 : e_tau_b_tau_b => dummy
955 2750 : e_rhoa_rhoa_rhoa => dummy
956 2750 : e_rhoa_rhoa_rhob => dummy
957 2750 : e_rhoa_rhob_rhob => dummy
958 2750 : e_rhob_rhob_rhob => dummy
959 :
960 2750 : IF (grad_deriv >= 0) THEN
961 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
962 2750 : allocate_deriv=.TRUE.)
963 2750 : CALL xc_derivative_get(deriv, deriv_data=e_0)
964 : END IF
965 2750 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
966 1376 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
967 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
968 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
969 1376 : allocate_deriv=.TRUE.)
970 1376 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
971 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
972 1376 : allocate_deriv=.TRUE.)
973 1376 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
974 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
975 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
976 386 : allocate_deriv=.TRUE.)
977 386 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
978 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
979 386 : allocate_deriv=.TRUE.)
980 386 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
981 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
982 386 : allocate_deriv=.TRUE.)
983 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
984 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
985 386 : allocate_deriv=.TRUE.)
986 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
987 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
988 386 : allocate_deriv=.TRUE.)
989 386 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
990 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
991 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
992 942 : allocate_deriv=.TRUE.)
993 942 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
994 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
995 942 : allocate_deriv=.TRUE.)
996 942 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
997 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
998 942 : allocate_deriv=.TRUE.)
999 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
1000 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
1001 942 : allocate_deriv=.TRUE.)
1002 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
1003 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
1004 942 : allocate_deriv=.TRUE.)
1005 942 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
1006 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
1007 942 : allocate_deriv=.TRUE.)
1008 942 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
1009 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
1010 942 : allocate_deriv=.TRUE.)
1011 942 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
1012 942 : IF (has_laplace) THEN
1013 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
1014 180 : allocate_deriv=.TRUE.)
1015 180 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
1016 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
1017 180 : allocate_deriv=.TRUE.)
1018 180 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
1019 : END IF
1020 : CASE default
1021 2704 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1022 : END SELECT
1023 : END IF
1024 2750 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1025 38 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1026 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1027 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1028 38 : allocate_deriv=.TRUE.)
1029 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1030 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1031 38 : allocate_deriv=.TRUE.)
1032 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1033 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1034 38 : allocate_deriv=.TRUE.)
1035 38 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1036 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1037 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1038 34 : allocate_deriv=.TRUE.)
1039 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1040 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1041 34 : allocate_deriv=.TRUE.)
1042 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1043 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1044 34 : allocate_deriv=.TRUE.)
1045 34 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1046 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1047 34 : allocate_deriv=.TRUE.)
1048 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1049 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1050 34 : allocate_deriv=.TRUE.)
1051 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1052 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1053 34 : allocate_deriv=.TRUE.)
1054 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1055 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1056 34 : allocate_deriv=.TRUE.)
1057 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1058 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1059 34 : allocate_deriv=.TRUE.)
1060 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1061 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1062 34 : allocate_deriv=.TRUE.)
1063 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1064 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1065 34 : allocate_deriv=.TRUE.)
1066 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1067 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1068 34 : allocate_deriv=.TRUE.)
1069 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1070 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1071 34 : allocate_deriv=.TRUE.)
1072 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1073 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1074 34 : allocate_deriv=.TRUE.)
1075 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1076 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1077 34 : allocate_deriv=.TRUE.)
1078 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1079 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1080 34 : allocate_deriv=.TRUE.)
1081 34 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1082 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1083 :
1084 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
1085 14 : allocate_deriv=.TRUE.)
1086 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa)
1087 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
1088 14 : allocate_deriv=.TRUE.)
1089 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob)
1090 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
1091 14 : allocate_deriv=.TRUE.)
1092 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob)
1093 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
1094 14 : allocate_deriv=.TRUE.)
1095 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa)
1096 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
1097 14 : allocate_deriv=.TRUE.)
1098 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob)
1099 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
1100 14 : allocate_deriv=.TRUE.)
1101 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa)
1102 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
1103 14 : allocate_deriv=.TRUE.)
1104 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob)
1105 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa], &
1106 14 : allocate_deriv=.TRUE.)
1107 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa)
1108 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
1109 14 : allocate_deriv=.TRUE.)
1110 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob)
1111 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho], &
1112 14 : allocate_deriv=.TRUE.)
1113 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
1114 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa], &
1115 14 : allocate_deriv=.TRUE.)
1116 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa)
1117 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob], &
1118 14 : allocate_deriv=.TRUE.)
1119 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob)
1120 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa], &
1121 14 : allocate_deriv=.TRUE.)
1122 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa)
1123 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob], &
1124 14 : allocate_deriv=.TRUE.)
1125 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob)
1126 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob], &
1127 14 : allocate_deriv=.TRUE.)
1128 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob)
1129 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_a], &
1130 14 : allocate_deriv=.TRUE.)
1131 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_a)
1132 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_tau_b], &
1133 14 : allocate_deriv=.TRUE.)
1134 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_tau_b)
1135 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_a], &
1136 14 : allocate_deriv=.TRUE.)
1137 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_a)
1138 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_tau_b], &
1139 14 : allocate_deriv=.TRUE.)
1140 14 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_tau_b)
1141 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_a], &
1142 14 : allocate_deriv=.TRUE.)
1143 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_a)
1144 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_tau_b], &
1145 14 : allocate_deriv=.TRUE.)
1146 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_tau_b)
1147 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_a], &
1148 14 : allocate_deriv=.TRUE.)
1149 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_a)
1150 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_tau_b], &
1151 14 : allocate_deriv=.TRUE.)
1152 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_tau_b)
1153 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_a], &
1154 14 : allocate_deriv=.TRUE.)
1155 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_a)
1156 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_tau_b], &
1157 14 : allocate_deriv=.TRUE.)
1158 14 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_tau_b)
1159 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_a], &
1160 14 : allocate_deriv=.TRUE.)
1161 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_a)
1162 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a, deriv_tau_b], &
1163 14 : allocate_deriv=.TRUE.)
1164 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a_tau_b)
1165 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b, deriv_tau_b], &
1166 14 : allocate_deriv=.TRUE.)
1167 14 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b_tau_b)
1168 14 : IF (has_laplace) THEN
1169 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhoa], &
1170 6 : allocate_deriv=.TRUE.)
1171 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhoa)
1172 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_laplace_rhob], &
1173 6 : allocate_deriv=.TRUE.)
1174 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_laplace_rhob)
1175 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhoa], &
1176 6 : allocate_deriv=.TRUE.)
1177 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhoa)
1178 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_laplace_rhob], &
1179 6 : allocate_deriv=.TRUE.)
1180 6 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_laplace_rhob)
1181 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhoa], &
1182 6 : allocate_deriv=.TRUE.)
1183 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhoa)
1184 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_laplace_rhob], &
1185 6 : allocate_deriv=.TRUE.)
1186 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_laplace_rhob)
1187 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhoa], &
1188 6 : allocate_deriv=.TRUE.)
1189 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhoa)
1190 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_laplace_rhob], &
1191 6 : allocate_deriv=.TRUE.)
1192 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_laplace_rhob)
1193 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhoa], &
1194 6 : allocate_deriv=.TRUE.)
1195 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhoa)
1196 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_laplace_rhob], &
1197 6 : allocate_deriv=.TRUE.)
1198 6 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_laplace_rhob)
1199 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhoa], &
1200 6 : allocate_deriv=.TRUE.)
1201 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhoa)
1202 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_laplace_rhob], &
1203 6 : allocate_deriv=.TRUE.)
1204 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_laplace_rhob)
1205 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_laplace_rhob], &
1206 6 : allocate_deriv=.TRUE.)
1207 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_laplace_rhob)
1208 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_a], &
1209 6 : allocate_deriv=.TRUE.)
1210 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_a)
1211 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa, deriv_tau_b], &
1212 6 : allocate_deriv=.TRUE.)
1213 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa_tau_b)
1214 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_a], &
1215 6 : allocate_deriv=.TRUE.)
1216 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_a)
1217 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob, deriv_tau_b], &
1218 6 : allocate_deriv=.TRUE.)
1219 6 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob_tau_b)
1220 : END IF
1221 : CASE default
1222 86 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1223 : END SELECT
1224 : END IF
1225 2750 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
1226 0 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1227 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1228 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
1229 0 : allocate_deriv=.TRUE.)
1230 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhoa)
1231 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
1232 0 : allocate_deriv=.TRUE.)
1233 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhob)
1234 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
1235 0 : allocate_deriv=.TRUE.)
1236 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob_rhob)
1237 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
1238 0 : allocate_deriv=.TRUE.)
1239 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob_rhob)
1240 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1241 0 : CPABORT("derivatives larger than 2 not implemented")
1242 : CASE default
1243 0 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1244 : END SELECT
1245 : END IF
1246 2750 : IF (grad_deriv >= 4 .OR. grad_deriv <= -4) THEN
1247 0 : CPABORT("derivatives larger than 3 not implemented")
1248 : END IF
1249 :
1250 : !$OMP PARALLEL DEFAULT(NONE), &
1251 : !$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob),&
1252 : !$OMP SHARED(laplace_rhoa,laplace_rhob,tau_a,tau_b),&
1253 : !$OMP SHARED(e_0,e_rhoa,e_rhob,e_ndrho,e_ndrhoa,e_ndrhob),&
1254 : !$OMP SHARED(e_laplace_rhoa,e_laplace_rhob,e_tau_a,e_tau_b),&
1255 : !$OMP SHARED(e_rhoa_rhoa,e_rhoa_rhob,e_rhob_rhob),&
1256 : !$OMP SHARED(e_ndrho_rhoa,e_ndrho_rhob),&
1257 : !$OMP SHARED(e_ndrhoa_rhoa,e_ndrhoa_rhob,e_ndrhob_rhoa,e_ndrhob_rhob),&
1258 : !$OMP SHARED(e_ndrho_ndrho,e_ndrho_ndrhoa,e_ndrho_ndrhob),&
1259 : !$OMP SHARED(e_ndrhoa_ndrhoa,e_ndrhoa_ndrhob,e_ndrhob_ndrhob),&
1260 : !$OMP SHARED(e_rhoa_laplace_rhoa,e_rhoa_laplace_rhob,e_rhob_laplace_rhoa,e_rhob_laplace_rhob),&
1261 : !$OMP SHARED(e_rhoa_tau_a,e_rhoa_tau_b,e_rhob_tau_a,e_rhob_tau_b),&
1262 : !$OMP SHARED(e_ndrho_laplace_rhoa,e_ndrho_laplace_rhob),&
1263 : !$OMP SHARED(e_ndrhoa_laplace_rhoa,e_ndrhoa_laplace_rhob,e_ndrhob_laplace_rhoa,e_ndrhob_laplace_rhob),&
1264 : !$OMP SHARED(e_ndrho_tau_a,e_ndrho_tau_b),&
1265 : !$OMP SHARED(e_ndrhoa_tau_a,e_ndrhoa_tau_b,e_ndrhob_tau_a,e_ndrhob_tau_b),&
1266 : !$OMP SHARED(e_laplace_rhoa_laplace_rhoa,e_laplace_rhoa_laplace_rhob,e_laplace_rhob_laplace_rhob),&
1267 : !$OMP SHARED(e_laplace_rhoa_tau_a,e_laplace_rhoa_tau_b,e_laplace_rhob_tau_a,e_laplace_rhob_tau_b),&
1268 : !$OMP SHARED(e_tau_a_tau_a,e_tau_a_tau_b,e_tau_b_tau_b),&
1269 : !$OMP SHARED(e_rhoa_rhoa_rhoa,e_rhoa_rhoa_rhob,e_rhoa_rhob_rhob,e_rhob_rhob_rhob),&
1270 : !$OMP SHARED(grad_deriv,npoints),&
1271 : !$OMP SHARED(epsilon_rho,epsilon_tau),&
1272 2750 : !$OMP SHARED(func_name,func_scale,xc_func,xc_info, no_exc, has_laplace)
1273 :
1274 : CALL libxc_spin_polarized_calc(rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1275 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, laplace_rhoa=laplace_rhoa, &
1276 : laplace_rhob=laplace_rhob, tau_a=tau_a, tau_b=tau_b, &
1277 : e_0=e_0, e_rhoa=e_rhoa, e_rhob=e_rhob, e_ndrho=e_ndrho, &
1278 : e_ndrhoa=e_ndrhoa, e_ndrhob=e_ndrhob, e_laplace_rhoa=e_laplace_rhoa, &
1279 : e_laplace_rhob=e_laplace_rhob, e_tau_a=e_tau_a, e_tau_b=e_tau_b, &
1280 : e_rhoa_rhoa=e_rhoa_rhoa, e_rhoa_rhob=e_rhoa_rhob, e_rhob_rhob=e_rhob_rhob, &
1281 : e_ndrho_rhoa=e_ndrho_rhoa, e_ndrho_rhob=e_ndrho_rhob, &
1282 : e_ndrhoa_rhoa=e_ndrhoa_rhoa, e_ndrhoa_rhob=e_ndrhoa_rhob, &
1283 : e_ndrhob_rhoa=e_ndrhob_rhoa, e_ndrhob_rhob=e_ndrhob_rhob, &
1284 : e_ndrho_ndrho=e_ndrho_ndrho, e_ndrho_ndrhoa=e_ndrho_ndrhoa, &
1285 : e_ndrho_ndrhob=e_ndrho_ndrhob, e_ndrhoa_ndrhoa=e_ndrhoa_ndrhoa, &
1286 : e_ndrhoa_ndrhob=e_ndrhoa_ndrhob, e_ndrhob_ndrhob=e_ndrhob_ndrhob, &
1287 : e_rhoa_laplace_rhoa=e_rhoa_laplace_rhoa, &
1288 : e_rhoa_laplace_rhob=e_rhoa_laplace_rhob, &
1289 : e_rhob_laplace_rhoa=e_rhob_laplace_rhoa, &
1290 : e_rhob_laplace_rhob=e_rhob_laplace_rhob, &
1291 : e_rhoa_tau_a=e_rhoa_tau_a, e_rhoa_tau_b=e_rhoa_tau_b, &
1292 : e_rhob_tau_a=e_rhob_tau_a, e_rhob_tau_b=e_rhob_tau_b, &
1293 : e_ndrho_laplace_rhoa=e_ndrho_laplace_rhoa, &
1294 : e_ndrho_laplace_rhob=e_ndrho_laplace_rhob, &
1295 : e_ndrhoa_laplace_rhoa=e_ndrhoa_laplace_rhoa, &
1296 : e_ndrhoa_laplace_rhob=e_ndrhoa_laplace_rhob, &
1297 : e_ndrhob_laplace_rhoa=e_ndrhob_laplace_rhoa, &
1298 : e_ndrhob_laplace_rhob=e_ndrhob_laplace_rhob, &
1299 : e_ndrho_tau_a=e_ndrho_tau_a, e_ndrho_tau_b=e_ndrho_tau_b, &
1300 : e_ndrhoa_tau_a=e_ndrhoa_tau_a, e_ndrhoa_tau_b=e_ndrhoa_tau_b, &
1301 : e_ndrhob_tau_a=e_ndrhob_tau_a, e_ndrhob_tau_b=e_ndrhob_tau_b, &
1302 : e_laplace_rhoa_laplace_rhoa=e_laplace_rhoa_laplace_rhoa, &
1303 : e_laplace_rhoa_laplace_rhob=e_laplace_rhoa_laplace_rhob, &
1304 : e_laplace_rhob_laplace_rhob=e_laplace_rhob_laplace_rhob, &
1305 : e_laplace_rhoa_tau_a=e_laplace_rhoa_tau_a, &
1306 : e_laplace_rhoa_tau_b=e_laplace_rhoa_tau_b, &
1307 : e_laplace_rhob_tau_a=e_laplace_rhob_tau_a, &
1308 : e_laplace_rhob_tau_b=e_laplace_rhob_tau_b, &
1309 : e_tau_a_tau_a=e_tau_a_tau_a, &
1310 : e_tau_a_tau_b=e_tau_a_tau_b, &
1311 : e_tau_b_tau_b=e_tau_b_tau_b, &
1312 : e_rhoa_rhoa_rhoa=e_rhoa_rhoa_rhoa, &
1313 : e_rhoa_rhoa_rhob=e_rhoa_rhoa_rhob, &
1314 : e_rhoa_rhob_rhob=e_rhoa_rhob_rhob, &
1315 : e_rhob_rhob_rhob=e_rhob_rhob_rhob, &
1316 : grad_deriv=grad_deriv, npoints=npoints, &
1317 : epsilon_rho=epsilon_rho, &
1318 : epsilon_tau=epsilon_tau, func_name=func_name, &
1319 : sc=func_scale, xc_func=xc_func, xc_info=xc_info, no_exc=no_exc, has_laplace=has_laplace)
1320 :
1321 : !$OMP END PARALLEL
1322 :
1323 2750 : NULLIFY (dummy)
1324 :
1325 2750 : CALL xc_f03_func_end(xc_func)
1326 :
1327 2750 : CALL timestop(handle)
1328 : #else
1329 : MARK_USED(rho_set)
1330 : MARK_USED(deriv_set)
1331 : MARK_USED(grad_deriv)
1332 : MARK_USED(libxc_params)
1333 : MARK_USED(func_name_override)
1334 :
1335 : CALL cp_abort(__LOCATION__, "Unknown functional! If you are asking "// &
1336 : "for a functional of the LibXC library, "// &
1337 : "you have to download and install the library!")
1338 : #endif
1339 2750 : END SUBROUTINE libxc_spin_polarized_eval
1340 :
1341 : ! **************************************************************************************************
1342 : !> \brief libxc exchange-correlation functionals
1343 : !> \param rho density
1344 : !> \param norm_drho norm of the gradient of the density
1345 : !> \param laplace_rho laplacian of the density
1346 : !> \param tau kinetic-energy density
1347 : !> \param e_0 energy density
1348 : !> \param e_rho derivative of the energy density with respect to rho
1349 : !> \param e_ndrho derivative of the energy density with respect to ndrho
1350 : !> \param e_laplace_rho derivative of the energy density with respect to laplace_rho
1351 : !> \param e_tau derivative of the energy density with respect to tau
1352 : !> \param e_rho_rho derivative of the energy density with respect to rho_rho
1353 : !> \param e_ndrho_rho derivative of the energy density with respect to ndrho_rho
1354 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
1355 : !> \param e_rho_laplace_rho derivative of the energy density with respect to rho_laplace_rho
1356 : !> \param e_rho_tau derivative of the energy density with respect to rho_tau
1357 : !> \param e_ndrho_laplace_rho derivative of the energy density with respect to ndrho_laplace_rho
1358 : !> \param e_ndrho_tau derivative of the energy density with respect to ndrho_tau
1359 : !> \param e_laplace_rho_laplace_rho derivative of the energy density with respect to laplace_rho_laplace_rho
1360 : !> \param e_laplace_rho_tau derivative of the energy density with respect to laplace_rho_tau
1361 : !> \param e_tau_tau derivative of the energy density with respect to tau_tau
1362 : !> \param e_rho_rho_rho derivative of the energy density with respect to rho_rho_rho
1363 : !> \param grad_deriv degree of the derivative that should be evaluated,
1364 : !> if positive all the derivatives up to the given degree are evaluated,
1365 : !> if negative only the given degree is calculated
1366 : !> \param npoints number of points on the grid
1367 : !> \param epsilon_rho ...
1368 : !> \param epsilon_tau ...
1369 : !> \param func_name name of the functional
1370 : !> \param sc scaling factor of the functional
1371 : !> \param xc_func libxc functional object
1372 : !> \param xc_info libxc functional info object
1373 : !> \param no_exc whether the EXC function is not available for the given functional
1374 : !> \param has_laplace ...
1375 : !> \author F. Tran
1376 : ! **************************************************************************************************
1377 : #if defined (__LIBXC)
1378 18110 : SUBROUTINE libxc_spin_unpolarized_calc(rho, norm_drho, laplace_rho, tau, &
1379 : e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, e_rho_rho, e_ndrho_rho, &
1380 : e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1381 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, &
1382 : e_tau_tau, e_rho_rho_rho, &
1383 : grad_deriv, npoints, epsilon_rho, &
1384 : epsilon_tau, func_name, sc, xc_func, xc_info, no_exc, has_laplace)
1385 :
1386 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho, laplace_rho, tau
1387 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, &
1388 : e_rho_rho, e_ndrho_rho, e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1389 : e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_tau_tau, e_rho_rho_rho
1390 : INTEGER, INTENT(in) :: grad_deriv, npoints
1391 : REAL(KIND=dp), INTENT(in) :: epsilon_rho, epsilon_tau
1392 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
1393 : REAL(KIND=dp), INTENT(in) :: sc
1394 : TYPE(xc_f03_func_t), INTENT(IN) :: xc_func
1395 : TYPE(xc_f03_func_info_t), INTENT(IN) :: xc_info
1396 : LOGICAL, INTENT(IN) :: no_exc, has_laplace
1397 :
1398 : INTEGER :: ii
1399 : REAL(KIND=dp), DIMENSION(1) :: exc, my_tau, sigma, v2lapl2, v2lapltau, v2rho2, v2rholapl, &
1400 : v2rhosigma, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2tau2, v3rho3, vlapl, vrho, &
1401 : vsigma, vtau
1402 :
1403 : ! init vlapl (prevent libxc-4.0.x bug)
1404 18110 : vlapl = 0.0_dp
1405 :
1406 10390 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1407 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1408 10390 : IF (grad_deriv == 0) THEN
1409 84 : !$OMP DO
1410 : DO ii = 1, npoints
1411 6410954 : IF (rho(ii) > epsilon_rho) THEN
1412 6408106 : CALL xc_f03_lda_exc(xc_func, one, rho(ii), exc)
1413 6408106 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1414 : END IF
1415 : END DO
1416 : !$OMP END DO
1417 : ELSE IF (grad_deriv == -1) THEN
1418 0 : !$OMP DO
1419 : DO ii = 1, npoints
1420 0 : IF (rho(ii) > epsilon_rho) THEN
1421 0 : CALL xc_f03_lda_vxc(xc_func, one, rho(ii), vrho)
1422 0 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1423 : END IF
1424 : END DO
1425 : !$OMP END DO
1426 : ELSE IF (grad_deriv == 1) THEN
1427 8722 : !$OMP DO
1428 : DO ii = 1, npoints
1429 230058306 : IF (rho(ii) > epsilon_rho) THEN
1430 220249911 : CALL xc_f03_lda_exc_vxc(xc_func, one, rho(ii), exc, vrho)
1431 220249911 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1432 220249911 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1433 : END IF
1434 : END DO
1435 : !$OMP END DO
1436 : ELSE IF (grad_deriv == -2) THEN
1437 0 : !$OMP DO
1438 : DO ii = 1, npoints
1439 0 : IF (rho(ii) > epsilon_rho) THEN
1440 0 : CALL xc_f03_lda_fxc(xc_func, one, rho(ii), v2rho2)
1441 0 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1442 : END IF
1443 : END DO
1444 : !$OMP END DO
1445 : ELSE IF (grad_deriv == 2) THEN
1446 1584 : !$OMP DO
1447 : DO ii = 1, npoints
1448 9897486 : IF (rho(ii) > epsilon_rho) THEN
1449 9365806 : CALL xc_f03_lda_exc_vxc_fxc(xc_func, one, rho(ii), exc, vrho, v2rho2)
1450 9365806 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1451 9365806 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1452 9365806 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1453 : END IF
1454 : END DO
1455 : !$OMP END DO
1456 : ELSE IF (grad_deriv == -3) THEN
1457 0 : !$OMP DO
1458 : DO ii = 1, npoints
1459 0 : IF (rho(ii) > epsilon_rho) THEN
1460 0 : CALL xc_f03_lda_kxc(xc_func, one, rho(ii), v3rho3)
1461 0 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1)
1462 : END IF
1463 : END DO
1464 : !$OMP END DO
1465 : ELSE IF (grad_deriv == 3) THEN
1466 0 : !$OMP DO
1467 : DO ii = 1, npoints
1468 0 : IF (rho(ii) > epsilon_rho) THEN
1469 0 : CALL xc_f03_lda(xc_func, one, rho(ii), exc, vrho, v2rho2, v3rho3)
1470 0 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1471 0 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1472 0 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1473 0 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1)
1474 : END IF
1475 : END DO
1476 : !$OMP END DO
1477 : END IF
1478 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1479 5482 : IF (grad_deriv == 0) THEN
1480 142 : !$OMP DO
1481 : DO ii = 1, npoints
1482 6104994 : IF (rho(ii) > epsilon_rho) THEN
1483 12037086 : sigma = norm_drho(ii)**2
1484 6018543 : CALL xc_f03_gga_exc(xc_func, one, rho(ii), sigma, exc)
1485 6018543 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1486 : END IF
1487 : END DO
1488 : !$OMP END DO
1489 : ELSE IF (grad_deriv == -1) THEN
1490 0 : !$OMP DO
1491 : DO ii = 1, npoints
1492 0 : IF (rho(ii) > epsilon_rho) THEN
1493 0 : sigma = norm_drho(ii)**2
1494 0 : CALL xc_f03_gga_vxc(xc_func, one, rho(ii), sigma, vrho, vsigma)
1495 0 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1496 0 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1497 : END IF
1498 : END DO
1499 : !$OMP END DO
1500 : ELSE IF (grad_deriv == 1) THEN
1501 4640 : !$OMP DO
1502 : DO ii = 1, npoints
1503 164167114 : IF (rho(ii) > epsilon_rho) THEN
1504 220792972 : sigma = norm_drho(ii)**2
1505 110396486 : IF (no_exc) THEN
1506 0 : CALL xc_f03_gga_vxc(xc_func, one, rho(ii), sigma, vrho, vsigma)
1507 0 : exc = 0.0_dp
1508 : ELSE
1509 : CALL xc_f03_gga_exc_vxc(xc_func, one, rho(ii), sigma, &
1510 110396486 : exc, vrho, vsigma)
1511 : END IF
1512 110396486 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1513 110396486 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1514 110396486 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1515 : END IF
1516 : END DO
1517 : !$OMP END DO
1518 : ELSE IF (grad_deriv == -2) THEN
1519 0 : !$OMP DO
1520 : DO ii = 1, npoints
1521 0 : IF (rho(ii) > epsilon_rho) THEN
1522 0 : sigma = norm_drho(ii)**2
1523 0 : IF (no_exc) THEN
1524 : CALL xc_f03_gga_vxc_fxc(xc_func, one, rho(ii), sigma, vrho, vsigma, &
1525 0 : v2rho2, v2rhosigma, v2sigma2)
1526 : ELSE
1527 : CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rho(ii), sigma, &
1528 : exc, vrho, vsigma, v2rho2, &
1529 0 : v2rhosigma, v2sigma2)
1530 : END IF
1531 0 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1532 0 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1533 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1534 0 : sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1535 : END IF
1536 : END DO
1537 : !$OMP END DO
1538 : ELSE IF (grad_deriv == 2) THEN
1539 700 : !$OMP DO
1540 : DO ii = 1, npoints
1541 4286530 : IF (rho(ii) > epsilon_rho) THEN
1542 8153720 : sigma = norm_drho(ii)**2
1543 4076860 : IF (no_exc) THEN
1544 : CALL xc_f03_gga_vxc_fxc(xc_func, one, rho(ii), sigma, vrho, vsigma, &
1545 0 : v2rho2, v2rhosigma, v2sigma2)
1546 0 : exc = 0.0_dp
1547 : ELSE
1548 : CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rho(ii), sigma, &
1549 : exc, vrho, vsigma, &
1550 4076860 : v2rho2, v2rhosigma, v2sigma2)
1551 : END IF
1552 4076860 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1553 4076860 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1554 4076860 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1555 4076860 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1556 4076860 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1557 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1558 4076860 : sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1559 : END IF
1560 : END DO
1561 : !$OMP END DO
1562 : END IF
1563 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
1564 2238 : IF (grad_deriv == 0) THEN
1565 28 : !$OMP DO
1566 : DO ii = 1, npoints
1567 526000 : IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN
1568 1007308 : sigma = norm_drho(ii)**2
1569 503654 : my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1570 : CALL xc_f03_mgga_exc(xc_func, one, rho(ii), sigma, &
1571 503654 : laplace_rho(ii), my_tau, exc)
1572 503654 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1573 : END IF
1574 : END DO
1575 : !$OMP END DO
1576 : ELSE IF (grad_deriv == -1) THEN
1577 0 : !$OMP DO
1578 : DO ii = 1, npoints
1579 0 : IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN
1580 0 : sigma = norm_drho(ii)**2
1581 0 : my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1582 : CALL xc_f03_mgga_vxc(xc_func, one, rho(ii), sigma, &
1583 0 : laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau)
1584 0 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1585 0 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1586 0 : IF (has_laplace) e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1587 0 : e_tau(ii) = e_tau(ii) + sc*vtau(1)
1588 : END IF
1589 : END DO
1590 : !$OMP END DO
1591 : ELSE IF (grad_deriv == 1) THEN
1592 1902 : !$OMP DO
1593 : DO ii = 1, npoints
1594 82174141 : IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN
1595 80973655 : sigma(1) = norm_drho(ii)**2
1596 80973655 : my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1597 80973655 : IF (no_exc) THEN
1598 : CALL xc_f03_mgga_vxc(xc_func, one, rho(ii), sigma, &
1599 0 : laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau)
1600 0 : exc = 0.0_dp
1601 : ELSE
1602 : CALL xc_f03_mgga_exc_vxc(xc_func, one, rho(ii), sigma, &
1603 80973655 : laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau)
1604 : END IF
1605 80973655 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1606 80973655 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1607 80973655 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1608 80973655 : IF (has_laplace) e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1609 80973655 : e_tau(ii) = e_tau(ii) + sc*vtau(1)
1610 : END IF
1611 : END DO
1612 : !$OMP END DO
1613 : ELSE IF (grad_deriv == -2) THEN
1614 0 : !$OMP DO
1615 : DO ii = 1, npoints
1616 0 : IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN
1617 0 : sigma = norm_drho(ii)**2
1618 0 : my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1619 0 : IF (no_exc) THEN
1620 : CALL xc_f03_mgga_vxc_fxc(xc_func, one, rho(ii), sigma, &
1621 : laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau, &
1622 : v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1623 0 : v2lapl2, v2lapltau, v2tau2)
1624 : ELSE
1625 : CALL xc_f03_mgga(xc_func, one, rho(ii), sigma, &
1626 : laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, &
1627 : v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1628 0 : v2lapl2, v2lapltau, v2tau2)
1629 : END IF
1630 0 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1631 0 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1632 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1633 0 : sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1634 0 : e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1)
1635 0 : e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii)
1636 0 : e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1)
1637 0 : IF (has_laplace) THEN
1638 0 : e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1)
1639 : e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
1640 0 : sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii)
1641 0 : e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1)
1642 0 : e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1)
1643 : END IF
1644 : END IF
1645 : END DO
1646 : !$OMP END DO
1647 : ELSE IF (grad_deriv == 2) THEN
1648 308 : !$OMP DO
1649 : DO ii = 1, npoints
1650 7209168 : IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN
1651 14150184 : sigma = norm_drho(ii)**2
1652 7075092 : my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1653 7075092 : IF (no_exc) THEN
1654 : CALL xc_f03_mgga_vxc_fxc(xc_func, one, rho(ii), sigma, &
1655 : laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau, &
1656 : v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1657 0 : v2lapl2, v2lapltau, v2tau2)
1658 0 : exc = 0.0_dp
1659 : ELSE
1660 : CALL xc_f03_mgga(xc_func, one, rho(ii), sigma, &
1661 : laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, &
1662 : v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1663 7075092 : v2lapl2, v2lapltau, v2tau2)
1664 : END IF
1665 7075092 : e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1666 7075092 : e_rho(ii) = e_rho(ii) + sc*vrho(1)
1667 7075092 : e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1668 7075092 : e_tau(ii) = e_tau(ii) + sc*vtau(1)
1669 7075092 : e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1670 7075092 : e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1671 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1672 7075092 : sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1673 7075092 : e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1)
1674 7075092 : e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii)
1675 7075092 : e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1)
1676 7075092 : IF (has_laplace) THEN
1677 2342952 : e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1678 2342952 : e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1)
1679 : e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
1680 2342952 : sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii)
1681 2342952 : e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1)
1682 2342952 : e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1)
1683 : END IF
1684 : END IF
1685 : END DO
1686 : !$OMP END DO
1687 : END IF
1688 : CASE default
1689 18110 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
1690 : END SELECT
1691 :
1692 18110 : END SUBROUTINE libxc_spin_unpolarized_calc
1693 : #endif
1694 :
1695 : ! **************************************************************************************************
1696 : !> \brief libxc exchange-correlation functionals
1697 : !> \param rhoa alpha density
1698 : !> \param rhob beta density
1699 : !> \param norm_drho ...
1700 : !> \param norm_drhoa norm of the gradient of the alpha density
1701 : !> \param norm_drhob norm of the gradient of the beta density
1702 : !> \param laplace_rhoa laplacian of the alpha density
1703 : !> \param laplace_rhob laplacian of the beta density
1704 : !> \param tau_a alpha kinetic-energy density
1705 : !> \param tau_b beta kinetic-energy density
1706 : !> \param e_0 energy density
1707 : !> \param e_rhoa derivative of the energy density with respect to rhoa
1708 : !> \param e_rhob derivative of the energy density with respect to rhob
1709 : !> \param e_ndrho derivative of the energy density with respect to ndrho
1710 : !> \param e_ndrhoa derivative of the energy density with respect to ndrhoa
1711 : !> \param e_ndrhob derivative of the energy density with respect to ndrhob
1712 : !> \param e_laplace_rhoa derivative of the energy density with respect to laplace_rhoa
1713 : !> \param e_laplace_rhob derivative of the energy density with respect to laplace_rhob
1714 : !> \param e_tau_a derivative of the energy density with respect to tau_a
1715 : !> \param e_tau_b derivative of the energy density with respect to tau_b
1716 : !> \param e_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa
1717 : !> \param e_rhoa_rhob derivative of the energy density with respect to rhoa_rhob
1718 : !> \param e_rhob_rhob derivative of the energy density with respect to rhob_rhob
1719 : !> \param e_ndrho_rhoa derivative of the energy density with respect to ndrho_rhoa
1720 : !> \param e_ndrho_rhob derivative of the energy density with respect to ndrho_rhob
1721 : !> \param e_ndrhoa_rhoa derivative of the energy density with respect to ndrhoa_rhoa
1722 : !> \param e_ndrhoa_rhob derivative of the energy density with respect to ndrhoa_rhob
1723 : !> \param e_ndrhob_rhoa derivative of the energy density with respect to ndrhob_rhoa
1724 : !> \param e_ndrhob_rhob derivative of the energy density with respect to ndrhob_rhob
1725 : !> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho
1726 : !> \param e_ndrho_ndrhoa derivative of the energy density with respect to ndrho_ndrhoa
1727 : !> \param e_ndrho_ndrhob derivative of the energy density with respect to ndrho_ndrhob
1728 : !> \param e_ndrhoa_ndrhoa derivative of the energy density with respect to ndrhoa_ndrhoa
1729 : !> \param e_ndrhoa_ndrhob derivative of the energy density with respect to ndrhoa_ndrhob
1730 : !> \param e_ndrhob_ndrhob derivative of the energy density with respect to ndrhob_ndrhob
1731 : !> \param e_rhoa_laplace_rhoa derivative of the energy density with respect to rhoa_laplace_rhoa
1732 : !> \param e_rhoa_laplace_rhob derivative of the energy density with respect to rhoa_laplace_rhob
1733 : !> \param e_rhob_laplace_rhoa derivative of the energy density with respect to rhob_laplace_rhoa
1734 : !> \param e_rhob_laplace_rhob derivative of the energy density with respect to rhob_laplace_rhob
1735 : !> \param e_rhoa_tau_a derivative of the energy density with respect to rhoa_tau_a
1736 : !> \param e_rhoa_tau_b derivative of the energy density with respect to rhoa_tau_b
1737 : !> \param e_rhob_tau_a derivative of the energy density with respect to rhob_tau_a
1738 : !> \param e_rhob_tau_b derivative of the energy density with respect to rhob_tau_b
1739 : !> \param e_ndrho_laplace_rhoa derivative of the energy density with respect to ndrho_laplace_rhoa
1740 : !> \param e_ndrho_laplace_rhob derivative of the energy density with respect to ndrho_laplace_rhob
1741 : !> \param e_ndrhoa_laplace_rhoa derivative of the energy density with respect to ndrhoa_laplace_rhoa
1742 : !> \param e_ndrhoa_laplace_rhob derivative of the energy density with respect to ndrhoa_laplace_rhob
1743 : !> \param e_ndrhob_laplace_rhoa derivative of the energy density with respect to ndrhob_laplace_rhoa
1744 : !> \param e_ndrhob_laplace_rhob derivative of the energy density with respect to ndrhob_laplace_rhob
1745 : !> \param e_ndrho_tau_a derivative of the energy density with respect to ndrho_tau_a
1746 : !> \param e_ndrho_tau_b derivative of the energy density with respect to ndrho_tau_b
1747 : !> \param e_ndrhoa_tau_a derivative of the energy density with respect to ndrhoa_tau_a
1748 : !> \param e_ndrhoa_tau_b derivative of the energy density with respect to ndrhoa_tau_b
1749 : !> \param e_ndrhob_tau_a derivative of the energy density with respect to ndrhob_tau_a
1750 : !> \param e_ndrhob_tau_b derivative of the energy density with respect to ndrhob_tau_b
1751 : !> \param e_laplace_rhoa_laplace_rhoa derivative of the energy density with respect to laplace_rhoa_laplace_rhoa
1752 : !> \param e_laplace_rhoa_laplace_rhob derivative of the energy density with respect to laplace_rhoa_laplace_rhob
1753 : !> \param e_laplace_rhob_laplace_rhob derivative of the energy density with respect to laplace_rhob_laplace_rhob
1754 : !> \param e_laplace_rhoa_tau_a derivative of the energy density with respect to laplace_rhoa_tau_a
1755 : !> \param e_laplace_rhoa_tau_b derivative of the energy density with respect to laplace_rhoa_tau_b
1756 : !> \param e_laplace_rhob_tau_a derivative of the energy density with respect to laplace_rhob_tau_a
1757 : !> \param e_laplace_rhob_tau_b derivative of the energy density with respect to laplace_rhob_tau_b
1758 : !> \param e_tau_a_tau_a derivative of the energy density with respect to tau_a_tau_a
1759 : !> \param e_tau_a_tau_b derivative of the energy density with respect to tau_a_tau_b
1760 : !> \param e_tau_b_tau_b derivative of the energy density with respect to tau_b_tau_b
1761 : !> \param e_rhoa_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa_rhoa
1762 : !> \param e_rhoa_rhoa_rhob derivative of the energy density with respect to rhoa_rhoa_rhob
1763 : !> \param e_rhoa_rhob_rhob derivative of the energy density with respect to rhoa_rhob_rhob
1764 : !> \param e_rhob_rhob_rhob derivative of the energy density with respect to rhob_rhob_rhob
1765 : !> \param grad_deriv degree of the derivative that should be evaluated,
1766 : !> if positive all the derivatives up to the given degree are evaluated,
1767 : !> if negative only the given degree is calculated
1768 : !> \param npoints number of points on the grid
1769 : !> \param epsilon_rho ...
1770 : !> \param epsilon_tau ...
1771 : !> \param func_name name of the functional
1772 : !> \param sc scaling factor of the functional
1773 : !> \param xc_func libxc functional object
1774 : !> \param xc_info libxc functional info object
1775 : !> \param no_exc whether the EXC function is not available for the given functional
1776 : !> \param has_laplace ...
1777 : !> \author F. Tran
1778 : ! **************************************************************************************************
1779 : #if defined (__LIBXC)
1780 2750 : SUBROUTINE libxc_spin_polarized_calc(rhoa, rhob, norm_drho, norm_drhoa, &
1781 : norm_drhob, laplace_rhoa, laplace_rhob, tau_a, tau_b, &
1782 : e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, e_ndrhob, &
1783 : e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, &
1784 : e_rhoa_rhoa, e_rhoa_rhob, e_rhob_rhob, &
1785 : e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, &
1786 : e_ndrhoa_rhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
1787 : e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
1788 : e_ndrhoa_ndrhoa, e_ndrhoa_ndrhob, e_ndrhob_ndrhob, &
1789 : e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
1790 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, &
1791 : e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, e_rhob_tau_b, &
1792 : e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, &
1793 : e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, &
1794 : e_ndrhob_laplace_rhoa, e_ndrhob_laplace_rhob, &
1795 : e_ndrho_tau_a, e_ndrho_tau_b, &
1796 : e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
1797 : e_ndrhob_tau_a, e_ndrhob_tau_b, &
1798 : e_laplace_rhoa_laplace_rhoa, &
1799 : e_laplace_rhoa_laplace_rhob, &
1800 : e_laplace_rhob_laplace_rhob, &
1801 : e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
1802 : e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, &
1803 : e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
1804 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
1805 : e_rhoa_rhob_rhob, e_rhob_rhob_rhob, &
1806 : grad_deriv, npoints, epsilon_rho, &
1807 : epsilon_tau, func_name, sc, xc_func, xc_info, no_exc, has_laplace)
1808 :
1809 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, norm_drho, norm_drhoa, &
1810 : norm_drhob, laplace_rhoa, &
1811 : laplace_rhob, tau_a, tau_b
1812 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, &
1813 : e_ndrhob, e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, e_rhoa_rhoa, e_rhoa_rhob, &
1814 : e_rhob_rhob, e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhob_rhoa, &
1815 : e_ndrhob_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, e_ndrhoa_ndrhoa, &
1816 : e_ndrhoa_ndrhob, e_ndrhob_ndrhob, e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
1817 : e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, &
1818 : e_rhob_tau_b, e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, e_ndrhoa_laplace_rhoa
1819 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_ndrhoa_laplace_rhob, e_ndrhob_laplace_rhoa, &
1820 : e_ndrhob_laplace_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
1821 : e_ndrhob_tau_a, e_ndrhob_tau_b, e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
1822 : e_laplace_rhob_laplace_rhob, e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
1823 : 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, &
1824 : e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, e_rhoa_rhob_rhob, e_rhob_rhob_rhob
1825 : INTEGER, INTENT(in) :: grad_deriv, npoints
1826 : REAL(KIND=dp), INTENT(in) :: epsilon_rho, epsilon_tau
1827 : CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name
1828 : REAL(KIND=dp), INTENT(in) :: sc
1829 : TYPE(xc_f03_func_t), INTENT(IN) :: xc_func
1830 : TYPE(xc_f03_func_info_t), INTENT(IN) :: xc_info
1831 : LOGICAL, INTENT(IN) :: no_exc, has_laplace
1832 :
1833 : INTEGER :: ii
1834 : REAL(KIND=dp) :: my_norm_drho, my_norm_drhoa, &
1835 : my_norm_drhob, my_rhoa, my_rhob, &
1836 : my_tau_a, my_tau_b
1837 : REAL(KIND=dp), DIMENSION(1) :: exc
1838 : REAL(KIND=dp), DIMENSION(2, 1) :: laplace_rhov, rhov, tauv, vlapl, vrho, &
1839 : vtau
1840 : REAL(KIND=dp), DIMENSION(3, 1) :: sigmav, v2lapl2, v2rho2, v2tau2, vsigma
1841 : REAL(KIND=dp), DIMENSION(4, 1) :: v2lapltau, v2rholapl, v2rhotau, v3rho3
1842 : REAL(KIND=dp), DIMENSION(6, 1) :: v2rhosigma, v2sigma2, v2sigmalapl, &
1843 : v2sigmatau
1844 :
1845 2750 : vlapl(1, 1) = 0.0_dp
1846 2750 : vlapl(2, 1) = 0.0_dp
1847 :
1848 1376 : SELECT CASE (xc_f03_func_info_get_family(xc_info))
1849 : CASE (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA)
1850 1376 : IF (grad_deriv == 0) THEN
1851 0 : !$OMP DO
1852 : DO ii = 1, npoints
1853 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1854 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1855 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1856 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1857 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1858 0 : CALL xc_f03_lda_exc(xc_func, one, rhov(1, 1), exc)
1859 0 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1860 : END IF
1861 : END DO
1862 : !$OMP END DO
1863 : ELSE IF (grad_deriv == -1) THEN
1864 0 : !$OMP DO
1865 : DO ii = 1, npoints
1866 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1867 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1868 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1869 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1870 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1871 0 : CALL xc_f03_lda_vxc(xc_func, one, rhov(1, 1), vrho(1, 1))
1872 0 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1873 0 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1874 : END IF
1875 : END DO
1876 : !$OMP END DO
1877 : ELSE IF (grad_deriv == 1) THEN
1878 1338 : !$OMP DO
1879 : DO ii = 1, npoints
1880 59646240 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1881 59646240 : my_rhob = MAX(rhob(ii), 0.0_dp)
1882 59646240 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1883 56904753 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1884 56904753 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1885 56904753 : CALL xc_f03_lda_exc_vxc(xc_func, one, rhov(1, 1), exc, vrho(1, 1))
1886 56904753 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1887 56904753 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1888 56904753 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1889 : END IF
1890 : END DO
1891 : !$OMP END DO
1892 : ELSE IF (grad_deriv == -2) THEN
1893 0 : !$OMP DO
1894 : DO ii = 1, npoints
1895 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1896 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1897 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1898 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1899 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1900 0 : CALL xc_f03_lda_fxc(xc_func, one, rhov(1, 1), v2rho2(1, 1))
1901 0 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1902 0 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1903 0 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1904 : END IF
1905 : END DO
1906 : !$OMP END DO
1907 : ELSE IF (grad_deriv == 2) THEN
1908 38 : !$OMP DO
1909 : DO ii = 1, npoints
1910 873348 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1911 873348 : my_rhob = MAX(rhob(ii), 0.0_dp)
1912 873348 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1913 848540 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1914 848540 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1915 848540 : CALL xc_f03_lda_exc_vxc_fxc(xc_func, one, rhov(1, 1), exc, vrho(1, 1), v2rho2(1, 1))
1916 848540 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1917 848540 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1918 848540 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1919 848540 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1920 848540 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1921 848540 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1922 : END IF
1923 : END DO
1924 : !$OMP END DO
1925 : ELSE IF (grad_deriv == -3) THEN
1926 0 : !$OMP DO
1927 : DO ii = 1, npoints
1928 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1929 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1930 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1931 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1932 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1933 0 : CALL xc_f03_lda_kxc(xc_func, one, rhov(1, 1), v3rho3(1, 1))
1934 0 : e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1)
1935 0 : e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1)
1936 0 : e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1)
1937 0 : e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1)
1938 : END IF
1939 : END DO
1940 : !$OMP END DO
1941 : ELSE IF (grad_deriv == 3) THEN
1942 0 : !$OMP DO
1943 : DO ii = 1, npoints
1944 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1945 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1946 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1947 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1948 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1949 0 : CALL xc_f03_lda(xc_func, one, rhov(1, 1), exc, vrho(1, 1), v2rho2(1, 1), v3rho3(1, 1))
1950 0 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1951 0 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1952 0 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1953 0 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1954 0 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1955 0 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1956 0 : e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1)
1957 0 : e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1)
1958 0 : e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1)
1959 0 : e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1)
1960 : END IF
1961 : END DO
1962 : !$OMP END DO
1963 : END IF
1964 : CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
1965 402 : IF (grad_deriv == 0) THEN
1966 16 : !$OMP DO
1967 : DO ii = 1, npoints
1968 262144 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1969 262144 : my_rhob = MAX(rhob(ii), 0.0_dp)
1970 262144 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1971 262144 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1972 262144 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1973 262144 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
1974 262144 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
1975 262144 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
1976 262144 : sigmav(1, 1) = my_norm_drhoa**2
1977 262144 : sigmav(3, 1) = my_norm_drhob**2
1978 262144 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
1979 262144 : CALL xc_f03_gga_exc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc)
1980 262144 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1981 : END IF
1982 : END DO
1983 : !$OMP END DO
1984 : ELSE IF (grad_deriv == -1) THEN
1985 0 : !$OMP DO
1986 : DO ii = 1, npoints
1987 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
1988 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
1989 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
1990 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
1991 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
1992 0 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
1993 0 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
1994 0 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
1995 0 : sigmav(1, 1) = my_norm_drhoa**2
1996 0 : sigmav(3, 1) = my_norm_drhob**2
1997 0 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
1998 0 : CALL xc_f03_gga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1))
1999 0 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2000 0 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2001 0 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2002 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2003 0 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2004 : e_ndrhob(ii) = e_ndrhob(ii) + &
2005 0 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2006 : END IF
2007 : END DO
2008 : !$OMP END DO
2009 : ELSE IF (grad_deriv == 1) THEN
2010 352 : !$OMP DO
2011 : DO ii = 1, npoints
2012 7639848 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2013 7639848 : my_rhob = MAX(rhob(ii), 0.0_dp)
2014 7639848 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
2015 7535189 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2016 7535189 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2017 7535189 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2018 7535189 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2019 7535189 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2020 7535189 : sigmav(1, 1) = my_norm_drhoa**2
2021 7535189 : sigmav(3, 1) = my_norm_drhob**2
2022 7535189 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2023 7535189 : IF (no_exc) THEN
2024 0 : CALL xc_f03_gga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1))
2025 0 : exc = 0.0_dp
2026 : ELSE
2027 7535189 : CALL xc_f03_gga_exc_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1))
2028 : END IF
2029 7535189 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2030 7535189 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2031 7535189 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2032 7535189 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2033 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2034 7535189 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2035 : e_ndrhob(ii) = e_ndrhob(ii) + &
2036 7535189 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2037 : END IF
2038 : END DO
2039 : !$OMP END DO
2040 : ELSE IF (grad_deriv == -2) THEN
2041 0 : !$OMP DO
2042 : DO ii = 1, npoints
2043 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2044 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
2045 0 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
2046 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2047 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2048 0 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2049 0 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2050 0 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2051 0 : sigmav(1, 1) = my_norm_drhoa**2
2052 0 : sigmav(3, 1) = my_norm_drhob**2
2053 0 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2054 0 : IF (no_exc) THEN
2055 : CALL xc_f03_gga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1), &
2056 0 : v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2057 : ELSE
2058 : CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2059 0 : v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2060 : END IF
2061 0 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2062 0 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2063 0 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2064 0 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2065 0 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2066 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2067 0 : sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2068 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2069 0 : sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2070 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2071 0 : sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2072 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2073 0 : sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2074 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2075 0 : sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2076 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2077 0 : sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2078 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2079 0 : sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2080 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2081 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2082 0 : 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2083 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2084 : sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2085 0 : 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2086 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2087 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2088 0 : 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2089 : END IF
2090 : END DO
2091 : !$OMP END DO
2092 : ELSE IF (grad_deriv == 2) THEN
2093 34 : !$OMP DO
2094 : DO ii = 1, npoints
2095 663036 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2096 663036 : my_rhob = MAX(rhob(ii), 0.0_dp)
2097 663036 : IF ((my_rhoa + my_rhob) > epsilon_rho) THEN
2098 627564 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2099 627564 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2100 627564 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2101 627564 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2102 627564 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2103 627564 : sigmav(1, 1) = my_norm_drhoa**2
2104 627564 : sigmav(3, 1) = my_norm_drhob**2
2105 627564 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2106 627564 : IF (no_exc) THEN
2107 : CALL xc_f03_gga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1), &
2108 0 : v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2109 0 : exc = 0.0_dp
2110 : ELSE
2111 : CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2112 627564 : v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2113 : END IF
2114 627564 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2115 627564 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2116 627564 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2117 627564 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2118 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2119 627564 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2120 : e_ndrhob(ii) = e_ndrhob(ii) + &
2121 627564 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2122 627564 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2123 627564 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2124 627564 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2125 627564 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2126 627564 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2127 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2128 627564 : sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2129 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2130 627564 : sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2131 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2132 627564 : sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2133 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2134 627564 : sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2135 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2136 627564 : sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2137 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2138 627564 : sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2139 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2140 627564 : sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2141 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2142 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2143 627564 : 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2144 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2145 : sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2146 627564 : 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2147 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2148 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2149 627564 : 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2150 : END IF
2151 : END DO
2152 : !$OMP END DO
2153 : END IF
2154 : CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA)
2155 972 : IF (grad_deriv == 0) THEN
2156 30 : !$OMP DO
2157 : DO ii = 1, npoints
2158 314916 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2159 314916 : my_rhob = MAX(rhob(ii), 0.0_dp)
2160 314916 : my_tau_a = MAX(tau_a(ii), 0.0_dp)
2161 314916 : my_tau_b = MAX(tau_b(ii), 0.0_dp)
2162 314916 : IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN
2163 314916 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2164 314916 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2165 314916 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2166 314916 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2167 314916 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2168 314916 : sigmav(1, 1) = my_norm_drhoa**2
2169 314916 : sigmav(3, 1) = my_norm_drhob**2
2170 314916 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2171 314916 : tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp)
2172 314916 : tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp)
2173 314916 : tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2174 314916 : tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2175 314916 : laplace_rhov(1, 1) = laplace_rhoa(ii)
2176 314916 : laplace_rhov(2, 1) = laplace_rhob(ii)
2177 : CALL xc_f03_mgga_exc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2178 314916 : laplace_rhov(1, 1), tauv(1, 1), exc)
2179 314916 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2180 : END IF
2181 : END DO
2182 : !$OMP END DO
2183 : ELSE IF (grad_deriv == -1) THEN
2184 0 : !$OMP DO
2185 : DO ii = 1, npoints
2186 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2187 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
2188 0 : my_tau_a = MAX(tau_a(ii), 0.0_dp)
2189 0 : my_tau_b = MAX(tau_b(ii), 0.0_dp)
2190 0 : IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN
2191 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2192 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2193 0 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2194 0 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2195 0 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2196 0 : sigmav(1, 1) = my_norm_drhoa**2
2197 0 : sigmav(3, 1) = my_norm_drhob**2
2198 0 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2199 0 : laplace_rhov(1, 1) = laplace_rhoa(ii)
2200 0 : laplace_rhov(2, 1) = laplace_rhob(ii)
2201 0 : tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp)
2202 0 : tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp)
2203 0 : tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2204 0 : tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2205 : CALL xc_f03_mgga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2206 0 : laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1))
2207 0 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2208 0 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2209 0 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2210 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2211 0 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2212 : e_ndrhob(ii) = e_ndrhob(ii) + &
2213 0 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2214 0 : e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2215 0 : e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2216 0 : IF (has_laplace) THEN
2217 0 : e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2218 0 : e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2219 : END IF
2220 : END IF
2221 : END DO
2222 : !$OMP END DO
2223 : ELSE IF (grad_deriv == 1) THEN
2224 928 : !$OMP DO
2225 : DO ii = 1, npoints
2226 9065364 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2227 9065364 : my_rhob = MAX(rhob(ii), 0.0_dp)
2228 9065364 : my_tau_a = MAX(tau_a(ii), 0.0_dp)
2229 9065364 : my_tau_b = MAX(tau_b(ii), 0.0_dp)
2230 9065364 : IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN
2231 9054536 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2232 9054536 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2233 9054536 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2234 9054536 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2235 9054536 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2236 9054536 : sigmav(1, 1) = my_norm_drhoa**2
2237 9054536 : sigmav(3, 1) = my_norm_drhob**2
2238 9054536 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2239 9054536 : laplace_rhov(1, 1) = laplace_rhoa(ii)
2240 9054536 : laplace_rhov(2, 1) = laplace_rhob(ii)
2241 9054536 : tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp)
2242 9054536 : tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp)
2243 9054536 : tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2244 9054536 : tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2245 9054536 : IF (no_exc) THEN
2246 : CALL xc_f03_mgga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2247 : laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2248 0 : vlapl(1, 1), vtau(1, 1))
2249 0 : exc = 0.0_dp
2250 : ELSE
2251 : CALL xc_f03_mgga_exc_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2252 : laplace_rhov(1, 1), tauv(1, 1), exc, &
2253 9054536 : vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1))
2254 : END IF
2255 9054536 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2256 9054536 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2257 9054536 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2258 9054536 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2259 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2260 9054536 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2261 : e_ndrhob(ii) = e_ndrhob(ii) + &
2262 9054536 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2263 9054536 : e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2264 9054536 : e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2265 9054536 : IF (has_laplace) THEN
2266 1202688 : e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2267 1202688 : e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2268 : END IF
2269 : END IF
2270 : END DO
2271 : !$OMP END DO
2272 : ELSE IF (grad_deriv == -2) THEN
2273 0 : !$OMP DO
2274 : DO ii = 1, npoints
2275 0 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2276 0 : my_rhob = MAX(rhob(ii), 0.0_dp)
2277 0 : my_tau_a = MAX(tau_a(ii), 0.0_dp)
2278 0 : my_tau_b = MAX(tau_b(ii), 0.0_dp)
2279 0 : IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN
2280 0 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2281 0 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2282 0 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2283 0 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2284 0 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2285 0 : sigmav(1, 1) = my_norm_drhoa**2
2286 0 : sigmav(3, 1) = my_norm_drhob**2
2287 0 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2288 0 : laplace_rhov(1, 1) = laplace_rhoa(ii)
2289 0 : laplace_rhov(2, 1) = laplace_rhob(ii)
2290 0 : tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp)
2291 0 : tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp)
2292 0 : tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2293 0 : tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2294 0 : IF (no_exc) THEN
2295 : CALL xc_f03_mgga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2296 : laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2297 : vlapl(1, 1), vtau(1, 1), &
2298 : v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), &
2299 : v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2300 0 : v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2301 : ELSE
2302 : CALL xc_f03_mgga(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2303 : laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2304 : vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), &
2305 : v2rhotau(1, 1), v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2306 0 : v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2307 : END IF
2308 0 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2309 0 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2310 0 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2311 0 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2312 0 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2313 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2314 0 : sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2315 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2316 0 : sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2317 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2318 0 : sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2319 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2320 0 : sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2321 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2322 0 : sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2323 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2324 0 : sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2325 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2326 0 : sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2327 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2328 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2329 0 : 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2330 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2331 : sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2332 0 : 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2333 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2334 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2335 0 : 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2336 0 : e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1)
2337 0 : e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1)
2338 0 : e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1)
2339 0 : e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1)
2340 0 : e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho
2341 0 : e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho
2342 : e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2343 0 : sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa
2344 : e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2345 0 : sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa
2346 : e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2347 0 : sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob
2348 : e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2349 0 : sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob
2350 0 : e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1)
2351 0 : e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1)
2352 0 : e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1)
2353 0 : IF (has_laplace) THEN
2354 0 : e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1)
2355 0 : e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1)
2356 0 : e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1)
2357 0 : e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1)
2358 0 : e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho
2359 0 : e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho
2360 : e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2361 0 : sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa
2362 : e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2363 0 : sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa
2364 : e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2365 0 : sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob
2366 : e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2367 0 : sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob
2368 0 : e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1)
2369 0 : e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1)
2370 0 : e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1)
2371 0 : e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1)
2372 0 : e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1)
2373 0 : e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1)
2374 0 : e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1)
2375 : END IF
2376 : END IF
2377 : END DO
2378 : !$OMP END DO
2379 : ELSE IF (grad_deriv == 2) THEN
2380 14 : !$OMP DO
2381 : DO ii = 1, npoints
2382 96768 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
2383 96768 : my_rhob = MAX(rhob(ii), 0.0_dp)
2384 96768 : my_tau_a = MAX(tau_a(ii), 0.0_dp)
2385 96768 : my_tau_b = MAX(tau_b(ii), 0.0_dp)
2386 96768 : IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN
2387 96768 : rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp)
2388 96768 : rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp)
2389 96768 : my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp)
2390 96768 : my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp)
2391 96768 : my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp)
2392 96768 : sigmav(1, 1) = my_norm_drhoa**2
2393 96768 : sigmav(3, 1) = my_norm_drhob**2
2394 96768 : sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2395 96768 : laplace_rhov(1, 1) = laplace_rhoa(ii)
2396 96768 : laplace_rhov(2, 1) = laplace_rhob(ii)
2397 96768 : tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp)
2398 96768 : tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp)
2399 96768 : tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2400 96768 : tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2401 96768 : IF (no_exc) THEN
2402 : CALL xc_f03_mgga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2403 : laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2404 : vlapl(1, 1), vtau(1, 1), &
2405 : v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), &
2406 : v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2407 0 : v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2408 0 : exc = 0.0_dp
2409 : ELSE
2410 : CALL xc_f03_mgga(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2411 : laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2412 : vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), &
2413 : v2rhotau(1, 1), v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2414 96768 : v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2415 : END IF
2416 96768 : e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2417 96768 : e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2418 96768 : e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2419 96768 : e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2420 : e_ndrhoa(ii) = e_ndrhoa(ii) + &
2421 96768 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2422 : e_ndrhob(ii) = e_ndrhob(ii) + &
2423 96768 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2424 96768 : e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2425 96768 : e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2426 96768 : e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2427 96768 : e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2428 96768 : e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2429 96768 : e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2430 96768 : e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2431 : e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2432 96768 : sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2433 : e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2434 96768 : sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2435 : e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2436 96768 : sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2437 : e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2438 96768 : sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2439 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2440 96768 : sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2441 : e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2442 96768 : sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2443 : e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2444 96768 : sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2445 : e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2446 : sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2447 96768 : 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2448 : e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2449 : sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2450 96768 : 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2451 : e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2452 : sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2453 96768 : 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2454 96768 : e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1)
2455 96768 : e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1)
2456 96768 : e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1)
2457 96768 : e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1)
2458 96768 : e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho
2459 96768 : e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho
2460 : e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2461 96768 : sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa
2462 : e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2463 96768 : sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa
2464 : e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2465 96768 : sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob
2466 : e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2467 96768 : sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob
2468 96768 : e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1)
2469 96768 : e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1)
2470 96768 : e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1)
2471 96768 : IF (has_laplace) THEN
2472 41472 : e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2473 41472 : e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2474 41472 : e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1)
2475 41472 : e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1)
2476 41472 : e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1)
2477 41472 : e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1)
2478 41472 : e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho
2479 41472 : e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho
2480 : e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2481 41472 : sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa
2482 : e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2483 41472 : sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa
2484 : e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2485 41472 : sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob
2486 : e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2487 41472 : sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob
2488 41472 : e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1)
2489 41472 : e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1)
2490 41472 : e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1)
2491 41472 : e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1)
2492 41472 : e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1)
2493 41472 : e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1)
2494 41472 : e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1)
2495 : END IF
2496 : END IF
2497 : END DO
2498 : !$OMP END DO
2499 : END IF
2500 : CASE default
2501 2750 : CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.")
2502 : END SELECT
2503 :
2504 2750 : END SUBROUTINE libxc_spin_polarized_calc
2505 : #endif
2506 :
2507 : END MODULE xc_libxc
|
