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