(git:08825b9)
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-2026 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. Zhang 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 Zhang 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 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
257 allocate_deriv=.true.)
258 CALL xc_derivative_get(deriv, deriv_data=e_0)
259
260 CALL x_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, e_0, &
261 npoints)
262 IF (order >= 1) THEN
263 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
264 allocate_deriv=.true.)
265 CALL xc_derivative_get(deriv, deriv_data=e_rho)
266 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
267 allocate_deriv=.true.)
268 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
269
270 CALL x_p_1(rho(ispin)%array, &
271 rho_1_3(ispin)%array, s, fs, e_rho, e_ndrho, npoints)
272 END IF
273 IF (order >= 2) THEN
274 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
275 rho_spin_name(ispin)], allocate_deriv=.true.)
276 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
277 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
278 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
279 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
280 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
281 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
282 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
283
284 CALL x_p_2(rho(ispin)%array, &
285 rho_1_3(ispin)%array, s, fs, e_rho_rho, e_rho_ndrho, &
286 e_ndrho_ndrho, npoints)
287 END IF
288 IF (order >= 3) THEN
289 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
290 rho_spin_name(ispin), rho_spin_name(ispin)], &
291 allocate_deriv=.true.)
292 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
293 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
294 rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
295 allocate_deriv=.true.)
296 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
297 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
298 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
299 allocate_deriv=.true.)
300 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
301 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
302 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
303 allocate_deriv=.true.)
304 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
305
306 CALL x_p_3(rho(ispin)%array, &
307 rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, &
308 e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
309 npoints)
310 END IF
311 IF (order > 3) THEN
312 cpabort("derivatives bigger than 3 not implemented")
313 END IF
314 END DO
315
316 DEALLOCATE (s)
317 DEALLOCATE (fs)
318
319 CALL timestop(handle)
320 END SUBROUTINE xgga_eval
321
322! **************************************************************************************************
323!> \brief ...
324!> \param cutoff ...
325! **************************************************************************************************
326 SUBROUTINE xgga_init(cutoff)
327
328 REAL(kind=dp), INTENT(IN) :: cutoff
329
330 eps_rho = cutoff
331 CALL set_util(cutoff)
332
333 cx = -0.75_dp*(3.0_dp/pi)**f13
334 t13 = 2.0_dp**f13
335 flda = cx
336 flsd = cx*t13
337
338 sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
339
340 END SUBROUTINE xgga_init
341
342! **************************************************************************************************
343!> \brief ...
344!> \param rho ...
345!> \param r13 ...
346!> \param fs ...
347!> \param e_0 ...
348!> \param npoints ...
349! **************************************************************************************************
350 SUBROUTINE x_p_0(rho, r13, fs, e_0, npoints)
351
352 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
353 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
354 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
355 INTEGER, INTENT(in) :: npoints
356
357 INTEGER :: ip
358
359!$OMP PARALLEL DO DEFAULT (NONE) &
360!$OMP SHARED (npoints, rho, eps_rho, fact, r13, fs, e_0) &
361!$OMP PRIVATE(ip)
362
363 DO ip = 1, npoints
364
365 IF (rho(ip) > eps_rho) THEN
366 e_0(ip) = e_0(ip) + fact*r13(ip)*rho(ip)*fs(ip, 1)
367 END IF
368
369 END DO
370
371!$OMP END PARALLEL DO
372
373 END SUBROUTINE x_p_0
374
375! **************************************************************************************************
376!> \brief ...
377!> \param rho ...
378!> \param r13 ...
379!> \param s ...
380!> \param fs ...
381!> \param e_rho ...
382!> \param e_ndrho ...
383!> \param npoints ...
384! **************************************************************************************************
385 SUBROUTINE x_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
386
387 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
388 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
389 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
390 INTEGER, INTENT(in) :: npoints
391
392 INTEGER :: ip
393 REAL(kind=dp) :: a0, a1, sx, sy
394
395!$OMP PARALLEL DO DEFAULT(NONE) &
396!$OMP SHARED(npoints, rho, eps_rho, fact, r13, tact, fs) &
397!$OMP SHARED(e_rho, e_ndrho, sfac, s) &
398!$OMP PRIVATE(ip,a0,a1,sx,sy)
399
400 DO ip = 1, npoints
401
402 IF (rho(ip) > eps_rho) THEN
403
404 a0 = fact*r13(ip)*rho(ip)
405 a1 = f43*fact*r13(ip)
406 sx = -f43*s(ip)/rho(ip)
407 sy = sfac*tact/(r13(ip)*rho(ip))
408 e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
409 e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
410
411 END IF
412
413 END DO
414
415!$OMP END PARALLEL DO
416
417 END SUBROUTINE x_p_1
418
419! **************************************************************************************************
420!> \brief ...
421!> \param rho ...
422!> \param r13 ...
423!> \param s ...
424!> \param fs ...
425!> \param e_rho_rho ...
426!> \param e_rho_ndrho ...
427!> \param e_ndrho_ndrho ...
428!> \param npoints ...
429! **************************************************************************************************
430 SUBROUTINE x_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, &
431 e_ndrho_ndrho, npoints)
432
433 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
434 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
435 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
436 INTEGER, INTENT(in) :: npoints
437
438 INTEGER :: ip
439 REAL(kind=dp) :: a0, a1, a2, sx, sxx, sxy, sy
440
441!$OMP PARALLEL DO DEFAULT(NONE) &
442!$OMP SHARED(npoints, rho, eps_rho, r13, fact, e_rho_rho) &
443!$OMP SHARED(e_rho_ndrho, e_ndrho_ndrho, fs, sfac, tact, s) &
444!$OMP PRIVATE(ip,a0,a1,a2,sx,sy,sxx,sxy)
445
446 DO ip = 1, npoints
447
448 IF (rho(ip) > eps_rho) THEN
449
450 a0 = fact*r13(ip)*rho(ip)
451 a1 = f43*fact*r13(ip)
452 a2 = f13*f43*fact/(r13(ip)*r13(ip))
453 sx = -f43*s(ip)/rho(ip)
454 sy = sfac*tact/(r13(ip)*rho(ip))
455 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
456 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
457 e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
458 a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
459 e_rho_ndrho(ip) = e_rho_ndrho(ip) &
460 + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + a0*fs(ip, 2)*sxy
461 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
462
463 END IF
464
465 END DO
466
467!$OMP END PARALLEL DO
468
469 END SUBROUTINE x_p_2
470
471! **************************************************************************************************
472!> \brief ...
473!> \param rho ...
474!> \param r13 ...
475!> \param s ...
476!> \param fs ...
477!> \param e_rho_rho_rho ...
478!> \param e_rho_rho_ndrho ...
479!> \param e_rho_ndrho_ndrho ...
480!> \param e_ndrho_ndrho_ndrho ...
481!> \param npoints ...
482! **************************************************************************************************
483 SUBROUTINE x_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
484 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
485
486 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
487 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: fs
488 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
489 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
490 INTEGER, INTENT(in) :: npoints
491
492 INTEGER :: ip
493 REAL(kind=dp) :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
494 sxy, sy
495
496!$OMP PARALLEL DO DEFAULT (NONE) &
497!$OMP SHARED(npoints, rho, eps_rho, r13, fact, fs) &
498!$OMP SHARED(e_rho_rho_rho, e_rho_rho_ndrho) &
499!$OMP SHARED(e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho) &
500!$OMP SHARED(sfac, tact, s) &
501!$OMP PRIVATE(ip,a0,a1,a2,a3,sx,sy,sxx,sxy,sxxx,sxxy)
502
503 DO ip = 1, npoints
504
505 IF (rho(ip) > eps_rho) THEN
506
507 a0 = fact*r13(ip)*rho(ip)
508 a1 = f43*fact*r13(ip)
509 a2 = f13*f43*fact/(r13(ip)*r13(ip))
510 a3 = -f23*f13*f43*fact/(r13(ip)*r13(ip)*rho(ip))
511 sx = -f43*s(ip)/rho(ip)
512 sy = sfac*tact/(r13(ip)*rho(ip))
513 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
514 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
515 sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
516 sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
517 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
518 + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
519 3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
520 a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
521 a0*fs(ip, 2)*sxxx
522 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
523 + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
524 2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
525 2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
526 a0*fs(ip, 2)*sxxy
527 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
528 + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
529 2.0_dp*a0*fs(ip, 3)*sxy*sy
530 e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) &
531 + a0*fs(ip, 4)*sy*sy*sy
532
533 END IF
534
535 END DO
536
537!$OMP END PARALLEL DO
538
539 END SUBROUTINE x_p_3
540
541! Enhancement Factors
542! **************************************************************************************************
543!> \brief ...
544!> \param s ...
545!> \param fs ...
546!> \param m ...
547! **************************************************************************************************
548 SUBROUTINE efactor_b88(s, fs, m)
549 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
550 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
551 INTEGER, INTENT(IN) :: m
552
553 INTEGER :: ip
554 REAL(kind=dp) :: as, asp, beta, bs, f0, p, q, sas, sbs, sbs3, t1, t10, t13, t15, t16, t19, &
555 t2, t22, t24, t25, t32, t34, t36, t39, t4, t40, t41, t44, t48, t49, t5, t6, t65, t8, t87, &
556 t9, x, ys
557
558 beta = 0.0042_dp
559 f0 = 1.0_dp/sfac
560 p = -beta/flsd
561 q = 6.0_dp*beta
562
563!$OMP PARALLEL DO DEFAULT(NONE) &
564!$OMP SHARED(s,fs,m,beta,f0,p,q) &
565!$OMP PRIVATE(ip,x,bs,sbs,as,sas,ys,asp,sbs3) &
566!$OMP PRIVATE(t1,t2,t4,t5,t6,t8,t9,t10,t13,t15,t16,t19,t22) &
567!$OMP PRIVATE(t24,t25,t32,t34) &
568!$OMP PRIVATE(t36,t39,t40,t41,t44,t48,t49,t65,t87)
569
570 DO ip = 1, SIZE(s)
571 x = s(ip)*f0
572 bs = beta*x
573 sbs = sqrt(x*x + 1.0_dp)
574 as = log(x + sbs)
575 sas = x*as
576 ys = 1.0_dp/(1.0_dp + q*sas)
577 SELECT CASE (m)
578 CASE (0)
579 fs(ip, 1) = 1.0_dp + p*x*x*ys
580 CASE (1)
581 asp = as + x/sbs
582 fs(ip, 1) = 1.0_dp + p*x*x*ys
583 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
584 CASE (2)
585 asp = as + x/sbs
586 sbs3 = 1.0_dp/(sbs*sbs*sbs)
587 fs(ip, 1) = 1.0_dp + p*x*x*ys
588 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
589 fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
590 - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
591 + 3.0_dp*q - sbs) - sbs))
592 CASE (3)
593 asp = as + x/sbs
594 sbs3 = 1.0_dp/(sbs*sbs*sbs)
595 fs(ip, 1) = 1.0_dp + p*x*x*ys
596 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
597 fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
598 - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
599 + 3.0_dp*q - sbs) - sbs))
600 t1 = q*x
601 t2 = x**2
602 t4 = sqrt(1 + t2)
603 t5 = x + t4
604 t6 = log(t5)
605 t8 = 1 + t1*t6
606 t9 = t8**2
607 t10 = 1/t9
608 t13 = 1/t4
609 t15 = 1 + t13*x
610 t16 = 1/t5
611 t19 = q*t6 + t1*t15*t16
612 t22 = p*x
613 t24 = 1/t9/t8
614 t25 = t19**2
615 t32 = t4**2
616 t34 = 1/t32/t4
617 t36 = -t34*t2 + t13
618 t39 = t15**2
619 t40 = t5**2
620 t41 = 1/t40
621 t44 = 2*q*t15*t16 + t1*t36*t16 - t1*t39*t41
622 t48 = p*t2
623 t49 = t9**2
624 t65 = t32**2
625 t87 = -6*p*t10*t19 + 12*t22*t24*t25 - 6*t22*t10*t44 - 6*t48/t49*t25*t19 + &
626 6*t48*t24*t19*t44 - t48*t10*(3*q*t36*t16 - 3*q*t39*t41 + 3*t1*(1/t65/t4* &
627 t2*x - t34*x)*t16 - 3*t1*t36*t41*t15 + 2*t1*t39*t15/t40/t5)
628
629 fs(ip, 4) = t87
630 fs(ip, 4) = f0*f0*f0*fs(ip, 4)
631
632 CASE DEFAULT
633 cpabort("Illegal order")
634 END SELECT
635 END DO
636
637!$OMP END PARALLEL DO
638
639 END SUBROUTINE efactor_b88
640! **************************************************************************************************
641!> \brief ...
642!> \param s ...
643!> \param fs ...
644!> \param m ...
645! **************************************************************************************************
646 SUBROUTINE efactor_pw86(s, fs, m)
647 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
648 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
649 INTEGER, INTENT(IN) :: m
650
651 INTEGER :: ip
652 REAL(kind=dp) :: f15, p0, p1, p15, s2, s4, s6, t1, t10, &
653 t12, t13, t14, t19, t2, t25, t3, t8, t9
654
655 t1 = 1.296_dp
656 t2 = 14.0_dp
657 t3 = 0.2_dp
658 f15 = 1.0_dp/15.0_dp
659
660!$OMP PARALLEL DO DEFAULT(NONE) &
661!$OMP SHARED(s, fs, m, t1, t2, t3, f15) &
662!$OMP PRIVATE(ip,s2,s4,s6,p0,p1,p15)&
663!$OMP PRIVATE(t8, t9, t10, t12, t13, t14, t19, t25)
664
665 DO ip = 1, SIZE(s)
666 s2 = s(ip)*s(ip)
667 s4 = s2*s2
668 s6 = s2*s4
669 SELECT CASE (m)
670 CASE (0)
671 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
672 fs(ip, 1) = p0**f15
673 CASE (1)
674 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
675 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
676 p15 = p0**f15
677 fs(ip, 1) = p15
678 fs(ip, 2) = f15*p1*p15/p0
679 CASE (2)
680 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
681 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
682 p15 = p0**f15
683 fs(ip, 1) = p15
684 fs(ip, 2) = f15*p1*p15/p0
685 t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
686 t25 = p1*p1
687 fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
688 1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
689 CASE (3)
690 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
691 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
692 p15 = p0**f15
693 fs(ip, 1) = p15
694 fs(ip, 2) = f15*p1*p15/p0
695 t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
696 t25 = p1*p1
697 fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
698 1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
699 t8 = p0**2; t9 = p0**f15; t14 = p0/t9; t19 = s2*s(ip)
700 fs(ip, 4) = 406.0_dp/3375.0_dp/t14/t8*p1*p1*p1 - 14.0_dp/ &
701 75.0_dp/t14/p0*p1*(2*t1 + 12*t2*s2 + 30*t3*s4) + &
702 1/t14*(24*t2*s(ip) + 120*t3*t19)*f15
703 CASE DEFAULT
704 cpabort("Illegal order")
705 END SELECT
706 END DO
707
708!$OMP END PARALLEL DO
709
710 END SUBROUTINE efactor_pw86
711! **************************************************************************************************
712!> \brief ...
713!> \param s ...
714!> \param fs ...
715!> \param m ...
716! **************************************************************************************************
717 SUBROUTINE efactor_ev93(s, fs, m)
718 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
719 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: fs
720 INTEGER, INTENT(IN) :: m
721
722 INTEGER :: ip
723 REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, d0, d1, d2, d3, &
724 f0, f1, f2, n0, n1, n2, n3, s2, s4, &
725 s6, scale_s, ss
726
727 a1 = 1.647127_dp
728 a2 = 0.980118_dp
729 a3 = 0.017399_dp
730 b1 = 1.523671_dp
731 b2 = 0.367229_dp
732 b3 = 0.011282_dp
733 scale_s = 1._dp/tact
734
735!$OMP PARALLEL DO DEFAULT(NONE) &
736!$OMP SHARED(s, fs, m, a1, a2, a3, b1, b2, b3, scale_s) &
737!$OMP PRIVATE(ip,ss,s2,s4,s6) &
738!$OMP PRIVATE(n0,n1,n2,n3,d0,d1,d2,d3,f0,f1,f2)
739
740 DO ip = 1, SIZE(s)
741! ss = s(ip)
742!
743 ss = scale_s*s(ip)
744 s2 = ss*ss
745 s4 = s2*s2
746 s6 = s2*s4
747 SELECT CASE (m)
748 CASE (0)
749 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
750 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
751 fs(ip, 1) = n0/d0
752 CASE (1)
753 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
754 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
755 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
756 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
757 f0 = n0/d0
758 fs(ip, 1) = f0
759 fs(ip, 2) = (n1 - f0*d1)/d0*scale_s
760 CASE (2)
761 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
762 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
763 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
764 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
765 n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
766 d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
767 f0 = n0/d0
768 f1 = (n1 - f0*d1)/d0
769 fs(ip, 1) = f0
770 fs(ip, 2) = f1*scale_s
771 fs(ip, 3) = (n2 - f0*d2 - 2._dp*f1*d1)/d0*scale_s*scale_s
772 CASE (3)
773 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
774 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
775 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
776 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
777 n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
778 d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
779 n3 = ss*(24._dp*a2 + 120._dp*a3*s2)
780 d3 = ss*(24._dp*b2 + 120._dp*b3*s2)
781 f0 = n0/d0
782 f1 = (n1 - f0*d1)/d0
783 f2 = (n2 - f0*d2 - 2._dp*f1*d1)/d0
784 fs(ip, 1) = f0
785 fs(ip, 2) = f1*scale_s
786 fs(ip, 3) = f2*scale_s*scale_s
787 fs(ip, 4) = (n3 - f0*d3 - 3._dp*f2*d1 - 3._dp*f1*d2)/d0*scale_s*scale_s*scale_s
788 CASE DEFAULT
789 cpabort("Illegal order")
790 END SELECT
791 END DO
792
793!$OMP END PARALLEL DO
794
795 END SUBROUTINE efactor_ev93
796! **************************************************************************************************
797!> \brief ...
798!> \param s ...
799!> \param fs ...
800!> \param m ...
801! **************************************************************************************************
802 SUBROUTINE efactor_optx(s, fs, m)
803 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
804 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
805 INTEGER, INTENT(IN) :: m
806
807 REAL(kind=dp), PARAMETER :: a1 = 1.05151_dp, a2 = 1.43169_dp, &
808 gamma_bo = 0.006_dp
809
810 INTEGER :: ip
811 REAL(kind=dp) :: a, b, f0, x, y
812
813 f0 = 1.0_dp/sfac
814 b = -a2/flsd
815
816!$OMP PARALLEL DO DEFAULT(NONE) &
817!$OMP SHARED(s, fs, m, f0, b) &
818!$OMP PRIVATE(ip,x,A,y)
819
820 DO ip = 1, SIZE(s)
821 x = s(ip)*f0
822 a = gamma_bo*x*x
823 y = 1.0_dp/(1.0_dp + a)
824 SELECT CASE (m)
825 CASE (0)
826 fs(ip, 1) = a1 + b*a*a*y*y
827 CASE (1)
828 fs(ip, 1) = a1 + b*a*a*y*y
829 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
830 CASE (2)
831 fs(ip, 1) = a1 + b*a*a*y*y
832 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
833 fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
834 CASE (3)
835 fs(ip, 1) = a1 + b*a*a*y*y
836 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
837 fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
838 fs(ip, 4) = 24.0_dp*b*f0*f0*f0*gamma_bo*gamma_bo*x* &
839 (1.0_dp - 5.0_dp*a + 2.0_dp*a*a)*y*y*y*y*y
840 CASE DEFAULT
841 cpabort("Illegal order")
842 END SELECT
843 END DO
844
845!$OMP END PARALLEL DO
846
847 END SUBROUTINE efactor_optx
848
849! **************************************************************************************************
850!> \brief ...
851!> \param s ...
852!> \param fs ...
853!> \param m ...
854! **************************************************************************************************
855 SUBROUTINE efactor_pw91(s, fs, m)
856 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
857 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
858 INTEGER, INTENT(IN) :: m
859
860 INTEGER :: ip
861 REAL(dp) :: t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
862 t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
863 t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
864 t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
865 t96, t98
866 REAL(kind=dp) :: a1, a2, a3, a4, a5, b, o, x
867
868 o = 1.0_dp
869 a1 = 0.19645_dp
870 a2 = 0.2743_dp
871 a3 = 0.1508_dp
872 a4 = 100.0_dp
873 b = 0.8145161_dp
874 a5 = 0.004_dp
875 IF (m >= 0) THEN
876
877!$OMP DO
878
879 DO ip = 1, SIZE(s)
880 x = s(ip)
881 t3 = b**2
882 t4 = x**2
883 t7 = sqrt(o + t3*t4)
884 t9 = log(b*x + t7)
885 t10 = a1*x*t9
886 t12 = exp(-a4*t4)
887 t17 = t4**2
888 fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
889 END DO
890
891!$OMP END DO
892
893 END IF
894 IF (m >= 1) THEN
895
896!$OMP DO
897
898 DO ip = 1, SIZE(s)
899 x = s(ip)
900 t2 = b**2
901 t3 = x**2
902 t6 = sqrt(o + t2*t3)
903 t7 = b*x + t6
904 t8 = log(t7)
905 t9 = a1*t8
906 t10 = a1*x
907 t17 = t10*(b + 1/t6*t2*x)/t7
908 t19 = t3*x
909 t21 = exp(-a4*t3)
910 t26 = a2 - a3*t21
911 t30 = t10*t8
912 t31 = t3**2
913 t33 = o + t30 + a5*t31
914 t38 = t33**2
915 fs(ip, 2) = &
916 (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
917 t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
918 END DO
919
920!$OMP END DO
921
922 END IF
923 IF (m >= 2) THEN
924
925!$OMP DO
926
927 DO ip = 1, SIZE(s)
928 x = s(ip)
929 t1 = b**2
930 t2 = x**2
931 t5 = sqrt(o + t1*t2)
932 t7 = o/t5*t1
933 t9 = b + t7*x
934 t12 = b*x + t5
935 t13 = o/t12
936 t15 = 2._dp*a1*t9*t13
937 t16 = a1*x
938 t17 = t5**2
939 t20 = t1**2
940 t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
941 t26 = t9**2
942 t27 = t12**2
943 t30 = t16*t26/t27
944 t31 = a3*a4
945 t33 = exp(-a4*t2)
946 t37 = a4**2
947 t39 = t2**2
948 t44 = a3*t33
949 t47 = log(t12)
950 t48 = t16*t47
951 t50 = o + t48 + a5*t39
952 t53 = a1*t47
953 t55 = t16*t9*t13
954 t56 = t2*x
955 t60 = a2 - t44
956 t64 = t50**2
957 t65 = o/t64
958 t69 = t53 + t55 + 4._dp*a5*t56
959 t73 = o + t48 + t60*t2
960 t77 = t69**2
961 fs(ip, 3) = &
962 (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
963 2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
964 (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
965 t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
966 END DO
967
968!$OMP END DO
969
970 END IF
971 IF (m >= 3) THEN
972
973!$OMP DO
974
975 DO ip = 1, SIZE(s)
976 x = s(ip)
977 t1 = b**2
978 t2 = x**2
979 t5 = sqrt(0.1e1_dp + t1*t2)
980 t6 = t5**2
981 t9 = t1**2
982 t10 = 1/t6/t5*t9
983 t13 = 1/t5*t1
984 t14 = -t10*t2 + t13
985 t17 = b*x + t5
986 t18 = 1/t17
987 t20 = 3*a1*t14*t18
988 t22 = b + t13*x
989 t23 = t22**2
990 t25 = t17**2
991 t26 = 1/t25
992 t28 = 3*a1*t23*t26
993 t29 = a1*x
994 t30 = t6**2
995 t35 = t2*x
996 t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
997 t44 = 3*t29*t14*t26*t22
998 t50 = 2*t29*t23*t22/t25/t17
999 t51 = a3*a4
1000 t53 = exp(-a4*t2)
1001 t57 = a4**2
1002 t58 = a3*t57
1003 t59 = t35*t53
1004 t64 = t2**2
1005 t70 = log(t17)
1006 t71 = t29*t70
1007 t73 = 0.1e1_dp + t71 + a5*t64
1008 t78 = 2*a1*t22*t18
1009 t80 = t29*t14*t18
1010 t82 = t29*t23*t26
1011 t90 = a3*t53
1012 t93 = t73**2
1013 t94 = 1/t93
1014 t96 = a1*t70
1015 t98 = t29*t18*t22
1016 t101 = t96 + t98 + 4*a5*t35
1017 t106 = a2 - t90
1018 t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1019 t111 = 1/t93/t73
1020 t113 = t101**2
1021 t119 = t78 + t80 - t82 + 12*a5*t2
1022 t123 = 0.1e1_dp + t71 + t106*t2
1023 t124 = t93**2
1024 fs(ip, 4) = &
1025 (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
1026 x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
1027 4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
1028 6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
1029 6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
1030 END DO
1031
1032!$OMP END DO
1033
1034 END IF
1035
1036 END SUBROUTINE efactor_pw91
1037
1038! **************************************************************************************************
1039!> \brief ...
1040!> \param s ...
1041!> \param fs ...
1042!> \param m ...
1043!> \param pset ...
1044! **************************************************************************************************
1045 SUBROUTINE efactor_pbex(s, fs, m, pset)
1046 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: s
1047 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1048 INTEGER, INTENT(IN) :: m, pset
1049
1050 REAL(kind=dp), PARAMETER :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
1051 mu = 0.2195149727645171_dp
1052
1053 INTEGER :: ip
1054 REAL(kind=dp) :: f0, mk, x, x2, y
1055
1056 IF (pset == 1) mk = mu/kappa1
1057 IF (pset == 2) mk = mu/kappa2
1058
1059 f0 = 1.0_dp/tact
1060
1061!$OMP PARALLEL DO DEFAULT(NONE) &
1062!$OMP SHARED (s, fs, m, mk, f0) &
1063!$OMP PRIVATE(ip,x,x2,y)
1064
1065 DO ip = 1, SIZE(s)
1066 x = s(ip)*f0
1067 x2 = x*x
1068 y = 1.0_dp/(1.0_dp + mk*x2)
1069 SELECT CASE (m)
1070 CASE (0)
1071 fs(ip, 1) = 1.0_dp + mu*x2*y
1072 CASE (1)
1073 fs(ip, 1) = 1.0_dp + mu*x2*y
1074 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1075 CASE (2)
1076 fs(ip, 1) = 1.0_dp + mu*x2*y
1077 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1078 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1079 CASE (3)
1080 fs(ip, 1) = 1.0_dp + mu*x2*y
1081 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1082 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1083 fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
1084 CASE DEFAULT
1085 cpabort("Illegal order")
1086 END SELECT
1087 END DO
1088
1089!$OMP END PARALLEL DO
1090
1091 END SUBROUTINE efactor_pbex
1092
1093END MODULE xc_exchange_gga
1094
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:65
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:45
xc_input_constants::xgga_opt
integer, parameter, public xgga_opt
Definition xc_input_constants.F:45
xc_input_constants::xgga_revpbe
integer, parameter, public xgga_revpbe
Definition xc_input_constants.F:45
xc_input_constants::xgga_pw86
integer, parameter, public xgga_pw86
Definition xc_input_constants.F:45
xc_input_constants::xgga_pbe
integer, parameter, public xgga_pbe
Definition xc_input_constants.F:45
xc_input_constants::xgga_pw91
integer, parameter, public xgga_pw91
Definition xc_input_constants.F:45
xc_input_constants::xgga_ev93
integer, parameter, public xgga_ev93
Definition xc_input_constants.F:45
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:48
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