(git:b279b6b)
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 ! **************************************************************************************************
16 MODULE 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 
63 CONTAINS
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 
1329 END 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