(git:e7e05ae)
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 ! **************************************************************************************************
28 MODULE xc_libxc
29  USE bibliography, ONLY: lehtola2018, &
30  marques2012, &
31  cite_reference
32  USE input_section_types, ONLY: section_add_keyword, &
33  section_add_subsection, &
34  section_create, &
35  section_release, &
36  section_type, &
37  section_vals_type, &
38  section_vals_val_get
39  USE kinds, ONLY: default_string_length, &
40  dp
41  USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
42  xc_dset_get_derivative
43  USE xc_derivative_types, ONLY: xc_derivative_get, &
44  xc_derivative_type
45  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
46  USE xc_rho_set_types, ONLY: xc_rho_set_get, &
47  xc_rho_set_type
48 #if defined (__LIBXC)
49  USE input_keyword_types, ONLY: keyword_create, &
50  keyword_release, &
51  keyword_type
52  USE iso_c_binding, ONLY: c_size_t, c_int, c_double
53  USE xc_derivative_desc, ONLY: &
54  deriv_rho, deriv_rhoa, deriv_rhob, &
55  deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
56  deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob
57  USE xc_libxc_wrap, ONLY: xc_f03_func_t, &
58  xc_f03_func_init, &
59  xc_f03_func_end, &
60  xc_f03_func_info_t, &
61  xc_f03_functional_get_name, &
62  xc_f03_func_get_info, &
63  xc_f03_func_info_get_family, &
64  xc_f03_func_info_get_kind, &
65  xc_f03_func_info_get_n_ext_params, &
66  xc_f03_func_info_get_name, &
67  xc_f03_available_functional_numbers, &
68  xc_f03_available_functional_names, &
69  xc_f03_maximum_name_length, &
70  xc_f03_number_of_functionals, &
71  xc_f03_func_info_get_ext_params_name, &
72  xc_f03_func_info_get_ext_params_description, &
73  xc_f03_func_info_get_ext_params_default_value, &
74  xc_f03_gga_exc, &
75  xc_f03_gga_exc_vxc, &
76  xc_f03_gga_exc_vxc_fxc, &
77  xc_f03_gga_fxc, &
78  xc_f03_gga_vxc, &
79  xc_f03_gga_vxc_fxc, &
80  xc_f03_lda, &
81  xc_f03_lda_exc, &
82  xc_f03_lda_exc_vxc, &
83  xc_f03_lda_exc_vxc_fxc, &
84  xc_f03_lda_fxc, &
85  xc_f03_lda_kxc, &
86  xc_f03_lda_vxc, &
87  xc_f03_mgga, &
88  xc_f03_mgga_exc, &
89  xc_f03_mgga_exc_vxc, &
90  xc_f03_mgga_fxc, &
91  xc_f03_mgga_vxc, &
92  xc_f03_mgga_vxc_fxc, &
93  xc_polarized, &
94  xc_unpolarized, &
95  xc_family_lda, &
96  xc_family_gga, &
97  xc_family_mgga, &
98  xc_family_hyb_lda, &
99  xc_family_hyb_gga, &
100  xc_family_hyb_mgga, &
101  xc_correlation, &
102  xc_exchange, &
103  xc_exchange_correlation, &
104  xc_kinetic, &
105  xc_libxc_wrap_info_refs, &
106  xc_libxc_wrap_version, &
107  xc_libxc_wrap_functional_get_number, &
108  xc_libxc_wrap_needs_laplace, &
109  xc_libxc_wrap_functional_set_params, &
110  xc_libxc_wrap_is_under_development, &
111  xc_libxc_get_reference_length, &
112  xc_libxc_check_functional
113 #endif
114 
115 #include "../base/base_uses.f90"
116 
117  IMPLICIT NONE
118  PRIVATE
119 
120  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_libxc'
121 
122  PUBLIC :: libxc_lda_info, libxc_lda_eval, libxc_lsd_info, libxc_lsd_eval, &
123  libxc_version_info, libxc_get_reference_length, libxc_add_sections, &
124  libxc_check_existence_in_libxc
125 
126 #if defined (__LIBXC)
127  INTEGER(C_SIZE_T), PARAMETER, PRIVATE :: one = 1
128 #endif
129 
130 CONTAINS
131 
132 ! **************************************************************************************************
133 !> \brief This function checks whether a functional name belongs to LibXC
134 !> \param libxc_params (possible) LibXC input section
135 !> \return exists whether the functional exists in LibXC
136 ! **************************************************************************************************
137  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 
2461 END 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