36 #include "../base/base_uses.f90"
44 REAL(KIND=
dp),
PARAMETER :: pi = 3.14159265358979323846264338_dp
45 REAL(KIND=
dp),
PARAMETER :: f13 = 1.0_dp/3.0_dp, &
51 REAL(KIND=
dp) :: xparam, flda, flsd
52 REAL(KIND=
dp) :: eps_rho
53 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_xalpha'
62 SUBROUTINE xalpha_init(cutoff, xalpha)
64 REAL(KIND=
dp),
INTENT(IN) :: cutoff
65 REAL(KIND=
dp),
INTENT(IN),
OPTIONAL :: xalpha
69 IF (
PRESENT(xalpha))
THEN
72 xparam = 2.0_dp/3.0_dp
75 flda = -9.0_dp/8.0_dp*xparam*(3.0_dp/pi)**f13
76 flsd = flda*2.0_dp**f13
78 END SUBROUTINE xalpha_init
90 SUBROUTINE xalpha_info(lsd, reference, shortform, needs, max_deriv, &
91 xa_parameter, scaling)
92 LOGICAL,
INTENT(in) :: lsd
93 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
94 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
95 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
96 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: xa_parameter, scaling
98 REAL(kind=
dp) :: my_scaling, my_xparam
100 my_xparam = 2.0_dp/3.0_dp
101 IF (
PRESENT(xa_parameter)) my_xparam = xa_parameter
103 IF (
PRESENT(scaling)) my_scaling = scaling
105 IF (
PRESENT(reference))
THEN
106 IF (my_scaling /= 1._dp)
THEN
107 WRITE (reference,
'(A,F8.4,A,F8.4)') &
108 "Dirac/Slater local exchange; parameter=", my_xparam,
" scaling=", my_scaling
110 WRITE (reference,
'(A,F8.4)') &
111 "Dirac/Slater local exchange; parameter=", my_xparam
114 IF (len_trim(reference) + 6 < len(reference))
THEN
115 reference(len_trim(reference):len_trim(reference) + 6) =
' {LDA}'
119 IF (
PRESENT(shortform))
THEN
120 IF (my_scaling /= 1._dp)
THEN
121 WRITE (shortform,
'(A,F8.4,F8.4)')
"Dirac/Slater exchange", my_xparam, my_scaling
123 WRITE (shortform,
'(A,F8.4)')
"Dirac/Slater exchange", my_xparam
126 IF (len_trim(shortform) + 6 < len(shortform))
THEN
127 shortform(len_trim(shortform):len_trim(shortform) + 6) =
' {LDA}'
131 IF (
PRESENT(needs))
THEN
133 needs%rho_spin = .true.
134 needs%rho_spin_1_3 = .true.
137 needs%rho_1_3 = .true.
140 IF (
PRESENT(max_deriv)) max_deriv = 3
153 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
154 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
155 INTEGER,
INTENT(in) :: order
156 TYPE(section_vals_type),
POINTER :: xa_params
157 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: xa_parameter
159 CHARACTER(len=*),
PARAMETER :: routinen =
'xalpha_lda_eval'
161 INTEGER :: handle, npoints
162 INTEGER,
DIMENSION(2, 3) :: bo
163 REAL(kind=
dp) :: epsilon_rho, sx
164 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
165 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
167 TYPE(xc_derivative_type),
POINTER :: deriv
169 CALL timeset(routinen, handle)
174 local_bounds=bo, rho_cutoff=epsilon_rho)
175 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
176 CALL xalpha_init(epsilon_rho, xa_parameter)
180 allocate_deriv=.true.)
183 CALL xalpha_lda_0(npoints, rho, r13, e_0, sx)
186 IF (order >= 1 .OR. order == -1)
THEN
188 allocate_deriv=.true.)
191 CALL xalpha_lda_1(npoints, rho, r13, e_rho, sx)
193 IF (order >= 2 .OR. order == -2)
THEN
195 allocate_deriv=.true.)
198 CALL xalpha_lda_2(npoints, rho, r13, e_rho_rho, sx)
200 IF (order >= 3 .OR. order == -3)
THEN
202 allocate_deriv=.true.)
205 CALL xalpha_lda_3(npoints, rho, r13, e_rho_rho_rho, sx)
207 IF (order > 3 .OR. order < -3)
THEN
208 cpabort(
"derivatives bigger than 3 not implemented")
210 CALL timestop(handle)
223 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
224 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
225 INTEGER,
INTENT(in) :: order
226 TYPE(section_vals_type),
POINTER :: xa_params
227 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: xa_parameter
229 CHARACTER(len=*),
PARAMETER :: routinen =
'xalpha_lsd_eval'
232 INTEGER :: handle, i, ispin, npoints
233 INTEGER,
DIMENSION(2, 3) :: bo
234 REAL(kind=
dp) :: epsilon_rho, sx
235 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
236 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
237 TYPE(cp_3d_r_cp_type),
DIMENSION(2) :: rho, rho_1_3
238 TYPE(xc_derivative_type),
POINTER :: deriv
240 CALL timeset(routinen, handle)
243 NULLIFY (rho(i)%array, rho_1_3(i)%array)
249 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
250 rhob=rho(2)%array, rho_cutoff=epsilon_rho, &
252 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
253 CALL xalpha_init(epsilon_rho, xa_parameter)
258 allocate_deriv=.true.)
261 CALL xalpha_lsd_0(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
264 IF (order >= 1 .OR. order == -1)
THEN
266 allocate_deriv=.true.)
269 CALL xalpha_lsd_1(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
272 IF (order >= 2 .OR. order == -2)
THEN
274 rho_spin_name(ispin)], allocate_deriv=.true.)
277 CALL xalpha_lsd_2(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
280 IF (order >= 3 .OR. order == -3)
THEN
282 rho_spin_name(ispin), rho_spin_name(ispin)], &
283 allocate_deriv=.true.)
286 CALL xalpha_lsd_3(npoints, rho(ispin)%array, rho_1_3(ispin)%array, &
289 IF (order > 3 .OR. order < -3)
THEN
290 cpabort(
"derivatives bigger than 3 not implemented")
293 CALL timestop(handle)
304 SUBROUTINE xalpha_lda_0(n, rho, r13, pot, sx)
306 INTEGER,
INTENT(IN) :: n
307 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
308 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pot
309 REAL(kind=
dp),
INTENT(IN) :: sx
319 IF (rho(ip) > eps_rho)
THEN
320 pot(ip) = pot(ip) + f*r13(ip)*rho(ip)
324 END SUBROUTINE xalpha_lda_0
334 SUBROUTINE xalpha_lda_1(n, rho, r13, pot, sx)
336 INTEGER,
INTENT(IN) :: n
337 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
338 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pot
348 IF (rho(ip) > eps_rho)
THEN
349 pot(ip) = pot(ip) + f*r13(ip)
353 END SUBROUTINE xalpha_lda_1
363 SUBROUTINE xalpha_lda_2(n, rho, r13, pot, sx)
365 INTEGER,
INTENT(IN) :: n
366 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
367 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pot
377 IF (rho(ip) > eps_rho)
THEN
378 pot(ip) = pot(ip) + f*r13(ip)/rho(ip)
382 END SUBROUTINE xalpha_lda_2
392 SUBROUTINE xalpha_lda_3(n, rho, r13, pot, sx)
394 INTEGER,
INTENT(IN) :: n
395 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
396 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pot
402 f = -f23*f13*f43*flda*sx
406 IF (rho(ip) > eps_rho)
THEN
407 pot(ip) = pot(ip) + f*r13(ip)/(rho(ip)*rho(ip))
411 END SUBROUTINE xalpha_lda_3
421 SUBROUTINE xalpha_lsd_0(n, rhoa, r13a, pot, sx)
423 INTEGER,
INTENT(IN) :: n
424 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
425 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pot
439 IF (rhoa(ip) > eps_rho)
THEN
440 pot(ip) = pot(ip) + f*r13a(ip)*rhoa(ip)
445 END SUBROUTINE xalpha_lsd_0
455 SUBROUTINE xalpha_lsd_1(n, rhoa, r13a, pota, sx)
457 INTEGER,
INTENT(IN) :: n
458 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
459 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: pota
473 IF (rhoa(ip) > eps_rho)
THEN
474 pota(ip) = pota(ip) + f*r13a(ip)
479 END SUBROUTINE xalpha_lsd_1
489 SUBROUTINE xalpha_lsd_2(n, rhoa, r13a, potaa, sx)
491 INTEGER,
INTENT(IN) :: n
492 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
493 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: potaa
507 IF (rhoa(ip) > eps_rho)
THEN
508 potaa(ip) = potaa(ip) + f*r13a(ip)/rhoa(ip)
513 END SUBROUTINE xalpha_lsd_2
523 SUBROUTINE xalpha_lsd_3(n, rhoa, r13a, potaaa, sx)
525 INTEGER,
INTENT(IN) :: n
526 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, r13a
527 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: potaaa
535 f = -f23*f13*f43*flsd*sx
541 IF (rhoa(ip) > eps_rho)
THEN
542 potaaa(ip) = potaaa(ip) + f*r13a(ip)/(rhoa(ip)*rhoa(ip))
547 END SUBROUTINE xalpha_lsd_3
559 TYPE(pw_r3d_rs_type),
INTENT(IN) :: rho_a, rho_b
560 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: fxc_aa, fxc_bb
561 REAL(kind=
dp),
INTENT(IN) :: scale_x, eps_rho
564 INTEGER,
DIMENSION(2, 3) :: bo
565 REAL(kind=
dp) :: eaa, ebb, f, flda, flsd, rhoa, rhob
567 flda = -3.0_dp/4.0_dp*(3.0_dp/pi)**f13
568 flsd = flda*2.0_dp**f13
569 f = f13*f43*flsd*scale_x
570 bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
573 DO k = bo(1, 3), bo(2, 3)
574 DO j = bo(1, 2), bo(2, 2)
575 DO i = bo(1, 1), bo(2, 1)
577 rhoa = rho_a%array(i, j, k)
578 IF (rhoa > eps_rho)
THEN
580 fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + f*eaa
582 rhob = rho_b%array(i, j, k)
583 IF (rhob > eps_rho)
THEN
585 fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + f*ebb
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
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 local exchange functional.
subroutine, public xalpha_fxc_eval(rho_a, rho_b, fxc_aa, fxc_bb, scale_x, eps_rho)
...
subroutine, public xalpha_lda_eval(rho_set, deriv_set, order, xa_params, xa_parameter)
...
subroutine, public xalpha_info(lsd, reference, shortform, needs, max_deriv, xa_parameter, scaling)
...
subroutine, public xalpha_lsd_eval(rho_set, deriv_set, order, xa_params, xa_parameter)
...