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