(git:ccc2433)
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 
29 MODULE 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
103 CONTAINS
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 
703 END 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
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_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_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_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
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