(git:0de0cc2)
xc_exchange_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 several different exchange 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 (27.02.2002)
15 ! **************************************************************************************************
17 
18  USE cp_array_utils, ONLY: cp_3d_r_cp_type
19  USE cp_log_handling, ONLY: cp_to_string
20  USE kinds, ONLY: dp
21  USE mathconstants, ONLY: pi
25  deriv_rho,&
26  deriv_rhoa,&
28  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
31  xc_derivative_type
33  set_util
34  USE xc_input_constants, ONLY: xgga_b88,&
35  xgga_ev93,&
36  xgga_opt,&
37  xgga_pbe,&
38  xgga_pw86,&
39  xgga_pw91,&
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  PUBLIC :: xgga_info, xgga_eval
51 
52  REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
53  f23 = 2.0_dp*f13, &
54  f43 = 4.0_dp*f13, &
55  f53 = 5.0_dp*f13
56 
57  REAL(KIND=dp) :: cx, flda, flsd, sfac, t13
58  REAL(KIND=dp) :: fact, tact
59  REAL(KIND=dp) :: eps_rho
60  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_exchange_gga'
61 
62 CONTAINS
63 
64 ! **************************************************************************************************
65 !> \brief return various information on the xgga functionals
66 !> \param functional integer selecting the xgga functional, it should be one of
67 !> the constants defined in this module: xgga_b88, xgga_pw86,...
68 !> \param lsd a logical that specifies if you are asking about the lsd or lda
69 !> version of the functional
70 !> \param reference string with the reference of the actual functional
71 !> \param shortform string with the shortform of the functional name
72 !> \param needs the components needed by this functional are set to
73 !> true (does not set the unneeded components to false)
74 !> \param max_deriv ...
75 !> \author fawzi
76 ! **************************************************************************************************
77  SUBROUTINE xgga_info(functional, lsd, reference, shortform, needs, max_deriv)
78  INTEGER, INTENT(in) :: functional
79  LOGICAL, INTENT(in) :: lsd
80  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
81  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
82  INTEGER, INTENT(out), OPTIONAL :: max_deriv
83 
84  IF (PRESENT(reference)) THEN
85  SELECT CASE (functional)
86  CASE (xgga_b88)
87  reference = "A. Becke, Phys. Rev. A 38, 3098 (1988)"
88  CASE (xgga_pw86)
89  reference = "J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
90  CASE (xgga_pw91)
91  reference = "J.P. Perdew et al., Phys. Rev. B, 46, 6671 (1992)"
92  CASE (xgga_pbe)
93  reference = "J.P. Perdew, K. Burke, M Ernzerhof, Phys. Rev. Lett, 77, 3865 (1996)"
94  CASE (xgga_revpbe)
95  reference = "Y. Zang et al., PRL, 80, 890 (1998) (Revised PBEX)"
96  CASE (xgga_opt)
97  reference = "Wee-Meng Hoe, A.J. Cohen, N.C. Handy, CPL, 341, 319 (2001)"
98  CASE (xgga_ev93)
99  reference = "E. Engel and S.H. Vosko, Phys. Rev. B, 47, 13164 (1993)"
100  CASE default
101  cpabort("Invalid functional requested ("//cp_to_string(functional)//")")
102  END SELECT
103  IF (.NOT. lsd) THEN
104  IF (len_trim(reference) + 6 < len(reference)) THEN
105  reference(len_trim(reference):len_trim(reference) + 6) = ' {LDA}'
106  END IF
107  END IF
108  END IF
109  IF (PRESENT(shortform)) THEN
110  SELECT CASE (functional)
111  CASE (xgga_b88)
112  shortform = "Becke 1988 Exchange Functional"
113  CASE (xgga_pw86)
114  shortform = "Perdew-Wang 1986 Functional (exchange energy)"
115  CASE (xgga_pw91)
116  shortform = "Perdew-Wang 1991 Functional (exchange energy)"
117  CASE (xgga_pbe)
118  shortform = "PBE exchange energy functional"
119  CASE (xgga_revpbe)
120  shortform = "Revised PBEX by Zang et al."
121  CASE (xgga_opt)
122  shortform = "OPTX exchange energy functional"
123  CASE (xgga_ev93)
124  shortform = "Engel-Vosko exchange energy from virial relation"
125  CASE default
126  cpabort("Invalid functional requested ("//cp_to_string(functional)//")")
127  END SELECT
128  IF (.NOT. lsd) THEN
129  IF (len_trim(shortform) + 6 < len(shortform)) THEN
130  shortform(len_trim(shortform):len_trim(shortform) + 6) = ' {LDA}'
131  END IF
132  END IF
133  END IF
134  IF (PRESENT(needs)) THEN
135  IF (lsd) THEN
136  needs%rho_spin = .true.
137  needs%rho_spin_1_3 = .true.
138  needs%norm_drho_spin = .true.
139  ELSE
140  needs%rho = .true.
141  needs%rho_1_3 = .true.
142  needs%norm_drho = .true.
143  END IF
144  END IF
145  IF (PRESENT(max_deriv)) max_deriv = 3
146 
147  END SUBROUTINE xgga_info
148 
149 ! **************************************************************************************************
150 !> \brief evaluates different exchange gga
151 !> \param functional integer to select the functional that should be evaluated
152 !> \param lsd if the lsd version of the functional should be used
153 !> \param rho_set the density where you want to evaluate the functional
154 !> \param deriv_set place where to store the functional derivatives (they are
155 !> added to the derivatives)
156 !> \param order ...
157 !> \author fawzi
158 ! **************************************************************************************************
159  SUBROUTINE xgga_eval(functional, lsd, rho_set, deriv_set, order)
160  INTEGER, INTENT(in) :: functional
161  LOGICAL, INTENT(in) :: lsd
162  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
163  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
164  INTEGER, INTENT(in) :: order
165 
166  CHARACTER(len=*), PARAMETER :: routinen = 'xgga_eval'
167 
168  INTEGER :: handle, ispin, m, npoints, nspin
169  INTEGER, DIMENSION(2) :: norm_drho_spin_name, rho_spin_name
170  INTEGER, DIMENSION(2, 3) :: bo
171  REAL(kind=dp) :: drho_cutoff, rho_cutoff
172  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
173  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
174  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
175  e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
176  e_rho_rho_rho
177  TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
178  TYPE(xc_derivative_type), POINTER :: deriv
179 
180  CALL timeset(routinen, handle)
181  NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
182  e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
183  DO ispin = 1, 2
184  NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
185  END DO
186 
187  IF (lsd) THEN
188  CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
189  rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
190  rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
191  norm_drhob=norm_drho(2)%array, rho_cutoff=rho_cutoff, &
192  drho_cutoff=drho_cutoff, local_bounds=bo)
193  nspin = 2
194  rho_spin_name = [deriv_rhoa, deriv_rhob]
195  norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob]
196  ELSE
197  CALL xc_rho_set_get(rho_set, rho=rho(1)%array, rho_1_3=rho_1_3(1)%array, &
198  norm_drho=norm_drho(1)%array, local_bounds=bo, rho_cutoff=rho_cutoff, &
199  drho_cutoff=drho_cutoff)
200  nspin = 1
201  rho_spin_name = [deriv_rho, deriv_rho]
202  norm_drho_spin_name = [deriv_norm_drho, deriv_norm_drho]
203  END IF
204  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
205  m = abs(order)
206  CALL xgga_init(rho_cutoff)
207 
208  ALLOCATE (s(npoints))
209  ALLOCATE (fs(npoints, m + 1))
210 
211  DO ispin = 1, nspin
212  IF (lsd) THEN
213  fact = flsd
214  tact = 1.0_dp
215  CALL calc_wave_vector("p", rho(ispin)%array, norm_drho(ispin)%array, s)
216  ELSE
217  fact = flda
218  tact = t13
219  CALL calc_wave_vector("u", rho(ispin)%array, &
220  norm_drho(ispin)%array, s)
221  END IF
222 
223  SELECT CASE (functional)
224  CASE (xgga_b88)
225  CALL efactor_b88(s, fs, m)
226  CASE (xgga_pw86)
227  CALL efactor_pw86(s, fs, m)
228  CASE (xgga_pw91)
229 
230  !! omp: note this is handled slightly differently to the
231  !! other cases to prevent sprawling scope declarations
232  !! in efactor_pw91()
233 
234 !$OMP PARALLEL DEFAULT (NONE) SHARED(s, fs, m)
235  CALL efactor_pw91(s, fs, m)
236 !$OMP END PARALLEL
237 
238  CASE (xgga_pbe)
239  tact = t13
240  CALL efactor_pbex(s, fs, m, 1)
241  IF (lsd) tact = 1.0_dp
242  CASE (xgga_revpbe)
243  tact = t13
244  CALL efactor_pbex(s, fs, m, 2)
245  IF (lsd) tact = 1.0_dp
246  CASE (xgga_opt)
247  CALL efactor_optx(s, fs, m)
248  CASE (xgga_ev93)
249  tact = t13
250  CALL efactor_ev93(s, fs, m)
251  IF (lsd) tact = 1.0_dp
252  CASE DEFAULT
253  cpabort("Unsupported functional")
254  END SELECT
255 
256  IF (order >= 0) THEN
257  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
258  allocate_deriv=.true.)
259  CALL xc_derivative_get(deriv, deriv_data=e_0)
260 
261  CALL x_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, e_0, &
262  npoints)
263  END IF
264  IF (order >= 1 .OR. order == -1) THEN
265  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
266  allocate_deriv=.true.)
267  CALL xc_derivative_get(deriv, deriv_data=e_rho)
268  deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
269  allocate_deriv=.true.)
270  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
271 
272  CALL x_p_1(rho(ispin)%array, &
273  rho_1_3(ispin)%array, s, fs, e_rho, e_ndrho, npoints)
274  END IF
275  IF (order >= 2 .OR. order == -2) THEN
276  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
277  rho_spin_name(ispin)], allocate_deriv=.true.)
278  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
279  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
280  norm_drho_spin_name(ispin)], allocate_deriv=.true.)
281  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
282  deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
283  norm_drho_spin_name(ispin)], allocate_deriv=.true.)
284  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
285 
286  CALL x_p_2(rho(ispin)%array, &
287  rho_1_3(ispin)%array, s, fs, e_rho_rho, e_rho_ndrho, &
288  e_ndrho_ndrho, npoints)
289  END IF
290  IF (order >= 3 .OR. order == -3) THEN
291  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
292  rho_spin_name(ispin), rho_spin_name(ispin)], &
293  allocate_deriv=.true.)
294  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
295  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
296  rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
297  allocate_deriv=.true.)
298  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
299  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
300  norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
301  allocate_deriv=.true.)
302  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
303  deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
304  norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
305  allocate_deriv=.true.)
306  CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
307 
308  CALL x_p_3(rho(ispin)%array, &
309  rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, &
310  e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
311  npoints)
312  END IF
313  IF (order > 3 .OR. order < -3) THEN
314  cpabort("derivatives bigger than 3 not implemented")
315  END IF
316  END DO
317 
318  DEALLOCATE (s)
319  DEALLOCATE (fs)
320 
321  CALL timestop(handle)
322  END SUBROUTINE xgga_eval
323 
324 ! **************************************************************************************************
325 !> \brief ...
326 !> \param cutoff ...
327 ! **************************************************************************************************
328  SUBROUTINE xgga_init(cutoff)
329 
330  REAL(kind=dp), INTENT(IN) :: cutoff
331 
332  eps_rho = cutoff
333  CALL set_util(cutoff)
334 
335  cx = -0.75_dp*(3.0_dp/pi)**f13
336  t13 = 2.0_dp**f13
337  flda = cx
338  flsd = cx*t13
339 
340  sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
341 
342  END SUBROUTINE xgga_init
343 
344 ! **************************************************************************************************
345 !> \brief ...
346 !> \param rho ...
347 !> \param r13 ...
348 !> \param fs ...
349 !> \param e_0 ...
350 !> \param npoints ...
351 ! **************************************************************************************************
352  SUBROUTINE x_p_0(rho, r13, fs, e_0, npoints)
353 
354  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
355  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
356  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
357  INTEGER, INTENT(in) :: npoints
358 
359  INTEGER :: ip
360 
361 !$OMP PARALLEL DO DEFAULT (NONE) &
362 !$OMP SHARED (npoints, rho, eps_rho, fact, r13, fs, e_0) &
363 !$OMP PRIVATE(ip)
364 
365  DO ip = 1, npoints
366 
367  IF (rho(ip) > eps_rho) THEN
368  e_0(ip) = e_0(ip) + fact*r13(ip)*rho(ip)*fs(ip, 1)
369  END IF
370 
371  END DO
372 
373 !$OMP END PARALLEL DO
374 
375  END SUBROUTINE x_p_0
376 
377 ! **************************************************************************************************
378 !> \brief ...
379 !> \param rho ...
380 !> \param r13 ...
381 !> \param s ...
382 !> \param fs ...
383 !> \param e_rho ...
384 !> \param e_ndrho ...
385 !> \param npoints ...
386 ! **************************************************************************************************
387  SUBROUTINE x_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
388 
389  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
390  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
391  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
392  INTEGER, INTENT(in) :: npoints
393 
394  INTEGER :: ip
395  REAL(kind=dp) :: a0, a1, sx, sy
396 
397 !$OMP PARALLEL DO DEFAULT(NONE) &
398 !$OMP SHARED(npoints, rho, eps_rho, fact, r13, tact, fs) &
399 !$OMP SHARED(e_rho, e_ndrho, sfac, s) &
400 !$OMP PRIVATE(ip,a0,a1,sx,sy)
401 
402  DO ip = 1, npoints
403 
404  IF (rho(ip) > eps_rho) THEN
405 
406  a0 = fact*r13(ip)*rho(ip)
407  a1 = f43*fact*r13(ip)
408  sx = -f43*s(ip)/rho(ip)
409  sy = sfac*tact/(r13(ip)*rho(ip))
410  e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
411  e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
412 
413  END IF
414 
415  END DO
416 
417 !$OMP END PARALLEL DO
418 
419  END SUBROUTINE x_p_1
420 
421 ! **************************************************************************************************
422 !> \brief ...
423 !> \param rho ...
424 !> \param r13 ...
425 !> \param s ...
426 !> \param fs ...
427 !> \param e_rho_rho ...
428 !> \param e_rho_ndrho ...
429 !> \param e_ndrho_ndrho ...
430 !> \param npoints ...
431 ! **************************************************************************************************
432  SUBROUTINE x_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, &
433  e_ndrho_ndrho, npoints)
434 
435  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
436  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
437  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
438  INTEGER, INTENT(in) :: npoints
439 
440  INTEGER :: ip
441  REAL(kind=dp) :: a0, a1, a2, sx, sxx, sxy, sy
442 
443 !$OMP PARALLEL DO DEFAULT(NONE) &
444 !$OMP SHARED(npoints, rho, eps_rho, r13, fact, e_rho_rho) &
445 !$OMP SHARED(e_rho_ndrho, e_ndrho_ndrho, fs, sfac, tact, s) &
446 !$OMP PRIVATE(ip,a0,a1,a2,sx,sy,sxx,sxy)
447 
448  DO ip = 1, npoints
449 
450  IF (rho(ip) > eps_rho) THEN
451 
452  a0 = fact*r13(ip)*rho(ip)
453  a1 = f43*fact*r13(ip)
454  a2 = f13*f43*fact/(r13(ip)*r13(ip))
455  sx = -f43*s(ip)/rho(ip)
456  sy = sfac*tact/(r13(ip)*rho(ip))
457  sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
458  sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
459  e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
460  a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
461  e_rho_ndrho(ip) = e_rho_ndrho(ip) &
462  + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + a0*fs(ip, 2)*sxy
463  e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
464 
465  END IF
466 
467  END DO
468 
469 !$OMP END PARALLEL DO
470 
471  END SUBROUTINE x_p_2
472 
473 ! **************************************************************************************************
474 !> \brief ...
475 !> \param rho ...
476 !> \param r13 ...
477 !> \param s ...
478 !> \param fs ...
479 !> \param e_rho_rho_rho ...
480 !> \param e_rho_rho_ndrho ...
481 !> \param e_rho_ndrho_ndrho ...
482 !> \param e_ndrho_ndrho_ndrho ...
483 !> \param npoints ...
484 ! **************************************************************************************************
485  SUBROUTINE x_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
486  e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
487 
488  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
489  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
490  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
491  e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
492  INTEGER, INTENT(in) :: npoints
493 
494  INTEGER :: ip
495  REAL(kind=dp) :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
496  sxy, sy
497 
498 !$OMP PARALLEL DO DEFAULT (NONE) &
499 !$OMP SHARED(npoints, rho, eps_rho, r13, fact, fs) &
500 !$OMP SHARED(e_rho_rho_rho, e_rho_rho_ndrho) &
501 !$OMP SHARED(e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho) &
502 !$OMP SHARED(sfac, tact, s) &
503 !$OMP PRIVATE(ip,a0,a1,a2,a3,sx,sy,sxx,sxy,sxxx,sxxy)
504 
505  DO ip = 1, npoints
506 
507  IF (rho(ip) > eps_rho) THEN
508 
509  a0 = fact*r13(ip)*rho(ip)
510  a1 = f43*fact*r13(ip)
511  a2 = f13*f43*fact/(r13(ip)*r13(ip))
512  a3 = -f23*f13*f43*fact/(r13(ip)*r13(ip)*rho(ip))
513  sx = -f43*s(ip)/rho(ip)
514  sy = sfac*tact/(r13(ip)*rho(ip))
515  sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
516  sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
517  sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
518  sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
519  e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
520  + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
521  3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
522  a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
523  a0*fs(ip, 2)*sxxx
524  e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
525  + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
526  2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
527  2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
528  a0*fs(ip, 2)*sxxy
529  e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
530  + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
531  2.0_dp*a0*fs(ip, 3)*sxy*sy
532  e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) &
533  + a0*fs(ip, 4)*sy*sy*sy
534 
535  END IF
536 
537  END DO
538 
539 !$OMP END PARALLEL DO
540 
541  END SUBROUTINE x_p_3
542 
543 ! Enhancement Factors
544 ! **************************************************************************************************
545 !> \brief ...
546 !> \param s ...
547 !> \param fs ...
548 !> \param m ...
549 ! **************************************************************************************************
550  SUBROUTINE efactor_b88(s, fs, m)
551  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
552  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
553  INTEGER, INTENT(IN) :: m
554 
555  INTEGER :: ip
556  REAL(kind=dp) :: as, asp, beta, bs, f0, p, q, sas, sbs, sbs3, t1, t10, t13, t15, t16, t19, &
557  t2, t22, t24, t25, t32, t34, t36, t39, t4, t40, t41, t44, t48, t49, t5, t6, t65, t8, t87, &
558  t9, x, ys
559 
560  beta = 0.0042_dp
561  f0 = 1.0_dp/sfac
562  p = -beta/flsd
563  q = 6.0_dp*beta
564 
565 !$OMP PARALLEL DO DEFAULT(NONE) &
566 !$OMP SHARED(s,fs,m,beta,f0,p,q) &
567 !$OMP PRIVATE(ip,x,bs,sbs,as,sas,ys,asp,sbs3) &
568 !$OMP PRIVATE(t1,t2,t4,t5,t6,t8,t9,t10,t13,t15,t16,t19,t22) &
569 !$OMP PRIVATE(t24,t25,t32,t34) &
570 !$OMP PRIVATE(t36,t39,t40,t41,t44,t48,t49,t65,t87)
571 
572  DO ip = 1, SIZE(s)
573  x = s(ip)*f0
574  bs = beta*x
575  sbs = sqrt(x*x + 1.0_dp)
576  as = log(x + sbs)
577  sas = x*as
578  ys = 1.0_dp/(1.0_dp + q*sas)
579  SELECT CASE (m)
580  CASE (0)
581  fs(ip, 1) = 1.0_dp + p*x*x*ys
582  CASE (1)
583  asp = as + x/sbs
584  fs(ip, 1) = 1.0_dp + p*x*x*ys
585  fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
586  CASE (2)
587  asp = as + x/sbs
588  sbs3 = 1.0_dp/(sbs*sbs*sbs)
589  fs(ip, 1) = 1.0_dp + p*x*x*ys
590  fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
591  fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
592  - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
593  + 3.0_dp*q - sbs) - sbs))
594  CASE (3)
595  asp = as + x/sbs
596  sbs3 = 1.0_dp/(sbs*sbs*sbs)
597  fs(ip, 1) = 1.0_dp + p*x*x*ys
598  fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
599  fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
600  - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
601  + 3.0_dp*q - sbs) - sbs))
602  t1 = q*x
603  t2 = x**2
604  t4 = sqrt(1 + t2)
605  t5 = x + t4
606  t6 = log(t5)
607  t8 = 1 + t1*t6
608  t9 = t8**2
609  t10 = 1/t9
610  t13 = 1/t4
611  t15 = 1 + t13*x
612  t16 = 1/t5
613  t19 = q*t6 + t1*t15*t16
614  t22 = p*x
615  t24 = 1/t9/t8
616  t25 = t19**2
617  t32 = t4**2
618  t34 = 1/t32/t4
619  t36 = -t34*t2 + t13
620  t39 = t15**2
621  t40 = t5**2
622  t41 = 1/t40
623  t44 = 2*q*t15*t16 + t1*t36*t16 - t1*t39*t41
624  t48 = p*t2
625  t49 = t9**2
626  t65 = t32**2
627  t87 = -6*p*t10*t19 + 12*t22*t24*t25 - 6*t22*t10*t44 - 6*t48/t49*t25*t19 + &
628  6*t48*t24*t19*t44 - t48*t10*(3*q*t36*t16 - 3*q*t39*t41 + 3*t1*(1/t65/t4* &
629  t2*x - t34*x)*t16 - 3*t1*t36*t41*t15 + 2*t1*t39*t15/t40/t5)
630 
631  fs(ip, 4) = t87
632  fs(ip, 4) = f0*f0*f0*fs(ip, 4)
633 
634  CASE DEFAULT
635  cpabort("Illegal order")
636  END SELECT
637  END DO
638 
639 !$OMP END PARALLEL DO
640 
641  END SUBROUTINE efactor_b88
642 ! **************************************************************************************************
643 !> \brief ...
644 !> \param s ...
645 !> \param fs ...
646 !> \param m ...
647 ! **************************************************************************************************
648  SUBROUTINE efactor_pw86(s, fs, m)
649  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
650  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
651  INTEGER, INTENT(IN) :: m
652 
653  INTEGER :: ip
654  REAL(kind=dp) :: f15, p0, p1, p15, s2, s4, s6, t1, t10, &
655  t12, t13, t14, t19, t2, t25, t3, t8, t9
656 
657  t1 = 1.296_dp
658  t2 = 14.0_dp
659  t3 = 0.2_dp
660  f15 = 1.0_dp/15.0_dp
661 
662 !$OMP PARALLEL DO DEFAULT(NONE) &
663 !$OMP SHARED(s, fs, m, t1, t2, t3, f15) &
664 !$OMP PRIVATE(ip,s2,s4,s6,p0,p1,p15)&
665 !$OMP PRIVATE(t8, t9, t10, t12, t13, t14, t19, t25)
666 
667  DO ip = 1, SIZE(s)
668  s2 = s(ip)*s(ip)
669  s4 = s2*s2
670  s6 = s2*s4
671  SELECT CASE (m)
672  CASE (0)
673  p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
674  fs(ip, 1) = p0**f15
675  CASE (1)
676  p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
677  p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
678  p15 = p0**f15
679  fs(ip, 1) = p15
680  fs(ip, 2) = f15*p1*p15/p0
681  CASE (2)
682  p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
683  p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
684  p15 = p0**f15
685  fs(ip, 1) = p15
686  fs(ip, 2) = f15*p1*p15/p0
687  t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
688  t25 = p1*p1
689  fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
690  1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
691  CASE (3)
692  p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
693  p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
694  p15 = p0**f15
695  fs(ip, 1) = p15
696  fs(ip, 2) = f15*p1*p15/p0
697  t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
698  t25 = p1*p1
699  fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
700  1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
701  t8 = p0**2; t9 = p0**f15; t14 = p0/t9; t19 = s2*s(ip)
702  fs(ip, 4) = 406.0_dp/3375.0_dp/t14/t8*p1*p1*p1 - 14.0_dp/ &
703  75.0_dp/t14/p0*p1*(2*t1 + 12*t2*s2 + 30*t3*s4) + &
704  1/t14*(24*t2*s(ip) + 120*t3*t19)*f15
705  CASE DEFAULT
706  cpabort("Illegal order")
707  END SELECT
708  END DO
709 
710 !$OMP END PARALLEL DO
711 
712  END SUBROUTINE efactor_pw86
713 ! **************************************************************************************************
714 !> \brief ...
715 !> \param s ...
716 !> \param fs ...
717 !> \param m ...
718 ! **************************************************************************************************
719  SUBROUTINE efactor_ev93(s, fs, m)
720  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
721  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
722  INTEGER, INTENT(IN) :: m
723 
724  INTEGER :: ip
725  REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, d0, d1, d2, d3, &
726  f0, f1, f2, n0, n1, n2, n3, s2, s4, &
727  s6, scale_s, ss
728 
729  a1 = 1.647127_dp
730  a2 = 0.980118_dp
731  a3 = 0.017399_dp
732  b1 = 1.523671_dp
733  b2 = 0.367229_dp
734  b3 = 0.011282_dp
735  scale_s = 1._dp/tact
736 
737 !$OMP PARALLEL DO DEFAULT(NONE) &
738 !$OMP SHARED(s, fs, m, a1, a2, a3, b1, b2, b3, scale_s) &
739 !$OMP PRIVATE(ip,ss,s2,s4,s6) &
740 !$OMP PRIVATE(n0,n1,n2,n3,d0,d1,d2,d3,f0,f1,f2)
741 
742  DO ip = 1, SIZE(s)
743 ! ss = s(ip)
744 !
745  ss = scale_s*s(ip)
746  s2 = ss*ss
747  s4 = s2*s2
748  s6 = s2*s4
749  SELECT CASE (m)
750  CASE (0)
751  n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
752  d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
753  fs(ip, 1) = n0/d0
754  CASE (1)
755  n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
756  d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
757  n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
758  d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
759  f0 = n0/d0
760  fs(ip, 1) = f0
761  fs(ip, 2) = (n1 - f0*d1)/d0*scale_s
762  CASE (2)
763  n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
764  d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
765  n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
766  d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
767  n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
768  d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
769  f0 = n0/d0
770  f1 = (n1 - f0*d1)/d0
771  fs(ip, 1) = f0
772  fs(ip, 2) = f1*scale_s
773  fs(ip, 3) = (n2 - f0*d2 - 2._dp*f1*d1)/d0*scale_s*scale_s
774  CASE (3)
775  n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
776  d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
777  n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
778  d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
779  n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
780  d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
781  n3 = ss*(24._dp*a2 + 120._dp*a3*s2)
782  d3 = ss*(24._dp*b2 + 120._dp*b3*s2)
783  f0 = n0/d0
784  f1 = (n1 - f0*d1)/d0
785  f2 = (n2 - f0*d2 - 2._dp*f1*d1)/d0
786  fs(ip, 1) = f0
787  fs(ip, 2) = f1*scale_s
788  fs(ip, 3) = f2*scale_s*scale_s
789  fs(ip, 4) = (n3 - f0*d3 - 3._dp*f2*d1 - 3._dp*f1*d2)/d0*scale_s*scale_s*scale_s
790  CASE DEFAULT
791  cpabort("Illegal order")
792  END SELECT
793  END DO
794 
795 !$OMP END PARALLEL DO
796 
797  END SUBROUTINE efactor_ev93
798 ! **************************************************************************************************
799 !> \brief ...
800 !> \param s ...
801 !> \param fs ...
802 !> \param m ...
803 ! **************************************************************************************************
804  SUBROUTINE efactor_optx(s, fs, m)
805  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
806  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
807  INTEGER, INTENT(IN) :: m
808 
809  REAL(kind=dp), PARAMETER :: a1 = 1.05151_dp, a2 = 1.43169_dp, &
810  gamma_bo = 0.006_dp
811 
812  INTEGER :: ip
813  REAL(kind=dp) :: a, b, f0, x, y
814 
815  f0 = 1.0_dp/sfac
816  b = -a2/flsd
817 
818 !$OMP PARALLEL DO DEFAULT(NONE) &
819 !$OMP SHARED(s, fs, m, f0, b) &
820 !$OMP PRIVATE(ip,x,A,y)
821 
822  DO ip = 1, SIZE(s)
823  x = s(ip)*f0
824  a = gamma_bo*x*x
825  y = 1.0_dp/(1.0_dp + a)
826  SELECT CASE (m)
827  CASE (0)
828  fs(ip, 1) = a1 + b*a*a*y*y
829  CASE (1)
830  fs(ip, 1) = a1 + b*a*a*y*y
831  fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
832  CASE (2)
833  fs(ip, 1) = a1 + b*a*a*y*y
834  fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
835  fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
836  CASE (3)
837  fs(ip, 1) = a1 + b*a*a*y*y
838  fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
839  fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
840  fs(ip, 4) = 24.0_dp*b*f0*f0*f0*gamma_bo*gamma_bo*x* &
841  (1.0_dp - 5.0_dp*a + 2.0_dp*a*a)*y*y*y*y*y
842  CASE DEFAULT
843  cpabort("Illegal order")
844  END SELECT
845  END DO
846 
847 !$OMP END PARALLEL DO
848 
849  END SUBROUTINE efactor_optx
850 
851 ! **************************************************************************************************
852 !> \brief ...
853 !> \param s ...
854 !> \param fs ...
855 !> \param m ...
856 ! **************************************************************************************************
857  SUBROUTINE efactor_pw91(s, fs, m)
858  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
859  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
860  INTEGER, INTENT(IN) :: m
861 
862  INTEGER :: ip
863  REAL(dp) :: t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
864  t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
865  t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
866  t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
867  t96, t98
868  REAL(kind=dp) :: a1, a2, a3, a4, a5, b, o, x
869 
870  o = 1.0_dp
871  a1 = 0.19645_dp
872  a2 = 0.2743_dp
873  a3 = 0.1508_dp
874  a4 = 100.0_dp
875  b = 0.8145161_dp
876  a5 = 0.004_dp
877  IF (m >= 0) THEN
878 
879 !$OMP DO
880 
881  DO ip = 1, SIZE(s)
882  x = s(ip)
883  t3 = b**2
884  t4 = x**2
885  t7 = sqrt(o + t3*t4)
886  t9 = log(b*x + t7)
887  t10 = a1*x*t9
888  t12 = exp(-a4*t4)
889  t17 = t4**2
890  fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
891  END DO
892 
893 !$OMP END DO
894 
895  END IF
896  IF (m >= 1) THEN
897 
898 !$OMP DO
899 
900  DO ip = 1, SIZE(s)
901  x = s(ip)
902  t2 = b**2
903  t3 = x**2
904  t6 = sqrt(o + t2*t3)
905  t7 = b*x + t6
906  t8 = log(t7)
907  t9 = a1*t8
908  t10 = a1*x
909  t17 = t10*(b + 1/t6*t2*x)/t7
910  t19 = t3*x
911  t21 = exp(-a4*t3)
912  t26 = a2 - a3*t21
913  t30 = t10*t8
914  t31 = t3**2
915  t33 = o + t30 + a5*t31
916  t38 = t33**2
917  fs(ip, 2) = &
918  (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
919  t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
920  END DO
921 
922 !$OMP END DO
923 
924  END IF
925  IF (m >= 2) THEN
926 
927 !$OMP DO
928 
929  DO ip = 1, SIZE(s)
930  x = s(ip)
931  t1 = b**2
932  t2 = x**2
933  t5 = sqrt(o + t1*t2)
934  t7 = o/t5*t1
935  t9 = b + t7*x
936  t12 = b*x + t5
937  t13 = o/t12
938  t15 = 2._dp*a1*t9*t13
939  t16 = a1*x
940  t17 = t5**2
941  t20 = t1**2
942  t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
943  t26 = t9**2
944  t27 = t12**2
945  t30 = t16*t26/t27
946  t31 = a3*a4
947  t33 = exp(-a4*t2)
948  t37 = a4**2
949  t39 = t2**2
950  t44 = a3*t33
951  t47 = log(t12)
952  t48 = t16*t47
953  t50 = o + t48 + a5*t39
954  t53 = a1*t47
955  t55 = t16*t9*t13
956  t56 = t2*x
957  t60 = a2 - t44
958  t64 = t50**2
959  t65 = o/t64
960  t69 = t53 + t55 + 4._dp*a5*t56
961  t73 = o + t48 + t60*t2
962  t77 = t69**2
963  fs(ip, 3) = &
964  (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
965  2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
966  (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
967  t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
968  END DO
969 
970 !$OMP END DO
971 
972  END IF
973  IF (m >= 3) THEN
974 
975 !$OMP DO
976 
977  DO ip = 1, SIZE(s)
978  x = s(ip)
979  t1 = b**2
980  t2 = x**2
981  t5 = sqrt(0.1e1_dp + t1*t2)
982  t6 = t5**2
983  t9 = t1**2
984  t10 = 1/t6/t5*t9
985  t13 = 1/t5*t1
986  t14 = -t10*t2 + t13
987  t17 = b*x + t5
988  t18 = 1/t17
989  t20 = 3*a1*t14*t18
990  t22 = b + t13*x
991  t23 = t22**2
992  t25 = t17**2
993  t26 = 1/t25
994  t28 = 3*a1*t23*t26
995  t29 = a1*x
996  t30 = t6**2
997  t35 = t2*x
998  t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
999  t44 = 3*t29*t14*t26*t22
1000  t50 = 2*t29*t23*t22/t25/t17
1001  t51 = a3*a4
1002  t53 = exp(-a4*t2)
1003  t57 = a4**2
1004  t58 = a3*t57
1005  t59 = t35*t53
1006  t64 = t2**2
1007  t70 = log(t17)
1008  t71 = t29*t70
1009  t73 = 0.1e1_dp + t71 + a5*t64
1010  t78 = 2*a1*t22*t18
1011  t80 = t29*t14*t18
1012  t82 = t29*t23*t26
1013  t90 = a3*t53
1014  t93 = t73**2
1015  t94 = 1/t93
1016  t96 = a1*t70
1017  t98 = t29*t18*t22
1018  t101 = t96 + t98 + 4*a5*t35
1019  t106 = a2 - t90
1020  t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1021  t111 = 1/t93/t73
1022  t113 = t101**2
1023  t119 = t78 + t80 - t82 + 12*a5*t2
1024  t123 = 0.1e1_dp + t71 + t106*t2
1025  t124 = t93**2
1026  fs(ip, 4) = &
1027  (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
1028  x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
1029  4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
1030  6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
1031  6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
1032  END DO
1033 
1034 !$OMP END DO
1035 
1036  END IF
1037 
1038  END SUBROUTINE efactor_pw91
1039 
1040 ! **************************************************************************************************
1041 !> \brief ...
1042 !> \param s ...
1043 !> \param fs ...
1044 !> \param m ...
1045 !> \param pset ...
1046 ! **************************************************************************************************
1047  SUBROUTINE efactor_pbex(s, fs, m, pset)
1048  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
1049  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1050  INTEGER, INTENT(IN) :: m, pset
1051 
1052  REAL(kind=dp), PARAMETER :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
1053  mu = 0.2195149727645171_dp
1054 
1055  INTEGER :: ip
1056  REAL(kind=dp) :: f0, mk, x, x2, y
1057 
1058  IF (pset == 1) mk = mu/kappa1
1059  IF (pset == 2) mk = mu/kappa2
1060 
1061  f0 = 1.0_dp/tact
1062 
1063 !$OMP PARALLEL DO DEFAULT(NONE) &
1064 !$OMP SHARED (s, fs, m, mk, f0) &
1065 !$OMP PRIVATE(ip,x,x2,y)
1066 
1067  DO ip = 1, SIZE(s)
1068  x = s(ip)*f0
1069  x2 = x*x
1070  y = 1.0_dp/(1.0_dp + mk*x2)
1071  SELECT CASE (m)
1072  CASE (0)
1073  fs(ip, 1) = 1.0_dp + mu*x2*y
1074  CASE (1)
1075  fs(ip, 1) = 1.0_dp + mu*x2*y
1076  fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1077  CASE (2)
1078  fs(ip, 1) = 1.0_dp + mu*x2*y
1079  fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1080  fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1081  CASE (3)
1082  fs(ip, 1) = 1.0_dp + mu*x2*y
1083  fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1084  fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1085  fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
1086  CASE DEFAULT
1087  cpabort("Illegal order")
1088  END SELECT
1089  END DO
1090 
1091 !$OMP END PARALLEL DO
1092 
1093  END SUBROUTINE efactor_pbex
1094 
1095 END MODULE xc_exchange_gga
1096 
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
various routines to log and control the output. The idea is that decisions about where to log should ...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_norm_drho
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Calculate several different exchange energy functionals with a GGA form.
subroutine, public xgga_eval(functional, lsd, rho_set, deriv_set, order)
evaluates different exchange gga
subroutine, public xgga_info(functional, lsd, reference, shortform, needs, max_deriv)
return various information on the xgga functionals
Utility routines for the functional calculations.
subroutine, public calc_wave_vector(tag, rho, grho, s)
...
subroutine, public set_util(cutoff)
...
input constants for xc
integer, parameter, public xgga_b88
integer, parameter, public xgga_opt
integer, parameter, public xgga_revpbe
integer, parameter, public xgga_pw86
integer, parameter, public xgga_pbe
integer, parameter, public xgga_pw91
integer, parameter, public xgga_ev93
contains the structure
contains the structure
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