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 the Thomas-Fermi kinetic energy functional
10 : !> plus the von Weizsaecker term
11 : !> \par History
12 : !> JGH (26.02.2003) : OpenMP enabled
13 : !> fawzi (04.2004) : adapted to the new xc interface
14 : !> \author JGH (18.02.2002)
15 : ! **************************************************************************************************
16 : MODULE xc_tfw
17 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: pi
20 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
21 : deriv_norm_drhoa,&
22 : deriv_norm_drhob,&
23 : deriv_rho,&
24 : deriv_rhoa,&
25 : deriv_rhob
26 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
27 : xc_dset_get_derivative
28 : USE xc_derivative_types, ONLY: xc_derivative_get,&
29 : xc_derivative_type
30 : USE xc_functionals_utilities, ONLY: set_util
31 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
32 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
33 : xc_rho_set_type
34 : #include "../base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : ! *** Global parameters ***
41 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
42 : f23 = 2.0_dp*f13, &
43 : f43 = 4.0_dp*f13, &
44 : f53 = 5.0_dp*f13
45 :
46 : PUBLIC :: tfw_lda_info, tfw_lda_eval, tfw_lsd_info, tfw_lsd_eval
47 :
48 : REAL(KIND=dp) :: cf, flda, flsd, fvw
49 : REAL(KIND=dp) :: eps_rho
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tfw'
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief ...
56 : !> \param cutoff ...
57 : ! **************************************************************************************************
58 0 : SUBROUTINE tfw_init(cutoff)
59 :
60 : REAL(KIND=dp), INTENT(IN) :: cutoff
61 :
62 0 : eps_rho = cutoff
63 0 : CALL set_util(cutoff)
64 :
65 0 : cf = 0.3_dp*(3.0_dp*pi*pi)**f23
66 0 : flda = cf
67 0 : flsd = flda*2.0_dp**f23
68 0 : fvw = 1.0_dp/72.0_dp
69 :
70 0 : END SUBROUTINE tfw_init
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param reference ...
75 : !> \param shortform ...
76 : !> \param needs ...
77 : !> \param max_deriv ...
78 : ! **************************************************************************************************
79 0 : SUBROUTINE tfw_lda_info(reference, shortform, needs, max_deriv)
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 0 : IF (PRESENT(reference)) THEN
85 0 : reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}"
86 : END IF
87 0 : IF (PRESENT(shortform)) THEN
88 0 : shortform = "TF+vW kinetic energy functional {LDA}"
89 : END IF
90 0 : IF (PRESENT(needs)) THEN
91 0 : needs%rho = .TRUE.
92 0 : needs%rho_1_3 = .TRUE.
93 0 : needs%norm_drho = .TRUE.
94 : END IF
95 0 : IF (PRESENT(max_deriv)) max_deriv = 3
96 :
97 0 : END SUBROUTINE tfw_lda_info
98 :
99 : ! **************************************************************************************************
100 : !> \brief ...
101 : !> \param reference ...
102 : !> \param shortform ...
103 : !> \param needs ...
104 : !> \param max_deriv ...
105 : ! **************************************************************************************************
106 0 : SUBROUTINE tfw_lsd_info(reference, shortform, needs, max_deriv)
107 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
108 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
109 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
110 :
111 0 : IF (PRESENT(reference)) THEN
112 0 : reference = "Thomas-Fermi-Weizsaecker kinetic energy functional"
113 : END IF
114 0 : IF (PRESENT(shortform)) THEN
115 0 : shortform = "TF+vW kinetic energy functional"
116 : END IF
117 0 : IF (PRESENT(needs)) THEN
118 0 : needs%rho_spin = .TRUE.
119 0 : needs%rho_spin_1_3 = .TRUE.
120 0 : needs%norm_drho = .TRUE.
121 : END IF
122 0 : IF (PRESENT(max_deriv)) max_deriv = 3
123 :
124 0 : END SUBROUTINE tfw_lsd_info
125 :
126 : ! **************************************************************************************************
127 : !> \brief ...
128 : !> \param rho_set ...
129 : !> \param deriv_set ...
130 : !> \param order ...
131 : ! **************************************************************************************************
132 0 : SUBROUTINE tfw_lda_eval(rho_set, deriv_set, order)
133 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
134 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
135 : INTEGER, INTENT(in) :: order
136 :
137 : CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lda_eval'
138 :
139 : INTEGER :: handle, npoints
140 : INTEGER, DIMENSION(2, 3) :: bo
141 : REAL(KIND=dp) :: epsilon_rho
142 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
143 0 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
144 0 : e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, e_rho_rho_rho, grho, &
145 0 : r13, rho
146 : TYPE(xc_derivative_type), POINTER :: deriv
147 :
148 0 : CALL timeset(routineN, handle)
149 :
150 : CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
151 0 : norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
152 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
153 0 : CALL tfw_init(epsilon_rho)
154 :
155 0 : ALLOCATE (s(npoints))
156 0 : CALL calc_s(rho, grho, s, npoints)
157 :
158 0 : IF (order >= 0) THEN
159 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
160 0 : allocate_deriv=.TRUE.)
161 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
162 :
163 0 : CALL tfw_u_0(rho, r13, s, e_0, npoints)
164 : END IF
165 0 : IF (order >= 1 .OR. order == -1) THEN
166 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
167 0 : allocate_deriv=.TRUE.)
168 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
169 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
170 0 : allocate_deriv=.TRUE.)
171 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
172 :
173 0 : CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
174 : END IF
175 0 : IF (order >= 2 .OR. order == -2) THEN
176 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
177 0 : allocate_deriv=.TRUE.)
178 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
179 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
180 0 : allocate_deriv=.TRUE.)
181 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
182 : deriv => xc_dset_get_derivative(deriv_set, &
183 0 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
184 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
185 :
186 : CALL tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, &
187 0 : e_ndrho_ndrho, npoints)
188 : END IF
189 0 : IF (order >= 3 .OR. order == -3) THEN
190 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
191 0 : allocate_deriv=.TRUE.)
192 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
193 : deriv => xc_dset_get_derivative(deriv_set, &
194 0 : [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
195 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
196 : deriv => xc_dset_get_derivative(deriv_set, &
197 0 : [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
198 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
199 :
200 : CALL tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
201 0 : e_rho_ndrho_ndrho, npoints)
202 : END IF
203 0 : IF (order > 3 .OR. order < -3) THEN
204 0 : CPABORT("derivatives bigger than 3 not implemented")
205 : END IF
206 :
207 0 : DEALLOCATE (s)
208 0 : CALL timestop(handle)
209 0 : END SUBROUTINE tfw_lda_eval
210 :
211 : ! **************************************************************************************************
212 : !> \brief ...
213 : !> \param rho ...
214 : !> \param grho ...
215 : !> \param s ...
216 : !> \param npoints ...
217 : ! **************************************************************************************************
218 0 : SUBROUTINE calc_s(rho, grho, s, npoints)
219 : REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rho, grho
220 : REAL(KIND=dp), DIMENSION(*), INTENT(out) :: s
221 : INTEGER, INTENT(in) :: npoints
222 :
223 : INTEGER :: ip
224 :
225 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
226 0 : !$OMP SHARED(npoints,rho,eps_rho,s,grho)
227 : DO ip = 1, npoints
228 : IF (rho(ip) < eps_rho) THEN
229 : s(ip) = 0.0_dp
230 : ELSE
231 : s(ip) = grho(ip)*grho(ip)/rho(ip)
232 : END IF
233 : END DO
234 0 : END SUBROUTINE calc_s
235 :
236 : ! **************************************************************************************************
237 : !> \brief ...
238 : !> \param rho_set ...
239 : !> \param deriv_set ...
240 : !> \param order ...
241 : ! **************************************************************************************************
242 0 : SUBROUTINE tfw_lsd_eval(rho_set, deriv_set, order)
243 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
244 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
245 : INTEGER, INTENT(in) :: order
246 :
247 : CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lsd_eval'
248 : INTEGER, DIMENSION(2), PARAMETER :: &
249 : norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
250 : rho_spin_name = [deriv_rhoa, deriv_rhob]
251 :
252 : INTEGER :: handle, i, ispin, npoints
253 : INTEGER, DIMENSION(2, 3) :: bo
254 : REAL(KIND=dp) :: epsilon_rho
255 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
256 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
257 0 : POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, &
258 0 : e_rho_ndrho, e_rho_ndrho_ndrho, &
259 0 : e_rho_rho, e_rho_rho_ndrho, &
260 0 : e_rho_rho_rho
261 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
262 : TYPE(xc_derivative_type), POINTER :: deriv
263 :
264 0 : CALL timeset(routineN, handle)
265 0 : NULLIFY (deriv)
266 0 : DO i = 1, 2
267 0 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
268 : END DO
269 :
270 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
271 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
272 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
273 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
274 0 : local_bounds=bo)
275 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
276 0 : CALL tfw_init(epsilon_rho)
277 :
278 0 : ALLOCATE (s(npoints))
279 :
280 0 : DO ispin = 1, 2
281 0 : CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints)
282 :
283 0 : IF (order >= 0) THEN
284 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
285 0 : allocate_deriv=.TRUE.)
286 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
287 :
288 : CALL tfw_p_0(rho(ispin)%array, &
289 0 : rho_1_3(ispin)%array, s, e_0, npoints)
290 : END IF
291 0 : IF (order >= 1 .OR. order == -1) THEN
292 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
293 0 : allocate_deriv=.TRUE.)
294 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
295 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
296 0 : allocate_deriv=.TRUE.)
297 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
298 :
299 : CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, &
300 0 : rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints)
301 : END IF
302 0 : IF (order >= 2 .OR. order == -2) THEN
303 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
304 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
305 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
306 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
307 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
308 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
309 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
310 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
311 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
312 :
313 : CALL tfw_p_2(rho(ispin)%array, norm_drho(ispin)%array, &
314 : rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho, &
315 0 : e_ndrho_ndrho, npoints)
316 : END IF
317 0 : IF (order >= 3 .OR. order == -3) THEN
318 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
319 : rho_spin_name(ispin), rho_spin_name(ispin)], &
320 0 : allocate_deriv=.TRUE.)
321 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
322 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
323 : rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
324 0 : allocate_deriv=.TRUE.)
325 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
326 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
327 : norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
328 0 : allocate_deriv=.TRUE.)
329 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
330 :
331 : CALL tfw_p_3(rho(ispin)%array, norm_drho(ispin)%array, &
332 : rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho, &
333 0 : e_rho_ndrho_ndrho, npoints)
334 : END IF
335 0 : IF (order > 3 .OR. order < -3) THEN
336 0 : CPABORT("derivatives bigger than 3 not implemented")
337 : END IF
338 : END DO
339 :
340 0 : DEALLOCATE (s)
341 0 : CALL timestop(handle)
342 0 : END SUBROUTINE tfw_lsd_eval
343 :
344 : ! **************************************************************************************************
345 : !> \brief ...
346 : !> \param rho ...
347 : !> \param r13 ...
348 : !> \param s ...
349 : !> \param e_0 ...
350 : !> \param npoints ...
351 : ! **************************************************************************************************
352 0 : SUBROUTINE tfw_u_0(rho, r13, s, e_0, npoints)
353 :
354 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
355 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
356 : INTEGER, INTENT(in) :: npoints
357 :
358 : INTEGER :: ip
359 :
360 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
361 0 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13,s,fvw)
362 : DO ip = 1, npoints
363 :
364 : IF (rho(ip) > eps_rho) THEN
365 :
366 : e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw*s(ip)
367 :
368 : END IF
369 :
370 : END DO
371 :
372 0 : END SUBROUTINE tfw_u_0
373 :
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param rho ...
377 : !> \param grho ...
378 : !> \param r13 ...
379 : !> \param s ...
380 : !> \param e_rho ...
381 : !> \param e_ndrho ...
382 : !> \param npoints ...
383 : ! **************************************************************************************************
384 0 : SUBROUTINE tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
385 :
386 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
387 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
388 : INTEGER, INTENT(in) :: npoints
389 :
390 : INTEGER :: ip
391 : REAL(KIND=dp) :: f
392 :
393 0 : f = f53*flda
394 :
395 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
396 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,e_ndrho,grho,s,r13,f,fvw)
397 : DO ip = 1, npoints
398 :
399 : IF (rho(ip) > eps_rho) THEN
400 :
401 : e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip) - fvw*s(ip)/rho(ip)
402 : e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grho(ip)/rho(ip)
403 :
404 : END IF
405 :
406 : END DO
407 :
408 0 : END SUBROUTINE tfw_u_1
409 :
410 : ! **************************************************************************************************
411 : !> \brief ...
412 : !> \param rho ...
413 : !> \param grho ...
414 : !> \param r13 ...
415 : !> \param s ...
416 : !> \param e_rho_rho ...
417 : !> \param e_rho_ndrho ...
418 : !> \param e_ndrho_ndrho ...
419 : !> \param npoints ...
420 : ! **************************************************************************************************
421 0 : SUBROUTINE tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
422 : npoints)
423 :
424 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
425 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
426 : INTEGER, INTENT(in) :: npoints
427 :
428 : INTEGER :: ip
429 : REAL(KIND=dp) :: f
430 :
431 0 : f = f23*f53*flda
432 :
433 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
434 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho,grho,f,fvw)
435 : DO ip = 1, npoints
436 :
437 : IF (rho(ip) > eps_rho) THEN
438 :
439 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip) + 2.0_dp*fvw*s(ip)/(rho(ip)*rho(ip))
440 : e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip))
441 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rho(ip)
442 :
443 : END IF
444 :
445 : END DO
446 :
447 0 : END SUBROUTINE tfw_u_2
448 :
449 : ! **************************************************************************************************
450 : !> \brief ...
451 : !> \param rho ...
452 : !> \param grho ...
453 : !> \param r13 ...
454 : !> \param s ...
455 : !> \param e_rho_rho_rho ...
456 : !> \param e_rho_rho_ndrho ...
457 : !> \param e_rho_ndrho_ndrho ...
458 : !> \param npoints ...
459 : ! **************************************************************************************************
460 0 : SUBROUTINE tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
461 : e_rho_ndrho_ndrho, npoints)
462 :
463 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
464 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
465 : e_rho_ndrho_ndrho
466 : INTEGER, INTENT(in) :: npoints
467 :
468 : INTEGER :: ip
469 : REAL(KIND=dp) :: f
470 :
471 0 : f = -f13*f23*f53*flda
472 :
473 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
474 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,r13,s,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw)
475 : DO ip = 1, npoints
476 :
477 : IF (rho(ip) > eps_rho) THEN
478 :
479 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip)) &
480 : - 6.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)*rho(ip))
481 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
482 : + 4.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)*rho(ip))
483 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
484 : - 2.0_dp*fvw/(rho(ip)*rho(ip))
485 : END IF
486 :
487 : END DO
488 :
489 0 : END SUBROUTINE tfw_u_3
490 :
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param rhoa ...
494 : !> \param r13a ...
495 : !> \param sa ...
496 : !> \param e_0 ...
497 : !> \param npoints ...
498 : ! **************************************************************************************************
499 0 : SUBROUTINE tfw_p_0(rhoa, r13a, sa, e_0, npoints)
500 :
501 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a, sa
502 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
503 : INTEGER, INTENT(in) :: npoints
504 :
505 : INTEGER :: ip
506 :
507 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
508 0 : !$OMP SHARED(npoints, rhoa,eps_rho,e_0,r13a,sa,flsd,fvw)
509 : DO ip = 1, npoints
510 :
511 : IF (rhoa(ip) > eps_rho) THEN
512 : e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip) + fvw*sa(ip)
513 : END IF
514 :
515 : END DO
516 :
517 0 : END SUBROUTINE tfw_p_0
518 :
519 : ! **************************************************************************************************
520 : !> \brief ...
521 : !> \param rhoa ...
522 : !> \param grhoa ...
523 : !> \param r13a ...
524 : !> \param sa ...
525 : !> \param e_rho ...
526 : !> \param e_ndrho ...
527 : !> \param npoints ...
528 : ! **************************************************************************************************
529 0 : SUBROUTINE tfw_p_1(rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints)
530 :
531 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
532 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
533 : INTEGER, INTENT(in) :: npoints
534 :
535 : INTEGER :: ip
536 : REAL(KIND=dp) :: f
537 :
538 0 : f = f53*flsd
539 :
540 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
541 0 : !$OMP SHARED(npoints,rhoa,eps_rho,r13a,sa,fvw,grhoa,e_rho,e_ndrho,f)
542 : DO ip = 1, npoints
543 :
544 : IF (rhoa(ip) > eps_rho) THEN
545 : e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip) - fvw*sa(ip)/rhoa(ip)
546 : e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grhoa(ip)/rhoa(ip)
547 : END IF
548 :
549 : END DO
550 :
551 0 : END SUBROUTINE tfw_p_1
552 :
553 : ! **************************************************************************************************
554 : !> \brief ...
555 : !> \param rhoa ...
556 : !> \param grhoa ...
557 : !> \param r13a ...
558 : !> \param sa ...
559 : !> \param e_rho_rho ...
560 : !> \param e_rho_ndrho ...
561 : !> \param e_ndrho_ndrho ...
562 : !> \param npoints ...
563 : ! **************************************************************************************************
564 0 : SUBROUTINE tfw_p_2(rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho, &
565 : e_ndrho_ndrho, npoints)
566 :
567 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
568 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
569 : INTEGER, INTENT(in) :: npoints
570 :
571 : INTEGER :: ip
572 : REAL(KIND=dp) :: f
573 :
574 0 : f = f23*f53*flsd
575 :
576 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
577 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,fvw,r13a,sa,e_rho_ndrho,e_ndrho_ndrho)
578 : DO ip = 1, npoints
579 :
580 : IF (rhoa(ip) > eps_rho) THEN
581 : e_rho_rho(ip) = e_rho_rho(ip) &
582 : + f/r13a(ip) + 2.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip))
583 : e_rho_ndrho(ip) = e_rho_ndrho(ip) &
584 : - 2.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip))
585 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rhoa(ip)
586 : END IF
587 :
588 : END DO
589 :
590 0 : END SUBROUTINE tfw_p_2
591 :
592 : ! **************************************************************************************************
593 : !> \brief ...
594 : !> \param rhoa ...
595 : !> \param grhoa ...
596 : !> \param r13a ...
597 : !> \param sa ...
598 : !> \param e_rho_rho_rho ...
599 : !> \param e_rho_rho_ndrho ...
600 : !> \param e_rho_ndrho_ndrho ...
601 : !> \param npoints ...
602 : ! **************************************************************************************************
603 0 : SUBROUTINE tfw_p_3(rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho, &
604 : e_rho_ndrho_ndrho, npoints)
605 :
606 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
607 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
608 : e_rho_ndrho_ndrho
609 : INTEGER, INTENT(in) :: npoints
610 :
611 : INTEGER :: ip
612 : REAL(KIND=dp) :: f
613 :
614 0 : f = -f13*f23*f53*flsd
615 :
616 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
617 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw,sa,grhoa)
618 : DO ip = 1, npoints
619 :
620 : IF (rhoa(ip) > eps_rho) THEN
621 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
622 : + f/(r13a(ip)*rhoa(ip)) &
623 : - 6.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
624 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
625 : + 4.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
626 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
627 : - 2.0_dp*fvw/(rhoa(ip)*rhoa(ip))
628 : END IF
629 :
630 : END DO
631 :
632 0 : END SUBROUTINE tfw_p_3
633 :
634 : END MODULE xc_tfw
635 :
|
