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