Line data Source code
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 : ! **************************************************************************************************
16 : MODULE 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 :
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 14 : 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 14 : IF (PRESENT(reference)) THEN
85 2 : SELECT CASE (functional)
86 : CASE (xgga_b88)
87 0 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988)"
88 : CASE (xgga_pw86)
89 0 : reference = "J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
90 : CASE (xgga_pw91)
91 0 : reference = "J.P. Perdew et al., Phys. Rev. B, 46, 6671 (1992)"
92 : CASE (xgga_pbe)
93 0 : reference = "J.P. Perdew, K. Burke, M Ernzerhof, Phys. Rev. Lett, 77, 3865 (1996)"
94 : CASE (xgga_revpbe)
95 0 : reference = "Y. Zhang et al., PRL, 80, 890 (1998) (Revised PBEX)"
96 : CASE (xgga_opt)
97 0 : reference = "Wee-Meng Hoe, A.J. Cohen, N.C. Handy, CPL, 341, 319 (2001)"
98 : CASE (xgga_ev93)
99 2 : reference = "E. Engel and S.H. Vosko, Phys. Rev. B, 47, 13164 (1993)"
100 : CASE default
101 2 : CPABORT("Invalid functional requested ("//cp_to_string(functional)//")")
102 : END SELECT
103 2 : IF (.NOT. lsd) THEN
104 2 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
105 2 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
106 : END IF
107 : END IF
108 : END IF
109 14 : IF (PRESENT(shortform)) THEN
110 2 : SELECT CASE (functional)
111 : CASE (xgga_b88)
112 0 : shortform = "Becke 1988 Exchange Functional"
113 : CASE (xgga_pw86)
114 0 : shortform = "Perdew-Wang 1986 Functional (exchange energy)"
115 : CASE (xgga_pw91)
116 0 : shortform = "Perdew-Wang 1991 Functional (exchange energy)"
117 : CASE (xgga_pbe)
118 0 : shortform = "PBE exchange energy functional"
119 : CASE (xgga_revpbe)
120 0 : shortform = "Revised PBEX by Zhang et al."
121 : CASE (xgga_opt)
122 0 : shortform = "OPTX exchange energy functional"
123 : CASE (xgga_ev93)
124 2 : shortform = "Engel-Vosko exchange energy from virial relation"
125 : CASE default
126 2 : CPABORT("Invalid functional requested ("//cp_to_string(functional)//")")
127 : END SELECT
128 2 : IF (.NOT. lsd) THEN
129 2 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
130 2 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
131 : END IF
132 : END IF
133 : END IF
134 14 : IF (PRESENT(needs)) THEN
135 12 : IF (lsd) THEN
136 0 : needs%rho_spin = .TRUE.
137 0 : needs%rho_spin_1_3 = .TRUE.
138 0 : needs%norm_drho_spin = .TRUE.
139 : ELSE
140 12 : needs%rho = .TRUE.
141 12 : needs%rho_1_3 = .TRUE.
142 12 : needs%norm_drho = .TRUE.
143 : END IF
144 : END IF
145 14 : IF (PRESENT(max_deriv)) max_deriv = 3
146 :
147 14 : 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 4 : 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 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
173 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
174 4 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
175 4 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
176 4 : e_rho_rho_rho
177 36 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
178 : TYPE(xc_derivative_type), POINTER :: deriv
179 :
180 4 : CALL timeset(routineN, handle)
181 4 : NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
182 4 : e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
183 12 : DO ispin = 1, 2
184 12 : NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
185 : END DO
186 :
187 4 : 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 0 : drho_cutoff=drho_cutoff, local_bounds=bo)
193 0 : nspin = 2
194 0 : rho_spin_name = [deriv_rhoa, deriv_rhob]
195 0 : 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 4 : drho_cutoff=drho_cutoff)
200 4 : nspin = 1
201 4 : rho_spin_name = [deriv_rho, deriv_rho]
202 4 : norm_drho_spin_name = [deriv_norm_drho, deriv_norm_drho]
203 : END IF
204 4 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
205 4 : m = ABS(order)
206 4 : CALL xgga_init(rho_cutoff)
207 :
208 12 : ALLOCATE (s(npoints))
209 16 : ALLOCATE (fs(npoints, m + 1))
210 :
211 8 : DO ispin = 1, nspin
212 4 : IF (lsd) THEN
213 0 : fact = flsd
214 0 : tact = 1.0_dp
215 0 : CALL calc_wave_vector("p", rho(ispin)%array, norm_drho(ispin)%array, s)
216 : ELSE
217 4 : fact = flda
218 4 : tact = t13
219 : CALL calc_wave_vector("u", rho(ispin)%array, &
220 4 : norm_drho(ispin)%array, s)
221 : END IF
222 :
223 4 : SELECT CASE (functional)
224 : CASE (xgga_b88)
225 0 : CALL efactor_b88(s, fs, m)
226 : CASE (xgga_pw86)
227 0 : 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 0 : !$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 0 : tact = t13
240 0 : CALL efactor_pbex(s, fs, m, 1)
241 0 : IF (lsd) tact = 1.0_dp
242 : CASE (xgga_revpbe)
243 0 : tact = t13
244 0 : CALL efactor_pbex(s, fs, m, 2)
245 0 : IF (lsd) tact = 1.0_dp
246 : CASE (xgga_opt)
247 0 : CALL efactor_optx(s, fs, m)
248 : CASE (xgga_ev93)
249 4 : tact = t13
250 4 : CALL efactor_ev93(s, fs, m)
251 4 : IF (lsd) tact = 1.0_dp
252 : CASE DEFAULT
253 4 : CPABORT("Unsupported functional")
254 : END SELECT
255 :
256 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
257 4 : allocate_deriv=.TRUE.)
258 4 : 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 4 : npoints)
262 4 : IF (order >= 1) THEN
263 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
264 8 : allocate_deriv=.TRUE.)
265 4 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
266 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
267 8 : allocate_deriv=.TRUE.)
268 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
269 :
270 : CALL x_p_1(rho(ispin)%array, &
271 4 : rho_1_3(ispin)%array, s, fs, e_rho, e_ndrho, npoints)
272 : END IF
273 4 : IF (order >= 2) THEN
274 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
275 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
276 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
277 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
278 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
279 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
280 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
281 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
282 0 : 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 0 : e_ndrho_ndrho, npoints)
287 : END IF
288 4 : 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 0 : allocate_deriv=.TRUE.)
292 0 : 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 0 : allocate_deriv=.TRUE.)
296 0 : 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 0 : allocate_deriv=.TRUE.)
300 0 : 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 0 : allocate_deriv=.TRUE.)
304 0 : 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 0 : npoints)
310 : END IF
311 8 : IF (order > 3) THEN
312 0 : CPABORT("derivatives bigger than 3 not implemented")
313 : END IF
314 : END DO
315 :
316 4 : DEALLOCATE (s)
317 4 : DEALLOCATE (fs)
318 :
319 4 : CALL timestop(handle)
320 12 : END SUBROUTINE xgga_eval
321 :
322 : ! **************************************************************************************************
323 : !> \brief ...
324 : !> \param cutoff ...
325 : ! **************************************************************************************************
326 4 : SUBROUTINE xgga_init(cutoff)
327 :
328 : REAL(KIND=dp), INTENT(IN) :: cutoff
329 :
330 4 : eps_rho = cutoff
331 4 : CALL set_util(cutoff)
332 :
333 4 : cx = -0.75_dp*(3.0_dp/pi)**f13
334 4 : t13 = 2.0_dp**f13
335 4 : flda = cx
336 4 : flsd = cx*t13
337 :
338 4 : sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
339 :
340 4 : 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 4 : 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 4 : !$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 4 : 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 4 : 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 4 : !$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 4 : 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 0 : 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 0 : !$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 0 : 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 0 : 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 0 : !$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 0 : END SUBROUTINE x_p_3
540 :
541 : ! Enhancement Factors
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param s ...
545 : !> \param fs ...
546 : !> \param m ...
547 : ! **************************************************************************************************
548 0 : 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 0 : beta = 0.0042_dp
559 0 : f0 = 1.0_dp/sfac
560 0 : p = -beta/flsd
561 0 : 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 0 : !$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 0 : END SUBROUTINE efactor_b88
640 : ! **************************************************************************************************
641 : !> \brief ...
642 : !> \param s ...
643 : !> \param fs ...
644 : !> \param m ...
645 : ! **************************************************************************************************
646 0 : 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 0 : t1 = 1.296_dp
656 0 : t2 = 14.0_dp
657 0 : t3 = 0.2_dp
658 0 : 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 0 : !$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 0 : END SUBROUTINE efactor_pw86
711 : ! **************************************************************************************************
712 : !> \brief ...
713 : !> \param s ...
714 : !> \param fs ...
715 : !> \param m ...
716 : ! **************************************************************************************************
717 4 : 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 4 : a1 = 1.647127_dp
728 4 : a2 = 0.980118_dp
729 4 : a3 = 0.017399_dp
730 4 : b1 = 1.523671_dp
731 4 : b2 = 0.367229_dp
732 4 : b3 = 0.011282_dp
733 4 : 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 4 : !$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 4 : END SUBROUTINE efactor_ev93
796 : ! **************************************************************************************************
797 : !> \brief ...
798 : !> \param s ...
799 : !> \param fs ...
800 : !> \param m ...
801 : ! **************************************************************************************************
802 0 : 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 0 : f0 = 1.0_dp/sfac
814 0 : b = -a2/flsd
815 :
816 : !$OMP PARALLEL DO DEFAULT(NONE) &
817 : !$OMP SHARED(s, fs, m, f0, b) &
818 0 : !$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 0 : END SUBROUTINE efactor_optx
848 :
849 : ! **************************************************************************************************
850 : !> \brief ...
851 : !> \param s ...
852 : !> \param fs ...
853 : !> \param m ...
854 : ! **************************************************************************************************
855 0 : 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 0 : o = 1.0_dp
869 0 : a1 = 0.19645_dp
870 0 : a2 = 0.2743_dp
871 0 : a3 = 0.1508_dp
872 0 : a4 = 100.0_dp
873 0 : b = 0.8145161_dp
874 0 : a5 = 0.004_dp
875 0 : IF (m >= 0) THEN
876 :
877 0 : !$OMP DO
878 :
879 : DO ip = 1, SIZE(s)
880 0 : x = s(ip)
881 0 : t3 = b**2
882 0 : t4 = x**2
883 0 : t7 = SQRT(o + t3*t4)
884 0 : t9 = LOG(b*x + t7)
885 0 : t10 = a1*x*t9
886 0 : t12 = EXP(-a4*t4)
887 0 : t17 = t4**2
888 0 : 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 0 : IF (m >= 1) THEN
895 :
896 0 : !$OMP DO
897 :
898 : DO ip = 1, SIZE(s)
899 0 : x = s(ip)
900 0 : t2 = b**2
901 0 : t3 = x**2
902 0 : t6 = SQRT(o + t2*t3)
903 0 : t7 = b*x + t6
904 0 : t8 = LOG(t7)
905 0 : t9 = a1*t8
906 0 : t10 = a1*x
907 0 : t17 = t10*(b + 1/t6*t2*x)/t7
908 0 : t19 = t3*x
909 0 : t21 = EXP(-a4*t3)
910 0 : t26 = a2 - a3*t21
911 0 : t30 = t10*t8
912 0 : t31 = t3**2
913 0 : t33 = o + t30 + a5*t31
914 0 : t38 = t33**2
915 : fs(ip, 2) = &
916 : (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
917 0 : 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 0 : IF (m >= 2) THEN
924 :
925 0 : !$OMP DO
926 :
927 : DO ip = 1, SIZE(s)
928 0 : x = s(ip)
929 0 : t1 = b**2
930 0 : t2 = x**2
931 0 : t5 = SQRT(o + t1*t2)
932 0 : t7 = o/t5*t1
933 0 : t9 = b + t7*x
934 0 : t12 = b*x + t5
935 0 : t13 = o/t12
936 0 : t15 = 2._dp*a1*t9*t13
937 0 : t16 = a1*x
938 0 : t17 = t5**2
939 0 : t20 = t1**2
940 0 : t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
941 0 : t26 = t9**2
942 0 : t27 = t12**2
943 0 : t30 = t16*t26/t27
944 0 : t31 = a3*a4
945 0 : t33 = EXP(-a4*t2)
946 0 : t37 = a4**2
947 0 : t39 = t2**2
948 0 : t44 = a3*t33
949 0 : t47 = LOG(t12)
950 0 : t48 = t16*t47
951 0 : t50 = o + t48 + a5*t39
952 0 : t53 = a1*t47
953 0 : t55 = t16*t9*t13
954 0 : t56 = t2*x
955 0 : t60 = a2 - t44
956 0 : t64 = t50**2
957 0 : t65 = o/t64
958 0 : t69 = t53 + t55 + 4._dp*a5*t56
959 0 : t73 = o + t48 + t60*t2
960 0 : 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 0 : 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 0 : IF (m >= 3) THEN
972 :
973 0 : !$OMP DO
974 :
975 : DO ip = 1, SIZE(s)
976 0 : x = s(ip)
977 0 : t1 = b**2
978 0 : t2 = x**2
979 0 : t5 = SQRT(0.1e1_dp + t1*t2)
980 0 : t6 = t5**2
981 0 : t9 = t1**2
982 0 : t10 = 1/t6/t5*t9
983 0 : t13 = 1/t5*t1
984 0 : t14 = -t10*t2 + t13
985 0 : t17 = b*x + t5
986 0 : t18 = 1/t17
987 0 : t20 = 3*a1*t14*t18
988 0 : t22 = b + t13*x
989 0 : t23 = t22**2
990 0 : t25 = t17**2
991 0 : t26 = 1/t25
992 0 : t28 = 3*a1*t23*t26
993 0 : t29 = a1*x
994 0 : t30 = t6**2
995 0 : t35 = t2*x
996 0 : t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
997 0 : t44 = 3*t29*t14*t26*t22
998 0 : t50 = 2*t29*t23*t22/t25/t17
999 0 : t51 = a3*a4
1000 0 : t53 = EXP(-a4*t2)
1001 0 : t57 = a4**2
1002 0 : t58 = a3*t57
1003 0 : t59 = t35*t53
1004 0 : t64 = t2**2
1005 0 : t70 = LOG(t17)
1006 0 : t71 = t29*t70
1007 0 : t73 = 0.1e1_dp + t71 + a5*t64
1008 0 : t78 = 2*a1*t22*t18
1009 0 : t80 = t29*t14*t18
1010 0 : t82 = t29*t23*t26
1011 0 : t90 = a3*t53
1012 0 : t93 = t73**2
1013 0 : t94 = 1/t93
1014 0 : t96 = a1*t70
1015 0 : t98 = t29*t18*t22
1016 0 : t101 = t96 + t98 + 4*a5*t35
1017 0 : t106 = a2 - t90
1018 0 : t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1019 0 : t111 = 1/t93/t73
1020 0 : t113 = t101**2
1021 0 : t119 = t78 + t80 - t82 + 12*a5*t2
1022 0 : t123 = 0.1e1_dp + t71 + t106*t2
1023 0 : 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 0 : 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 0 : END SUBROUTINE efactor_pw91
1037 :
1038 : ! **************************************************************************************************
1039 : !> \brief ...
1040 : !> \param s ...
1041 : !> \param fs ...
1042 : !> \param m ...
1043 : !> \param pset ...
1044 : ! **************************************************************************************************
1045 0 : 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 0 : IF (pset == 2) mk = mu/kappa2
1058 :
1059 0 : f0 = 1.0_dp/tact
1060 :
1061 : !$OMP PARALLEL DO DEFAULT(NONE) &
1062 : !$OMP SHARED (s, fs, m, mk, f0) &
1063 0 : !$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 0 : END SUBROUTINE efactor_pbex
1092 :
1093 : END MODULE xc_exchange_gga
1094 :
|
