Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 28 : 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 28 : IF (PRESENT(reference)) THEN
85 4 : 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. Zang 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 4 : reference = "E. Engel and S.H. Vosko, Phys. Rev. B, 47, 13164 (1993)"
100 : CASE default
101 4 : CPABORT("Invalid functional requested ("//cp_to_string(functional)//")")
102 : END SELECT
103 4 : IF (.NOT. lsd) THEN
104 4 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
105 4 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
106 : END IF
107 : END IF
108 : END IF
109 28 : IF (PRESENT(shortform)) THEN
110 4 : 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 Zang et al."
121 : CASE (xgga_opt)
122 0 : shortform = "OPTX exchange energy functional"
123 : CASE (xgga_ev93)
124 4 : shortform = "Engel-Vosko exchange energy from virial relation"
125 : CASE default
126 4 : CPABORT("Invalid functional requested ("//cp_to_string(functional)//")")
127 : END SELECT
128 4 : IF (.NOT. lsd) THEN
129 4 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
130 4 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
131 : END IF
132 : END IF
133 : END IF
134 28 : IF (PRESENT(needs)) THEN
135 24 : 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 24 : needs%rho = .TRUE.
141 24 : needs%rho_1_3 = .TRUE.
142 24 : needs%norm_drho = .TRUE.
143 : END IF
144 : END IF
145 28 : IF (PRESENT(max_deriv)) max_deriv = 3
146 :
147 28 : 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 8 : 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 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
173 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
174 8 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
175 8 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
176 8 : e_rho_rho_rho
177 72 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
178 : TYPE(xc_derivative_type), POINTER :: deriv
179 :
180 8 : CALL timeset(routineN, handle)
181 8 : NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
182 8 : e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
183 24 : DO ispin = 1, 2
184 24 : NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
185 : END DO
186 :
187 8 : 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 8 : drho_cutoff=drho_cutoff)
200 8 : nspin = 1
201 8 : rho_spin_name = [deriv_rho, deriv_rho]
202 8 : norm_drho_spin_name = [deriv_norm_drho, deriv_norm_drho]
203 : END IF
204 8 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
205 8 : m = ABS(order)
206 8 : CALL xgga_init(rho_cutoff)
207 :
208 24 : ALLOCATE (s(npoints))
209 32 : ALLOCATE (fs(npoints, m + 1))
210 :
211 16 : DO ispin = 1, nspin
212 8 : 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 8 : fact = flda
218 8 : tact = t13
219 : CALL calc_wave_vector("u", rho(ispin)%array, &
220 8 : norm_drho(ispin)%array, s)
221 : END IF
222 :
223 8 : 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 8 : tact = t13
250 8 : CALL efactor_ev93(s, fs, m)
251 8 : IF (lsd) tact = 1.0_dp
252 : CASE DEFAULT
253 8 : CPABORT("Unsupported functional")
254 : END SELECT
255 :
256 8 : IF (order >= 0) THEN
257 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
258 8 : allocate_deriv=.TRUE.)
259 8 : 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 8 : npoints)
263 : END IF
264 8 : IF (order >= 1 .OR. order == -1) THEN
265 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
266 16 : allocate_deriv=.TRUE.)
267 8 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
268 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
269 16 : allocate_deriv=.TRUE.)
270 8 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
271 :
272 : CALL x_p_1(rho(ispin)%array, &
273 8 : rho_1_3(ispin)%array, s, fs, e_rho, e_ndrho, npoints)
274 : END IF
275 8 : IF (order >= 2 .OR. order == -2) THEN
276 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
277 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
278 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
279 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
280 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
281 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
282 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
283 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
284 0 : 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 0 : e_ndrho_ndrho, npoints)
289 : END IF
290 8 : 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 0 : allocate_deriv=.TRUE.)
294 0 : 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 0 : allocate_deriv=.TRUE.)
298 0 : 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 0 : allocate_deriv=.TRUE.)
302 0 : 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 0 : allocate_deriv=.TRUE.)
306 0 : 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 0 : npoints)
312 : END IF
313 16 : IF (order > 3 .OR. order < -3) THEN
314 0 : CPABORT("derivatives bigger than 3 not implemented")
315 : END IF
316 : END DO
317 :
318 8 : DEALLOCATE (s)
319 8 : DEALLOCATE (fs)
320 :
321 8 : CALL timestop(handle)
322 24 : END SUBROUTINE xgga_eval
323 :
324 : ! **************************************************************************************************
325 : !> \brief ...
326 : !> \param cutoff ...
327 : ! **************************************************************************************************
328 8 : SUBROUTINE xgga_init(cutoff)
329 :
330 : REAL(KIND=dp), INTENT(IN) :: cutoff
331 :
332 8 : eps_rho = cutoff
333 8 : CALL set_util(cutoff)
334 :
335 8 : cx = -0.75_dp*(3.0_dp/pi)**f13
336 8 : t13 = 2.0_dp**f13
337 8 : flda = cx
338 8 : flsd = cx*t13
339 :
340 8 : sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
341 :
342 8 : 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 8 : 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 8 : !$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 8 : 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 8 : 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 8 : !$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 8 : 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 0 : 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 0 : !$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 0 : 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 0 : 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 0 : !$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 0 : END SUBROUTINE x_p_3
542 :
543 : ! Enhancement Factors
544 : ! **************************************************************************************************
545 : !> \brief ...
546 : !> \param s ...
547 : !> \param fs ...
548 : !> \param m ...
549 : ! **************************************************************************************************
550 0 : 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 0 : beta = 0.0042_dp
561 0 : f0 = 1.0_dp/sfac
562 0 : p = -beta/flsd
563 0 : 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 0 : !$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 0 : END SUBROUTINE efactor_b88
642 : ! **************************************************************************************************
643 : !> \brief ...
644 : !> \param s ...
645 : !> \param fs ...
646 : !> \param m ...
647 : ! **************************************************************************************************
648 0 : 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 0 : t1 = 1.296_dp
658 0 : t2 = 14.0_dp
659 0 : t3 = 0.2_dp
660 0 : 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 0 : !$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 0 : END SUBROUTINE efactor_pw86
713 : ! **************************************************************************************************
714 : !> \brief ...
715 : !> \param s ...
716 : !> \param fs ...
717 : !> \param m ...
718 : ! **************************************************************************************************
719 8 : 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 8 : a1 = 1.647127_dp
730 8 : a2 = 0.980118_dp
731 8 : a3 = 0.017399_dp
732 8 : b1 = 1.523671_dp
733 8 : b2 = 0.367229_dp
734 8 : b3 = 0.011282_dp
735 8 : 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 8 : !$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 8 : END SUBROUTINE efactor_ev93
798 : ! **************************************************************************************************
799 : !> \brief ...
800 : !> \param s ...
801 : !> \param fs ...
802 : !> \param m ...
803 : ! **************************************************************************************************
804 0 : 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 0 : f0 = 1.0_dp/sfac
816 0 : b = -a2/flsd
817 :
818 : !$OMP PARALLEL DO DEFAULT(NONE) &
819 : !$OMP SHARED(s, fs, m, f0, b) &
820 0 : !$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 0 : END SUBROUTINE efactor_optx
850 :
851 : ! **************************************************************************************************
852 : !> \brief ...
853 : !> \param s ...
854 : !> \param fs ...
855 : !> \param m ...
856 : ! **************************************************************************************************
857 0 : 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 0 : o = 1.0_dp
871 0 : a1 = 0.19645_dp
872 0 : a2 = 0.2743_dp
873 0 : a3 = 0.1508_dp
874 0 : a4 = 100.0_dp
875 0 : b = 0.8145161_dp
876 0 : a5 = 0.004_dp
877 0 : IF (m >= 0) THEN
878 :
879 0 : !$OMP DO
880 :
881 : DO ip = 1, SIZE(s)
882 0 : x = s(ip)
883 0 : t3 = b**2
884 0 : t4 = x**2
885 0 : t7 = SQRT(o + t3*t4)
886 0 : t9 = LOG(b*x + t7)
887 0 : t10 = a1*x*t9
888 0 : t12 = EXP(-a4*t4)
889 0 : t17 = t4**2
890 0 : 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 0 : IF (m >= 1) THEN
897 :
898 0 : !$OMP DO
899 :
900 : DO ip = 1, SIZE(s)
901 0 : x = s(ip)
902 0 : t2 = b**2
903 0 : t3 = x**2
904 0 : t6 = SQRT(o + t2*t3)
905 0 : t7 = b*x + t6
906 0 : t8 = LOG(t7)
907 0 : t9 = a1*t8
908 0 : t10 = a1*x
909 0 : t17 = t10*(b + 1/t6*t2*x)/t7
910 0 : t19 = t3*x
911 0 : t21 = EXP(-a4*t3)
912 0 : t26 = a2 - a3*t21
913 0 : t30 = t10*t8
914 0 : t31 = t3**2
915 0 : t33 = o + t30 + a5*t31
916 0 : t38 = t33**2
917 : fs(ip, 2) = &
918 : (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
919 0 : 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 0 : IF (m >= 2) THEN
926 :
927 0 : !$OMP DO
928 :
929 : DO ip = 1, SIZE(s)
930 0 : x = s(ip)
931 0 : t1 = b**2
932 0 : t2 = x**2
933 0 : t5 = SQRT(o + t1*t2)
934 0 : t7 = o/t5*t1
935 0 : t9 = b + t7*x
936 0 : t12 = b*x + t5
937 0 : t13 = o/t12
938 0 : t15 = 2._dp*a1*t9*t13
939 0 : t16 = a1*x
940 0 : t17 = t5**2
941 0 : t20 = t1**2
942 0 : t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
943 0 : t26 = t9**2
944 0 : t27 = t12**2
945 0 : t30 = t16*t26/t27
946 0 : t31 = a3*a4
947 0 : t33 = EXP(-a4*t2)
948 0 : t37 = a4**2
949 0 : t39 = t2**2
950 0 : t44 = a3*t33
951 0 : t47 = LOG(t12)
952 0 : t48 = t16*t47
953 0 : t50 = o + t48 + a5*t39
954 0 : t53 = a1*t47
955 0 : t55 = t16*t9*t13
956 0 : t56 = t2*x
957 0 : t60 = a2 - t44
958 0 : t64 = t50**2
959 0 : t65 = o/t64
960 0 : t69 = t53 + t55 + 4._dp*a5*t56
961 0 : t73 = o + t48 + t60*t2
962 0 : 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 0 : 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 0 : IF (m >= 3) THEN
974 :
975 0 : !$OMP DO
976 :
977 : DO ip = 1, SIZE(s)
978 0 : x = s(ip)
979 0 : t1 = b**2
980 0 : t2 = x**2
981 0 : t5 = SQRT(0.1e1_dp + t1*t2)
982 0 : t6 = t5**2
983 0 : t9 = t1**2
984 0 : t10 = 1/t6/t5*t9
985 0 : t13 = 1/t5*t1
986 0 : t14 = -t10*t2 + t13
987 0 : t17 = b*x + t5
988 0 : t18 = 1/t17
989 0 : t20 = 3*a1*t14*t18
990 0 : t22 = b + t13*x
991 0 : t23 = t22**2
992 0 : t25 = t17**2
993 0 : t26 = 1/t25
994 0 : t28 = 3*a1*t23*t26
995 0 : t29 = a1*x
996 0 : t30 = t6**2
997 0 : t35 = t2*x
998 0 : t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
999 0 : t44 = 3*t29*t14*t26*t22
1000 0 : t50 = 2*t29*t23*t22/t25/t17
1001 0 : t51 = a3*a4
1002 0 : t53 = EXP(-a4*t2)
1003 0 : t57 = a4**2
1004 0 : t58 = a3*t57
1005 0 : t59 = t35*t53
1006 0 : t64 = t2**2
1007 0 : t70 = LOG(t17)
1008 0 : t71 = t29*t70
1009 0 : t73 = 0.1e1_dp + t71 + a5*t64
1010 0 : t78 = 2*a1*t22*t18
1011 0 : t80 = t29*t14*t18
1012 0 : t82 = t29*t23*t26
1013 0 : t90 = a3*t53
1014 0 : t93 = t73**2
1015 0 : t94 = 1/t93
1016 0 : t96 = a1*t70
1017 0 : t98 = t29*t18*t22
1018 0 : t101 = t96 + t98 + 4*a5*t35
1019 0 : t106 = a2 - t90
1020 0 : t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1021 0 : t111 = 1/t93/t73
1022 0 : t113 = t101**2
1023 0 : t119 = t78 + t80 - t82 + 12*a5*t2
1024 0 : t123 = 0.1e1_dp + t71 + t106*t2
1025 0 : 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 0 : 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 0 : END SUBROUTINE efactor_pw91
1039 :
1040 : ! **************************************************************************************************
1041 : !> \brief ...
1042 : !> \param s ...
1043 : !> \param fs ...
1044 : !> \param m ...
1045 : !> \param pset ...
1046 : ! **************************************************************************************************
1047 0 : 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 0 : IF (pset == 2) mk = mu/kappa2
1060 :
1061 0 : f0 = 1.0_dp/tact
1062 :
1063 : !$OMP PARALLEL DO DEFAULT(NONE) &
1064 : !$OMP SHARED (s, fs, m, mk, f0) &
1065 0 : !$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 0 : END SUBROUTINE efactor_pbex
1094 :
1095 : END MODULE xc_exchange_gga
1096 :
|
