37 #include "../base/base_uses.f90"
42 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
43 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_xbeef'
48 REAL(kind=
dp),
PARAMETER :: a(0:29) = (/1.516501714304992365356_dp, 0.441353209874497942611_dp, -0.091821352411060291887_dp, &
49 -0.023527543314744041314_dp, 0.034188284548603550816_dp, 0.002411870075717384172_dp, &
50 -0.014163813515916020766_dp, 0.000697589558149178113_dp, 0.009859205136982565273_dp, &
51 -0.006737855050935187551_dp, -0.001573330824338589097_dp, 0.005036146253345903309_dp, &
52 -0.002569472452841069059_dp, -0.000987495397608761146_dp, 0.002033722894696920677_dp, &
53 -0.000801871884834044583_dp, -0.000668807872347525591_dp, 0.001030936331268264214_dp, &
54 -0.000367383865990214423_dp, -0.000421363539352619543_dp, 0.000576160799160517858_dp, &
55 -0.000083465037349510408_dp, -0.000445844758523195788_dp, 0.000460129009232047457_dp, &
56 -0.000005231775398304339_dp, -0.000423957047149510404_dp, 0.000375019067938866537_dp, &
57 0.000021149381251344578_dp, -0.000190491156503997170_dp, 0.000073843624209823442_dp/)
73 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
74 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
75 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
77 IF (
PRESENT(reference))
THEN
78 reference =
"Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LDA}"
80 IF (
PRESENT(shortform))
THEN
81 shortform =
"Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LDA}"
83 IF (
PRESENT(needs))
THEN
85 needs%rho_1_3 = .true.
86 needs%norm_drho = .true.
88 IF (
PRESENT(max_deriv)) max_deriv = 1
104 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
105 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
106 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
108 IF (
PRESENT(reference))
THEN
109 reference =
"Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LSD}"
111 IF (
PRESENT(shortform))
THEN
112 shortform =
"Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LSD}"
114 IF (
PRESENT(needs))
THEN
115 needs%rho_spin = .true.
116 needs%rho_spin_1_3 = .true.
117 needs%norm_drho_spin = .true.
119 IF (
PRESENT(max_deriv)) max_deriv = 1
137 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
138 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
139 INTEGER,
INTENT(in) :: grad_deriv
140 TYPE(section_vals_type),
POINTER :: xbeef_params
142 CHARACTER(len=*),
PARAMETER :: routinen =
'xbeef_lda_eval'
144 INTEGER :: handle, npoints
145 INTEGER,
DIMENSION(2, 3) :: bo
146 REAL(kind=
dp) :: epsilon_rho, sx
147 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
148 POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
150 TYPE(xc_derivative_type),
POINTER :: deriv
152 CALL timeset(routinen, handle)
159 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
160 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
168 IF (grad_deriv >= 0)
THEN
170 allocate_deriv=.true.)
173 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
175 allocate_deriv=.true.)
178 allocate_deriv=.true.)
181 IF (grad_deriv > 1 .OR. grad_deriv < -1)
THEN
182 cpabort(
"derivatives greater than 1 not implemented")
190 CALL xbeef_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
191 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
192 grad_deriv=grad_deriv, &
193 npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
196 CALL timestop(handle)
217 SUBROUTINE xbeef_lda_calc(rho, rho_1_3, norm_drho, &
218 e_0, e_rho, e_ndrho, &
219 grad_deriv, npoints, epsilon_rho, sx)
220 INTEGER,
INTENT(in) :: npoints, grad_deriv
221 REAL(kind=
dp),
DIMENSION(1:npoints),
INTENT(inout) :: e_ndrho, e_rho, e_0
222 REAL(kind=
dp),
DIMENSION(1:npoints),
INTENT(in) :: norm_drho, rho_1_3, rho
223 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, sx
225 INTEGER,
PARAMETER :: m = 30
228 REAL(kind=
dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
229 epsilon_rho43, k_lda, kf, my_rho, &
230 my_rho_1_3, s, s2, t, t3
231 REAL(kind=
dp),
DIMENSION(0:29) :: de_leg, e_leg
235 kf = (3.0_dp*
pi**2)**(1.0_dp/3.0_dp)
236 k_lda = -(2.0_dp/(2.0_dp**(4._dp/3._dp)))*(3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*
pi))**(1.0_dp/3.0_dp)
239 epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
246 IF (my_rho > epsilon_rho)
THEN
247 my_rho_1_3 = rho_1_3(ii)
249 e_ueg = k_lda*my_rho*my_rho_1_3
250 e_ueg_drho = (4.0_dp/3.0_dp)*k_lda*my_rho_1_3
252 t3 = my_rho_1_3*my_rho*2*kf
254 s = norm_drho(ii)/max(t3, epsilon_rho43)
256 t = 2.0_dp*s2/(4.0_dp + s2) - 1.0_dp
258 IF (grad_deriv >= 0)
THEN
263 IF ((grad_deriv >= 1) .OR. (grad_deriv == -1))
THEN
266 dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
267 ds_rho = -(4.0_dp*s)/(3.0_dp*max(my_rho, epsilon_rho))
268 ds_ndrho = 1.0_dp/(max(t3, epsilon_rho43))
272 e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(real(i, kind=
dp))
275 IF (abs(grad_deriv) >= 1)
THEN
277 de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
282 IF (grad_deriv >= 0)
THEN
284 e_0(ii) = e_0(ii) + sum(e_leg*a)*e_ueg*sx
288 IF ((grad_deriv >= 1) .OR. (grad_deriv == -1))
THEN
289 e_rho(ii) = e_rho(ii) + (sum(e_leg*a)*e_ueg_drho + sum(de_leg*a)*dt*ds_rho*e_ueg)*sx
290 e_ndrho(ii) = e_ndrho(ii) + (sum(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
298 END SUBROUTINE xbeef_lda_calc
311 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
312 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
313 INTEGER,
INTENT(in) :: grad_deriv
314 TYPE(section_vals_type),
POINTER :: xbeef_params
316 CHARACTER(len=*),
PARAMETER :: routinen =
'xbeef_lsd_eval'
318 INTEGER :: handle, i, ispin, npoints
319 INTEGER,
DIMENSION(2, 3) :: bo
320 REAL(kind=
dp) :: epsilon_rho, sx
321 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :), &
322 POINTER :: dummy, e_0
323 TYPE(cp_3d_r_cp_type),
DIMENSION(2) :: e_ndrho, e_rho, norm_drho, rho, rho_1_3
324 TYPE(xc_derivative_type),
POINTER :: deriv
326 CALL timeset(routinen, handle)
332 NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
337 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
338 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
339 norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
341 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
343 dummy => rho(1)%array
347 e_rho(i)%array => dummy
348 e_ndrho(i)%array => dummy
351 IF (grad_deriv >= 0)
THEN
353 allocate_deriv=.true.)
356 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
358 allocate_deriv=.true.)
361 allocate_deriv=.true.)
364 allocate_deriv=.true.)
367 allocate_deriv=.true.)
370 IF (grad_deriv > 1 .OR. grad_deriv < -1)
THEN
371 cpabort(
"derivatives greater than 1 not implemented")
382 CALL xbeef_lsd_calc( &
383 rho_spin=rho(ispin)%array, &
384 rho_1_3_spin=rho_1_3(ispin)%array, &
385 norm_drho_spin=norm_drho(ispin)%array, &
386 e_0=e_0, e_rho_spin=e_rho(ispin)%array, &
387 e_ndrho_spin=e_ndrho(ispin)%array, &
388 grad_deriv=grad_deriv, npoints=npoints, &
389 epsilon_rho=epsilon_rho, sx=sx)
395 CALL timestop(handle)
415 SUBROUTINE xbeef_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
416 e_rho_spin, e_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
417 REAL(kind=
dp),
DIMENSION(*),
INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
418 REAL(kind=
dp),
DIMENSION(*),
INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin
419 INTEGER,
INTENT(in) :: grad_deriv, npoints
420 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, sx
422 INTEGER,
PARAMETER :: m = 30
425 REAL(kind=
dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
426 epsilon_rho43, k_lsd, kf, &
427 my_epsilon_rho, my_rho, my_rho_1_3, s, &
429 REAL(kind=
dp),
DIMENSION(0:29) :: de_leg, e_leg
433 kf = (3.0_dp*
pi**2)**(1.0_dp/3.0_dp)
434 k_lsd = (3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*
pi))**(1.0_dp/3.0_dp)
437 my_epsilon_rho = 0.5_dp*epsilon_rho
438 epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
442 my_rho = rho_spin(ii)
444 IF (my_rho > epsilon_rho)
THEN
445 my_rho_1_3 = rho_1_3_spin(ii)
447 e_ueg = k_lsd*my_rho*my_rho_1_3
448 e_ueg_drho = (4.0_dp/3.0_dp)*k_lsd*my_rho_1_3
450 t3 = my_rho_1_3*my_rho*2*kf
452 s = norm_drho_spin(ii)/max(t3, epsilon_rho43)
454 t = 2.0_dp*s**2/(4.0_dp + s**2) - 1.0_dp
456 IF (grad_deriv >= 0)
THEN
461 IF ((grad_deriv >= 1) .OR. (grad_deriv == -1))
THEN
464 dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
465 ds_rho = -(4.0_dp*s)/(3.0_dp*max(my_rho, epsilon_rho))
466 ds_ndrho = 1.0_dp/(max(t3, epsilon_rho43))
470 e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(real(i, kind=
dp))
473 IF (abs(grad_deriv) >= 1)
THEN
475 de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
481 IF (grad_deriv >= 0)
THEN
483 e_0(ii) = e_0(ii) + sum(e_leg*a)*e_ueg*sx
487 IF ((grad_deriv >= 1) .OR. (grad_deriv == -1))
THEN
488 e_rho_spin(ii) = e_rho_spin(ii) + (sum(e_leg*a)*e_ueg_drho + sum(de_leg*a)*dt*ds_rho*e_ueg)*sx
489 e_ndrho_spin(ii) = e_ndrho_spin(ii) + (sum(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
495 END SUBROUTINE xbeef_lsd_calc
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public wellendorff2012
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_norm_drho
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
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
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
calculates the Exchange contribution in the BEEF-vdW functional
subroutine, public xbeef_lda_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
evaluates the beef exchange functional for lda
subroutine, public xbeef_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
subroutine, public xbeef_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
subroutine, public xbeef_lsd_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
evaluates the beef 88 exchange functional for lsd