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 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
159 0 : allocate_deriv=.TRUE.)
160 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
161 :
162 0 : CALL tfw_u_0(rho, r13, s, e_0, npoints)
163 0 : IF (order >= 1) THEN
164 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
165 0 : allocate_deriv=.TRUE.)
166 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
167 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
168 0 : allocate_deriv=.TRUE.)
169 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
170 :
171 0 : CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
172 : END IF
173 0 : IF (order >= 2) THEN
174 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
175 0 : allocate_deriv=.TRUE.)
176 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
177 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
178 0 : allocate_deriv=.TRUE.)
179 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
180 : deriv => xc_dset_get_derivative(deriv_set, &
181 0 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
182 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
183 :
184 : CALL tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, &
185 0 : e_ndrho_ndrho, npoints)
186 : END IF
187 0 : IF (order >= 3) THEN
188 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
189 0 : allocate_deriv=.TRUE.)
190 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
191 : deriv => xc_dset_get_derivative(deriv_set, &
192 0 : [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
193 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
194 : deriv => xc_dset_get_derivative(deriv_set, &
195 0 : [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
196 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
197 :
198 : CALL tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
199 0 : e_rho_ndrho_ndrho, npoints)
200 : END IF
201 0 : IF (order > 3) THEN
202 0 : CPABORT("derivatives bigger than 3 not implemented")
203 : END IF
204 :
205 0 : DEALLOCATE (s)
206 0 : CALL timestop(handle)
207 0 : END SUBROUTINE tfw_lda_eval
208 :
209 : ! **************************************************************************************************
210 : !> \brief ...
211 : !> \param rho ...
212 : !> \param grho ...
213 : !> \param s ...
214 : !> \param npoints ...
215 : ! **************************************************************************************************
216 0 : SUBROUTINE calc_s(rho, grho, s, npoints)
217 : REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rho, grho
218 : REAL(KIND=dp), DIMENSION(*), INTENT(out) :: s
219 : INTEGER, INTENT(in) :: npoints
220 :
221 : INTEGER :: ip
222 :
223 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
224 0 : !$OMP SHARED(npoints,rho,eps_rho,s,grho)
225 : DO ip = 1, npoints
226 : IF (rho(ip) < eps_rho) THEN
227 : s(ip) = 0.0_dp
228 : ELSE
229 : s(ip) = grho(ip)*grho(ip)/rho(ip)
230 : END IF
231 : END DO
232 0 : END SUBROUTINE calc_s
233 :
234 : ! **************************************************************************************************
235 : !> \brief ...
236 : !> \param rho_set ...
237 : !> \param deriv_set ...
238 : !> \param order ...
239 : ! **************************************************************************************************
240 0 : SUBROUTINE tfw_lsd_eval(rho_set, deriv_set, order)
241 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
242 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
243 : INTEGER, INTENT(in) :: order
244 :
245 : CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lsd_eval'
246 : INTEGER, DIMENSION(2), PARAMETER :: &
247 : norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
248 : rho_spin_name = [deriv_rhoa, deriv_rhob]
249 :
250 : INTEGER :: handle, i, ispin, npoints
251 : INTEGER, DIMENSION(2, 3) :: bo
252 : REAL(KIND=dp) :: epsilon_rho
253 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s
254 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
255 0 : POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, &
256 0 : e_rho_ndrho, e_rho_ndrho_ndrho, &
257 0 : e_rho_rho, e_rho_rho_ndrho, &
258 0 : e_rho_rho_rho
259 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
260 : TYPE(xc_derivative_type), POINTER :: deriv
261 :
262 0 : CALL timeset(routineN, handle)
263 0 : NULLIFY (deriv)
264 0 : DO i = 1, 2
265 0 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
266 : END DO
267 :
268 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
269 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
270 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
271 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
272 0 : local_bounds=bo)
273 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
274 0 : CALL tfw_init(epsilon_rho)
275 :
276 0 : ALLOCATE (s(npoints))
277 :
278 0 : DO ispin = 1, 2
279 0 : CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints)
280 :
281 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
282 0 : allocate_deriv=.TRUE.)
283 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
284 :
285 : CALL tfw_p_0(rho(ispin)%array, &
286 0 : rho_1_3(ispin)%array, s, e_0, npoints)
287 0 : IF (order >= 1) THEN
288 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
289 0 : allocate_deriv=.TRUE.)
290 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
291 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
292 0 : allocate_deriv=.TRUE.)
293 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
294 :
295 : CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, &
296 0 : rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints)
297 : END IF
298 0 : IF (order >= 2) THEN
299 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
300 0 : rho_spin_name(ispin)], allocate_deriv=.TRUE.)
301 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
302 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
303 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
304 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
305 : deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
306 0 : norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
307 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
308 :
309 : CALL tfw_p_2(rho(ispin)%array, norm_drho(ispin)%array, &
310 : rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho, &
311 0 : e_ndrho_ndrho, npoints)
312 : END IF
313 0 : IF (order >= 3) THEN
314 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
315 : rho_spin_name(ispin), rho_spin_name(ispin)], &
316 0 : allocate_deriv=.TRUE.)
317 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
318 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
319 : rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
320 0 : allocate_deriv=.TRUE.)
321 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
322 : deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
323 : norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
324 0 : allocate_deriv=.TRUE.)
325 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
326 :
327 : CALL tfw_p_3(rho(ispin)%array, norm_drho(ispin)%array, &
328 : rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho, &
329 0 : e_rho_ndrho_ndrho, npoints)
330 : END IF
331 0 : IF (order > 3) THEN
332 0 : CPABORT("derivatives bigger than 3 not implemented")
333 : END IF
334 : END DO
335 :
336 0 : DEALLOCATE (s)
337 0 : CALL timestop(handle)
338 0 : END SUBROUTINE tfw_lsd_eval
339 :
340 : ! **************************************************************************************************
341 : !> \brief ...
342 : !> \param rho ...
343 : !> \param r13 ...
344 : !> \param s ...
345 : !> \param e_0 ...
346 : !> \param npoints ...
347 : ! **************************************************************************************************
348 0 : SUBROUTINE tfw_u_0(rho, r13, s, e_0, npoints)
349 :
350 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
351 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
352 : INTEGER, INTENT(in) :: npoints
353 :
354 : INTEGER :: ip
355 :
356 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
357 0 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13,s,fvw)
358 : DO ip = 1, npoints
359 :
360 : IF (rho(ip) > eps_rho) THEN
361 :
362 : e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw*s(ip)
363 :
364 : END IF
365 :
366 : END DO
367 :
368 0 : END SUBROUTINE tfw_u_0
369 :
370 : ! **************************************************************************************************
371 : !> \brief ...
372 : !> \param rho ...
373 : !> \param grho ...
374 : !> \param r13 ...
375 : !> \param s ...
376 : !> \param e_rho ...
377 : !> \param e_ndrho ...
378 : !> \param npoints ...
379 : ! **************************************************************************************************
380 0 : SUBROUTINE tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
381 :
382 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
383 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
384 : INTEGER, INTENT(in) :: npoints
385 :
386 : INTEGER :: ip
387 : REAL(KIND=dp) :: f
388 :
389 0 : f = f53*flda
390 :
391 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
392 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,e_ndrho,grho,s,r13,f,fvw)
393 : DO ip = 1, npoints
394 :
395 : IF (rho(ip) > eps_rho) THEN
396 :
397 : e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip) - fvw*s(ip)/rho(ip)
398 : e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grho(ip)/rho(ip)
399 :
400 : END IF
401 :
402 : END DO
403 :
404 0 : END SUBROUTINE tfw_u_1
405 :
406 : ! **************************************************************************************************
407 : !> \brief ...
408 : !> \param rho ...
409 : !> \param grho ...
410 : !> \param r13 ...
411 : !> \param s ...
412 : !> \param e_rho_rho ...
413 : !> \param e_rho_ndrho ...
414 : !> \param e_ndrho_ndrho ...
415 : !> \param npoints ...
416 : ! **************************************************************************************************
417 0 : SUBROUTINE tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
418 : npoints)
419 :
420 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
421 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
422 : INTEGER, INTENT(in) :: npoints
423 :
424 : INTEGER :: ip
425 : REAL(KIND=dp) :: f
426 :
427 0 : f = f23*f53*flda
428 :
429 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
430 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho,grho,f,fvw)
431 : DO ip = 1, npoints
432 :
433 : IF (rho(ip) > eps_rho) THEN
434 :
435 : e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip) + 2.0_dp*fvw*s(ip)/(rho(ip)*rho(ip))
436 : e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip))
437 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rho(ip)
438 :
439 : END IF
440 :
441 : END DO
442 :
443 0 : END SUBROUTINE tfw_u_2
444 :
445 : ! **************************************************************************************************
446 : !> \brief ...
447 : !> \param rho ...
448 : !> \param grho ...
449 : !> \param r13 ...
450 : !> \param s ...
451 : !> \param e_rho_rho_rho ...
452 : !> \param e_rho_rho_ndrho ...
453 : !> \param e_rho_ndrho_ndrho ...
454 : !> \param npoints ...
455 : ! **************************************************************************************************
456 0 : SUBROUTINE tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
457 : e_rho_ndrho_ndrho, npoints)
458 :
459 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
460 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
461 : e_rho_ndrho_ndrho
462 : INTEGER, INTENT(in) :: npoints
463 :
464 : INTEGER :: ip
465 : REAL(KIND=dp) :: f
466 :
467 0 : f = -f13*f23*f53*flda
468 :
469 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
470 0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,r13,s,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw)
471 : DO ip = 1, npoints
472 :
473 : IF (rho(ip) > eps_rho) THEN
474 :
475 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip)) &
476 : - 6.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)*rho(ip))
477 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
478 : + 4.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)*rho(ip))
479 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
480 : - 2.0_dp*fvw/(rho(ip)*rho(ip))
481 : END IF
482 :
483 : END DO
484 :
485 0 : END SUBROUTINE tfw_u_3
486 :
487 : ! **************************************************************************************************
488 : !> \brief ...
489 : !> \param rhoa ...
490 : !> \param r13a ...
491 : !> \param sa ...
492 : !> \param e_0 ...
493 : !> \param npoints ...
494 : ! **************************************************************************************************
495 0 : SUBROUTINE tfw_p_0(rhoa, r13a, sa, e_0, npoints)
496 :
497 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a, sa
498 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
499 : INTEGER, INTENT(in) :: npoints
500 :
501 : INTEGER :: ip
502 :
503 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
504 0 : !$OMP SHARED(npoints, rhoa,eps_rho,e_0,r13a,sa,flsd,fvw)
505 : DO ip = 1, npoints
506 :
507 : IF (rhoa(ip) > eps_rho) THEN
508 : e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip) + fvw*sa(ip)
509 : END IF
510 :
511 : END DO
512 :
513 0 : END SUBROUTINE tfw_p_0
514 :
515 : ! **************************************************************************************************
516 : !> \brief ...
517 : !> \param rhoa ...
518 : !> \param grhoa ...
519 : !> \param r13a ...
520 : !> \param sa ...
521 : !> \param e_rho ...
522 : !> \param e_ndrho ...
523 : !> \param npoints ...
524 : ! **************************************************************************************************
525 0 : SUBROUTINE tfw_p_1(rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints)
526 :
527 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
528 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
529 : INTEGER, INTENT(in) :: npoints
530 :
531 : INTEGER :: ip
532 : REAL(KIND=dp) :: f
533 :
534 0 : f = f53*flsd
535 :
536 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
537 0 : !$OMP SHARED(npoints,rhoa,eps_rho,r13a,sa,fvw,grhoa,e_rho,e_ndrho,f)
538 : DO ip = 1, npoints
539 :
540 : IF (rhoa(ip) > eps_rho) THEN
541 : e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip) - fvw*sa(ip)/rhoa(ip)
542 : e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grhoa(ip)/rhoa(ip)
543 : END IF
544 :
545 : END DO
546 :
547 0 : END SUBROUTINE tfw_p_1
548 :
549 : ! **************************************************************************************************
550 : !> \brief ...
551 : !> \param rhoa ...
552 : !> \param grhoa ...
553 : !> \param r13a ...
554 : !> \param sa ...
555 : !> \param e_rho_rho ...
556 : !> \param e_rho_ndrho ...
557 : !> \param e_ndrho_ndrho ...
558 : !> \param npoints ...
559 : ! **************************************************************************************************
560 0 : SUBROUTINE tfw_p_2(rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho, &
561 : e_ndrho_ndrho, npoints)
562 :
563 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
564 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
565 : INTEGER, INTENT(in) :: npoints
566 :
567 : INTEGER :: ip
568 : REAL(KIND=dp) :: f
569 :
570 0 : f = f23*f53*flsd
571 :
572 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
573 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,fvw,r13a,sa,e_rho_ndrho,e_ndrho_ndrho)
574 : DO ip = 1, npoints
575 :
576 : IF (rhoa(ip) > eps_rho) THEN
577 : e_rho_rho(ip) = e_rho_rho(ip) &
578 : + f/r13a(ip) + 2.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip))
579 : e_rho_ndrho(ip) = e_rho_ndrho(ip) &
580 : - 2.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip))
581 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rhoa(ip)
582 : END IF
583 :
584 : END DO
585 :
586 0 : END SUBROUTINE tfw_p_2
587 :
588 : ! **************************************************************************************************
589 : !> \brief ...
590 : !> \param rhoa ...
591 : !> \param grhoa ...
592 : !> \param r13a ...
593 : !> \param sa ...
594 : !> \param e_rho_rho_rho ...
595 : !> \param e_rho_rho_ndrho ...
596 : !> \param e_rho_ndrho_ndrho ...
597 : !> \param npoints ...
598 : ! **************************************************************************************************
599 0 : SUBROUTINE tfw_p_3(rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho, &
600 : e_rho_ndrho_ndrho, npoints)
601 :
602 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
603 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
604 : e_rho_ndrho_ndrho
605 : INTEGER, INTENT(in) :: npoints
606 :
607 : INTEGER :: ip
608 : REAL(KIND=dp) :: f
609 :
610 0 : f = -f13*f23*f53*flsd
611 :
612 : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
613 0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw,sa,grhoa)
614 : DO ip = 1, npoints
615 :
616 : IF (rhoa(ip) > eps_rho) THEN
617 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
618 : + f/(r13a(ip)*rhoa(ip)) &
619 : - 6.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
620 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
621 : + 4.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
622 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
623 : - 2.0_dp*fvw/(rhoa(ip)*rhoa(ip))
624 : END IF
625 :
626 : END DO
627 :
628 0 : END SUBROUTINE tfw_p_3
629 :
630 : END MODULE xc_tfw
631 :
|
