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