(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
16MODULE xc_exchange_gga
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
22 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
23 deriv_norm_drhoa,&
24 deriv_norm_drhob,&
25 deriv_rho,&
26 deriv_rhoa,&
27 deriv_rhob
28 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
29 xc_dset_get_derivative
30 USE xc_derivative_types, ONLY: xc_derivative_get,&
31 xc_derivative_type
32 USE xc_functionals_utilities, ONLY: calc_wave_vector,&
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,&
40 xgga_revpbe
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
62CONTAINS
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
1095END MODULE xc_exchange_gga
1096
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
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_exchange_gga
Calculate several different exchange energy functionals with a GGA form.
Definition xc_exchange_gga.F:16
xc_exchange_gga::xgga_eval
subroutine, public xgga_eval(functional, lsd, rho_set, deriv_set, order)
evaluates different exchange gga
Definition xc_exchange_gga.F:160
xc_exchange_gga::xgga_info
subroutine, public xgga_info(functional, lsd, reference, shortform, needs, max_deriv)
return various information on the xgga functionals
Definition xc_exchange_gga.F:78
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::xgga_b88
integer, parameter, public xgga_b88
Definition xc_input_constants.F:41
xc_input_constants::xgga_opt
integer, parameter, public xgga_opt
Definition xc_input_constants.F:41
xc_input_constants::xgga_revpbe
integer, parameter, public xgga_revpbe
Definition xc_input_constants.F:41
xc_input_constants::xgga_pw86
integer, parameter, public xgga_pw86
Definition xc_input_constants.F:41
xc_input_constants::xgga_pbe
integer, parameter, public xgga_pbe
Definition xc_input_constants.F:41
xc_input_constants::xgga_pw91
integer, parameter, public xgga_pw91
Definition xc_input_constants.F:41
xc_input_constants::xgga_ev93
integer, parameter, public xgga_ev93
Definition xc_input_constants.F:41
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