(git:e7e05ae)
xc_tfw.F
Go to the documentation of this file.
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  SUBROUTINE tfw_init(cutoff)
60 
61  REAL(KIND=dp), INTENT(IN) :: cutoff
62 
63  eps_rho = cutoff
64  CALL set_util(cutoff)
65 
66  cf = 0.3_dp*(3.0_dp*pi*pi)**f23
67  flda = cf
68  flsd = flda*2.0_dp**f23
69  fvw = 1.0_dp/72.0_dp
70 
71  END SUBROUTINE tfw_init
72 
73 ! **************************************************************************************************
74 !> \brief ...
75 !> \param reference ...
76 !> \param shortform ...
77 !> \param needs ...
78 !> \param max_deriv ...
79 ! **************************************************************************************************
80  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  IF (PRESENT(reference)) THEN
86  reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}"
87  END IF
88  IF (PRESENT(shortform)) THEN
89  shortform = "TF+vW kinetic energy functional {LDA}"
90  END IF
91  IF (PRESENT(needs)) THEN
92  needs%rho = .true.
93  needs%rho_1_3 = .true.
94  needs%norm_drho = .true.
95  END IF
96  IF (PRESENT(max_deriv)) max_deriv = 3
97 
98  END SUBROUTINE tfw_lda_info
99 
100 ! **************************************************************************************************
101 !> \brief ...
102 !> \param reference ...
103 !> \param shortform ...
104 !> \param needs ...
105 !> \param max_deriv ...
106 ! **************************************************************************************************
107  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  IF (PRESENT(reference)) THEN
113  reference = "Thomas-Fermi-Weizsaecker kinetic energy functional"
114  END IF
115  IF (PRESENT(shortform)) THEN
116  shortform = "TF+vW kinetic energy functional"
117  END IF
118  IF (PRESENT(needs)) THEN
119  needs%rho_spin = .true.
120  needs%rho_spin_1_3 = .true.
121  needs%norm_drho = .true.
122  END IF
123  IF (PRESENT(max_deriv)) max_deriv = 3
124 
125  END SUBROUTINE tfw_lsd_info
126 
127 ! **************************************************************************************************
128 !> \brief ...
129 !> \param rho_set ...
130 !> \param deriv_set ...
131 !> \param order ...
132 ! **************************************************************************************************
133  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
144  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
145  e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, e_rho_rho_rho, grho, &
146  r13, rho
147  TYPE(xc_derivative_type), POINTER :: deriv
148 
149  CALL timeset(routinen, handle)
150 
151  CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
152  norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
153  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
154  CALL tfw_init(epsilon_rho)
155 
156  ALLOCATE (s(npoints))
157  CALL calc_s(rho, grho, s, npoints)
158 
159  IF (order >= 0) THEN
160  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
161  allocate_deriv=.true.)
162  CALL xc_derivative_get(deriv, deriv_data=e_0)
163 
164  CALL tfw_u_0(rho, r13, s, e_0, npoints)
165  END IF
166  IF (order >= 1 .OR. order == -1) THEN
167  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
168  allocate_deriv=.true.)
169  CALL xc_derivative_get(deriv, deriv_data=e_rho)
170  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
171  allocate_deriv=.true.)
172  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
173 
174  CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
175  END IF
176  IF (order >= 2 .OR. order == -2) THEN
177  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
178  allocate_deriv=.true.)
179  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
180  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
181  allocate_deriv=.true.)
182  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
183  deriv => xc_dset_get_derivative(deriv_set, &
184  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
185  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  e_ndrho_ndrho, npoints)
189  END IF
190  IF (order >= 3 .OR. order == -3) THEN
191  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
192  allocate_deriv=.true.)
193  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
194  deriv => xc_dset_get_derivative(deriv_set, &
195  [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.true.)
196  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
197  deriv => xc_dset_get_derivative(deriv_set, &
198  [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
199  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  e_rho_ndrho_ndrho, npoints)
203  END IF
204  IF (order > 3 .OR. order < -3) THEN
205  cpabort("derivatives bigger than 3 not implemented")
206  END IF
207 
208  DEALLOCATE (s)
209  CALL timestop(handle)
210  END SUBROUTINE tfw_lda_eval
211 
212 ! **************************************************************************************************
213 !> \brief ...
214 !> \param rho ...
215 !> \param grho ...
216 !> \param s ...
217 !> \param npoints ...
218 ! **************************************************************************************************
219  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 !$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  END SUBROUTINE calc_s
236 
237 ! **************************************************************************************************
238 !> \brief ...
239 !> \param rho_set ...
240 !> \param deriv_set ...
241 !> \param order ...
242 ! **************************************************************************************************
243  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
257  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
258  POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, &
259  e_rho_ndrho, e_rho_ndrho_ndrho, &
260  e_rho_rho, e_rho_rho_ndrho, &
261  e_rho_rho_rho
262  TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
263  TYPE(xc_derivative_type), POINTER :: deriv
264 
265  CALL timeset(routinen, handle)
266  NULLIFY (deriv)
267  DO i = 1, 2
268  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  local_bounds=bo)
276  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
277  CALL tfw_init(epsilon_rho)
278 
279  ALLOCATE (s(npoints))
280 
281  DO ispin = 1, 2
282  CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints)
283 
284  IF (order >= 0) THEN
285  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
286  allocate_deriv=.true.)
287  CALL xc_derivative_get(deriv, deriv_data=e_0)
288 
289  CALL tfw_p_0(rho(ispin)%array, &
290  rho_1_3(ispin)%array, s, e_0, npoints)
291  END IF
292  IF (order >= 1 .OR. order == -1) THEN
293  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
294  allocate_deriv=.true.)
295  CALL xc_derivative_get(deriv, deriv_data=e_rho)
296  deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
297  allocate_deriv=.true.)
298  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
299 
300  CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, &
301  rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints)
302  END IF
303  IF (order >= 2 .OR. order == -2) THEN
304  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
305  rho_spin_name(ispin)], allocate_deriv=.true.)
306  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
307  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
308  norm_drho_spin_name(ispin)], allocate_deriv=.true.)
309  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
310  deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
311  norm_drho_spin_name(ispin)], allocate_deriv=.true.)
312  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  e_ndrho_ndrho, npoints)
317  END IF
318  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  allocate_deriv=.true.)
322  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  allocate_deriv=.true.)
326  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  allocate_deriv=.true.)
330  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  e_rho_ndrho_ndrho, npoints)
335  END IF
336  IF (order > 3 .OR. order < -3) THEN
337  cpabort("derivatives bigger than 3 not implemented")
338  END IF
339  END DO
340 
341  DEALLOCATE (s)
342  CALL timestop(handle)
343  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  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 !$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  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  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  f = f53*flda
395 
396 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
397 !$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  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  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  f = f23*f53*flda
433 
434 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
435 !$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  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  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  f = -f13*f23*f53*flda
473 
474 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
475 !$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  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  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 !$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  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  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  f = f53*flsd
540 
541 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
542 !$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  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  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  f = f23*f53*flsd
576 
577 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
578 !$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  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  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  f = -f13*f23*f53*flsd
616 
617 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
618 !$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  END SUBROUTINE tfw_p_3
634 
635 END MODULE xc_tfw
636 
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition: xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition: xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition: xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition: xc_derivative_set_types.F:66
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition: xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition: xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition: xc_functionals_utilities.F:15
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition: xc_functionals_utilities.F:54
xc_rho_cflags_types
contains the structure
Definition: xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition: xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition: xc_rho_set_types.F:281
xc_tfw
Calculate the Thomas-Fermi kinetic energy functional plus the von Weizsaecker term.
Definition: xc_tfw.F:16
xc_tfw::tfw_lda_info
subroutine, public tfw_lda_info(reference, shortform, needs, max_deriv)
...
Definition: xc_tfw.F:81
xc_tfw::tfw_lda_eval
subroutine, public tfw_lda_eval(rho_set, deriv_set, order)
...
Definition: xc_tfw.F:134
xc_tfw::tfw_lsd_info
subroutine, public tfw_lsd_info(reference, shortform, needs, max_deriv)
...
Definition: xc_tfw.F:108
xc_tfw::tfw_lsd_eval
subroutine, public tfw_lsd_eval(rho_set, deriv_set, order)
...
Definition: xc_tfw.F:244