Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate the several different kinetic 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 (20.02.2002)
15 : ! **************************************************************************************************
16 : MODULE xc_ke_gga
17 :
18 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: pi
21 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
22 : deriv_norm_drhoa,&
23 : deriv_norm_drhob,&
24 : deriv_rho,&
25 : deriv_rhoa,&
26 : deriv_rhob
27 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
28 : xc_dset_get_derivative
29 : USE xc_derivative_types, ONLY: xc_derivative_get,&
30 : xc_derivative_type
31 : USE xc_functionals_utilities, ONLY: calc_wave_vector,&
32 : set_util
33 : USE xc_input_constants, ONLY: ke_lc,&
34 : ke_llp,&
35 : ke_ol1,&
36 : ke_ol2,&
37 : ke_pbe,&
38 : ke_pw86,&
39 : ke_pw91,&
40 : ke_t92
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 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
51 : f23 = 2.0_dp*f13, &
52 : f43 = 4.0_dp*f13, &
53 : f53 = 5.0_dp*f13
54 :
55 : PUBLIC :: ke_gga_info, ke_gga_lda_eval, ke_gga_lsd_eval
56 :
57 : REAL(KIND=dp) :: cf, b, b_lda, b_lsd, flda, flsd, sfac, t13
58 : REAL(KIND=dp) :: fact, tact
59 : REAL(KIND=dp) :: eps_rho
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_ke_gga'
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief ...
67 : !> \param cutoff ...
68 : ! **************************************************************************************************
69 1648 : SUBROUTINE ke_gga_init(cutoff)
70 :
71 : REAL(KIND=dp), INTENT(IN) :: cutoff
72 :
73 1648 : eps_rho = cutoff
74 1648 : CALL set_util(cutoff)
75 :
76 1648 : cf = 0.3_dp*(3.0_dp*pi*pi)**f23
77 1648 : flda = cf
78 1648 : flsd = flda*2.0_dp**f23
79 : ! the_factor 2^(1/3) for LDA is here
80 1648 : b_lda = 2.0_dp**f43*(3.0_dp*pi*pi)**(f13)
81 1648 : b_lsd = 2.0_dp*(3.0_dp*pi*pi)**(f13)
82 1648 : sfac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
83 1648 : t13 = 2.0_dp**f13
84 :
85 1648 : END SUBROUTINE ke_gga_init
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param functional ...
90 : !> \param lsd ...
91 : !> \param reference ...
92 : !> \param shortform ...
93 : !> \param needs ...
94 : !> \param max_deriv ...
95 : ! **************************************************************************************************
96 1786 : SUBROUTINE ke_gga_info(functional, lsd, reference, shortform, needs, max_deriv)
97 : INTEGER, INTENT(in) :: functional
98 : LOGICAL, INTENT(in) :: lsd
99 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
100 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
101 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
102 :
103 1786 : IF (PRESENT(reference)) THEN
104 0 : SELECT CASE (functional)
105 : CASE (ke_ol1)
106 : reference = "H. Ou-Yang and M. Levy, "// &
107 0 : "Intl. J. Quant. Chem. 40, 379 (1991); Functional 1"
108 : CASE (ke_ol2)
109 : reference = "H. Ou-Yang and M. Levy, "// &
110 0 : "Intl. J. Quant. Chem. 40, 379 (1991); Functional 2"
111 : CASE (ke_llp)
112 0 : reference = "H. Lee, C. Lee, R.G. Parr, Phys. Rev. A, 44, 768 (1991)"
113 : CASE (ke_pw86)
114 0 : reference = "J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
115 : CASE (ke_pw91)
116 0 : reference = "J.P. Perdew and Y. Wang, Electronic Structure of Solids 91"
117 : CASE (ke_lc)
118 0 : reference = "A. Lembarki and H. Chermette, Phys. Rev. A, 50, 5328 (1994)"
119 : CASE (ke_t92)
120 0 : reference = "A.J. Thakkar, Phys. Rev. A, 46, 6920 (1992)"
121 : CASE (ke_pbe)
122 0 : reference = "J.P.Perdew, K.Burke, M.Ernzerhof, Phys. Rev. Letter, 77, 3865 (1996)"
123 : END SELECT
124 0 : IF (.NOT. lsd) THEN
125 0 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
126 0 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {spin unpolarized}'
127 : END IF
128 : END IF
129 : END IF
130 1786 : IF (PRESENT(shortform)) THEN
131 0 : SELECT CASE (functional)
132 : CASE (ke_ol1)
133 0 : shortform = "Ou-Yang-Levy Functional 1"
134 : CASE (ke_ol2)
135 0 : shortform = "Ou-Yang-Levy Functional 2"
136 : CASE (ke_llp)
137 0 : shortform = "Lee-Lee-Parr Functional"
138 : CASE (ke_pw86)
139 0 : shortform = "Perdew-Wang 1986 Functional (kinetic energy)"
140 : CASE (ke_pw91)
141 0 : shortform = "Perdew-Wang 1991 Functional (kinetic energy)"
142 : CASE (ke_lc)
143 0 : shortform = "Lembarki-Chermette kinetic energy functional"
144 : CASE (ke_t92)
145 0 : shortform = "Thakkar 1992 Functional"
146 : CASE (ke_pbe)
147 0 : shortform = "Perdew-Burke-Ernzerhof Functional (kinetic energy)"
148 : END SELECT
149 0 : IF (.NOT. lsd) THEN
150 0 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
151 0 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {spin unpolarized}'
152 : END IF
153 : END IF
154 : END IF
155 1786 : IF (PRESENT(needs)) THEN
156 1786 : IF (lsd) THEN
157 0 : needs%rho_spin = .TRUE.
158 0 : needs%rho_spin_1_3 = .TRUE.
159 0 : needs%norm_drho_spin = .TRUE.
160 : ELSE
161 1786 : needs%rho = .TRUE.
162 1786 : needs%rho_1_3 = .TRUE.
163 1786 : needs%norm_drho = .TRUE.
164 : END IF
165 : END IF
166 1786 : IF (PRESENT(max_deriv)) max_deriv = 3
167 :
168 1786 : END SUBROUTINE ke_gga_info
169 :
170 : ! **************************************************************************************************
171 : !> \brief ...
172 : !> \param functional ...
173 : !> \param rho_set ...
174 : !> \param deriv_set ...
175 : !> \param order ...
176 : ! **************************************************************************************************
177 1648 : SUBROUTINE ke_gga_lda_eval(functional, rho_set, deriv_set, order)
178 :
179 : INTEGER, INTENT(IN) :: functional
180 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
181 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
182 : INTEGER, INTENT(IN) :: order
183 :
184 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ke_gga_lda_eval'
185 :
186 : INTEGER :: handle, m, npoints
187 : INTEGER, DIMENSION(2, 3) :: bo
188 : REAL(KIND=dp) :: drho_cutoff, rho_cutoff
189 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
190 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
191 1648 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
192 1648 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
193 1648 : e_rho_rho_rho, grho, rho, rho13
194 : TYPE(xc_derivative_type), POINTER :: deriv
195 :
196 1648 : CALL timeset(routineN, handle)
197 1648 : NULLIFY (rho, rho13, e_0, e_rho, e_ndrho, &
198 1648 : e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
199 1648 : e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho)
200 1648 : m = ABS(order)
201 :
202 : CALL xc_rho_set_get(rho_set, rho_1_3=rho13, rho=rho, &
203 : norm_drho=grho, local_bounds=bo, rho_cutoff=rho_cutoff, &
204 1648 : drho_cutoff=drho_cutoff)
205 1648 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
206 1648 : CALL ke_gga_init(rho_cutoff)
207 :
208 4944 : ALLOCATE (s(npoints))
209 6592 : ALLOCATE (fs(npoints, m + 1))
210 :
211 : ! s = norm_drho/(rho^(4/3)*2*(pi*pi*3)^(1/3))
212 1648 : CALL calc_wave_vector("p", rho, grho, s)
213 :
214 1648 : fact = flda
215 : ! Definition of s has changed
216 1648 : b = b_lda
217 : ! tact = t13
218 1648 : tact = 1.0_dp
219 :
220 1648 : SELECT CASE (functional)
221 : CASE (ke_ol1)
222 0 : CALL efactor_ol1(s, fs, m)
223 0 : CPABORT("OL1 functional currently not working properly")
224 : CASE (ke_ol2)
225 0 : CALL efactor_ol2(s, fs, m)
226 0 : CPABORT("OL2 functional currently not working properly")
227 : CASE (ke_llp)
228 736 : CALL efactor_llp(s, fs, m)
229 : CASE (ke_pw86)
230 54 : CALL efactor_pw86(s, fs, m)
231 : CASE (ke_pw91)
232 54 : CALL efactor_pw91(s, fs, m, 1)
233 : CASE (ke_lc)
234 54 : CALL efactor_pw91(s, fs, m, 2)
235 : CASE (ke_t92)
236 174 : CALL efactor_t92(s, fs, m)
237 : CASE (ke_pbe)
238 576 : CALL efactor_pbex(s, fs, m, 1)
239 : CASE DEFAULT
240 1648 : CPABORT("Unsupported functional")
241 : END SELECT
242 :
243 1648 : IF (order >= 0) THEN
244 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
245 1648 : allocate_deriv=.TRUE.)
246 1648 : CALL xc_derivative_get(deriv, deriv_data=e_0)
247 :
248 1648 : CALL kex_p_0(rho, rho13, fs, e_0, npoints)
249 : END IF
250 :
251 1648 : IF (order >= 1 .OR. order == -1) THEN
252 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
253 1648 : allocate_deriv=.TRUE.)
254 1648 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
255 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
256 1648 : allocate_deriv=.TRUE.)
257 1648 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
258 :
259 : CALL kex_p_1(rho, rho13, s, fs, e_rho=e_rho, e_ndrho=e_ndrho, &
260 1648 : npoints=npoints)
261 : END IF
262 1648 : IF (order >= 2 .OR. order == -2) THEN
263 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
264 138 : allocate_deriv=.TRUE.)
265 138 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
266 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
267 138 : allocate_deriv=.TRUE.)
268 138 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
269 : deriv => xc_dset_get_derivative(deriv_set, &
270 138 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
271 138 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
272 :
273 : CALL kex_p_2(rho, rho13, s, fs, e_rho_rho=e_rho_rho, &
274 138 : e_rho_ndrho=e_rho_ndrho, e_ndrho_ndrho=e_ndrho_ndrho, npoints=npoints)
275 : END IF
276 1648 : IF (order >= 3 .OR. order == -3) THEN
277 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
278 0 : allocate_deriv=.TRUE.)
279 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
280 : deriv => xc_dset_get_derivative(deriv_set, &
281 0 : [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
282 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
283 : deriv => xc_dset_get_derivative(deriv_set, &
284 0 : [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
285 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
286 : deriv => xc_dset_get_derivative(deriv_set, &
287 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
288 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
289 :
290 : CALL kex_p_3(rho, rho13, s, fs, e_rho_rho_rho=e_rho_rho_rho, &
291 : e_rho_rho_ndrho=e_rho_rho_ndrho, e_rho_ndrho_ndrho=e_rho_ndrho_ndrho, &
292 0 : e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, npoints=npoints)
293 : END IF
294 1648 : IF (order > 3 .OR. order < -3) THEN
295 0 : CPABORT("derivatives bigger than 3 not implemented")
296 : END IF
297 :
298 1648 : DEALLOCATE (s)
299 1648 : DEALLOCATE (fs)
300 :
301 1648 : CALL timestop(handle)
302 :
303 1648 : END SUBROUTINE ke_gga_lda_eval
304 :
305 : ! **************************************************************************************************
306 : !> \brief ...
307 : !> \param functional ...
308 : !> \param rho_set ...
309 : !> \param deriv_set ...
310 : !> \param order ...
311 : ! **************************************************************************************************
312 0 : SUBROUTINE ke_gga_lsd_eval(functional, rho_set, deriv_set, order)
313 :
314 : INTEGER, INTENT(IN) :: functional
315 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
316 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
317 : INTEGER, INTENT(IN), OPTIONAL :: order
318 :
319 : CHARACTER(len=*), PARAMETER :: routineN = 'ke_gga_lsd_eval'
320 : INTEGER, DIMENSION(2), PARAMETER :: &
321 : norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
322 : rho_spin_name = [deriv_rhoa, deriv_rhob]
323 :
324 : INTEGER :: handle, ispin, m, npoints
325 : INTEGER, DIMENSION(2, 3) :: bo
326 : REAL(KIND=dp) :: rho_cutoff
327 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
328 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fs
329 0 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
330 0 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
331 0 : e_rho_rho_rho
332 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
333 : TYPE(xc_derivative_type), POINTER :: deriv
334 :
335 0 : CALL timeset(routineN, handle)
336 0 : NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
337 0 : e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
338 0 : DO ispin = 1, 2
339 0 : NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
340 : END DO
341 :
342 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
343 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
344 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
345 : norm_drhob=norm_drho(2)%array, rho_cutoff=rho_cutoff, &
346 0 : local_bounds=bo)
347 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
348 0 : m = ABS(order)
349 0 : CALL ke_gga_init(rho_cutoff)
350 :
351 0 : ALLOCATE (s(npoints))
352 0 : ALLOCATE (fs(npoints, m + 1))
353 :
354 0 : fact = flsd
355 0 : b = b_lsd
356 0 : tact = 1.0_dp
357 :
358 0 : DO ispin = 1, 2
359 :
360 0 : CALL calc_wave_vector("p", rho(ispin)%array, norm_drho(ispin)%array, s)
361 :
362 0 : SELECT CASE (functional)
363 : CASE (ke_ol1)
364 0 : CALL efactor_ol1(s, fs, m)
365 : CASE (ke_ol2)
366 0 : CALL efactor_ol2(s, fs, m)
367 : CASE (ke_llp)
368 0 : CALL efactor_llp(s, fs, m)
369 : CASE (ke_pw86)
370 0 : tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
371 0 : CALL efactor_pw86(s, fs, m, f2_lsd=tact)
372 0 : tact = 1.0_dp
373 : CASE (ke_pbe)
374 0 : tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
375 0 : CALL efactor_pbex(s, fs, m, 1, f2_lsd=tact)
376 0 : tact = 1.0_dp
377 : CASE (ke_pw91)
378 0 : tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
379 0 : CALL efactor_pw91(s, fs, m, 1, f2_lsd=tact)
380 0 : tact = 1.0_dp
381 : CASE (ke_lc)
382 0 : tact = (1.0_dp/2.0_dp)**(1.0_dp/3.0_dp)
383 0 : CALL efactor_pw91(s, fs, m, 2, f2_lsd=tact)
384 0 : tact = 1.0_dp
385 : CASE (ke_t92)
386 0 : CALL efactor_t92(s, fs, m)
387 : CASE DEFAULT
388 0 : CPABORT("Unsupported functional")
389 : END SELECT
390 :
391 0 : IF (order >= 0) THEN
392 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
393 0 : allocate_deriv=.TRUE.)
394 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
395 :
396 : CALL kex_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, &
397 0 : e_0=e_0, npoints=npoints)
398 : END IF
399 0 : IF (order >= 1 .OR. order == -1) THEN
400 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
401 0 : allocate_deriv=.TRUE.)
402 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
403 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
404 0 : allocate_deriv=.TRUE.)
405 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
406 :
407 : CALL kex_p_1(rho=rho(ispin)%array, &
408 : r13=rho_1_3(ispin)%array, s=s, fs=fs, e_rho=e_rho, &
409 0 : e_ndrho=e_ndrho, npoints=npoints)
410 : END IF
411 0 : IF (order >= 2 .OR. order == -2) THEN
412 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
413 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
414 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
415 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
416 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
417 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
418 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
419 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
420 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
421 :
422 : CALL kex_p_2(rho(ispin)%array, rho_1_3(ispin)%array, &
423 0 : s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, npoints)
424 : END IF
425 0 : IF (order >= 3 .OR. order == -3) THEN
426 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
427 : rho_spin_name(ispin), rho_spin_name(ispin)], &
428 0 : allocate_deriv=.TRUE.)
429 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
430 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
431 : rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
432 0 : allocate_deriv=.TRUE.)
433 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
434 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
435 : norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
436 0 : allocate_deriv=.TRUE.)
437 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
438 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
439 : norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
440 0 : allocate_deriv=.TRUE.)
441 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
442 :
443 : CALL kex_p_3(rho(ispin)%array, &
444 : rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
445 0 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
446 : END IF
447 0 : IF (order > 3 .OR. order < -3) THEN
448 0 : CPABORT("derivatives bigger than 3 not implemented")
449 : END IF
450 :
451 : END DO
452 :
453 0 : DEALLOCATE (s)
454 0 : DEALLOCATE (fs)
455 0 : CALL timestop(handle)
456 0 : END SUBROUTINE ke_gga_lsd_eval
457 :
458 : ! **************************************************************************************************
459 : !> \brief ...
460 : !> \param rho ...
461 : !> \param r13 ...
462 : !> \param fs ...
463 : !> \param e_0 ...
464 : !> \param npoints ...
465 : ! **************************************************************************************************
466 1648 : SUBROUTINE kex_p_0(rho, r13, fs, e_0, npoints)
467 :
468 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13
469 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: fs
470 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
471 : INTEGER, INTENT(in) :: npoints
472 :
473 : INTEGER :: ip
474 :
475 : !$OMP PARALLEL DO DEFAULT(NONE) &
476 : !$OMP SHARED(npoints, rho, e_0, fact, r13, fs, eps_rho) &
477 1648 : !$OMP PRIVATE(ip)
478 :
479 : DO ip = 1, npoints
480 :
481 : IF (rho(ip) > eps_rho) THEN
482 : e_0(ip) = e_0(ip) + fact*r13(ip)*r13(ip)*rho(ip)*fs(ip, 1)
483 : END IF
484 :
485 : END DO
486 :
487 : !$OMP END PARALLEL DO
488 :
489 1648 : END SUBROUTINE kex_p_0
490 :
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param rho ...
494 : !> \param r13 ...
495 : !> \param s ...
496 : !> \param fs ...
497 : !> \param e_rho ...
498 : !> \param e_ndrho ...
499 : !> \param npoints ...
500 : ! **************************************************************************************************
501 1648 : SUBROUTINE kex_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
502 :
503 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
504 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: fs
505 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
506 : INTEGER, INTENT(in) :: npoints
507 :
508 : INTEGER :: ip
509 : REAL(KIND=dp) :: a0, a1, sx, sy
510 :
511 : !$OMP PARALLEL DO DEFAULT(NONE) &
512 : !$OMP SHARED(npoints, rho, eps_rho, fact, r13, sfac, tact) &
513 : !$OMP SHARED(fs, e_rho, e_ndrho, s) &
514 1648 : !$OMP PRIVATE(ip,a0,a1,sx,sy)
515 :
516 : DO ip = 1, npoints
517 :
518 : IF (rho(ip) > eps_rho) THEN
519 :
520 : a0 = fact*r13(ip)*r13(ip)*rho(ip)
521 : a1 = f53*fact*r13(ip)*r13(ip)
522 : sx = -f43*s(ip)/rho(ip)
523 : sy = sfac*tact/(r13(ip)*rho(ip))
524 : e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
525 : e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
526 :
527 : END IF
528 :
529 : END DO
530 :
531 : !$OMP END PARALLEL DO
532 :
533 1648 : END SUBROUTINE kex_p_1
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param rho ...
538 : !> \param r13 ...
539 : !> \param s ...
540 : !> \param fs ...
541 : !> \param e_rho_rho ...
542 : !> \param e_rho_ndrho ...
543 : !> \param e_ndrho_ndrho ...
544 : !> \param npoints ...
545 : ! **************************************************************************************************
546 138 : SUBROUTINE kex_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
547 : npoints)
548 :
549 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
550 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: fs
551 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
552 : INTEGER, INTENT(in) :: npoints
553 :
554 : INTEGER :: ip
555 : REAL(KIND=dp) :: a0, a1, a2, sx, sxx, sxy, sy
556 :
557 : !$OMP PARALLEL DO DEFAULT(NONE) &
558 : !$OMP SHARED (npoints, rho, eps_rho, fact, r13) &
559 : !$OMP SHARED (e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, fs) &
560 : !$OMP SHARED(s, sfac, tact) &
561 138 : !$OMP PRIVATE(ip,a0,a1,a2,sx,sy,sxx,sxy)
562 :
563 : DO ip = 1, npoints
564 :
565 : IF (rho(ip) > eps_rho) THEN
566 :
567 : a0 = fact*r13(ip)*r13(ip)*rho(ip)
568 : a1 = f53*fact*r13(ip)*r13(ip)
569 : a2 = f23*f53*fact/r13(ip)
570 : sx = -f43*s(ip)/rho(ip)
571 : sy = sfac*tact/(r13(ip)*rho(ip))
572 : sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
573 : sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
574 : e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
575 : a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
576 : e_rho_ndrho(ip) = e_rho_ndrho(ip) + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + &
577 : a0*fs(ip, 2)*sxy
578 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
579 :
580 : END IF
581 :
582 : END DO
583 :
584 : !$OMP END PARALLEL DO
585 :
586 138 : END SUBROUTINE kex_p_2
587 :
588 : ! **************************************************************************************************
589 : !> \brief ...
590 : !> \param rho ...
591 : !> \param r13 ...
592 : !> \param s ...
593 : !> \param fs ...
594 : !> \param e_rho_rho_rho ...
595 : !> \param e_rho_rho_ndrho ...
596 : !> \param e_rho_ndrho_ndrho ...
597 : !> \param e_ndrho_ndrho_ndrho ...
598 : !> \param npoints ...
599 : ! **************************************************************************************************
600 0 : SUBROUTINE kex_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
601 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
602 :
603 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
604 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: fs
605 : REAL(KIND=dp), DIMENSION(*), INTENT(inout) :: e_rho_rho_rho, e_rho_rho_ndrho, &
606 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
607 : INTEGER, INTENT(in) :: npoints
608 :
609 : INTEGER :: ip
610 : REAL(KIND=dp) :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
611 : sxy, sy
612 :
613 : !$OMP PARALLEL DO DEFAULT(NONE) &
614 : !$OMP SHARED(npoints, rho, eps_rho, fact, r13) &
615 : !$OMP SHARED(s, sfac, tact, fs, e_rho_rho_rho) &
616 : !$OMP SHARED(e_rho_rho_ndrho, e_rho_ndrho_ndrho) &
617 : !$OMP SHARED(e_ndrho_ndrho_ndrho) &
618 0 : !$OMP PRIVATE(ip,a0,a1,a2,a3,sx,sy,sxx,sxy,sxxx,sxxy)
619 :
620 : DO ip = 1, npoints
621 :
622 : IF (rho(ip) > eps_rho) THEN
623 :
624 : a0 = fact*r13(ip)*r13(ip)*rho(ip)
625 : a1 = f53*fact*r13(ip)*r13(ip)
626 : a2 = f23*f53*fact/r13(ip)
627 : a3 = -f13*f23*f53*fact/(r13(ip)*rho(ip))
628 : sx = -f43*s(ip)/rho(ip)
629 : sy = sfac*tact/(r13(ip)*rho(ip))
630 : sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
631 : sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
632 : sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
633 : sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
634 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
635 : 3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
636 : a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
637 : a0*fs(ip, 2)*sxxx
638 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
639 : 2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
640 : 2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
641 : a0*fs(ip, 2)*sxxy
642 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
643 : 2.0_dp*a0*fs(ip, 3)*sxy*sy
644 : e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) + a0*fs(ip, 4)*sy*sy*sy
645 :
646 : END IF
647 :
648 : END DO
649 :
650 : !$OMP END PARALLEL DO
651 :
652 0 : END SUBROUTINE kex_p_3
653 :
654 : ! Enhancement Factors
655 : ! **************************************************************************************************
656 : !> \brief ...
657 : !> \param s ...
658 : !> \param fs ...
659 : !> \param m ...
660 : ! **************************************************************************************************
661 0 : SUBROUTINE efactor_ol1(s, fs, m)
662 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
663 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
664 : INTEGER, INTENT(IN) :: m
665 :
666 : INTEGER :: ip
667 : REAL(KIND=dp) :: t1, t2
668 :
669 0 : t1 = b*b/(72.0_dp*cf)
670 0 : t2 = 0.001878_dp*b
671 :
672 : !$OMP PARALLEL DO DEFAULT(NONE) &
673 : !$OMP SHARED(s, m, fs, t1, t2) &
674 0 : !$OMP PRIVATE(ip)
675 :
676 : DO ip = 1, SIZE(s)
677 : SELECT CASE (m)
678 : CASE (0)
679 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
680 : CASE (1)
681 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
682 : fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
683 : CASE (2:3)
684 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)
685 : fs(ip, 2) = 2.0_dp*t1*s(ip) + t2
686 : fs(ip, 3) = 2.0_dp*t1
687 : CASE DEFAULT
688 : CPABORT("Illegal order.")
689 : END SELECT
690 : END DO
691 :
692 : !$OMP END PARALLEL DO
693 :
694 0 : IF (m == 3) fs(:, 4) = 0.0_dp
695 :
696 0 : END SUBROUTINE efactor_ol1
697 : ! **************************************************************************************************
698 : !> \brief ...
699 : !> \param s ...
700 : !> \param fs ...
701 : !> \param m ...
702 : ! **************************************************************************************************
703 0 : SUBROUTINE efactor_ol2(s, fs, m)
704 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
705 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
706 : INTEGER, INTENT(IN) :: m
707 :
708 : INTEGER :: ip
709 : REAL(KIND=dp) :: t1, t2, t3, y
710 :
711 0 : t1 = b*b/(72.0_dp*cf)
712 0 : t2 = 0.0245_dp*b
713 0 : t3 = 2.0_dp**f53*b
714 :
715 : !$OMP PARALLEL DO DEFAULT(NONE) &
716 : !$OMP SHARED(s, m, t1, t2, t3, fs) &
717 0 : !$OMP PRIVATE(ip,y)
718 :
719 : DO ip = 1, SIZE(s)
720 : y = 1.0_dp/(1.0_dp + t3*s(ip))
721 : SELECT CASE (m)
722 : CASE (0)
723 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
724 : CASE (1)
725 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
726 : fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
727 : CASE (2)
728 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
729 : fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
730 : fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
731 : CASE (3)
732 : fs(ip, 1) = 1.0_dp + t1*s(ip)*s(ip) + t2*s(ip)*y
733 : fs(ip, 2) = 2.0_dp*t1*s(ip) + t2*y*y
734 : fs(ip, 3) = 2.0_dp*(t1 - t2*t3*y*y*y)
735 : fs(ip, 4) = 6.0_dp*t2*t3*t3*y*y*y*y
736 : CASE DEFAULT
737 : CPABORT("Illegal order.")
738 : END SELECT
739 : END DO
740 :
741 : !$OMP END PARALLEL DO
742 :
743 0 : END SUBROUTINE efactor_ol2
744 : ! **************************************************************************************************
745 : !> \brief ...
746 : !> \param s ...
747 : !> \param fs ...
748 : !> \param m ...
749 : ! **************************************************************************************************
750 736 : SUBROUTINE efactor_llp(s, fs, m)
751 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
752 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
753 : INTEGER, INTENT(IN) :: m
754 :
755 : INTEGER :: ip
756 : REAL(KIND=dp) :: as, bs, p, q, sas, sbs, t1, t10, t11, t12, t133, t16, t17, t19, t2, t20, &
757 : t22, t23, t26, t28, t29, t3, t30, t33, t36, t39, t4, t40, t42, t43, t45, t46, t47, t49, &
758 : t5, t50, t54, t55, t7, t71, t8, t9, x, ys
759 :
760 736 : p = 0.0044188_dp*b*b
761 736 : q = 0.0253_dp*b
762 :
763 : !$OMP PARALLEL DO DEFAULT(NONE) &
764 : !$OMP SHARED(s, m, fs, p, q, b) &
765 : !$OMP PRIVATE(ip,x,bs,sbs,as,sas,ys,t2,t4,t5,t8,t9,t10,t12) &
766 : !$OMP PRIVATE(t1,t3,t7,t11,t16,t17,t20,t22,t23,t26,t30,t33) &
767 : !$OMP PRIVATE(t40,t42,t43,t45,t46,t47,t49,t50,t71,t133) &
768 736 : !$OMP PRIVATE(t19,t28,t29,t36,t39,t54,t55)
769 :
770 : DO ip = 1, SIZE(s)
771 : x = s(ip)
772 : bs = b*x
773 : sbs = SQRT(bs*bs + 1.0_dp)
774 : as = LOG(bs + sbs)
775 : sas = x*as
776 : ys = 1.0_dp/(1.0_dp + q*sas)
777 : SELECT CASE (m)
778 : CASE (0)
779 : fs(ip, 1) = 1.0_dp + p*x*x*ys
780 : CASE (1)
781 : fs(ip, 1) = 1.0_dp + p*x*x*ys
782 : t2 = q*x
783 : t4 = b**2
784 : t5 = x**2
785 : t8 = SQRT(1.0_dp + t4*t5)
786 : t9 = b*x + t8
787 : t10 = LOG(t9)
788 : t12 = 1.0_dp + t2*t10
789 : t17 = t12**2
790 : fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
791 :
792 : CASE (2)
793 : fs(ip, 1) = 1.0_dp + p*x*x*ys
794 : ! first der
795 : t2 = q*x
796 : t4 = b**2
797 : t5 = x**2
798 : t8 = SQRT(1.0_dp + t4*t5)
799 : t9 = b*x + t8
800 : t10 = LOG(t9)
801 : t12 = 1.0_dp + t2*t10
802 : t17 = t12**2
803 : fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
804 :
805 : ! second der
806 : t1 = q*x
807 : t3 = b**2
808 : t4 = x**2
809 : t7 = SQRT(1.0_dp + t3*t4)
810 : t8 = b*x + t7
811 : t9 = LOG(t8)
812 : t11 = 1.0_dp + t1*t9
813 : t16 = t11**2
814 : t17 = 1.0_dp/t16
815 : t20 = 1.0_dp/t7*t3
816 : t22 = b + t20*x
817 : t23 = 1/t8
818 : t26 = q*t9 + t1*t22*t23
819 : t30 = p*t4
820 : t33 = t26**2
821 : t40 = t7**2
822 : t43 = t3**2
823 : t49 = t22**2
824 : t50 = t8**2
825 : fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
826 : t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
827 : (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
828 :
829 : CASE (3)
830 :
831 : fs(ip, 1) = 1.0_dp + p*x*x*ys
832 : ! first der
833 : t2 = q*x
834 : t4 = b**2
835 : t5 = x**2
836 : t8 = SQRT(1.0_dp + t4*t5)
837 : t9 = b*x + t8
838 : t10 = LOG(t9)
839 : t12 = 1.0_dp + t2*t10
840 : t17 = t12**2
841 : fs(ip, 2) = 2.0_dp*p*x/t12 - p*t5/t17*(q*t10 + t2*(b + 1.0_dp/t8*t4*x)/t9)
842 :
843 : ! second der
844 : t1 = q*x
845 : t3 = b**2
846 : t4 = x**2
847 : t7 = SQRT(1.0_dp + t3*t4)
848 : t8 = b*x + t7
849 : t9 = LOG(t8)
850 : t11 = 1.0_dp + t1*t9
851 : t16 = t11**2
852 : t17 = 1.0_dp/t16
853 : t20 = 1.0_dp/t7*t3
854 : t22 = b + t20*x
855 : t23 = 1/t8
856 : t26 = q*t9 + t1*t22*t23
857 : t30 = p*t4
858 : t33 = t26**2
859 : t40 = t7**2
860 : t43 = t3**2
861 : t49 = t22**2
862 : t50 = t8**2
863 : fs(ip, 3) = 2.0_dp*p/t11 - 4.0_dp*p*x*t17*t26 + 2.0_dp*t30/t16/ &
864 : t11*t33 - t30*t17*(2.0_dp*q*t22*t23 + t1* &
865 : (-1.0_dp/t40/t7*t43*t4 + t20)*t23 - t1*t49/t50)
866 :
867 : t1 = q*x
868 : t3 = b**2
869 : t4 = x**2
870 : t7 = SQRT(1 + t3*t4)
871 : t8 = b*x + t7
872 : t9 = LOG(t8)
873 : t11 = 1.0_dp + t1*t9
874 : t12 = t11**2
875 : t133 = 1.0_dp/t12
876 : t17 = 1.0_dp/t7*t3
877 : t19 = b + t17*x
878 : t20 = 1.0_dp/t8
879 : t23 = q*t9 + t1*t19*t20
880 : t26 = p*x
881 : t28 = 1.0_dp/t12/t11
882 : t29 = t23**2
883 : t36 = t7**2
884 : t39 = t3**2
885 : t40 = 1.0_dp/t36/t7*t39
886 : t42 = -t40*t4 + t17
887 : t45 = t19**2
888 : t46 = t8**2
889 : t47 = 1.0_dp/t46
890 : t50 = 2.0_dp*q*t19*t20 + t1*t42*t20 - t1*t45*t47
891 : t54 = p*t4
892 : t55 = t12**2
893 : t71 = t36**2
894 : fs(ip, 4) = &
895 : -6.0_dp*p*t133*t23 + 12.0_dp*t26*t28*t29 - &
896 : 6.0_dp*t26*t133*t50 - 6.0_dp*t54/t55*t29*t23 + &
897 : 6.0_dp*t54*t28*t23*t50 - t54*t133* &
898 : (3.0_dp*q*t42*t20 - 3.0_dp*q*t45*t47 + 3.0_dp*t1* &
899 : (1.0_dp/t71/t7*t39*t3*t4*x - t40*x)*t20 - &
900 : 3.0_dp*t1*t42*t47*t19 + 2.0_dp*t1*t45*t19/t46/t8)
901 :
902 : CASE DEFAULT
903 : CPABORT("Illegal order.")
904 : END SELECT
905 : END DO
906 :
907 : !$OMP END PARALLEL DO
908 :
909 736 : END SUBROUTINE efactor_llp
910 : ! **************************************************************************************************
911 : !> \brief ...
912 : !> \param s ...
913 : !> \param fs ...
914 : !> \param m ...
915 : !> \param f2_lsd ...
916 : ! **************************************************************************************************
917 54 : SUBROUTINE efactor_pw86(s, fs, m, f2_lsd)
918 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
919 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
920 : INTEGER, INTENT(IN) :: m
921 : REAL(dp), OPTIONAL :: f2_lsd
922 :
923 : INTEGER :: ip
924 : REAL(KIND=dp) :: f15, ff, p0, p1, p15, p2, p3, s1, s2, &
925 : s4, s6, t1, t2, t3
926 :
927 54 : t1 = 1.296_dp
928 54 : t2 = 14.0_dp
929 54 : t3 = 0.2_dp
930 54 : f15 = 1.0_dp/15.0_dp
931 54 : ff = 1.0_dp
932 54 : IF (PRESENT(f2_lsd)) ff = f2_lsd
933 :
934 : !$OMP PARALLEL DO DEFAULT(NONE) &
935 : !$OMP SHARED(s, fs, m, t1, t2, t3, f15, ff) &
936 54 : !$OMP PRIVATE(ip, s1, s2, s4, s6, p0, p1, p2, p3, p15)
937 :
938 : DO ip = 1, SIZE(s)
939 : s1 = s(ip)*ff
940 : s2 = s1*s1
941 : s4 = s2*s2
942 : s6 = s2*s4
943 : SELECT CASE (m)
944 : CASE (0)
945 : p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
946 : fs(ip, 1) = p0**f15
947 : CASE (1)
948 : p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
949 : p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
950 : p15 = p0**f15
951 : fs(ip, 1) = p15
952 : fs(ip, 2) = f15*p1*p15/p0
953 : CASE (2)
954 : p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
955 : p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
956 : p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
957 : p15 = p0**f15
958 : fs(ip, 1) = p15
959 : fs(ip, 2) = f15*p1*p15/p0
960 : fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
961 : CASE (3)
962 : p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
963 : p1 = s1*ff*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
964 : p2 = ff*ff*(2.0_dp*t1 + 12.0_dp*t2*s2 + 30.0_dp*t3*s4)
965 : p3 = s1*ff*ff*ff*(24.0_dp*t2 + 120.0_dp*t3*s2)
966 : p15 = p0**f15
967 : fs(ip, 1) = p15
968 : fs(ip, 2) = f15*p1*p15/p0
969 : fs(ip, 3) = f15*p15/p0*(p2 - 14.0_dp/15.0_dp*p1*p1/p0)
970 : fs(ip, 4) = f15*p15/p0*(-14.0_dp*f15*p1*p1/p0 + 14.0_dp*14.0_dp*f15*p1*p1*p1/p0/p0 + &
971 : p3 - 14.0_dp*p2*p1/p0 + 14.0_dp*p1*p1/p0/p0)
972 : CASE DEFAULT
973 : CPABORT("Illegal order.")
974 : END SELECT
975 : END DO
976 :
977 : !$OMP END PARALLEL DO
978 :
979 54 : END SUBROUTINE efactor_pw86
980 : ! **************************************************************************************************
981 : !> \brief ...
982 : !> \param s ...
983 : !> \param fs ...
984 : !> \param m ...
985 : ! **************************************************************************************************
986 174 : SUBROUTINE efactor_t92(s, fs, m)
987 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
988 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
989 : INTEGER, INTENT(IN) :: m
990 :
991 : INTEGER :: ip
992 : REAL(KIND=dp) :: a1, a2, as, asp, asp2, asp3, bs, p, q, &
993 : s1, s2, sas, sbs, sbs3, sbs5, t0, w1, &
994 : x, ys
995 :
996 174 : p = 0.0055_dp*b*b
997 174 : q = 0.0253_dp*b
998 174 : a1 = 0.072_dp*b
999 174 : a2 = 2.0_dp**f53*b
1000 :
1001 : !$OMP PARALLEL DO DEFAULT(NONE) &
1002 : !$OMP SHARED(s, fs, m, p, q, a1, a2, b) &
1003 : !$OMP PRIVATE(ip, x, bs, sbs, sas, ys, asp, sbs3, asp2, sbs5) &
1004 174 : !$OMP PRIVATE(asp3, as, s2, s1, t0, w1)
1005 :
1006 : DO ip = 1, SIZE(s)
1007 : x = s(ip)
1008 : bs = b*x
1009 : sbs = SQRT(bs*bs + 1.0_dp)
1010 : as = LOG(bs + sbs)
1011 : sas = x*as
1012 : ys = 1.0_dp/(1.0_dp + q*sas)
1013 : SELECT CASE (m)
1014 : CASE (0)
1015 : fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1016 : CASE (1)
1017 : asp = as + bs/sbs
1018 : fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1019 : fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1020 : CASE (2)
1021 : asp = as + bs/sbs
1022 : sbs3 = sbs*sbs*sbs
1023 : asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
1024 : fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1025 : fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1026 : fs(ip, 3) = 2.0_dp*p*ys - p*q*x*(4.0_dp*asp + x*asp2)*ys*ys + &
1027 : 2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
1028 : CASE (3)
1029 : asp = as + bs/sbs
1030 : sbs3 = sbs*sbs*sbs
1031 : sbs5 = sbs3*sbs*sbs
1032 : asp2 = 2.0_dp*b/sbs - b*bs*bs/sbs3
1033 : asp3 = -4.0_dp*b*b*bs/sbs3 + 3.0_dp*b*b*bs*bs*bs/sbs5
1034 : w1 = (4.0_dp*asp + x*asp2)
1035 : fs(ip, 1) = 1.0_dp + p*x*x*ys - a1*x/(1 + a2*x)
1036 : fs(ip, 2) = 2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys - a1/(1 + a2*x)**2
1037 : fs(ip, 3) = 2.0_dp*p*ys - p*q*x*w1*ys*ys + &
1038 : 2.0_dp*p*q*q*x*x*asp*asp*ys*ys*ys + 2.0_dp*a1*a2/(1 + a2*x)**3
1039 :
1040 : s2 = -6*p/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
1041 : q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2))) + 12*p*x/ &
1042 : (1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**3*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
1043 : q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))**2
1044 : s1 = s2 - 6*p*x/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(2*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/ &
1045 : (b*x + SQRT(1 + b**2*x**2)) + q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)*b**2)/ &
1046 : (b*x + SQRT(1 + b**2*x**2)) - q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/ &
1047 : (b*x + SQRT(1 + b**2*x**2))**2) - 6*p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**4 &
1048 : *(q*LOG(b*x + SQRT(1 + b**2*x**2)) + q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))**3
1049 : s2 = s1 + 6*p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**3*(q*LOG(b*x + SQRT(1 + b**2*x**2)) + &
1050 : q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)/(b*x + SQRT(1 + b**2*x**2)))*(2*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x) &
1051 : /(b*x + SQRT(1 + b**2*x**2)) + q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)* &
1052 : b**2)/(b*x + SQRT(1 + b**2*x**2)) - q*x*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/(b*x + SQRT(1 + b**2*x**2))**2)
1053 : t0 = s2 - p*x**2/(1 + q*x*LOG(b*x + SQRT(1 + b**2*x**2)))**2*(3*q*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + &
1054 : 1/SQRT(1 + b**2*x**2)*b**2)/(b*x + SQRT(1 + b**2*x**2)) - 3*q*(b + 1/SQRT(1 + b**2*x**2)*b**2*x)**2/ &
1055 : (b*x + SQRT(1 + b**2*x**2))**2 + q*x*(3/SQRT(1 + b**2*x**2)**5*b**6*x**3 - 3/SQRT(1 + b**2*x**2)**3*b**4*x)/ &
1056 : (b*x + SQRT(1 + b**2*x**2)) - 3*q*x*(-1/SQRT(1 + b**2*x**2)**3*b**4*x**2 + 1/SQRT(1 + b**2*x**2)*b**2)/ &
1057 : (b*x + SQRT(1 + b**2*x**2))**2*(b + 1/SQRT(1 + b**2*x**2)*b**2*x) + 2*q*x*(b + 1/SQRT(1 + b**2*x**2)* &
1058 : b**2*x)**3/(b*x + SQRT(1 + b**2*x**2))**3) - 6*a1/(1 + a2*x)**3*a2**2 + 6*a1*x/(1 + a2*x)**4*a2**3
1059 :
1060 : fs(ip, 4) = t0
1061 :
1062 : CASE DEFAULT
1063 : CPABORT("Illegal order")
1064 : END SELECT
1065 : END DO
1066 :
1067 : !$OMP END PARALLEL DO
1068 :
1069 174 : END SUBROUTINE efactor_t92
1070 : ! **************************************************************************************************
1071 : !> \brief ...
1072 : !> \param s ...
1073 : !> \param fs ...
1074 : !> \param m ...
1075 : !> \param pset ...
1076 : !> \param f2_lsd ...
1077 : ! **************************************************************************************************
1078 576 : SUBROUTINE efactor_pbex(s, fs, m, pset, f2_lsd)
1079 :
1080 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
1081 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1082 : INTEGER, INTENT(IN) :: m, pset
1083 : REAL(dp), OPTIONAL :: f2_lsd
1084 :
1085 : REAL(KIND=dp), PARAMETER :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
1086 : mu = 0.2195149727645171_dp
1087 :
1088 : INTEGER :: ip
1089 : REAL(KIND=dp) :: f0, mk, x, x2, y
1090 :
1091 576 : IF (pset == 1) mk = mu/kappa1
1092 576 : IF (pset == 2) mk = mu/kappa2
1093 :
1094 576 : f0 = 1.0_dp/tact
1095 576 : IF (PRESENT(f2_lsd)) f0 = f2_lsd
1096 :
1097 : !$OMP PARALLEL DO DEFAULT(NONE) &
1098 : !$OMP SHARED(s, m, fs, f0, mk) &
1099 576 : !$OMP PRIVATE(ip,x,x2,y)
1100 :
1101 : DO ip = 1, SIZE(s)
1102 : x = s(ip)*f0
1103 : x2 = x*x
1104 : y = 1.0_dp/(1.0_dp + mk*x2)
1105 : SELECT CASE (m)
1106 : CASE (0)
1107 : fs(ip, 1) = 1.0_dp + mu*x2*y
1108 : CASE (1)
1109 : fs(ip, 1) = 1.0_dp + mu*x2*y
1110 : fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1111 : CASE (2)
1112 : fs(ip, 1) = 1.0_dp + mu*x2*y
1113 : fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1114 : fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1115 : CASE (3)
1116 : fs(ip, 1) = 1.0_dp + mu*x2*y
1117 : fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1118 : fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1119 : fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
1120 : CASE DEFAULT
1121 : CPABORT("Illegal order.")
1122 : END SELECT
1123 : END DO
1124 :
1125 : !$OMP END PARALLEL DO
1126 :
1127 576 : END SUBROUTINE efactor_pbex
1128 :
1129 : ! **************************************************************************************************
1130 : !> \brief ...
1131 : !> \param s ...
1132 : !> \param fs ...
1133 : !> \param m ...
1134 : !> \param pset ...
1135 : !> \param f2_lsd ...
1136 : ! **************************************************************************************************
1137 108 : SUBROUTINE efactor_pw91(s, fs, m, pset, f2_lsd)
1138 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: s
1139 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fs
1140 : INTEGER, INTENT(IN) :: m, pset
1141 : REAL(dp), OPTIONAL :: f2_lsd
1142 :
1143 : INTEGER :: ip
1144 : REAL(dp) :: ff, t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
1145 : t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
1146 : t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
1147 : t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
1148 : t96, t98
1149 : REAL(KIND=dp) :: a1, a2, a3, a4, a5, bb, o, pa(6, 2), x
1150 :
1151 : ! parameter set 1: Perdew-Wang
1152 : ! parameter set 2: Lembarki-Chermette
1153 :
1154 : pa(1:6, 1) = (/0.19645_dp, 0.2743_dp, &
1155 : 0.1508_dp, 100.0_dp, &
1156 756 : 7.7956_dp, 0.004_dp/)
1157 : pa(1:6, 2) = (/0.093907_dp, 0.26608_dp, &
1158 : 0.0809615_dp, 100.0_dp, &
1159 756 : 76.320_dp, 0.57767e-4_dp/)
1160 108 : o = 1.0_dp
1161 108 : ff = 1.0_dp
1162 108 : IF (PRESENT(f2_lsd)) ff = f2_lsd
1163 108 : a1 = pa(1, pset)*FF
1164 108 : a2 = pa(2, pset)*ff*ff
1165 108 : a3 = pa(3, pset)*ff*ff
1166 108 : a4 = pa(4, pset)*ff*ff
1167 108 : bb = pa(5, pset)*ff
1168 : ! it should be valid also for lsd
1169 108 : a5 = pa(6, pset)*ff*ff*ff*ff
1170 :
1171 : !$OMP PARALLEL DEFAULT(NONE) &
1172 : !$OMP SHARED(s, fs, m, a1, a2, a3, a4, a5, bb, pa, o, ff) &
1173 : !$OMP PRIVATE(x,t1,t10,t101,t106,t109,t111,t113,t119,t12,t123) &
1174 : !$OMP PRIVATE(t124,t13,t14,t15,t16,t17,t18,t19,t2,t20,t21,t22) &
1175 : !$OMP PRIVATE(t23,t25,t26,t27,t28,t29,t3,t30,t31,t33,t35,t37) &
1176 : !$OMP PRIVATE(t38,t39,t4,t40,t44,t47,t48,t5,t50,t51,t53,t55) &
1177 : !$OMP PRIVATE(t56,t57,t58,t59,t6,t60,t64,t65,t69,t7,t70,t71) &
1178 108 : !$OMP PRIVATE(t73,t77,t78,t8,t80,t82,t9,t90,t93,t94,t96,t98,ip)
1179 :
1180 : IF (m >= 0) THEN
1181 : !$OMP DO
1182 : DO ip = 1, SIZE(s)
1183 : x = s(ip)
1184 : t3 = bb**2
1185 : t4 = x**2
1186 : t7 = SQRT(o + t3*t4)
1187 : t9 = LOG(bb*x + t7)
1188 : t10 = a1*x*t9
1189 : t12 = EXP(-a4*t4)
1190 : t17 = t4**2
1191 : fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
1192 : END DO
1193 : !$OMP END DO
1194 : END IF
1195 : IF (m >= 1) THEN
1196 : !$OMP DO
1197 : DO ip = 1, SIZE(s)
1198 : x = s(ip)
1199 : t2 = bb**2
1200 : t3 = x**2
1201 : t6 = SQRT(o + t2*t3)
1202 : t7 = bb*x + t6
1203 : t8 = LOG(t7)
1204 : t9 = a1*t8
1205 : t10 = a1*x
1206 : t17 = t10*(bb + 1/t6*t2*x)/t7
1207 : t19 = t3*x
1208 : t21 = EXP(-a4*t3)
1209 : t26 = a2 - a3*t21
1210 : t30 = t10*t8
1211 : t31 = t3**2
1212 : t33 = o + t30 + a5*t31
1213 : t38 = t33**2
1214 : fs(ip, 2) = &
1215 : (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
1216 : t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
1217 : END DO
1218 : !$OMP END DO
1219 : END IF
1220 : IF (m >= 2) THEN
1221 : !$OMP DO
1222 : DO ip = 1, SIZE(s)
1223 : x = s(ip)
1224 : t1 = bb**2
1225 : t2 = x**2
1226 : t5 = SQRT(o + t1*t2)
1227 : t7 = o/t5*t1
1228 : t9 = bb + t7*x
1229 : t12 = bb*x + t5
1230 : t13 = o/t12
1231 : t15 = 2._dp*a1*t9*t13
1232 : t16 = a1*x
1233 : t17 = t5**2
1234 : t20 = t1**2
1235 : t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
1236 : t26 = t9**2
1237 : t27 = t12**2
1238 : t30 = t16*t26/t27
1239 : t31 = a3*a4
1240 : t33 = EXP(-a4*t2)
1241 : t37 = a4**2
1242 : t39 = t2**2
1243 : t44 = a3*t33
1244 : t47 = LOG(t12)
1245 : t48 = t16*t47
1246 : t50 = o + t48 + a5*t39
1247 : t53 = a1*t47
1248 : t55 = t16*t9*t13
1249 : t56 = t2*x
1250 : t60 = a2 - t44
1251 : t64 = t50**2
1252 : t65 = o/t64
1253 : t69 = t53 + t55 + 4._dp*a5*t56
1254 : t73 = o + t48 + t60*t2
1255 : t77 = t69**2
1256 : fs(ip, 3) = &
1257 : (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
1258 : 2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
1259 : (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
1260 : t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
1261 : END DO
1262 : !$OMP END DO
1263 : END IF
1264 : IF (m >= 3) THEN
1265 : !$OMP DO
1266 : DO ip = 1, SIZE(s)
1267 : x = s(ip)
1268 : t1 = bb**2
1269 : t2 = x**2
1270 : t5 = SQRT(0.1e1_dp + t1*t2)
1271 : t6 = t5**2
1272 : t9 = t1**2
1273 : t10 = 1/t6/t5*t9
1274 : t13 = 1/t5*t1
1275 : t14 = -t10*t2 + t13
1276 : t17 = bb*x + t5
1277 : t18 = 1/t17
1278 : t20 = 3*a1*t14*t18
1279 : t22 = bb + t13*x
1280 : t23 = t22**2
1281 : t25 = t17**2
1282 : t26 = 1/t25
1283 : t28 = 3*a1*t23*t26
1284 : t29 = a1*x
1285 : t30 = t6**2
1286 : t35 = t2*x
1287 : t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
1288 : t44 = 3*t29*t14*t26*t22
1289 : t50 = 2*t29*t23*t22/t25/t17
1290 : t51 = a3*a4
1291 : t53 = EXP(-a4*t2)
1292 : t57 = a4**2
1293 : t58 = a3*t57
1294 : t59 = t35*t53
1295 : t64 = t2**2
1296 : t70 = LOG(t17)
1297 : t71 = t29*t70
1298 : t73 = 0.1e1_dp + t71 + a5*t64
1299 : t78 = 2*a1*t22*t18
1300 : t80 = t29*t14*t18
1301 : t82 = t29*t23*t26
1302 : t90 = a3*t53
1303 : t93 = t73**2
1304 : t94 = 1/t93
1305 : t96 = a1*t70
1306 : t98 = t29*t18*t22
1307 : t101 = t96 + t98 + 4*a5*t35
1308 : t106 = a2 - t90
1309 : t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1310 : t111 = 1/t93/t73
1311 : t113 = t101**2
1312 : t119 = t78 + t80 - t82 + 12*a5*t2
1313 : t123 = 0.1e1_dp + t71 + t106*t2
1314 : t124 = t93**2
1315 : fs(ip, 4) = &
1316 : (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
1317 : x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
1318 : 4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
1319 : 6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
1320 : 6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
1321 : END DO
1322 : !$OMP END DO
1323 : END IF
1324 :
1325 : !$OMP END PARALLEL
1326 :
1327 108 : END SUBROUTINE efactor_pw91
1328 :
1329 : END MODULE xc_ke_gga
1330 :
|
