(git:374b731)
Loading...
Searching...
No Matches
xc_ke_gga.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the several different kinetic energy functionals
10!> with a GGA form
11!> \par History
12!> JGH (26.02.2003) : OpenMP enabled
13!> fawzi (04.2004) : adapted to the new xc interface
14!> \author JGH (20.02.2002)
15! **************************************************************************************************
16MODULE xc_ke_gga
17
18 USE cp_array_utils, ONLY: cp_3d_r_cp_type
19 USE kinds, ONLY: dp
20 USE mathconstants, ONLY: pi
21 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
22 deriv_norm_drhoa,&
23 deriv_norm_drhob,&
24 deriv_rho,&
25 deriv_rhoa,&
26 deriv_rhob
27 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
28 xc_dset_get_derivative
29 USE xc_derivative_types, ONLY: xc_derivative_get,&
30 xc_derivative_type
31 USE xc_functionals_utilities, ONLY: calc_wave_vector,&
32 set_util
33 USE xc_input_constants, ONLY: ke_lc,&
34 ke_llp,&
35 ke_ol1,&
36 ke_ol2,&
37 ke_pbe,&
38 ke_pw86,&
39 ke_pw91,&
40 ke_t92
41 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
42 USE xc_rho_set_types, ONLY: xc_rho_set_get,&
43 xc_rho_set_type
44#include "../base/base_uses.f90"
45
46 IMPLICIT NONE
47
48 PRIVATE
49
50 REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
51 f23 = 2.0_dp*f13, &
52 f43 = 4.0_dp*f13, &
53 f53 = 5.0_dp*f13
54
55 PUBLIC :: ke_gga_info, ke_gga_lda_eval, ke_gga_lsd_eval
56
57 REAL(KIND=dp) :: cf, b, b_lda, b_lsd, flda, flsd, sfac, t13
58 REAL(KIND=dp) :: fact, tact
59 REAL(KIND=dp) :: eps_rho
60
61 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_ke_gga'
62
63CONTAINS
64
65! **************************************************************************************************
66!> \brief ...
67!> \param cutoff ...
68! **************************************************************************************************
69 SUBROUTINE ke_gga_init(cutoff)
70
71 REAL(KIND=dp), INTENT(IN) :: cutoff
72
73 eps_rho = cutoff
74 CALL set_util(cutoff)
75
76 cf = 0.3_dp*(3.0_dp*pi*pi)**f23
77 flda = cf
78 flsd = flda*2.0_dp**f23
79! the_factor 2^(1/3) for LDA is here
80 b_lda = 2.0_dp**f43*(3.0_dp*pi*pi)**(f13)
81 b_lsd = 2.0_dp*(3.0_dp*pi*pi)**(f13)
82 sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
83 t13 = 2.0_dp**f13
84
85 END SUBROUTINE ke_gga_init
86
87! **************************************************************************************************
88!> \brief ...
89!> \param functional ...
90!> \param lsd ...
91!> \param reference ...
92!> \param shortform ...
93!> \param needs ...
94!> \param max_deriv ...
95! **************************************************************************************************
96 SUBROUTINE ke_gga_info(functional, lsd, reference, shortform, needs, max_deriv)
97 INTEGER, INTENT(in) :: functional
98 LOGICAL, INTENT(in) :: lsd
99 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
100 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
101 INTEGER, INTENT(out), OPTIONAL :: max_deriv
102
103 IF (PRESENT(reference)) THEN
104 SELECT CASE (functional)
105 CASE (ke_ol1)
106 reference = "H. Ou-Yang and M. Levy, "// &
107 "Intl. J. Quant. Chem. 40, 379 (1991); Functional 1"
108 CASE (ke_ol2)
109 reference = "H. Ou-Yang and M. Levy, "// &
110 "Intl. J. Quant. Chem. 40, 379 (1991); Functional 2"
111 CASE (ke_llp)
112 reference = "H. Lee, C. Lee, R.G. Parr, Phys. Rev. A, 44, 768 (1991)"
113 CASE (ke_pw86)
114 reference = "J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
115 CASE (ke_pw91)
116 reference = "J.P. Perdew and Y. Wang, Electronic Structure of Solids 91"
117 CASE (ke_lc)
118 reference = "A. Lembarki and H. Chermette, Phys. Rev. A, 50, 5328 (1994)"
119 CASE (ke_t92)
120 reference = "A.J. Thakkar, Phys. Rev. A, 46, 6920 (1992)"
121 CASE (ke_pbe)
122 reference = "J.P.Perdew, K.Burke, M.Ernzerhof, Phys. Rev. Letter, 77, 3865 (1996)"
123 END SELECT
124 IF (.NOT. lsd) THEN
125 IF (len_trim(reference) + 6 < len(reference)) THEN
126 reference(len_trim(reference):len_trim(reference) + 6) = ' {spin unpolarized}'
127 END IF
128 END IF
129 END IF
130 IF (PRESENT(shortform)) THEN
131 SELECT CASE (functional)
132 CASE (ke_ol1)
133 shortform = "Ou-Yang-Levy Functional 1"
134 CASE (ke_ol2)
135 shortform = "Ou-Yang-Levy Functional 2"
136 CASE (ke_llp)
137 shortform = "Lee-Lee-Parr Functional"
138 CASE (ke_pw86)
139 shortform = "Perdew-Wang 1986 Functional (kinetic energy)"
140 CASE (ke_pw91)
141 shortform = "Perdew-Wang 1991 Functional (kinetic energy)"
142 CASE (ke_lc)
143 shortform = "Lembarki-Chermette kinetic energy functional"
144 CASE (ke_t92)
145 shortform = "Thakkar 1992 Functional"
146 CASE (ke_pbe)
147 shortform = "Perdew-Burke-Ernzerhof Functional (kinetic energy)"
148 END SELECT
149 IF (.NOT. lsd) THEN
150 IF (len_trim(shortform) + 6 < len(shortform)) THEN
151 shortform(len_trim(shortform):len_trim(shortform) + 6) = ' {spin unpolarized}'
152 END IF
153 END IF
154 END IF
155 IF (PRESENT(needs)) THEN
156 IF (lsd) THEN
157 needs%rho_spin = .true.
158 needs%rho_spin_1_3 = .true.
159 needs%norm_drho_spin = .true.
160 ELSE
161 needs%rho = .true.
162 needs%rho_1_3 = .true.
163 needs%norm_drho = .true.
164 END IF
165 END IF
166 IF (PRESENT(max_deriv)) max_deriv = 3
167
168 END SUBROUTINE ke_gga_info
169
170! **************************************************************************************************
171!> \brief ...
172!> \param functional ...
173!> \param rho_set ...
174!> \param deriv_set ...
175!> \param order ...
176! **************************************************************************************************
177 SUBROUTINE ke_gga_lda_eval(functional, rho_set, deriv_set, order)
178
179 INTEGER, INTENT(IN) :: functional
180 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
181 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
182 INTEGER, INTENT(IN) :: order
183
184 CHARACTER(LEN=*), PARAMETER :: routinen = 'ke_gga_lda_eval'
185
186 INTEGER :: handle, m, npoints
187 INTEGER, DIMENSION(2, 3) :: bo
188 REAL(kind=dp) :: drho_cutoff, rho_cutoff
189 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
190 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
191 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
192 e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
193 e_rho_rho_rho, grho, rho, rho13
194 TYPE(xc_derivative_type), POINTER :: deriv
195
196 CALL timeset(routinen, handle)
197 NULLIFY (rho, rho13, e_0, e_rho, e_ndrho, &
198 e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
199 e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho)
200 m = abs(order)
201
202 CALL xc_rho_set_get(rho_set, rho_1_3=rho13, rho=rho, &
203 norm_drho=grho, local_bounds=bo, rho_cutoff=rho_cutoff, &
204 drho_cutoff=drho_cutoff)
205 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
206 CALL ke_gga_init(rho_cutoff)
207
208 ALLOCATE (s(npoints))
209 ALLOCATE (fs(npoints, m + 1))
210
211! s = norm_drho/(rho^(4/3)*2*(pi*pi*3)^(1/3))
212 CALL calc_wave_vector("p", rho, grho, s)
213
214 fact = flda
215! Definition of s has changed
216 b = b_lda
217! tact = t13
218 tact = 1.0_dp
219
220 SELECT CASE (functional)
221 CASE (ke_ol1)
222 CALL efactor_ol1(s, fs, m)
223 cpabort("OL1 functional currently not working properly")
224 CASE (ke_ol2)
225 CALL efactor_ol2(s, fs, m)
226 cpabort("OL2 functional currently not working properly")
227 CASE (ke_llp)
228 CALL efactor_llp(s, fs, m)
229 CASE (ke_pw86)
230 CALL efactor_pw86(s, fs, m)
231 CASE (ke_pw91)
232 CALL efactor_pw91(s, fs, m, 1)
233 CASE (ke_lc)
234 CALL efactor_pw91(s, fs, m, 2)
235 CASE (ke_t92)
236 CALL efactor_t92(s, fs, m)
237 CASE (ke_pbe)
238 CALL efactor_pbex(s, fs, m, 1)
239 CASE DEFAULT
240 cpabort("Unsupported functional")
241 END SELECT
242
243 IF (order >= 0) THEN
244 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
245 allocate_deriv=.true.)
246 CALL xc_derivative_get(deriv, deriv_data=e_0)
247
248 CALL kex_p_0(rho, rho13, fs, e_0, npoints)
249 END IF
250
251 IF (order >= 1 .OR. order == -1) THEN
252 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
253 allocate_deriv=.true.)
254 CALL xc_derivative_get(deriv, deriv_data=e_rho)
255 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
256 allocate_deriv=.true.)
257 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
258
259 CALL kex_p_1(rho, rho13, s, fs, e_rho=e_rho, e_ndrho=e_ndrho, &
260 npoints=npoints)
261 END IF
262 IF (order >= 2 .OR. order == -2) THEN
263 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
264 allocate_deriv=.true.)
265 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
266 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
267 allocate_deriv=.true.)
268 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
269 deriv => xc_dset_get_derivative(deriv_set, &
270 [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
271 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
272
273 CALL kex_p_2(rho, rho13, s, fs, e_rho_rho=e_rho_rho, &
274 e_rho_ndrho=e_rho_ndrho, e_ndrho_ndrho=e_ndrho_ndrho, npoints=npoints)
275 END IF
276 IF (order >= 3 .OR. order == -3) THEN
277 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
278 allocate_deriv=.true.)
279 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
280 deriv => xc_dset_get_derivative(deriv_set, &
281 [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.true.)
282 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
283 deriv => xc_dset_get_derivative(deriv_set, &
284 [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
285 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
286 deriv => xc_dset_get_derivative(deriv_set, &
287 [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
288 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
289
290 CALL kex_p_3(rho, rho13, s, fs, e_rho_rho_rho=e_rho_rho_rho, &
291 e_rho_rho_ndrho=e_rho_rho_ndrho, e_rho_ndrho_ndrho=e_rho_ndrho_ndrho, &
292 e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, npoints=npoints)
293 END IF
294 IF (order > 3 .OR. order < -3) THEN
295 cpabort("derivatives bigger than 3 not implemented")
296 END IF
297
298 DEALLOCATE (s)
299 DEALLOCATE (fs)
300
301 CALL timestop(handle)
302
303 END SUBROUTINE ke_gga_lda_eval
304
305! **************************************************************************************************
306!> \brief ...
307!> \param functional ...
308!> \param rho_set ...
309!> \param deriv_set ...
310!> \param order ...
311! **************************************************************************************************
312 SUBROUTINE ke_gga_lsd_eval(functional, rho_set, deriv_set, order)
313
314 INTEGER, INTENT(IN) :: functional
315 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
316 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
317 INTEGER, INTENT(IN), OPTIONAL :: order
318
319 CHARACTER(len=*), PARAMETER :: routinen = 'ke_gga_lsd_eval'
320 INTEGER, DIMENSION(2), PARAMETER :: &
321 norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
322 rho_spin_name = [deriv_rhoa, deriv_rhob]
323
324 INTEGER :: handle, ispin, m, npoints
325 INTEGER, DIMENSION(2, 3) :: bo
326 REAL(kind=dp) :: rho_cutoff
327 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
328 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
329 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
330 e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
331 e_rho_rho_rho
332 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
333 TYPE(xc_derivative_type), POINTER :: deriv
334
335 CALL timeset(routinen, handle)
336 NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
337 e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
338 DO ispin = 1, 2
339 NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
340 END DO
341
342 CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
343 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
344 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
345 norm_drhob=norm_drho(2)%array, rho_cutoff=rho_cutoff, &
346 local_bounds=bo)
347 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
348 m = abs(order)
349 CALL ke_gga_init(rho_cutoff)
350
351 ALLOCATE (s(npoints))
352 ALLOCATE (fs(npoints, m + 1))
353
354 fact = flsd
355 b = b_lsd
356 tact = 1.0_dp
357
358 DO ispin = 1, 2
359
360 CALL calc_wave_vector("p", rho(ispin)%array, norm_drho(ispin)%array, s)
361
362 SELECT CASE (functional)
363 CASE (ke_ol1)
364 CALL efactor_ol1(s, fs, m)
365 CASE (ke_ol2)
366 CALL efactor_ol2(s, fs, m)
367 CASE (ke_llp)
368 CALL efactor_llp(s, fs, m)
369 CASE (ke_pw86)
370 tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
371 CALL efactor_pw86(s, fs, m, f2_lsd=tact)
372 tact = 1.0_dp
373 CASE (ke_pbe)
374 tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
375 CALL efactor_pbex(s, fs, m, 1, f2_lsd=tact)
376 tact = 1.0_dp
377 CASE (ke_pw91)
378 tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
379 CALL efactor_pw91(s, fs, m, 1, f2_lsd=tact)
380 tact = 1.0_dp
381 CASE (ke_lc)
382 tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
383 CALL efactor_pw91(s, fs, m, 2, f2_lsd=tact)
384 tact = 1.0_dp
385 CASE (ke_t92)
386 CALL efactor_t92(s, fs, m)
387 CASE DEFAULT
388 cpabort("Unsupported functional")
389 END SELECT
390
391 IF (order >= 0) THEN
392 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
393 allocate_deriv=.true.)
394 CALL xc_derivative_get(deriv, deriv_data=e_0)
395
396 CALL kex_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, &
397 e_0=e_0, npoints=npoints)
398 END IF
399 IF (order >= 1 .OR. order == -1) THEN
400 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
401 allocate_deriv=.true.)
402 CALL xc_derivative_get(deriv, deriv_data=e_rho)
403 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
404 allocate_deriv=.true.)
405 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
406
407 CALL kex_p_1(rho=rho(ispin)%array, &
408 r13=rho_1_3(ispin)%array, s=s, fs=fs, e_rho=e_rho, &
409 e_ndrho=e_ndrho, npoints=npoints)
410 END IF
411 IF (order >= 2 .OR. order == -2) THEN
412 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
413 rho_spin_name(ispin)], allocate_deriv=.true.)
414 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
415 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
416 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
417 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
418 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
419 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
420 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
421
422 CALL kex_p_2(rho(ispin)%array, rho_1_3(ispin)%array, &
423 s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, npoints)
424 END IF
425 IF (order >= 3 .OR. order == -3) THEN
426 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
427 rho_spin_name(ispin), rho_spin_name(ispin)], &
428 allocate_deriv=.true.)
429 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
430 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
431 rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
432 allocate_deriv=.true.)
433 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
434 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
435 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
436 allocate_deriv=.true.)
437 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
438 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
439 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
440 allocate_deriv=.true.)
441 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
442
443 CALL kex_p_3(rho(ispin)%array, &
444 rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
445 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
446 END IF
447 IF (order > 3 .OR. order < -3) THEN
448 cpabort("derivatives bigger than 3 not implemented")
449 END IF
450
451 END DO
452
453 DEALLOCATE (s)
454 DEALLOCATE (fs)
455 CALL timestop(handle)
456 END SUBROUTINE ke_gga_lsd_eval
457
458! **************************************************************************************************
459!> \brief ...
460!> \param rho ...
461!> \param r13 ...
462!> \param fs ...
463!> \param e_0 ...
464!> \param npoints ...
465! **************************************************************************************************
466 SUBROUTINE kex_p_0(rho, r13, fs, e_0, npoints)
467
468 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
469 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
470 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
471 INTEGER, INTENT(in) :: npoints
472
473 INTEGER :: ip
474
475!$OMP PARALLEL DO DEFAULT(NONE) &
476!$OMP SHARED(npoints, rho, e_0, fact, r13, fs, eps_rho) &
477!$OMP PRIVATE(ip)
478
479 DO ip = 1, npoints
480
481 IF (rho(ip) > eps_rho) THEN
482 e_0(ip) = e_0(ip) + fact*r13(ip)*r13(ip)*rho(ip)*fs(ip, 1)
483 END IF
484
485 END DO
486
487!$OMP END PARALLEL DO
488
489 END SUBROUTINE kex_p_0
490
491! **************************************************************************************************
492!> \brief ...
493!> \param rho ...
494!> \param r13 ...
495!> \param s ...
496!> \param fs ...
497!> \param e_rho ...
498!> \param e_ndrho ...
499!> \param npoints ...
500! **************************************************************************************************
501 SUBROUTINE kex_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
502
503 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
504 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
505 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
506 INTEGER, INTENT(in) :: npoints
507
508 INTEGER :: ip
509 REAL(kind=dp) :: a0, a1, sx, sy
510
511!$OMP PARALLEL DO DEFAULT(NONE) &
512!$OMP SHARED(npoints, rho, eps_rho, fact, r13, sfac, tact) &
513!$OMP SHARED(fs, e_rho, e_ndrho, s) &
514!$OMP PRIVATE(ip,a0,a1,sx,sy)
515
516 DO ip = 1, npoints
517
518 IF (rho(ip) > eps_rho) THEN
519
520 a0 = fact*r13(ip)*r13(ip)*rho(ip)
521 a1 = f53*fact*r13(ip)*r13(ip)
522 sx = -f43*s(ip)/rho(ip)
523 sy = sfac*tact/(r13(ip)*rho(ip))
524 e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
525 e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
526
527 END IF
528
529 END DO
530
531!$OMP END PARALLEL DO
532
533 END SUBROUTINE kex_p_1
534
535! **************************************************************************************************
536!> \brief ...
537!> \param rho ...
538!> \param r13 ...
539!> \param s ...
540!> \param fs ...
541!> \param e_rho_rho ...
542!> \param e_rho_ndrho ...
543!> \param e_ndrho_ndrho ...
544!> \param npoints ...
545! **************************************************************************************************
546 SUBROUTINE kex_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
547 npoints)
548
549 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
550 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
551 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
552 INTEGER, INTENT(in) :: npoints
553
554 INTEGER :: ip
555 REAL(kind=dp) :: a0, a1, a2, sx, sxx, sxy, sy
556
557!$OMP PARALLEL DO DEFAULT(NONE) &
558!$OMP SHARED (npoints, rho, eps_rho, fact, r13) &
559!$OMP SHARED (e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, fs) &
560!$OMP SHARED(s, sfac, tact) &
561!$OMP PRIVATE(ip,a0,a1,a2,sx,sy,sxx,sxy)
562
563 DO ip = 1, npoints
564
565 IF (rho(ip) > eps_rho) THEN
566
567 a0 = fact*r13(ip)*r13(ip)*rho(ip)
568 a1 = f53*fact*r13(ip)*r13(ip)
569 a2 = f23*f53*fact/r13(ip)
570 sx = -f43*s(ip)/rho(ip)
571 sy = sfac*tact/(r13(ip)*rho(ip))
572 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
573 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
574 e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
575 a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
576 e_rho_ndrho(ip) = e_rho_ndrho(ip) + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + &
577 a0*fs(ip, 2)*sxy
578 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
579
580 END IF
581
582 END DO
583
584!$OMP END PARALLEL DO
585
586 END SUBROUTINE kex_p_2
587
588! **************************************************************************************************
589!> \brief ...
590!> \param rho ...
591!> \param r13 ...
592!> \param s ...
593!> \param fs ...
594!> \param e_rho_rho_rho ...
595!> \param e_rho_rho_ndrho ...
596!> \param e_rho_ndrho_ndrho ...
597!> \param e_ndrho_ndrho_ndrho ...
598!> \param npoints ...
599! **************************************************************************************************
600 SUBROUTINE kex_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
601 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
602
603 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
604 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
605 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_rho_rho_rho, e_rho_rho_ndrho, &
606 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
607 INTEGER, INTENT(in) :: npoints
608
609 INTEGER :: ip
610 REAL(kind=dp) :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
611 sxy, sy
612
613!$OMP PARALLEL DO DEFAULT(NONE) &
614!$OMP SHARED(npoints, rho, eps_rho, fact, r13) &
615!$OMP SHARED(s, sfac, tact, fs, e_rho_rho_rho) &
616!$OMP SHARED(e_rho_rho_ndrho, e_rho_ndrho_ndrho) &
617!$OMP SHARED(e_ndrho_ndrho_ndrho) &
618!$OMP PRIVATE(ip,a0,a1,a2,a3,sx,sy,sxx,sxy,sxxx,sxxy)
619
620 DO ip = 1, npoints
621
622 IF (rho(ip) > eps_rho) THEN
623
624 a0 = fact*r13(ip)*r13(ip)*rho(ip)
625 a1 = f53*fact*r13(ip)*r13(ip)
626 a2 = f23*f53*fact/r13(ip)
627 a3 = -f13*f23*f53*fact/(r13(ip)*rho(ip))
628 sx = -f43*s(ip)/rho(ip)
629 sy = sfac*tact/(r13(ip)*rho(ip))
630 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
631 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
632 sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
633 sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
634 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
635 3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
636 a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
637 a0*fs(ip, 2)*sxxx
638 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
639 2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
640 2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
641 a0*fs(ip, 2)*sxxy
642 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
643 2.0_dp*a0*fs(ip, 3)*sxy*sy
644 e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) + a0*fs(ip, 4)*sy*sy*sy
645
646 END IF
647
648 END DO
649
650!$OMP END PARALLEL DO
651
652 END SUBROUTINE kex_p_3
653
654! Enhancement Factors
655! **************************************************************************************************
656!> \brief ...
657!> \param s ...
658!> \param fs ...
659!> \param m ...
660! **************************************************************************************************
661 SUBROUTINE efactor_ol1(s, fs, m)
662 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
663 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
664 INTEGER, INTENT(IN) :: m
665
666 INTEGER :: ip
667 REAL(kind=dp) :: t1, t2
668
669 t1 = b*b/(72.0_dp*cf)
670 t2 = 0.001878_dp*b
671
672!$OMP PARALLEL DO DEFAULT(NONE) &
673!$OMP SHARED(s, m, fs, t1, t2) &
674!$OMP PRIVATE(ip)
675
676 DO ip = 1, SIZE(s)
677 SELECT CASE (m)
678 CASE (0)
679 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
680 CASE (1)
681 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
682 fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
683 CASE (2:3)
684 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
685 fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
686 fs(ip, 3) = 2.0_dp*t1
687 CASE DEFAULT
688 cpabort("Illegal order.")
689 END SELECT
690 END DO
691
692!$OMP END PARALLEL DO
693
694 IF (m == 3) fs(:, 4) = 0.0_dp
695
696 END SUBROUTINE efactor_ol1
697! **************************************************************************************************
698!> \brief ...
699!> \param s ...
700!> \param fs ...
701!> \param m ...
702! **************************************************************************************************
703 SUBROUTINE efactor_ol2(s, fs, m)
704 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
705 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
706 INTEGER, INTENT(IN) :: m
707
708 INTEGER :: ip
709 REAL(kind=dp) :: t1, t2, t3, y
710
711 t1 = b*b/(72.0_dp*cf)
712 t2 = 0.0245_dp*b
713 t3 = 2.0_dp**f53*b
714
715!$OMP PARALLEL DO DEFAULT(NONE) &
716!$OMP SHARED(s, m, t1, t2, t3, fs) &
717!$OMP PRIVATE(ip,y)
718
719 DO ip = 1, SIZE(s)
720 y = 1.0_dp/(1.0_dp + t3*s(ip))
721 SELECT CASE (m)
722 CASE (0)
723 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
724 CASE (1)
725 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
726 fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
727 CASE (2)
728 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
729 fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
730 fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
731 CASE (3)
732 fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
733 fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
734 fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
735 fs(ip, 4) = 6.0_dp*t2*t3*t3*y*y*y*y
736 CASE DEFAULT
737 cpabort("Illegal order.")
738 END SELECT
739 END DO
740
741!$OMP END PARALLEL DO
742
743 END SUBROUTINE efactor_ol2
744! **************************************************************************************************
745!> \brief ...
746!> \param s ...
747!> \param fs ...
748!> \param m ...
749! **************************************************************************************************
750 SUBROUTINE efactor_llp(s, fs, m)
751 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
752 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
753 INTEGER, INTENT(IN) :: m
754
755 INTEGER :: ip
756 REAL(kind=dp) :: as, bs, p, q, sas, sbs, t1, t10, t11, t12, t133, t16, t17, t19, t2, t20, &
757 t22, t23, t26, t28, t29, t3, t30, t33, t36, t39, t4, t40, t42, t43, t45, t46, t47, t49, &
758 t5, t50, t54, t55, t7, t71, t8, t9, x, ys
759
760 p = 0.0044188_dp*b*b
761 q = 0.0253_dp*b
762
763!$OMP PARALLEL DO DEFAULT(NONE) &
764!$OMP SHARED(s, m, fs, p, q, b) &
765!$OMP PRIVATE(ip,x,bs,sbs,as,sas,ys,t2,t4,t5,t8,t9,t10,t12) &
766!$OMP PRIVATE(t1,t3,t7,t11,t16,t17,t20,t22,t23,t26,t30,t33) &
767!$OMP PRIVATE(t40,t42,t43,t45,t46,t47,t49,t50,t71,t133) &
768!$OMP PRIVATE(t19,t28,t29,t36,t39,t54,t55)
769
770 DO ip = 1, SIZE(s)
771 x = s(ip)
772 bs = b*x
773 sbs = sqrt(bs*bs + 1.0_dp)
774 as = log(bs + sbs)
775 sas = x*as
776 ys = 1.0_dp/(1.0_dp + q*sas)
777 SELECT CASE (m)
778 CASE (0)
779 fs(ip, 1) = 1.0_dp + p*x*x*ys
780 CASE (1)
781 fs(ip, 1) = 1.0_dp + p*x*x*ys
782 t2 = q*x
783 t4 = b**2
784 t5 = x**2
785 t8 = sqrt(1.0_dp + t4*t5)
786 t9 = b*x + t8
787 t10 = log(t9)
788 t12 = 1.0_dp + t2*t10
789 t17 = t12**2
790 fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
791
792 CASE (2)
793 fs(ip, 1) = 1.0_dp + p*x*x*ys
794 ! first der
795 t2 = q*x
796 t4 = b**2
797 t5 = x**2
798 t8 = sqrt(1.0_dp + t4*t5)
799 t9 = b*x + t8
800 t10 = log(t9)
801 t12 = 1.0_dp + t2*t10
802 t17 = t12**2
803 fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
804
805 ! second der
806 t1 = q*x
807 t3 = b**2
808 t4 = x**2
809 t7 = sqrt(1.0_dp + t3*t4)
810 t8 = b*x + t7
811 t9 = log(t8)
812 t11 = 1.0_dp + t1*t9
813 t16 = t11**2
814 t17 = 1.0_dp/t16
815 t20 = 1.0_dp/t7*t3
816 t22 = b + t20*x
817 t23 = 1/t8
818 t26 = q*t9 + t1*t22*t23
819 t30 = p*t4
820 t33 = t26**2
821 t40 = t7**2
822 t43 = t3**2
823 t49 = t22**2
824 t50 = t8**2
825 fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
826 t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
827 (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
828
829 CASE (3)
830
831 fs(ip, 1) = 1.0_dp + p*x*x*ys
832 ! first der
833 t2 = q*x
834 t4 = b**2
835 t5 = x**2
836 t8 = sqrt(1.0_dp + t4*t5)
837 t9 = b*x + t8
838 t10 = log(t9)
839 t12 = 1.0_dp + t2*t10
840 t17 = t12**2
841 fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
842
843 ! second der
844 t1 = q*x
845 t3 = b**2
846 t4 = x**2
847 t7 = sqrt(1.0_dp + t3*t4)
848 t8 = b*x + t7
849 t9 = log(t8)
850 t11 = 1.0_dp + t1*t9
851 t16 = t11**2
852 t17 = 1.0_dp/t16
853 t20 = 1.0_dp/t7*t3
854 t22 = b + t20*x
855 t23 = 1/t8
856 t26 = q*t9 + t1*t22*t23
857 t30 = p*t4
858 t33 = t26**2
859 t40 = t7**2
860 t43 = t3**2
861 t49 = t22**2
862 t50 = t8**2
863 fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
864 t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
865 (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
866
867 t1 = q*x
868 t3 = b**2
869 t4 = x**2
870 t7 = sqrt(1 + t3*t4)
871 t8 = b*x + t7
872 t9 = log(t8)
873 t11 = 1.0_dp + t1*t9
874 t12 = t11**2
875 t133 = 1.0_dp/t12
876 t17 = 1.0_dp/t7*t3
877 t19 = b + t17*x
878 t20 = 1.0_dp/t8
879 t23 = q*t9 + t1*t19*t20
880 t26 = p*x
881 t28 = 1.0_dp/t12/t11
882 t29 = t23**2
883 t36 = t7**2
884 t39 = t3**2
885 t40 = 1.0_dp/t36/t7*t39
886 t42 = -t40*t4 + t17
887 t45 = t19**2
888 t46 = t8**2
889 t47 = 1.0_dp/t46
890 t50 = 2.0_dp*q*t19*t20 + t1*t42*t20 - t1*t45*t47
891 t54 = p*t4
892 t55 = t12**2
893 t71 = t36**2
894 fs(ip, 4) = &
895 -6.0_dp*p*t133*t23 + 12.0_dp*t26*t28*t29 - &
896 6.0_dp*t26*t133*t50 - 6.0_dp*t54/t55*t29*t23 + &
897 6.0_dp*t54*t28*t23*t50 - t54*t133* &
898 (3.0_dp*q*t42*t20 - 3.0_dp*q*t45*t47 + 3.0_dp*t1* &
899 (1.0_dp/t71/t7*t39*t3*t4*x - t40*x)*t20 - &
900 3.0_dp*t1*t42*t47*t19 + 2.0_dp*t1*t45*t19/t46/t8)
901
902 CASE DEFAULT
903 cpabort("Illegal order.")
904 END SELECT
905 END DO
906
907!$OMP END PARALLEL DO
908
909 END SUBROUTINE efactor_llp
910! **************************************************************************************************
911!> \brief ...
912!> \param s ...
913!> \param fs ...
914!> \param m ...
915!> \param f2_lsd ...
916! **************************************************************************************************
917 SUBROUTINE efactor_pw86(s, fs, m, f2_lsd)
918 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
919 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
920 INTEGER, INTENT(IN) :: m
921 REAL(dp), OPTIONAL :: f2_lsd
922
923 INTEGER :: ip
924 REAL(kind=dp) :: f15, ff, p0, p1, p15, p2, p3, s1, s2, &
925 s4, s6, t1, t2, t3
926
927 t1 = 1.296_dp
928 t2 = 14.0_dp
929 t3 = 0.2_dp
930 f15 = 1.0_dp/15.0_dp
931 ff = 1.0_dp
932 IF (PRESENT(f2_lsd)) ff = f2_lsd
933
934!$OMP PARALLEL DO DEFAULT(NONE) &
935!$OMP SHARED(s, fs, m, t1, t2, t3, f15, ff) &
936!$OMP PRIVATE(ip, s1, s2, s4, s6, p0, p1, p2, p3, p15)
937
938 DO ip = 1, SIZE(s)
939 s1 = s(ip)*ff
940 s2 = s1*s1
941 s4 = s2*s2
942 s6 = s2*s4
943 SELECT CASE (m)
944 CASE (0)
945 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
946 fs(ip, 1) = p0**f15
947 CASE (1)
948 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
949 p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
950 p15 = p0**f15
951 fs(ip, 1) = p15
952 fs(ip, 2) = f15*p1*p15/p0
953 CASE (2)
954 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
955 p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
956 p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
957 p15 = p0**f15
958 fs(ip, 1) = p15
959 fs(ip, 2) = f15*p1*p15/p0
960 fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
961 CASE (3)
962 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
963 p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
964 p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
965 p3 = s1*ff*ff*ff*(24.0_dp*t2 + 120.0_dp*t3*s2)
966 p15 = p0**f15
967 fs(ip, 1) = p15
968 fs(ip, 2) = f15*p1*p15/p0
969 fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
970 fs(ip, 4) = f15*p15/p0*(-14.0_dp*f15*p1*p1/p0 + 14.0_dp*14.0_dp*f15*p1*p1*p1/p0/p0 + &
971 p3 - 14.0_dp*p2*p1/p0 + 14.0_dp*p1*p1/p0/p0)
972 CASE DEFAULT
973 cpabort("Illegal order.")
974 END SELECT
975 END DO
976
977!$OMP END PARALLEL DO
978
979 END SUBROUTINE efactor_pw86
980! **************************************************************************************************
981!> \brief ...
982!> \param s ...
983!> \param fs ...
984!> \param m ...
985! **************************************************************************************************
986 SUBROUTINE efactor_t92(s, fs, m)
987 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
988 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
989 INTEGER, INTENT(IN) :: m
990
991 INTEGER :: ip
992 REAL(kind=dp) :: a1, a2, as, asp, asp2, asp3, bs, p, q, &
993 s1, s2, sas, sbs, sbs3, sbs5, t0, w1, &
994 x, ys
995
996 p = 0.0055_dp*b*b
997 q = 0.0253_dp*b
998 a1 = 0.072_dp*b
999 a2 = 2.0_dp**f53*b
1000
1001!$OMP PARALLEL DO DEFAULT(NONE) &
1002!$OMP SHARED(s, fs, m, p, q, a1, a2, b) &
1003!$OMP PRIVATE(ip, x, bs, sbs, sas, ys, asp, sbs3, asp2, sbs5) &
1004!$OMP PRIVATE(asp3, as, s2, s1, t0, w1)
1005
1006 DO ip = 1, SIZE(s)
1007 x = s(ip)
1008 bs = b*x
1009 sbs = sqrt(bs*bs + 1.0_dp)
1010 as = log(bs + sbs)
1011 sas = x*as
1012 ys = 1.0_dp/(1.0_dp + q*sas)
1013 SELECT CASE (m)
1014 CASE (0)
1015 fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1016 CASE (1)
1017 asp = as + bs/sbs
1018 fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1019 fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1020 CASE (2)
1021 asp = as + bs/sbs
1022 sbs3 = sbs*sbs*sbs
1023 asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
1024 fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1025 fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1026 fs(ip, 3) = 2.0_dp*p*ys - p*q*x*(4.0_dp*asp + x*asp2)*ys*ys + &
1027 2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
1028 CASE (3)
1029 asp = as + bs/sbs
1030 sbs3 = sbs*sbs*sbs
1031 sbs5 = sbs3*sbs*sbs
1032 asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
1033 asp3 = -4.0_dp*b*b*bs/sbs3 + 3.0_dp*b*b*bs*bs*bs/sbs5
1034 w1 = (4.0_dp*asp + x*asp2)
1035 fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1036 fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1037 fs(ip, 3) = 2.0_dp*p*ys - p*q*x*w1*ys*ys + &
1038 2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
1039
1040 s2 = -6*p/(1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**2*(q*log(b*x + sqrt(1 + b**2*x**2)) + &
1041 q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)/(b*x + sqrt(1 + b**2*x**2))) + 12*p*x/ &
1042 (1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**3*(q*log(b*x + sqrt(1 + b**2*x**2)) + &
1043 q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)/(b*x + sqrt(1 + b**2*x**2)))**2
1044 s1 = s2 - 6*p*x/(1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**2*(2*q*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)/ &
1045 (b*x + sqrt(1 + b**2*x**2)) + q*x*(-1/sqrt(1 + b**2*x**2)**3*b**4*x**2 + 1/sqrt(1 + b**2*x**2)*b**2)/ &
1046 (b*x + sqrt(1 + b**2*x**2)) - q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)**2/ &
1047 (b*x + sqrt(1 + b**2*x**2))**2) - 6*p*x**2/(1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**4 &
1048 *(q*log(b*x + sqrt(1 + b**2*x**2)) + q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)/(b*x + sqrt(1 + b**2*x**2)))**3
1049 s2 = s1 + 6*p*x**2/(1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**3*(q*log(b*x + sqrt(1 + b**2*x**2)) + &
1050 q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)/(b*x + sqrt(1 + b**2*x**2)))*(2*q*(b + 1/sqrt(1 + b**2*x**2)*b**2*x) &
1051 /(b*x + sqrt(1 + b**2*x**2)) + q*x*(-1/sqrt(1 + b**2*x**2)**3*b**4*x**2 + 1/sqrt(1 + b**2*x**2)* &
1052 b**2)/(b*x + sqrt(1 + b**2*x**2)) - q*x*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)**2/(b*x + sqrt(1 + b**2*x**2))**2)
1053 t0 = s2 - p*x**2/(1 + q*x*log(b*x + sqrt(1 + b**2*x**2)))**2*(3*q*(-1/sqrt(1 + b**2*x**2)**3*b**4*x**2 + &
1054 1/sqrt(1 + b**2*x**2)*b**2)/(b*x + sqrt(1 + b**2*x**2)) - 3*q*(b + 1/sqrt(1 + b**2*x**2)*b**2*x)**2/ &
1055 (b*x + sqrt(1 + b**2*x**2))**2 + q*x*(3/sqrt(1 + b**2*x**2)**5*b**6*x**3 - 3/sqrt(1 + b**2*x**2)**3*b**4*x)/ &
1056 (b*x + sqrt(1 + b**2*x**2)) - 3*q*x*(-1/sqrt(1 + b**2*x**2)**3*b**4*x**2 + 1/sqrt(1 + b**2*x**2)*b**2)/ &
1057 (b*x + sqrt(1 + b**2*x**2))**2*(b + 1/sqrt(1 + b**2*x**2)*b**2*x) + 2*q*x*(b + 1/sqrt(1 + b**2*x**2)* &
1058 b**2*x)**3/(b*x + sqrt(1 + b**2*x**2))**3) - 6*a1/(1 + a2*x)**3*a2**2 + 6*a1*x/(1 + a2*x)**4*a2**3
1059
1060 fs(ip, 4) = t0
1061
1062 CASE DEFAULT
1063 cpabort("Illegal order")
1064 END SELECT
1065 END DO
1066
1067!$OMP END PARALLEL DO
1068
1069 END SUBROUTINE efactor_t92
1070! **************************************************************************************************
1071!> \brief ...
1072!> \param s ...
1073!> \param fs ...
1074!> \param m ...
1075!> \param pset ...
1076!> \param f2_lsd ...
1077! **************************************************************************************************
1078 SUBROUTINE efactor_pbex(s, fs, m, pset, f2_lsd)
1079
1080 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
1081 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1082 INTEGER, INTENT(IN) :: m, pset
1083 REAL(dp), OPTIONAL :: f2_lsd
1084
1085 REAL(kind=dp), PARAMETER :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
1086 mu = 0.2195149727645171_dp
1087
1088 INTEGER :: ip
1089 REAL(kind=dp) :: f0, mk, x, x2, y
1090
1091 IF (pset == 1) mk = mu/kappa1
1092 IF (pset == 2) mk = mu/kappa2
1093
1094 f0 = 1.0_dp/tact
1095 IF (PRESENT(f2_lsd)) f0 = f2_lsd
1096
1097!$OMP PARALLEL DO DEFAULT(NONE) &
1098!$OMP SHARED(s, m, fs, f0, mk) &
1099!$OMP PRIVATE(ip,x,x2,y)
1100
1101 DO ip = 1, SIZE(s)
1102 x = s(ip)*f0
1103 x2 = x*x
1104 y = 1.0_dp/(1.0_dp + mk*x2)
1105 SELECT CASE (m)
1106 CASE (0)
1107 fs(ip, 1) = 1.0_dp + mu*x2*y
1108 CASE (1)
1109 fs(ip, 1) = 1.0_dp + mu*x2*y
1110 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1111 CASE (2)
1112 fs(ip, 1) = 1.0_dp + mu*x2*y
1113 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1114 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1115 CASE (3)
1116 fs(ip, 1) = 1.0_dp + mu*x2*y
1117 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1118 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1119 fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
1120 CASE DEFAULT
1121 cpabort("Illegal order.")
1122 END SELECT
1123 END DO
1124
1125!$OMP END PARALLEL DO
1126
1127 END SUBROUTINE efactor_pbex
1128
1129! **************************************************************************************************
1130!> \brief ...
1131!> \param s ...
1132!> \param fs ...
1133!> \param m ...
1134!> \param pset ...
1135!> \param f2_lsd ...
1136! **************************************************************************************************
1137 SUBROUTINE efactor_pw91(s, fs, m, pset, f2_lsd)
1138 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
1139 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1140 INTEGER, INTENT(IN) :: m, pset
1141 REAL(dp), OPTIONAL :: f2_lsd
1142
1143 INTEGER :: ip
1144 REAL(dp) :: ff, t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
1145 t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
1146 t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
1147 t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
1148 t96, t98
1149 REAL(kind=dp) :: a1, a2, a3, a4, a5, bb, o, pa(6, 2), x
1150
1151! parameter set 1: Perdew-Wang
1152! parameter set 2: Lembarki-Chermette
1153
1154 pa(1:6, 1) = (/0.19645_dp, 0.2743_dp, &
1155 0.1508_dp, 100.0_dp, &
1156 7.7956_dp, 0.004_dp/)
1157 pa(1:6, 2) = (/0.093907_dp, 0.26608_dp, &
1158 0.0809615_dp, 100.0_dp, &
1159 76.320_dp, 0.57767e-4_dp/)
1160 o = 1.0_dp
1161 ff = 1.0_dp
1162 IF (PRESENT(f2_lsd)) ff = f2_lsd
1163 a1 = pa(1, pset)*ff
1164 a2 = pa(2, pset)*ff*ff
1165 a3 = pa(3, pset)*ff*ff
1166 a4 = pa(4, pset)*ff*ff
1167 bb = pa(5, pset)*ff
1168! it should be valid also for lsd
1169 a5 = pa(6, pset)*ff*ff*ff*ff
1170
1171!$OMP PARALLEL DEFAULT(NONE) &
1172!$OMP SHARED(s, fs, m, a1, a2, a3, a4, a5, bb, pa, o, ff) &
1173!$OMP PRIVATE(x,t1,t10,t101,t106,t109,t111,t113,t119,t12,t123) &
1174!$OMP PRIVATE(t124,t13,t14,t15,t16,t17,t18,t19,t2,t20,t21,t22) &
1175!$OMP PRIVATE(t23,t25,t26,t27,t28,t29,t3,t30,t31,t33,t35,t37) &
1176!$OMP PRIVATE(t38,t39,t4,t40,t44,t47,t48,t5,t50,t51,t53,t55) &
1177!$OMP PRIVATE(t56,t57,t58,t59,t6,t60,t64,t65,t69,t7,t70,t71) &
1178!$OMP PRIVATE(t73,t77,t78,t8,t80,t82,t9,t90,t93,t94,t96,t98,ip)
1179
1180 IF (m >= 0) THEN
1181!$OMP DO
1182 DO ip = 1, SIZE(s)
1183 x = s(ip)
1184 t3 = bb**2
1185 t4 = x**2
1186 t7 = sqrt(o + t3*t4)
1187 t9 = log(bb*x + t7)
1188 t10 = a1*x*t9
1189 t12 = exp(-a4*t4)
1190 t17 = t4**2
1191 fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
1192 END DO
1193!$OMP END DO
1194 END IF
1195 IF (m >= 1) THEN
1196!$OMP DO
1197 DO ip = 1, SIZE(s)
1198 x = s(ip)
1199 t2 = bb**2
1200 t3 = x**2
1201 t6 = sqrt(o + t2*t3)
1202 t7 = bb*x + t6
1203 t8 = log(t7)
1204 t9 = a1*t8
1205 t10 = a1*x
1206 t17 = t10*(bb + 1/t6*t2*x)/t7
1207 t19 = t3*x
1208 t21 = exp(-a4*t3)
1209 t26 = a2 - a3*t21
1210 t30 = t10*t8
1211 t31 = t3**2
1212 t33 = o + t30 + a5*t31
1213 t38 = t33**2
1214 fs(ip, 2) = &
1215 (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
1216 t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
1217 END DO
1218!$OMP END DO
1219 END IF
1220 IF (m >= 2) THEN
1221!$OMP DO
1222 DO ip = 1, SIZE(s)
1223 x = s(ip)
1224 t1 = bb**2
1225 t2 = x**2
1226 t5 = sqrt(o + t1*t2)
1227 t7 = o/t5*t1
1228 t9 = bb + t7*x
1229 t12 = bb*x + t5
1230 t13 = o/t12
1231 t15 = 2._dp*a1*t9*t13
1232 t16 = a1*x
1233 t17 = t5**2
1234 t20 = t1**2
1235 t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
1236 t26 = t9**2
1237 t27 = t12**2
1238 t30 = t16*t26/t27
1239 t31 = a3*a4
1240 t33 = exp(-a4*t2)
1241 t37 = a4**2
1242 t39 = t2**2
1243 t44 = a3*t33
1244 t47 = log(t12)
1245 t48 = t16*t47
1246 t50 = o + t48 + a5*t39
1247 t53 = a1*t47
1248 t55 = t16*t9*t13
1249 t56 = t2*x
1250 t60 = a2 - t44
1251 t64 = t50**2
1252 t65 = o/t64
1253 t69 = t53 + t55 + 4._dp*a5*t56
1254 t73 = o + t48 + t60*t2
1255 t77 = t69**2
1256 fs(ip, 3) = &
1257 (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
1258 2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
1259 (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
1260 t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
1261 END DO
1262!$OMP END DO
1263 END IF
1264 IF (m >= 3) THEN
1265!$OMP DO
1266 DO ip = 1, SIZE(s)
1267 x = s(ip)
1268 t1 = bb**2
1269 t2 = x**2
1270 t5 = sqrt(0.1e1_dp + t1*t2)
1271 t6 = t5**2
1272 t9 = t1**2
1273 t10 = 1/t6/t5*t9
1274 t13 = 1/t5*t1
1275 t14 = -t10*t2 + t13
1276 t17 = bb*x + t5
1277 t18 = 1/t17
1278 t20 = 3*a1*t14*t18
1279 t22 = bb + t13*x
1280 t23 = t22**2
1281 t25 = t17**2
1282 t26 = 1/t25
1283 t28 = 3*a1*t23*t26
1284 t29 = a1*x
1285 t30 = t6**2
1286 t35 = t2*x
1287 t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
1288 t44 = 3*t29*t14*t26*t22
1289 t50 = 2*t29*t23*t22/t25/t17
1290 t51 = a3*a4
1291 t53 = exp(-a4*t2)
1292 t57 = a4**2
1293 t58 = a3*t57
1294 t59 = t35*t53
1295 t64 = t2**2
1296 t70 = log(t17)
1297 t71 = t29*t70
1298 t73 = 0.1e1_dp + t71 + a5*t64
1299 t78 = 2*a1*t22*t18
1300 t80 = t29*t14*t18
1301 t82 = t29*t23*t26
1302 t90 = a3*t53
1303 t93 = t73**2
1304 t94 = 1/t93
1305 t96 = a1*t70
1306 t98 = t29*t18*t22
1307 t101 = t96 + t98 + 4*a5*t35
1308 t106 = a2 - t90
1309 t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1310 t111 = 1/t93/t73
1311 t113 = t101**2
1312 t119 = t78 + t80 - t82 + 12*a5*t2
1313 t123 = 0.1e1_dp + t71 + t106*t2
1314 t124 = t93**2
1315 fs(ip, 4) = &
1316 (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
1317 x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
1318 4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
1319 6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
1320 6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
1321 END DO
1322!$OMP END DO
1323 END IF
1324
1325!$OMP END PARALLEL
1326
1327 END SUBROUTINE efactor_pw91
1328
1329END MODULE xc_ke_gga
1330
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:66
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition xc_functionals_utilities.F:15
xc_functionals_utilities::calc_wave_vector
subroutine, public calc_wave_vector(tag, rho, grho, s)
...
Definition xc_functionals_utilities.F:171
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition xc_functionals_utilities.F:54
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::ke_ol1
integer, parameter, public ke_ol1
Definition xc_input_constants.F:72
xc_input_constants::ke_ol2
integer, parameter, public ke_ol2
Definition xc_input_constants.F:72
xc_input_constants::ke_pw91
integer, parameter, public ke_pw91
Definition xc_input_constants.F:72
xc_input_constants::ke_pbe
integer, parameter, public ke_pbe
Definition xc_input_constants.F:72
xc_input_constants::ke_pw86
integer, parameter, public ke_pw86
Definition xc_input_constants.F:72
xc_input_constants::ke_lc
integer, parameter, public ke_lc
Definition xc_input_constants.F:72
xc_input_constants::ke_llp
integer, parameter, public ke_llp
Definition xc_input_constants.F:72
xc_input_constants::ke_t92
integer, parameter, public ke_t92
Definition xc_input_constants.F:72
xc_ke_gga
Calculate the several different kinetic energy functionals with a GGA form.
Definition xc_ke_gga.F:16
xc_ke_gga::ke_gga_lsd_eval
subroutine, public ke_gga_lsd_eval(functional, rho_set, deriv_set, order)
...
Definition xc_ke_gga.F:313
xc_ke_gga::ke_gga_info
subroutine, public ke_gga_info(functional, lsd, reference, shortform, needs, max_deriv)
...
Definition xc_ke_gga.F:97
xc_ke_gga::ke_gga_lda_eval
subroutine, public ke_gga_lda_eval(functional, rho_set, deriv_set, order)
...
Definition xc_ke_gga.F:178
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
cp_array_utils::cp_3d_r_cp_type
represent a pointer to a contiguous 3d array
Definition cp_array_utils.F:149
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78