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