33#include "../base/base_uses.f90"
40 REAL(KIND=
dp),
PARAMETER :: f13 = 1.0_dp/3.0_dp, &
47 REAL(KIND=
dp) :: cf, flda, flsd
48 REAL(KIND=
dp) :: eps_rho
50 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_thomas_fermi'
58 SUBROUTINE thomas_fermi_init(cutoff)
60 REAL(KIND=
dp),
INTENT(IN) :: cutoff
65 cf = 0.3_dp*(3.0_dp*
pi*
pi)**f23
67 flsd = flda*2.0_dp**f23
69 END SUBROUTINE thomas_fermi_init
80 LOGICAL,
INTENT(in) :: lsd
81 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
83 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
85 IF (
PRESENT(reference))
THEN
86 reference =
"Thomas-Fermi kinetic energy functional: see Parr and Yang"
88 IF (len_trim(reference) + 6 < len(reference))
THEN
89 reference(len_trim(reference):len_trim(reference) + 6) =
' {LDA}'
93 IF (
PRESENT(shortform))
THEN
94 shortform =
"Thomas-Fermi kinetic energy functional"
96 IF (len_trim(shortform) + 6 < len(shortform))
THEN
97 shortform(len_trim(shortform):len_trim(shortform) + 6) =
' {LDA}'
101 IF (
PRESENT(needs))
THEN
103 needs%rho_spin = .true.
104 needs%rho_spin_1_3 = .true.
107 needs%rho_1_3 = .true.
110 IF (
PRESENT(max_deriv)) max_deriv = 3
123 INTEGER,
INTENT(in) :: order
125 CHARACTER(len=*),
PARAMETER :: routinen =
'thomas_fermi_lda_eval'
127 INTEGER :: handle, npoints
128 INTEGER,
DIMENSION(2, 3) :: bo
129 REAL(kind=
dp) :: epsilon_rho
130 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
131 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
135 CALL timeset(routinen, handle)
138 local_bounds=bo, rho_cutoff=epsilon_rho)
139 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
140 CALL thomas_fermi_init(epsilon_rho)
144 allocate_deriv=.true.)
147 CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
149 IF (order >= 1 .OR. order == -1)
THEN
151 allocate_deriv=.true.)
154 CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
156 IF (order >= 2 .OR. order == -2)
THEN
158 allocate_deriv=.true.)
161 CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
163 IF (order >= 3 .OR. order == -3)
THEN
165 allocate_deriv=.true.)
168 CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
170 IF (order > 3 .OR. order < -3)
THEN
171 cpabort(
"derivatives bigger than 3 not implemented")
173 CALL timestop(handle)
185 INTEGER,
INTENT(in) :: order
187 CHARACTER(len=*),
PARAMETER :: routinen =
'thomas_fermi_lsd_eval'
190 INTEGER :: handle, i, ispin, npoints
191 INTEGER,
DIMENSION(2, 3) :: bo
192 REAL(kind=
dp) :: epsilon_rho
193 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
194 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
198 CALL timeset(routinen, handle)
201 NULLIFY (rho(i)%array, rho_1_3(i)%array)
205 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
207 rho_cutoff=epsilon_rho, &
209 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
210 CALL thomas_fermi_init(epsilon_rho)
215 allocate_deriv=.true.)
218 CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
221 IF (order >= 1 .OR. order == -1)
THEN
223 allocate_deriv=.true.)
226 CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
229 IF (order >= 2 .OR. order == -2)
THEN
231 rho_spin_name(ispin)], allocate_deriv=.true.)
234 CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
237 IF (order >= 3 .OR. order == -3)
THEN
239 rho_spin_name(ispin), rho_spin_name(ispin)], &
240 allocate_deriv=.true.)
243 CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
244 e_rho_rho_rho, npoints)
246 IF (order > 3 .OR. order < -3)
THEN
247 cpabort(
"derivatives bigger than 3 not implemented")
250 CALL timestop(handle)
260 SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
262 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
263 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_0
264 INTEGER,
INTENT(in) :: npoints
272 IF (rho(ip) > eps_rho)
THEN
274 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
280 END SUBROUTINE thomas_fermi_lda_0
289 SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
291 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
292 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho
293 INTEGER,
INTENT(in) :: npoints
304 IF (rho(ip) > eps_rho)
THEN
306 e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
312 END SUBROUTINE thomas_fermi_lda_1
321 SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
323 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
324 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho
325 INTEGER,
INTENT(in) :: npoints
336 IF (rho(ip) > eps_rho)
THEN
338 e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
344 END SUBROUTINE thomas_fermi_lda_2
353 SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
355 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
356 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho_rho
357 INTEGER,
INTENT(in) :: npoints
362 f = -f13*f23*f53*flda
368 IF (rho(ip) > eps_rho)
THEN
370 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
376 END SUBROUTINE thomas_fermi_lda_3
385 SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
387 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
388 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_0
389 INTEGER,
INTENT(in) :: npoints
397 IF (rhoa(ip) > eps_rho)
THEN
398 e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
403 END SUBROUTINE thomas_fermi_lsd_0
412 SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
414 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
415 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho
416 INTEGER,
INTENT(in) :: npoints
427 IF (rhoa(ip) > eps_rho)
THEN
428 e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
433 END SUBROUTINE thomas_fermi_lsd_1
442 SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
444 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
445 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho
446 INTEGER,
INTENT(in) :: npoints
458 IF (rhoa(ip) > eps_rho)
THEN
459 e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
464 END SUBROUTINE thomas_fermi_lsd_2
473 SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
475 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
476 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho_rho
477 INTEGER,
INTENT(in) :: npoints
482 f = -f13*f23*f53*flsd
488 IF (rhoa(ip) > eps_rho)
THEN
489 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
494 END SUBROUTINE thomas_fermi_lsd_3
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Utility routines for the functional calculations.
subroutine, public set_util(cutoff)
...
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
Calculate the Thomas-Fermi kinetic energy functional.
subroutine, public thomas_fermi_lsd_eval(rho_set, deriv_set, order)
...
subroutine, public thomas_fermi_lda_eval(rho_set, deriv_set, order)
...
subroutine, public thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
...
represent a pointer to a contiguous 3d array
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation