(git:374b731)
Loading...
Searching...
No Matches
xc_xbr_pbe_lda_hole_t_c_lr.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 This functional is a combination of three different exchange hole
10!> models. The ingredients are:
11!>
12!> 1. Becke Roussel exchange hole
13!> 2. PBE exchange hole
14!> 3. LDA exchange hole
15!>
16!> The full functionals is given as follows
17!>
18!> Fx = eps_lr_lda/eps_lr_br
19!> Fcorr = alpha/( exp( (Fx-mu)/N ) + 1)
20!> rhox = Fcorr * eps_lr_pbe + (1-Fcorr)*eps_lr_br
21!> eps = int_{R}^{\infty} rhox*s*ds
22!>
23!> with alpha, mu and N fitting parameters
24!> \par History
25!> 01.2009 created [mguidon]
26!> \author mguidon
27! **************************************************************************************************
28
29MODULE xc_xbr_pbe_lda_hole_t_c_lr
30
31 USE input_section_types, ONLY: section_vals_type,&
32 section_vals_val_get
33 USE kinds, ONLY: dp
34 USE mathconstants, ONLY: pi
35 USE xc_derivative_desc, ONLY: &
36 deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, deriv_norm_drho, &
37 deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, deriv_tau, &
38 deriv_tau_a, deriv_tau_b
39 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
40 xc_dset_get_derivative
41 USE xc_derivative_types, ONLY: xc_derivative_get,&
42 xc_derivative_type
43 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
44 USE xc_rho_set_types, ONLY: xc_rho_set_get,&
45 xc_rho_set_type
46 USE xc_xbecke_roussel, ONLY: x_br_lsd_y_gt_0,&
47 x_br_lsd_y_gt_0_cutoff,&
48 x_br_lsd_y_lte_0,&
49 x_br_lsd_y_lte_0_cutoff
50 USE xc_xlda_hole_t_c_lr, ONLY: xlda_hole_t_c_lr_lda_calc_0
51 USE xc_xpbe_hole_t_c_lr, ONLY: xpbe_hole_t_c_lr_lda_calc_1,&
52 xpbe_hole_t_c_lr_lda_calc_2
53#include "../base/base_uses.f90"
54
55 IMPLICIT NONE
56 PRIVATE
57
58 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbr_pbe_lda_hole_t_c_lr'
60
61 REAL(dp), PARAMETER, PRIVATE :: br_a1 = 1.5255251812009530_dp, &
62 br_a2 = 0.4576575543602858_dp, &
63 br_a3 = 0.4292036732051034_dp, &
64 br_c0 = 0.7566445420735584_dp, &
65 br_c1 = -2.6363977871370960_dp, &
66 br_c2 = 5.4745159964232880_dp, &
67 br_c3 = -12.657308127108290_dp, &
68 br_c4 = 4.1250584725121360_dp, &
69 br_c5 = -30.425133957163840_dp, &
70 br_b0 = 0.4771976183772063_dp, &
71 br_b1 = -1.7799813494556270_dp, &
72 br_b2 = 3.8433841862302150_dp, &
73 br_b3 = -9.5912050880518490_dp, &
74 br_b4 = 2.1730180285916720_dp, &
75 br_b5 = -30.425133851603660_dp, &
76 br_d0 = 0.00004435009886795587_dp, &
77 br_d1 = 0.58128653604457910_dp, &
78 br_d2 = 66.742764515940610_dp, &
79 br_d3 = 434.26780897229770_dp, &
80 br_d4 = 824.7765766052239000_dp, &
81 br_d5 = 1657.9652731582120_dp, &
82 br_e0 = 0.00003347285060926091_dp, &
83 br_e1 = 0.47917931023971350_dp, &
84 br_e2 = 62.392268338574240_dp, &
85 br_e3 = 463.14816427938120_dp, &
86 br_e4 = 785.2360350104029000_dp, &
87 br_e5 = 1657.962968223273000000_dp, &
88 br_bb = 2.085749716493756_dp
89
90 REAL(dp), PARAMETER, PRIVATE :: smax = 8.572844_dp, &
91 scutoff = 8.3_dp, &
92 sconst = 18.79622316_dp, &
93 gcutoff = 0.08_dp
94
95 REAL(dp), PARAMETER, PRIVATE :: alpha = 0.3956891_dp, &
96 n = -0.0009800242_dp, &
97 mu = 0.00118684_dp
98
99 PUBLIC :: xbr_pbe_lda_hole_tc_lr_lda_info, &
100 xbr_pbe_lda_hole_tc_lr_lsd_info, &
101 xbr_pbe_lda_hole_tc_lr_lda_eval, &
102 xbr_pbe_lda_hole_tc_lr_lsd_eval
103CONTAINS
104
105! **************************************************************************************************
106!> \brief return various information on the functional
107!> \param reference string with the reference of the actual functional
108!> \param shortform string with the shortform of the functional name
109!> \param needs the components needed by this functional are set to
110!> true (does not set the unneeded components to false)
111!> \param max_deriv ...
112!> \author mguidon (01.2009)
113! **************************************************************************************************
114 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info(reference, shortform, needs, max_deriv)
115 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
116 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
117 INTEGER, INTENT(out), OPTIONAL :: max_deriv
118
119 IF (PRESENT(reference)) THEN
120 reference = "{LDA version}"
121 END IF
122 IF (PRESENT(shortform)) THEN
123 shortform = "{LDA}"
124 END IF
125
126 IF (PRESENT(needs)) THEN
127 needs%rho = .true.
128 needs%norm_drho = .true.
129 needs%tau = .true.
130 needs%laplace_rho = .true.
131 END IF
132
133 IF (PRESENT(max_deriv)) max_deriv = 1
134
135 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info
136
137! **************************************************************************************************
138!> \brief return various information on the functional
139!> \param reference string with the reference of the actual functional
140!> \param shortform string with the shortform of the functional name
141!> \param needs the components needed by this functional are set to
142!> true (does not set the unneeded components to false)
143!> \param max_deriv ...
144!> \author mguidon (01.2009)
145! **************************************************************************************************
146 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info(reference, shortform, needs, max_deriv)
147 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
148 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
149 INTEGER, INTENT(out), OPTIONAL :: max_deriv
150
151 IF (PRESENT(reference)) THEN
152 reference = "{LDA version}"
153 END IF
154 IF (PRESENT(shortform)) THEN
155 shortform = "{LDA}"
156 END IF
157
158 IF (PRESENT(needs)) THEN
159 needs%rho_spin = .true.
160 needs%norm_drho_spin = .true.
161 needs%tau_spin = .true.
162 needs%laplace_rho_spin = .true.
163 END IF
164 IF (PRESENT(max_deriv)) max_deriv = 1
165
166 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info
167
168! **************************************************************************************************
169!> \brief Intermediate routine that gets grids, derivatives and some params
170!> \param rho_set the density where you want to evaluate the functional
171!> \param deriv_set place where to store the functional derivatives (they are
172!> added to the derivatives)
173!> \param grad_deriv degree of the derivative that should be evaluated,
174!> if positive all the derivatives up to the given degree are evaluated,
175!> if negative only the given degree is calculated
176!> \param params parameters for functional
177!> \author mguidon (01.2009)
178! **************************************************************************************************
179 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval(rho_set, deriv_set, grad_deriv, params)
180 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
181 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
182 INTEGER, INTENT(in) :: grad_deriv
183 TYPE(section_vals_type), POINTER :: params
184
185 CHARACTER(len=*), PARAMETER :: routinen = 'xbr_pbe_lda_hole_tc_lr_lda_eval'
186
187 INTEGER :: handle, npoints
188 INTEGER, DIMENSION(2, 3) :: bo
189 REAL(dp) :: gamma, r, sx
190 REAL(kind=dp) :: epsilon_rho
191 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
192 POINTER :: dummy, e_0, e_laplace_rho, e_ndrho, &
193 e_rho, e_tau, laplace_rho, norm_drho, &
194 rho, tau
195 TYPE(xc_derivative_type), POINTER :: deriv
196
197 CALL timeset(routinen, handle)
198
199 CALL xc_rho_set_get(rho_set, rho=rho, norm_drho=norm_drho, &
200 tau=tau, laplace_rho=laplace_rho, local_bounds=bo, &
201 rho_cutoff=epsilon_rho)
202 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
203
204 dummy => rho
205
206 e_0 => dummy
207 e_rho => dummy
208 e_ndrho => dummy
209 e_tau => dummy
210 e_laplace_rho => dummy
211
212 IF (grad_deriv >= 0) THEN
213 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
214 allocate_deriv=.true.)
215 CALL xc_derivative_get(deriv, deriv_data=e_0)
216 END IF
217 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
218 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
219 allocate_deriv=.true.)
220 CALL xc_derivative_get(deriv, deriv_data=e_rho)
221 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
222 allocate_deriv=.true.)
223 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
224 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
225 allocate_deriv=.true.)
226 CALL xc_derivative_get(deriv, deriv_data=e_tau)
227 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
228 allocate_deriv=.true.)
229 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
230 END IF
231 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
232 cpabort("derivatives bigger than 1 not implemented")
233 END IF
234
235 CALL section_vals_val_get(params, "scale_x", r_val=sx)
236 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
237 CALL section_vals_val_get(params, "GAMMA", r_val=gamma)
238
239 IF (r == 0.0_dp) THEN
240 cpabort("Cutoff_Radius 0.0 not implemented")
241 END IF
242
243!$OMP PARALLEL DEFAULT(NONE) &
244!$OMP SHARED(rho, norm_drho, laplace_rho, tau, e_0, e_rho) &
245!$OMP SHARED(e_ndrho, e_tau, e_laplace_rho, grad_deriv) &
246!$OMP SHARED(npoints, epsilon_rho) &
247!$OMP SHARED(sx, r, gamma)
248
249 CALL xbr_pbe_lda_hole_tc_lr_lda_calc(rho=rho, norm_drho=norm_drho, &
250 laplace_rho=laplace_rho, tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
251 e_tau=e_tau, e_laplace_rho=e_laplace_rho, grad_deriv=grad_deriv, &
252 npoints=npoints, epsilon_rho=epsilon_rho, sx=sx, r=r, gamma=gamma)
253
254!$OMP END PARALLEL
255
256 CALL timestop(handle)
257 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval
258
259! **************************************************************************************************
260!> \brief Low level routine that calls the three involved holes and puts them
261!> together
262!> \param rho values on the grid
263!> \param norm_drho values on the grid
264!> \param laplace_rho values on the grid
265!> \param tau values on the grid
266!> \param e_0 derivatives on the grid
267!> \param e_rho derivatives on the grid
268!> \param e_ndrho derivatives on the grid
269!> \param e_tau derivatives on the grid
270!> \param e_laplace_rho derivatives on the grid
271!> \param grad_deriv degree of the derivative that should be evaluated,
272!> if positive all the derivatives up to the given degree are evaluated,
273!> if negative only the given degree is calculated
274!> \param npoints number of gridpoints
275!> \param epsilon_rho cutoffs
276!> \param sx parameters for functional
277!> \param R parameters for functional
278!> \param gamma parameters for functional
279!> \author mguidon (01.2009)
280! **************************************************************************************************
281 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
282 e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
283 epsilon_rho, sx, R, gamma)
284
285 INTEGER, INTENT(in) :: npoints, grad_deriv
286 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
287 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
288 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r, gamma
289
290 INTEGER :: ip
291 REAL(dp) :: dfermi_dlaplace_rho, dfermi_dndrho, dfermi_drho, dfermi_dtau, e_0_br, e_0_lda, &
292 e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, e_rho_br, e_rho_lda, e_rho_pbe, &
293 e_tau_br, fermi, fx, my_laplace_rho, my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, &
294 t16, t2, t3, t4, t5, t6, t7, t8, t9, yval
295
296!$OMP DO
297
298 DO ip = 1, npoints
299 my_rho = 0.5_dp*max(rho(ip), 0.0_dp)
300 IF (my_rho > epsilon_rho) THEN
301 my_ndrho = 0.5_dp*max(norm_drho(ip), epsilon(0.0_dp)*1.e4_dp)
302 my_tau = 0.5_dp*max(epsilon(0.0_dp)*1.e4_dp, tau(ip))
303 my_laplace_rho = 0.5_dp*laplace_rho(ip)
304
305 ! ** We calculate first the Becke-Roussel part, saving everything in local variables
306 t1 = pi**(0.1e1_dp/0.3e1_dp)
307 t2 = t1**2
308 t3 = my_rho**(0.1e1_dp/0.3e1_dp)
309 t4 = t3**2
310 t5 = t4*my_rho
311 t8 = my_ndrho**2
312 t9 = 0.1e1_dp/my_rho
313 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
314 t16 = 0.1e1_dp/t15
315 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
316
317 e_0_br = 0.0_dp
318 e_rho_br = 0.0_dp
319 e_ndrho_br = 0.0_dp
320 e_tau_br = 0.0_dp
321 e_laplace_rho_br = 0.0_dp
322
323 IF (r == 0.0_dp) THEN
324 IF (yval <= 0.0_dp) THEN
325 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
326 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
327 sx, gamma, grad_deriv)
328 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
329 e_0_br = 2.0_dp*e_0_br
330 ELSE
331 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
332 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
333 sx, gamma, grad_deriv)
334 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
335 e_0_br = 2.0_dp*e_0_br
336 END IF
337 ELSE
338 IF (yval <= 0.0_dp) THEN
339 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
340 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
341 sx, r, gamma, grad_deriv)
342 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
343 e_0_br = 2.0_dp*e_0_br
344 ELSE
345 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
346 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
347 sx, r, gamma, grad_deriv)
348 ! VERY UGLY HACK e_0 has to multiplied by the factor 2
349 e_0_br = 2.0_dp*e_0_br
350 END IF
351 END IF
352
353 ! ** Now we calculate the pbe cutoff part
354 ! ** Attention we need to scale rho, ndrho first
355 my_rho = my_rho*2.0_dp
356 my_ndrho = my_ndrho*2.0_dp
357
358 ! ** Do some precalculation in order to catch the correct branch afterwards
359 sscale = 1.0_dp
360 t1 = pi**2
361 t2 = t1*my_rho
362 t3 = t2**(0.1e1_dp/0.3e1_dp)
363 t4 = 0.1e1_dp/t3
364 t6 = my_ndrho*t4
365 t7 = 0.1e1_dp/my_rho
366 t8 = t7*sscale
367 ss = 0.3466806371753173524216762e0_dp*t6*t8
368 IF (ss > scutoff) THEN
369 ss2 = ss*ss
370 sscale = (smax*ss2 - sconst)/(ss2*ss)
371 END IF
372 e_0_pbe = 0.0_dp
373 e_rho_pbe = 0.0_dp
374 e_ndrho_pbe = 0.0_dp
375 IF (ss*sscale > gcutoff) THEN
376 CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
377 my_rho, &
378 my_ndrho, sscale, sx, r, grad_deriv)
379 ELSE
380 CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
381 my_rho, &
382 my_ndrho, sscale, sx, r, grad_deriv)
383 END IF
384
385 ! ** Finally we get the LDA part
386
387 e_0_lda = 0.0_dp
388 e_rho_lda = 0.0_dp
389 CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda, &
390 sx, r)
391
392 fx = e_0_br/e_0_lda
393
394 fermi = alpha/(exp((fx - mu)/n) + 1.0_dp)
395
396 dfermi_drho = -fermi**2/alpha/n*(e_rho_br/e_0_lda - e_0_br*e_rho_lda/e_0_lda**2)*exp((fx - mu)/n)
397 dfermi_dndrho = -fermi**2/alpha/n*(e_ndrho_br/e_0_lda)*exp((fx - mu)/n)
398 dfermi_dtau = -fermi**2/alpha/n*(e_tau_br/e_0_lda)*exp((fx - mu)/n)
399 dfermi_dlaplace_rho = -fermi**2/alpha/n*(e_laplace_rho_br/e_0_lda)*exp((fx - mu)/n)
400
401 e_0(ip) = e_0(ip) + (fermi*e_0_pbe + (1.0_dp - fermi)*e_0_br)*sx
402
403 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
404
405 e_rho(ip) = e_rho(ip) + (fermi*e_rho_pbe + dfermi_drho*e_0_pbe + &
406 (1.0_dp - fermi)*e_rho_br - dfermi_drho*e_0_br)*sx
407
408 e_ndrho(ip) = e_ndrho(ip) + (fermi*e_ndrho_pbe + dfermi_dndrho*e_0_pbe + &
409 (1.0_dp - fermi)*e_ndrho_br - dfermi_dndrho*e_0_br)*sx
410
411 e_tau(ip) = e_tau(ip) + (dfermi_dtau*e_0_pbe + &
412 (1.0_dp - fermi)*e_tau_br - dfermi_dtau*e_0_br)*sx
413
414 e_laplace_rho(ip) = e_laplace_rho(ip) + (dfermi_dlaplace_rho*e_0_pbe + &
415 (1.0_dp - fermi)*e_laplace_rho_br - dfermi_dlaplace_rho*e_0_br)*sx
416 END IF
417
418 END IF
419 END DO
420
421!$OMP END DO
422
423 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc
424
425! **************************************************************************************************
426!> \brief Intermediate routine that gets grids, derivatives and some params
427!> \param rho_set the density where you want to evaluate the functional
428!> \param deriv_set place where to store the functional derivatives (they are
429!> added to the derivatives)
430!> \param grad_deriv degree of the derivative that should be evaluated,
431!> if positive all the derivatives up to the given degree are evaluated,
432!> if negative only the given degree is calculated
433!> \param params parameters for functional
434!> \author mguidon (01.2009)
435! **************************************************************************************************
436 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval(rho_set, deriv_set, grad_deriv, params)
437 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
438 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
439 INTEGER, INTENT(in) :: grad_deriv
440 TYPE(section_vals_type), POINTER :: params
441
442 CHARACTER(len=*), PARAMETER :: routinen = 'xbr_pbe_lda_hole_tc_lr_lsd_eval'
443
444 INTEGER :: handle, npoints
445 INTEGER, DIMENSION(2, 3) :: bo
446 REAL(dp) :: gamma, r, sx
447 REAL(kind=dp) :: epsilon_rho
448 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
449 e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, &
450 laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
451 TYPE(xc_derivative_type), POINTER :: deriv
452
453 CALL timeset(routinen, handle)
454
455 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
456 norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, laplace_rhoa=laplace_rhoa, &
457 laplace_rhob=laplace_rhob, local_bounds=bo, &
458 rho_cutoff=epsilon_rho)
459 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
460
461 dummy => rhoa
462
463 e_0 => dummy
464 e_rhoa => dummy
465 e_rhob => dummy
466 e_ndrhoa => dummy
467 e_ndrhob => dummy
468 e_tau_a => dummy
469 e_tau_b => dummy
470 e_laplace_rhoa => dummy
471 e_laplace_rhob => dummy
472
473 IF (grad_deriv >= 0) THEN
474 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
475 allocate_deriv=.true.)
476 CALL xc_derivative_get(deriv, deriv_data=e_0)
477 END IF
478 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
479 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
480 allocate_deriv=.true.)
481 CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
482 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
483 allocate_deriv=.true.)
484 CALL xc_derivative_get(deriv, deriv_data=e_rhob)
485
486 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
487 allocate_deriv=.true.)
488 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
489 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
490 allocate_deriv=.true.)
491 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
492
493 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
494 allocate_deriv=.true.)
495 CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
496 deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
497 allocate_deriv=.true.)
498 CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
499
500 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
501 allocate_deriv=.true.)
502 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
503 deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
504 allocate_deriv=.true.)
505 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
506 END IF
507 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
508 cpabort("derivatives bigger than 1 not implemented")
509 END IF
510
511 CALL section_vals_val_get(params, "scale_x", r_val=sx)
512 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
513 CALL section_vals_val_get(params, "GAMMA", r_val=gamma)
514
515 IF (r == 0.0_dp) THEN
516 cpabort("Cutoff_Radius 0.0 not implemented")
517 END IF
518
519!$OMP PARALLEL DEFAULT(NONE) &
520!$OMP SHARED(rhoa, norm_drhoa, laplace_rhoa, tau_a, e_0) &
521!$OMP SHARED(e_rhoa, e_ndrhoa, e_tau_a, e_laplace_rhoa) &
522!$OMP SHARED(grad_deriv, npoints, epsilon_rho) &
523!$OMP SHARED(sx, r, gamma) &
524!$OMP SHARED(rhob, norm_drhob, laplace_rhob, tau_b, e_rhob) &
525!$OMP SHARED(e_ndrhob, e_tau_b, e_laplace_rhob)
526
527 CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhoa, norm_drho=norm_drhoa, &
528 laplace_rho=laplace_rhoa, tau=tau_a, e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, &
529 e_tau=e_tau_a, e_laplace_rho=e_laplace_rhoa, grad_deriv=grad_deriv, &
530 npoints=npoints, epsilon_rho=epsilon_rho, &
531 sx=sx, r=r, gamma=gamma)
532
533 CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhob, norm_drho=norm_drhob, &
534 laplace_rho=laplace_rhob, tau=tau_b, e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, &
535 e_tau=e_tau_b, e_laplace_rho=e_laplace_rhob, grad_deriv=grad_deriv, &
536 npoints=npoints, epsilon_rho=epsilon_rho, &
537 sx=sx, r=r, gamma=gamma)
538
539!$OMP END PARALLEL
540
541 CALL timestop(handle)
542 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval
543! **************************************************************************************************
544!> \brief Low level routine that calls the three involved holes and puts them
545!> together
546!> \param rho values on the grid
547!> \param norm_drho values on the grid
548!> \param laplace_rho values on the grid
549!> \param tau values on the grid
550!> \param e_0 derivatives on the grid
551!> \param e_rho derivatives on the grid
552!> \param e_ndrho derivatives on the grid
553!> \param e_tau derivatives on the grid
554!> \param e_laplace_rho derivatives on the grid
555!> \param grad_deriv degree of the derivative that should be evaluated,
556!> if positive all the derivatives up to the given degree are evaluated,
557!> if negative only the given degree is calculated
558!> \param npoints number of gridpoints
559!> \param epsilon_rho cutoffs
560!> \param sx parameters for functional
561!> \param R parameters for functional
562!> \param gamma parameters for functional
563!> \author mguidon (01.2009)
564! **************************************************************************************************
565 SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
566 e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
567 epsilon_rho, sx, R, gamma)
568
569 INTEGER, INTENT(in) :: npoints, grad_deriv
570 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
571 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
572 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r, gamma
573
574 INTEGER :: ip
575 REAL(dp) :: dfermi_dlaplace_rho, dfermi_dndrho, dfermi_drho, dfermi_dtau, e_0_br, e_0_lda, &
576 e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, e_rho_br, e_rho_lda, e_rho_pbe, &
577 e_tau_br, fermi, fx, my_laplace_rho, my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, &
578 t16, t2, t3, t4, t5, t6, t7, t8, t9, yval
579
580!$OMP DO
581
582 DO ip = 1, npoints
583 my_rho = max(rho(ip), 0.0_dp)
584 IF (my_rho > epsilon_rho) THEN
585 my_ndrho = max(norm_drho(ip), epsilon(0.0_dp)*1.e4_dp)
586 my_tau = 1.0_dp*max(epsilon(0.0_dp)*1.e4_dp, tau(ip))
587 my_laplace_rho = 1.0_dp*laplace_rho(ip)
588
589 t1 = pi**(0.1e1_dp/0.3e1_dp)
590 t2 = t1**2
591 t3 = my_rho**(0.1e1_dp/0.3e1_dp)
592 t4 = t3**2
593 t5 = t4*my_rho
594 t8 = my_ndrho**2
595 t9 = 0.1e1_dp/my_rho
596 t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
597 t16 = 0.1e1_dp/t15
598 yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
599
600 e_0_br = 0.0_dp
601 e_rho_br = 0.0_dp
602 e_ndrho_br = 0.0_dp
603 e_tau_br = 0.0_dp
604 e_laplace_rho_br = 0.0_dp
605
606 IF (r == 0.0_dp) THEN
607 IF (yval <= 0.0_dp) THEN
608 CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
609 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
610 sx, gamma, grad_deriv)
611 ELSE
612 CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
613 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
614 sx, gamma, grad_deriv)
615 END IF
616 ELSE
617 IF (yval <= 0.0_dp) THEN
618 CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
619 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
620 sx, r, gamma, grad_deriv)
621 ELSE
622 CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
623 e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
624 sx, r, gamma, grad_deriv)
625 END IF
626 END IF
627
628 ! ** Now we calculate the pbe cutoff part
629 ! ** Attention we need to scale rho, ndrho first
630 my_rho = my_rho*2.0_dp
631 my_ndrho = my_ndrho*2.0_dp
632
633 ! ** Do some precalculation in order to catch the correct branch afterwards
634 sscale = 1.0_dp
635 t1 = pi**2
636 t2 = t1*my_rho
637 t3 = t2**(0.1e1_dp/0.3e1_dp)
638 t4 = 0.1e1_dp/t3
639 t6 = my_ndrho*t4
640 t7 = 0.1e1_dp/my_rho
641 t8 = t7*sscale
642 ss = 0.3466806371753173524216762e0_dp*t6*t8
643 IF (ss > scutoff) THEN
644 ss2 = ss*ss
645 sscale = (smax*ss2 - sconst)/(ss2*ss)
646 END IF
647 e_0_pbe = 0.0_dp
648 e_rho_pbe = 0.0_dp
649 e_ndrho_pbe = 0.0_dp
650 IF (ss*sscale > gcutoff) THEN
651 CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
652 my_rho, &
653 my_ndrho, sscale, sx, r, grad_deriv)
654 ELSE
655 CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
656 my_rho, &
657 my_ndrho, sscale, sx, r, grad_deriv)
658 END IF
659
660 e_0_pbe = 0.5_dp*e_0_pbe
661
662 ! ** Finally we get the LDA part
663
664 e_0_lda = 0.0_dp
665 e_rho_lda = 0.0_dp
666 CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda, &
667 sx, r)
668 e_0_lda = 0.5_dp*e_0_lda
669
670 fx = e_0_br/e_0_lda
671
672 fermi = alpha/(exp((fx - mu)/n) + 1.0_dp)
673
674 dfermi_drho = -fermi**2/alpha/n*(e_rho_br/e_0_lda - e_0_br*e_rho_lda/e_0_lda**2)*exp((fx - mu)/n)
675 dfermi_dndrho = -fermi**2/alpha/n*(e_ndrho_br/e_0_lda)*exp((fx - mu)/n)
676 dfermi_dtau = -fermi**2/alpha/n*(e_tau_br/e_0_lda)*exp((fx - mu)/n)
677 dfermi_dlaplace_rho = -fermi**2/alpha/n*(e_laplace_rho_br/e_0_lda)*exp((fx - mu)/n)
678
679 e_0(ip) = e_0(ip) + (fermi*e_0_pbe + (1.0_dp - fermi)*e_0_br)*sx
680
681 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
682
683 e_rho(ip) = e_rho(ip) + (fermi*e_rho_pbe + dfermi_drho*e_0_pbe + &
684 (1.0_dp - fermi)*e_rho_br - dfermi_drho*e_0_br)*sx
685
686 e_ndrho(ip) = e_ndrho(ip) + (fermi*e_ndrho_pbe + dfermi_dndrho*e_0_pbe + &
687 (1.0_dp - fermi)*e_ndrho_br - dfermi_dndrho*e_0_br)*sx
688
689 e_tau(ip) = e_tau(ip) + (dfermi_dtau*e_0_pbe + &
690 (1.0_dp - fermi)*e_tau_br - dfermi_dtau*e_0_br)*sx
691
692 e_laplace_rho(ip) = e_laplace_rho(ip) + (dfermi_dlaplace_rho*e_0_pbe + &
693 (1.0_dp - fermi)*e_laplace_rho_br - dfermi_dlaplace_rho*e_0_br)*sx
694 END IF
695
696 END IF
697 END DO
698
699!$OMP END DO
700
701 END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc
702
703END MODULE xc_xbr_pbe_lda_hole_t_c_lr
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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_laplace_rhob
integer, parameter, public deriv_laplace_rhob
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_tau
integer, parameter, public deriv_tau
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_b
integer, parameter, public deriv_tau_b
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_a
integer, parameter, public deriv_tau_a
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_laplace_rhoa
integer, parameter, public deriv_laplace_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_desc::deriv_laplace_rho
integer, parameter, public deriv_laplace_rho
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_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_xbecke_roussel
Calculates the exchange energy based on the Becke-Roussel exchange hole. Takes advantage of an analyt...
Definition xc_xbecke_roussel.F:18
xc_xbecke_roussel::x_br_lsd_y_gt_0
subroutine, public x_br_lsd_y_gt_0(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, gamma, grad_deriv)
Full range evaluation for y > 0.
Definition xc_xbecke_roussel.F:818
xc_xbecke_roussel::x_br_lsd_y_gt_0_cutoff
subroutine, public x_br_lsd_y_gt_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, r, gamma, grad_deriv)
Truncated long range part for y > 0.
Definition xc_xbecke_roussel.F:1834
xc_xbecke_roussel::x_br_lsd_y_lte_0_cutoff
subroutine, public x_br_lsd_y_lte_0_cutoff(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, r, gamma, grad_deriv)
Truncated long range part for y <= 0.
Definition xc_xbecke_roussel.F:1076
xc_xbecke_roussel::x_br_lsd_y_lte_0
subroutine, public x_br_lsd_y_lte_0(rho, ndrho, tau, laplace_rho, e_0, e_rho, e_ndrho, e_tau, e_laplace_rho, sx, gamma, grad_deriv)
full range evaluation for y <= 0
Definition xc_xbecke_roussel.F:574
xc_xbr_pbe_lda_hole_t_c_lr
This functional is a combination of three different exchange hole models. The ingredients are:
Definition xc_xbr_pbe_lda_hole_t_c_lr.F:29
xc_xbr_pbe_lda_hole_t_c_lr::xbr_pbe_lda_hole_tc_lr_lsd_eval
subroutine, public xbr_pbe_lda_hole_tc_lr_lsd_eval(rho_set, deriv_set, grad_deriv, params)
Intermediate routine that gets grids, derivatives and some params.
Definition xc_xbr_pbe_lda_hole_t_c_lr.F:437
xc_xbr_pbe_lda_hole_t_c_lr::xbr_pbe_lda_hole_tc_lr_lda_eval
subroutine, public xbr_pbe_lda_hole_tc_lr_lda_eval(rho_set, deriv_set, grad_deriv, params)
Intermediate routine that gets grids, derivatives and some params.
Definition xc_xbr_pbe_lda_hole_t_c_lr.F:180
xc_xbr_pbe_lda_hole_t_c_lr::xbr_pbe_lda_hole_tc_lr_lda_info
subroutine, public xbr_pbe_lda_hole_tc_lr_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbr_pbe_lda_hole_t_c_lr.F:115
xc_xbr_pbe_lda_hole_t_c_lr::xbr_pbe_lda_hole_tc_lr_lsd_info
subroutine, public xbr_pbe_lda_hole_tc_lr_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbr_pbe_lda_hole_t_c_lr.F:147
xc_xlda_hole_t_c_lr
Calculates the lda exchange hole in a truncated coulomb potential. Can be used as longrange correctio...
Definition xc_xlda_hole_t_c_lr.F:16
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lda_calc_0
subroutine, public xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, sx, r)
low level routine
Definition xc_xlda_hole_t_c_lr.F:237
xc_xpbe_hole_t_c_lr
Calculates the exchange energy for the pbe hole model in a truncated coulomb potential,...
Definition xc_xpbe_hole_t_c_lr.F:17
xc_xpbe_hole_t_c_lr::xpbe_hole_t_c_lr_lda_calc_1
elemental subroutine, public xpbe_hole_t_c_lr_lda_calc_1(e_0, e_rho, e_ndrho, rho, ndrho, sscale, sx, r, order)
low level routine that calculates the energy derivatives in one point
Definition xc_xpbe_hole_t_c_lr.F:471
xc_xpbe_hole_t_c_lr::xpbe_hole_t_c_lr_lda_calc_2
elemental subroutine, public xpbe_hole_t_c_lr_lda_calc_2(e_0, e_rho, e_ndrho, rho, ndrho, sscale, sx, r, order)
low level routine that calculates the energy derivatives in one point
Definition xc_xpbe_hole_t_c_lr.F:838
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78