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